README revision 112158
1112158SdasThis directory contains source for a library of binary -> decimal 2112158Sdasand decimal -> binary conversion routines, for single-, double-, 3112158Sdasand extended-precision IEEE binary floating-point arithmetic, and 4112158Sdasother IEEE-like binary floating-point, including "double double", 5112158Sdasas in 6112158Sdas 7112158Sdas T. J. Dekker, "A Floating-Point Technique for Extending the 8112158Sdas Available Precision", Numer. Math. 18 (1971), pp. 224-242 9112158Sdas 10112158Sdasand 11112158Sdas 12112158Sdas "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994 13112158Sdas 14112158SdasThe conversion routines use double-precision floating-point arithmetic 15112158Sdasand, where necessary, high precision integer arithmetic. The routines 16112158Sdasare generalizations of the strtod and dtoa routines described in 17112158Sdas 18112158Sdas David M. Gay, "Correctly Rounded Binary-Decimal and 19112158Sdas Decimal-Binary Conversions", Numerical Analysis Manuscript 20112158Sdas No. 90-10, Bell Labs, Murray Hill, 1990; 21112158Sdas http://cm.bell-labs.com/cm/cs/what/ampl/REFS/rounding.ps.gz 22112158Sdas 23112158Sdas(based in part on papers by Clinger and Steele & White: see the 24112158Sdasreferences in the above paper). 25112158Sdas 26112158SdasThe present conversion routines should be able to use any of IEEE binary, 27112158SdasVAX, or IBM-mainframe double-precision arithmetic internally, but I (dmg) 28112158Sdashave so far only had a chance to test them with IEEE double precision 29112158Sdasarithmetic. 30112158Sdas 31112158SdasThe core conversion routines are strtodg for decimal -> binary conversions 32112158Sdasand gdtoa for binary -> decimal conversions. These routines operate 33112158Sdason arrays of unsigned 32-bit integers of type ULong, a signed 32-bit 34112158Sdasexponent of type Long, and arithmetic characteristics described in 35112158Sdasstruct FPI; FPI, Long, and ULong are defined in gdtoa.h. File arith.h 36112158Sdasis supposed to provide #defines that cause gdtoa.h to define its 37112158Sdastypes correctly. File arithchk.c is source for a program that 38112158Sdasgenerates a suitable arith.h on all systems where I've been able to 39112158Sdastest it. 40112158Sdas 41112158SdasThe core conversion routines are meant to be called by helper routines 42112158Sdasthat know details of the particular binary arithmetic of interest and 43112158Sdasconvert. The present directory provides helper routines for 5 variants 44112158Sdasof IEEE binary floating-point arithmetic, each indicated by one or 45112158Sdastwo letters: 46112158Sdas 47112158Sdas f IEEE single precision 48112158Sdas d IEEE double precision 49112158Sdas x IEEE extended precision, as on Intel 80x87 50112158Sdas and software emulations of Motorola 68xxx chips 51112158Sdas that do not pad the way the 68xxx does, but 52112158Sdas only store 80 bits 53112158Sdas xL IEEE extended precision, as on Motorola 68xxx chips 54112158Sdas Q quad precision, as on Sun Sparc chips 55112158Sdas dd double double, pairs of IEEE double numbers 56112158Sdas whose sum is the desired value 57112158Sdas 58112158SdasFor decimal -> binary conversions, there are three families of 59112158Sdashelper routines: one for round-nearest: 60112158Sdas 61112158Sdas strtof 62112158Sdas strtod 63112158Sdas strtodd 64112158Sdas strtopd 65112158Sdas strtopf 66112158Sdas strtopx 67112158Sdas strtopxL 68112158Sdas strtopQ 69112158Sdas 70112158Sdasone with rounding direction specified: 71112158Sdas 72112158Sdas strtorf 73112158Sdas strtord 74112158Sdas strtordd 75112158Sdas strtorx 76112158Sdas strtorxL 77112158Sdas strtorQ 78112158Sdas 79112158Sdasand one for computing an interval (at most one bit wide) that contains 80112158Sdasthe decimal number: 81112158Sdas 82112158Sdas strtoIf 83112158Sdas strtoId 84112158Sdas strtoIdd 85112158Sdas strtoIx 86112158Sdas strtoIxL 87112158Sdas strtoIQ 88112158Sdas 89112158SdasThe latter call strtoIg, which makes one call on strtodg and adjusts 90112158Sdasthe result to provide the desired interval. On systems where native 91112158Sdasarithmetic can easily make one-ulp adjustments on values in the 92112158Sdasdesired floating-point format, it might be more efficient to use the 93112158Sdasnative arithmetic. Routine strtodI is a variant of strtoId that 94112158Sdasillustrates one way to do this for IEEE binary double-precision 95112158Sdasarithmetic -- but whether this is more efficient remains to be seen. 96112158Sdas 97112158SdasFunctions strtod and strtof have "natural" return types, float and 98112158Sdasdouble -- strtod is specified by the C standard, and strtof appears 99112158Sdasin the stdlib.h of some systems, such as (at least some) Linux systems. 100112158SdasThe other functions write their results to their final argument(s): 101112158Sdasto the final two argument for the strtoI... (interval) functions, 102112158Sdasand to the final argument for the others (strtop... and strtor...). 103112158SdasWhere possible, these arguments have "natural" return types (double* 104112158Sdasor float*), to permit at least some type checking. In reality, they 105112158Sdasare viewed as arrays of ULong (or, for the "x" functions, UShort) 106112158Sdasvalues. On systems where long double is the appropriate type, one can 107112158Sdaspass long double* final argument(s) to these routines. The int value 108112158Sdasthat these routines return is the return value from the call they make 109112158Sdason strtodg; see the enum of possible return values in gdtoa.h. 110112158Sdas 111112158SdasSource files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c 112112158Sdasshould use true IEEE double arithmetic (not, e.g., double extended), 113112158Sdasat least for storing (and viewing the bits of) the variables declared 114112158Sdas"double" within them. 115112158Sdas 116112158SdasOne detail indicated in struct FPI is whether the target binary 117112158Sdasarithmetic departs from the IEEE standard by flushing denormalized 118112158Sdasnumbers to 0. On systems that do this, the helper routines for 119112158Sdasconversion to double-double format (when compiled with 120112158SdasSudden_Underflow #defined) penalize the bottom of the exponent 121112158Sdasrange so that they return a nonzero result only when the least 122112158Sdassignificant bit of the less significant member of the pair of 123112158Sdasdouble values returned can be expressed as a normalized double 124112158Sdasvalue. An alternative would be to drop to 53-bit precision near 125112158Sdasthe bottom of the exponent range. To get correct rounding, this 126112158Sdaswould (in general) require two calls on strtodg (one specifying 127112158Sdas126-bit arithmetic, then, if necessary, one specifying 53-bit 128112158Sdasarithmetic). 129112158Sdas 130112158SdasBy default, the core routine strtodg and strtod set errno to ERANGE 131112158Sdasif the result overflows to +Infinity or underflows to 0. Compile 132112158Sdasthese routines with NO_ERRNO #defined to inhibit errno assignments. 133112158Sdas 134112158SdasRoutine strtod is based on netlib's "dtoa.c from fp", and 135112158Sdas(f = strtod(s,se)) is more efficient for some conversions than, say, 136112158Sdasstrtord(s,se,1,&f). Parts of strtod require true IEEE double 137112158Sdasarithmetic with the default rounding mode (round-to-nearest) and, on 138112158Sdassystems with IEEE extended-precision registers, double-precision 139112158Sdas(53-bit) rounding precision. If the machine uses (the equivalent of) 140112158SdasIntel 80x87 arithmetic, the call 141112158Sdas _control87(PC_53, MCW_PC); 142112158Sdasdoes this with many compilers. Whether this or another call is 143112158Sdasappropriate depends on the compiler; for this to work, it may be 144112158Sdasnecessary to #include "float.h" or another system-dependent header 145112158Sdasfile. 146112158Sdas 147112158SdasThe values returned for NaNs may be signaling NaNs on some systems, 148112158Sdassince the rules for distinguishing signaling from quiet NaNs are 149112158Sdassystem-dependent. You can easily fix this by suitably modifying the 150112158SdasULto* routines in strtor*.c. 151112158Sdas 152112158SdasC99's hexadecimal floating-point constants are recognized by the 153112158Sdasstrto* routines (but this feature has not yet been heavily tested). 154112158SdasCompiling with NO_HEX_FP #defined disables this feature. 155112158Sdas 156112158SdasThe strto* routines do not (yet) recognize C99's NaN(...) syntax; the 157112158Sdasstrto* routines simply regard '(' as the first unprocessed input 158112158Sdascharacter. 159112158Sdas 160112158SdasFor binary -> decimal conversions, I've provided just one family 161112158Sdasof helper routines: 162112158Sdas 163112158Sdas g_ffmt 164112158Sdas g_dfmt 165112158Sdas g_ddfmt 166112158Sdas g_xfmt 167112158Sdas g_xLfmt 168112158Sdas g_Qfmt 169112158Sdas 170112158Sdaswhich do a "%g" style conversion either to a specified number of decimal 171112158Sdasplaces (if their ndig argument is positive), or to the shortest 172112158Sdasdecimal string that rounds to the given binary floating-point value 173112158Sdas(if ndig <= 0). They write into a buffer supplied as an argument 174112158Sdasand return either a pointer to the end of the string (a null character) 175112158Sdasin the buffer, if the buffer was long enough, or 0. Other forms of 176112158Sdasconversion are easily done with the help of gdtoa(), such as %e or %f 177112158Sdasstyle and conversions with direction of rounding specified (so that, if 178112158Sdasdesired, the decimal value is either >= or <= the binary value). 179112158Sdas 180112158SdasFor an example of more general conversions based on dtoa(), see 181112158Sdasnetlib's "printf.c from ampl/solvers". 182112158Sdas 183112158SdasFor double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic 184112158Sdasof precision max(126, #bits(input)) bits, where #bits(input) is the 185112158Sdasnumber of mantissa bits needed to represent the sum of the two double 186112158Sdasvalues in the input. 187112158Sdas 188112158SdasThe makefile creates a library, gdtoa.a. To use the helper 189112158Sdasroutines, a program only needs to include gdtoa.h. All the 190112158Sdassource files for gdtoa.a include a more extensive gdtoaimp.h; 191112158Sdasamong other things, gdtoaimp.h has #defines that make "internal" 192112158Sdasnames end in _D2A. To make a "system" library, one could modify 193112158Sdasthese #defines to make the names start with __. 194112158Sdas 195112158SdasVarious comments about possible #defines appear in gdtoaimp.h, 196112158Sdasbut for most purposes, arith.h should set suitable #defines. 197112158Sdas 198112158SdasSystems with preemptive scheduling of multiple threads require some 199112158Sdasmanual intervention. On such systems, it's necessary to compile 200112158Sdasdmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined, 201112158Sdasand to provide (or suitably #define) two locks, acquired by 202112158SdasACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1. 203112158Sdas(The second lock, accessed in pow5mult, ensures lazy evaluation of 204112158Sdasonly one copy of high powers of 5; omitting this lock would introduce 205112158Sdasa small probability of wasting memory, but would otherwise be harmless.) 206112158SdasRoutines that call dtoa or gdtoa directly must also invoke freedtoa(s) 207112158Sdasto free the value s returned by dtoa or gdtoa. It's OK to do so whether 208112158Sdasor not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines 209112158Sdaslisted above all do this indirectly (in gfmt_D2A(), which they all call). 210112158Sdas 211112158SdasBy default, there is a private pool of memory of length 2000 bytes 212112158Sdasfor intermediate quantities, and MALLOC (see gdtoaimp.h) is called only 213112158Sdasif the private pool does not suffice. 2000 is large enough that MALLOC 214112158Sdasis called only under very unusual circumstances (decimal -> binary 215112158Sdasconversion of very long strings) for conversions to and from double 216112158Sdasprecision. For systems with preemptivaly scheduled multiple threads 217112158Sdasor for conversions to extended or quad, it may be appropriate to 218112158Sdas#define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000. 219112158SdasFor extended and quad precisions, -DPRIVATE_MEM=20000 is probably 220112158Sdasplenty even for many digits at the ends of the exponent range. 221112158SdasUse of the private pool avoids some overhead. 222112158Sdas 223112158SdasDirectory test provides some test routines. See its README. 224112158SdasI've also tested this stuff (except double double conversions) 225112158Sdaswith Vern Paxson's testbase program: see 226112158Sdas 227112158Sdas V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal 228112158Sdas Conversion", manuscript, May 1991, 229112158Sdas ftp://ftp.ee.lbl.gov/testbase-report.ps.Z . 230112158Sdas 231112158Sdas(The same ftp directory has source for testbase.) 232112158Sdas 233112158SdasSome system-dependent additions to CFLAGS in the makefile: 234112158Sdas 235112158Sdas HU-UX: -Aa -Ae 236112158Sdas OSF (DEC Unix): -ieee_with_no_inexact 237112158Sdas SunOS 4.1x: -DKR_headers -DBad_float_h 238112158Sdas 239112158SdasIf you want to put this stuff into a shared library and your 240112158Sdasoperating system requires export lists for shared libraries, 241112158Sdasthe following would be an appropriate export list: 242112158Sdas 243112158Sdas dtoa 244112158Sdas freedtoa 245112158Sdas g_Qfmt 246112158Sdas g_ddfmt 247112158Sdas g_dfmt 248112158Sdas g_ffmt 249112158Sdas g_xLfmt 250112158Sdas g_xfmt 251112158Sdas gdtoa 252112158Sdas strtoIQ 253112158Sdas strtoId 254112158Sdas strtoIdd 255112158Sdas strtoIf 256112158Sdas strtoIx 257112158Sdas strtoIxL 258112158Sdas strtod 259112158Sdas strtodI 260112158Sdas strtodg 261112158Sdas strtof 262112158Sdas strtopQ 263112158Sdas strtopd 264112158Sdas strtopdd 265112158Sdas strtopf 266112158Sdas strtopx 267112158Sdas strtopxL 268112158Sdas strtorQ 269112158Sdas strtord 270112158Sdas strtordd 271112158Sdas strtorf 272112158Sdas strtorx 273112158Sdas strtorxL 274112158Sdas 275112158SdasWhen time permits, I (dmg) hope to write in more detail about the 276112158Sdaspresent conversion routines; for now, this README file must suffice. 277112158SdasMeanwhile, if you wish to write helper functions for other kinds of 278112158SdasIEEE-like arithmetic, some explanation of struct FPI and the bits 279112158Sdasarray may be helpful. Both gdtoa and strtodg operate on a bits array 280112158Sdasdescribed by FPI *fpi. The bits array is of type ULong, a 32-bit 281112158Sdasunsigned integer type. Floating-point numbers have fpi->nbits bits, 282112158Sdaswith the least significant 32 bits in bits[0], the next 32 bits in 283112158Sdasbits[1], etc. These numbers are regarded as integers multiplied by 284112158Sdas2^e (i.e., 2 to the power of the exponent e), where e is the second 285112158Sdasargument (be) to gdtoa and is stored in *exp by strtodg. The minimum 286112158Sdasand maximum exponent values fpi->emin and fpi->emax for normalized 287112158Sdasfloating-point numbers reflect this arrangement. For example, the 288112158SdasP754 standard for binary IEEE arithmetic specifies doubles as having 289112158Sdas53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023), 290112158Sdaswith 52 bits (the x's) and the biased exponent b represented explicitly; 291112158Sdasb is an unsigned integer in the range 1 <= b <= 2046 for normalized 292112158Sdasfinite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs. 293112158SdasTo turn an IEEE double into the representation used by strtodg and gdtoa, 294112158Sdaswe multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the 295112158Sdasexponent e = (b-1023) by 52: 296112158Sdas 297112158Sdas fpi->emin = 1 - 1023 - 52 298112158Sdas fpi->emax = 1046 - 1023 - 52 299112158Sdas 300112158SdasIn various wrappers for IEEE double, we actually write -53 + 1 rather 301112158Sdasthan -52, to emphasize that there are 53 bits including one implicit bit. 302112158SdasField fpi->rounding indicates the desired rounding direction, with 303112158Sdaspossible values 304112158Sdas FPI_Round_zero = toward 0, 305112158Sdas FPI_Round_near = unbiased rounding -- the IEEE default, 306112158Sdas FPI_Round_up = toward +Infinity, and 307112158Sdas FPI_Round_down = toward -Infinity 308112158Sdasgiven in gdtoa.h. 309112158Sdas 310112158SdasField fpi->sudden_underflow indicates whether strtodg should return 311112158Sdasdenormals or flush them to zero. Normal floating-point numbers have 312112158Sdasbit fpi->nbits in the bits array on. Denormals have it off, with 313112158Sdasexponent = fpi->emin. Strtodg provides distinct return values for normals 314112158Sdasand denormals; see gdtoa.h. 315112158Sdas 316112158SdasPlease send comments to 317112158Sdas 318112158Sdas David M. Gay 319112158Sdas Bell Labs, Room 2C-463 320112158Sdas 600 Mountain Avenue 321112158Sdas Murray Hill, NJ 07974-0636, U.S.A. 322112158Sdas dmg@research.bell-labs.com 323