contrib/gdtoa/README

   1 This directory contains source for a library of binary -> decimal
   2 and decimal -> binary conversion routines, for single-, double-,
   3 and extended-precision IEEE binary floating-point arithmetic, and
   4 other IEEE-like binary floating-point, including "double double",
   5 as in
   6
   7         T. J. Dekker, "A Floating-Point Technique for Extending the
   8         Available Precision", Numer. Math. 18 (1971), pp. 224-242
   9
  10 and
  11
  12         "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994
  13
  14 The conversion routines use double-precision floating-point arithmetic
  15 and, where necessary, high precision integer arithmetic.  The routines
  16 are generalizations of the strtod and dtoa routines described in
  17
  18         David M. Gay, "Correctly Rounded Binary-Decimal and
  19         Decimal-Binary Conversions", Numerical Analysis Manuscript
  20         No. 90-10, Bell Labs, Murray Hill, 1990;
  21         http://cm.bell-labs.com/cm/cs/what/ampl/REFS/rounding.ps.gz
  22
  23 (based in part on papers by Clinger and Steele & White: see the
  24 references in the above paper).
  25
  26 The present conversion routines should be able to use any of IEEE binary,
  27 VAX, or IBM-mainframe double-precision arithmetic internally, but I (dmg)
  28 have so far only had a chance to test them with IEEE double precision
  29 arithmetic.
  30
  31 The core conversion routines are strtodg for decimal -> binary conversions
  32 and gdtoa for binary -> decimal conversions.  These routines operate
  33 on arrays of unsigned 32-bit integers of type ULong, a signed 32-bit
  34 exponent of type Long, and arithmetic characteristics described in
  35 struct FPI; FPI, Long, and ULong are defined in gdtoa.h.  File arith.h
  36 is supposed to provide #defines that cause gdtoa.h to define its
  37 types correctly.  File arithchk.c is source for a program that
  38 generates a suitable arith.h on all systems where I've been able to
  39 test it.
  40
  41 The core conversion routines are meant to be called by helper routines
  42 that know details of the particular binary arithmetic of interest and
  43 convert.  The present directory provides helper routines for 5 variants
  44 of IEEE binary floating-point arithmetic, each indicated by one or
  45 two letters:
  46
  47         f       IEEE single precision
  48         d       IEEE double precision
  49         x       IEEE extended precision, as on Intel 80x87
  50                 and software emulations of Motorola 68xxx chips
  51                 that do not pad the way the 68xxx does, but
  52                 only store 80 bits
  53         xL      IEEE extended precision, as on Motorola 68xxx chips
  54         Q       quad precision, as on Sun Sparc chips
  55         dd      double double, pairs of IEEE double numbers
  56                 whose sum is the desired value
  57
  58 For decimal -> binary conversions, there are three families of
  59 helper routines: one for round-nearest:
  60
  61         strtof
  62         strtod
  63         strtodd
  64         strtopd
  65         strtopf
  66         strtopx
  67         strtopxL
  68         strtopQ
  69
  70 one with rounding direction specified:
  71
  72         strtorf
  73         strtord
  74         strtordd
  75         strtorx
  76         strtorxL
  77         strtorQ
  78
  79 and one for computing an interval (at most one bit wide) that contains
  80 the decimal number:
  81
  82         strtoIf
  83         strtoId
  84         strtoIdd
  85         strtoIx
  86         strtoIxL
  87         strtoIQ
  88
  89 The latter call strtoIg, which makes one call on strtodg and adjusts
  90 the result to provide the desired interval.  On systems where native
  91 arithmetic can easily make one-ulp adjustments on values in the
  92 desired floating-point format, it might be more efficient to use the
  93 native arithmetic.  Routine strtodI is a variant of strtoId that
  94 illustrates one way to do this for IEEE binary double-precision
  95 arithmetic -- but whether this is more efficient remains to be seen.
  96
  97 Functions strtod and strtof have "natural" return types, float and
  98 double -- strtod is specified by the C standard, and strtof appears
  99 in the stdlib.h of some systems, such as (at least some) Linux systems.
 100 The other functions write their results to their final argument(s):
 101 to the final two argument for the strtoI... (interval) functions,
 102 and to the final argument for the others (strtop... and strtor...).
 103 Where possible, these arguments have "natural" return types (double*
 104 or float*), to permit at least some type checking.  In reality, they
 105 are viewed as arrays of ULong (or, for the "x" functions, UShort)
 106 values. On systems where long double is the appropriate type, one can
 107 pass long double* final argument(s) to these routines.  The int value
 108 that these routines return is the return value from the call they make
 109 on strtodg; see the enum of possible return values in gdtoa.h.
 110
 111 Source files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c
 112 should use true IEEE double arithmetic (not, e.g., double extended),
 113 at least for storing (and viewing the bits of) the variables declared
 114 "double" within them.
 115
 116 One detail indicated in struct FPI is whether the target binary
 117 arithmetic departs from the IEEE standard by flushing denormalized
 118 numbers to 0.  On systems that do this, the helper routines for
 119 conversion to double-double format (when compiled with
 120 Sudden_Underflow #defined) penalize the bottom of the exponent
 121 range so that they return a nonzero result only when the least
 122 significant bit of the less significant member of the pair of
 123 double values returned can be expressed as a normalized double
 124 value.  An alternative would be to drop to 53-bit precision near
 125 the bottom of the exponent range.  To get correct rounding, this
 126 would (in general) require two calls on strtodg (one specifying
 127 126-bit arithmetic, then, if necessary, one specifying 53-bit
 128 arithmetic).
 129
 130 By default, the core routine strtodg and strtod set errno to ERANGE
 131 if the result overflows to +Infinity or underflows to 0.  Compile
 132 these routines with NO_ERRNO #defined to inhibit errno assignments.
 133
 134 Routine strtod is based on netlib's "dtoa.c from fp", and
 135 (f = strtod(s,se)) is more efficient for some conversions than, say,
 136 strtord(s,se,1,&f).  Parts of strtod require true IEEE double
 137 arithmetic with the default rounding mode (round-to-nearest) and, on
 138 systems with IEEE extended-precision registers, double-precision
 139 (53-bit) rounding precision.  If the machine uses (the equivalent of)
 140 Intel 80x87 arithmetic, the call
 141         _control87(PC_53, MCW_PC);
 142 does this with many compilers.  Whether this or another call is
 143 appropriate depends on the compiler; for this to work, it may be
 144 necessary to #include "float.h" or another system-dependent header
 145 file.
 146
 147 The values returned for NaNs may be signaling NaNs on some systems,
 148 since the rules for distinguishing signaling from quiet NaNs are
 149 system-dependent.  You can easily fix this by suitably modifying the
 150 ULto* routines in strtor*.c.
 151
 152 C99's hexadecimal floating-point constants are recognized by the
 153 strto* routines (but this feature has not yet been heavily tested).
 154 Compiling with NO_HEX_FP #defined disables this feature.
 155
 156 The strto* routines do not (yet) recognize C99's NaN(...) syntax; the
 157 strto* routines simply regard '(' as the first unprocessed input
 158 character.
 159
 160 For binary -> decimal conversions, I've provided just one family
 161 of helper routines:
 162
 163         g_ffmt
 164         g_dfmt
 165         g_ddfmt
 166         g_xfmt
 167         g_xLfmt
 168         g_Qfmt
 169
 170 which do a "%g" style conversion either to a specified number of decimal
 171 places (if their ndig argument is positive), or to the shortest
 172 decimal string that rounds to the given binary floating-point value
 173 (if ndig <= 0).  They write into a buffer supplied as an argument
 174 and return either a pointer to the end of the string (a null character)
 175 in the buffer, if the buffer was long enough, or 0.  Other forms of
 176 conversion are easily done with the help of gdtoa(), such as %e or %f
 177 style and conversions with direction of rounding specified (so that, if
 178 desired, the decimal value is either >= or <= the binary value).
 179
 180 For an example of more general conversions based on dtoa(), see
 181 netlib's "printf.c from ampl/solvers".
 182
 183 For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic
 184 of precision max(126, #bits(input)) bits, where #bits(input) is the
 185 number of mantissa bits needed to represent the sum of the two double
 186 values in the input.
 187
 188 The makefile creates a library, gdtoa.a.  To use the helper
 189 routines, a program only needs to include gdtoa.h.  All the
 190 source files for gdtoa.a include a more extensive gdtoaimp.h;
 191 among other things, gdtoaimp.h has #defines that make "internal"
 192 names end in _D2A.  To make a "system" library, one could modify
 193 these #defines to make the names start with __.
 194
 195 Various comments about possible #defines appear in gdtoaimp.h,
 196 but for most purposes, arith.h should set suitable #defines.
 197
 198 Systems with preemptive scheduling of multiple threads require some
 199 manual intervention.  On such systems, it's necessary to compile
 200 dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined,
 201 and to provide (or suitably #define) two locks, acquired by
 202 ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1.
 203 (The second lock, accessed in pow5mult, ensures lazy evaluation of
 204 only one copy of high powers of 5; omitting this lock would introduce
 205 a small probability of wasting memory, but would otherwise be harmless.)
 206 Routines that call dtoa or gdtoa directly must also invoke freedtoa(s)
 207 to free the value s returned by dtoa or gdtoa.  It's OK to do so whether
 208 or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines
 209 listed above all do this indirectly (in gfmt_D2A(), which they all call).
 210
 211 By default, there is a private pool of memory of length 2000 bytes
 212 for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only
 213 if the private pool does not suffice.   2000 is large enough that MALLOC
 214 is called only under very unusual circumstances (decimal -> binary
 215 conversion of very long strings) for conversions to and from double
 216 precision.  For systems with preemptivaly scheduled multiple threads
 217 or for conversions to extended or quad, it may be appropriate to
 218 #define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000.
 219 For extended and quad precisions, -DPRIVATE_MEM=20000 is probably
 220 plenty even for many digits at the ends of the exponent range.
 221 Use of the private pool avoids some overhead.
 222
 223 Directory test provides some test routines.  See its README.
 224 I've also tested this stuff (except double double conversions)
 225 with Vern Paxson's testbase program: see
 226
 227         V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal
 228         Conversion", manuscript, May 1991,
 229         ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
 230
 231 (The same ftp directory has source for testbase.)
 232
 233 Some system-dependent additions to CFLAGS in the makefile:
 234
 235         HU-UX: -Aa -Ae
 236         OSF (DEC Unix): -ieee_with_no_inexact
 237         SunOS 4.1x: -DKR_headers -DBad_float_h
 238
 239 If you want to put this stuff into a shared library and your
 240 operating system requires export lists for shared libraries,
 241 the following would be an appropriate export list:
 242
 243         dtoa
 244         freedtoa
 245         g_Qfmt
 246         g_ddfmt
 247         g_dfmt
 248         g_ffmt
 249         g_xLfmt
 250         g_xfmt
 251         gdtoa
 252         strtoIQ
 253         strtoId
 254         strtoIdd
 255         strtoIf
 256         strtoIx
 257         strtoIxL
 258         strtod
 259         strtodI
 260         strtodg
 261         strtof
 262         strtopQ
 263         strtopd
 264         strtopdd
 265         strtopf
 266         strtopx
 267         strtopxL
 268         strtorQ
 269         strtord
 270         strtordd
 271         strtorf
 272         strtorx
 273         strtorxL
 274
 275 When time permits, I (dmg) hope to write in more detail about the
 276 present conversion routines; for now, this README file must suffice.
 277 Meanwhile, if you wish to write helper functions for other kinds of
 278 IEEE-like arithmetic, some explanation of struct FPI and the bits
 279 array may be helpful.  Both gdtoa and strtodg operate on a bits array
 280 described by FPI *fpi.  The bits array is of type ULong, a 32-bit
 281 unsigned integer type.  Floating-point numbers have fpi->nbits bits,
 282 with the least significant 32 bits in bits[0], the next 32 bits in
 283 bits[1], etc.  These numbers are regarded as integers multiplied by
 284 2^e (i.e., 2 to the power of the exponent e), where e is the second
 285 argument (be) to gdtoa and is stored in *exp by strtodg.  The minimum
 286 and maximum exponent values fpi->emin and fpi->emax for normalized
 287 floating-point numbers reflect this arrangement.  For example, the
 288 P754 standard for binary IEEE arithmetic specifies doubles as having
 289 53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023),
 290 with 52 bits (the x's) and the biased exponent b represented explicitly;
 291 b is an unsigned integer in the range 1 <= b <= 2046 for normalized
 292 finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs.
 293 To turn an IEEE double into the representation used by strtodg and gdtoa,
 294 we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the
 295 exponent e = (b-1023) by 52:
 296
 297         fpi->emin = 1 - 1023 - 52
 298         fpi->emax = 1046 - 1023 - 52
 299
 300 In various wrappers for IEEE double, we actually write -53 + 1 rather
 301 than -52, to emphasize that there are 53 bits including one implicit bit.
 302 Field fpi->rounding indicates the desired rounding direction, with
 303 possible values
 304         FPI_Round_zero = toward 0,
 305         FPI_Round_near = unbiased rounding -- the IEEE default,
 306         FPI_Round_up = toward +Infinity, and
 307         FPI_Round_down = toward -Infinity
 308 given in gdtoa.h.
 309
 310 Field fpi->sudden_underflow indicates whether strtodg should return
 311 denormals or flush them to zero.  Normal floating-point numbers have
 312 bit fpi->nbits in the bits array on.  Denormals have it off, with
 313 exponent = fpi->emin.  Strtodg provides distinct return values for normals
 314 and denormals; see gdtoa.h.
 315
 316 Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
 317 the decimal-point character to be taken from the current locale; otherwise
 318 it is '.'.
 319
 320 Please send comments to
 321
 322         David M. Gay
 323         dmg@acm.org