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 Source file strtodnrp.c gives a strtod that does not require 53-bit
 148 rounding precision on systems (such as Intel IA32 systems) that may
 149 suffer double rounding due to use of extended-precision registers.
 150 For some conversions this variant of strtod is less efficient than the
 151 one in strtod.c when the latter is run with 53-bit rounding precision.
 152
 153 The values that the strto* routines return for NaNs are determined by
 154 gd_qnan.h, which the makefile generates by running the program whose
 155 source is qnan.c.  Note that the rules for distinguishing signaling
 156 from quiet NaNs are system-dependent.  For cross-compilation, you need
 157 to determine arith.h and gd_qnan.h suitably, e.g., using the
 158 arithmetic of the target machine.
 159
 160 C99's hexadecimal floating-point constants are recognized by the
 161 strto* routines (but this feature has not yet been heavily tested).
 162 Compiling with NO_HEX_FP #defined disables this feature.
 163
 164 When compiled with -DINFNAN_CHECK, the strto* routines recognize C99's
 165 NaN and Infinity syntax.  Moreover, unless No_Hex_NaN is #defined, the
 166 strto* routines also recognize C99's NaN(...) syntax: they accept
 167 (case insensitively) strings of the form NaN(x), where x is a string
 168 of hexadecimal digits and spaces; if there is only one string of
 169 hexadecimal digits, it is taken for the fraction bits of the resulting
 170 NaN; if there are two or more strings of hexadecimal digits, each
 171 string is assigned to the next available sequence of 32-bit words of
 172 fractions bits (starting with the most significant), right-aligned in
 173 each sequence.
 174
 175 For binary -> decimal conversions, I've provided just one family
 176 of helper routines:
 177
 178         g_ffmt
 179         g_dfmt
 180         g_ddfmt
 181         g_xfmt
 182         g_xLfmt
 183         g_Qfmt
 184
 185 which do a "%g" style conversion either to a specified number of decimal
 186 places (if their ndig argument is positive), or to the shortest
 187 decimal string that rounds to the given binary floating-point value
 188 (if ndig <= 0).  They write into a buffer supplied as an argument
 189 and return either a pointer to the end of the string (a null character)
 190 in the buffer, if the buffer was long enough, or 0.  Other forms of
 191 conversion are easily done with the help of gdtoa(), such as %e or %f
 192 style and conversions with direction of rounding specified (so that, if
 193 desired, the decimal value is either >= or <= the binary value).
 194
 195 For an example of more general conversions based on dtoa(), see
 196 netlib's "printf.c from ampl/solvers".
 197
 198 For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic
 199 of precision max(126, #bits(input)) bits, where #bits(input) is the
 200 number of mantissa bits needed to represent the sum of the two double
 201 values in the input.
 202
 203 The makefile creates a library, gdtoa.a.  To use the helper
 204 routines, a program only needs to include gdtoa.h.  All the
 205 source files for gdtoa.a include a more extensive gdtoaimp.h;
 206 among other things, gdtoaimp.h has #defines that make "internal"
 207 names end in _D2A.  To make a "system" library, one could modify
 208 these #defines to make the names start with __.
 209
 210 Various comments about possible #defines appear in gdtoaimp.h,
 211 but for most purposes, arith.h should set suitable #defines.
 212
 213 Systems with preemptive scheduling of multiple threads require some
 214 manual intervention.  On such systems, it's necessary to compile
 215 dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined,
 216 and to provide (or suitably #define) two locks, acquired by
 217 ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1.
 218 (The second lock, accessed in pow5mult, ensures lazy evaluation of
 219 only one copy of high powers of 5; omitting this lock would introduce
 220 a small probability of wasting memory, but would otherwise be harmless.)
 221 Routines that call dtoa or gdtoa directly must also invoke freedtoa(s)
 222 to free the value s returned by dtoa or gdtoa.  It's OK to do so whether
 223 or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines
 224 listed above all do this indirectly (in gfmt_D2A(), which they all call).
 225
 226 By default, there is a private pool of memory of length 2000 bytes
 227 for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only
 228 if the private pool does not suffice.   2000 is large enough that MALLOC
 229 is called only under very unusual circumstances (decimal -> binary
 230 conversion of very long strings) for conversions to and from double
 231 precision.  For systems with preemptively scheduled multiple threads
 232 or for conversions to extended or quad, it may be appropriate to
 233 #define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000.
 234 For extended and quad precisions, -DPRIVATE_MEM=20000 is probably
 235 plenty even for many digits at the ends of the exponent range.
 236 Use of the private pool avoids some overhead.
 237
 238 Directory test provides some test routines.  See its README.
 239 I've also tested this stuff (except double double conversions)
 240 with Vern Paxson's testbase program: see
 241
 242         V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal
 243         Conversion", manuscript, May 1991,
 244         ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
 245
 246 (The same ftp directory has source for testbase.)
 247
 248 Some system-dependent additions to CFLAGS in the makefile:
 249
 250         HU-UX: -Aa -Ae
 251         OSF (DEC Unix): -ieee_with_no_inexact
 252         SunOS 4.1x: -DKR_headers -DBad_float_h
 253
 254 If you want to put this stuff into a shared library and your
 255 operating system requires export lists for shared libraries,
 256 the following would be an appropriate export list:
 257
 258         dtoa
 259         freedtoa
 260         g_Qfmt
 261         g_ddfmt
 262         g_dfmt
 263         g_ffmt
 264         g_xLfmt
 265         g_xfmt
 266         gdtoa
 267         strtoIQ
 268         strtoId
 269         strtoIdd
 270         strtoIf
 271         strtoIx
 272         strtoIxL
 273         strtod
 274         strtodI
 275         strtodg
 276         strtof
 277         strtopQ
 278         strtopd
 279         strtopdd
 280         strtopf
 281         strtopx
 282         strtopxL
 283         strtorQ
 284         strtord
 285         strtordd
 286         strtorf
 287         strtorx
 288         strtorxL
 289
 290 When time permits, I (dmg) hope to write in more detail about the
 291 present conversion routines; for now, this README file must suffice.
 292 Meanwhile, if you wish to write helper functions for other kinds of
 293 IEEE-like arithmetic, some explanation of struct FPI and the bits
 294 array may be helpful.  Both gdtoa and strtodg operate on a bits array
 295 described by FPI *fpi.  The bits array is of type ULong, a 32-bit
 296 unsigned integer type.  Floating-point numbers have fpi->nbits bits,
 297 with the least significant 32 bits in bits[0], the next 32 bits in
 298 bits[1], etc.  These numbers are regarded as integers multiplied by
 299 2^e (i.e., 2 to the power of the exponent e), where e is the second
 300 argument (be) to gdtoa and is stored in *exp by strtodg.  The minimum
 301 and maximum exponent values fpi->emin and fpi->emax for normalized
 302 floating-point numbers reflect this arrangement.  For example, the
 303 P754 standard for binary IEEE arithmetic specifies doubles as having
 304 53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023),
 305 with 52 bits (the x's) and the biased exponent b represented explicitly;
 306 b is an unsigned integer in the range 1 <= b <= 2046 for normalized
 307 finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs.
 308 To turn an IEEE double into the representation used by strtodg and gdtoa,
 309 we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the
 310 exponent e = (b-1023) by 52:
 311
 312         fpi->emin = 1 - 1023 - 52
 313         fpi->emax = 1046 - 1023 - 52
 314
 315 In various wrappers for IEEE double, we actually write -53 + 1 rather
 316 than -52, to emphasize that there are 53 bits including one implicit bit.
 317 Field fpi->rounding indicates the desired rounding direction, with
 318 possible values
 319         FPI_Round_zero = toward 0,
 320         FPI_Round_near = unbiased rounding -- the IEEE default,
 321         FPI_Round_up = toward +Infinity, and
 322         FPI_Round_down = toward -Infinity
 323 given in gdtoa.h.
 324
 325 Field fpi->sudden_underflow indicates whether strtodg should return
 326 denormals or flush them to zero.  Normal floating-point numbers have
 327 bit fpi->nbits in the bits array on.  Denormals have it off, with
 328 exponent = fpi->emin.  Strtodg provides distinct return values for normals
 329 and denormals; see gdtoa.h.
 330
 331 Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
 332 the decimal-point character to be taken from the current locale; otherwise
 333 it is '.'.
 334
 335 Please send comments to David M. Gay (dmg at acm dot org, with " at "
 336 changed at "@" and " dot " changed to ".").