README revision 182709
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 147165743SdasSource file strtodnrp.c gives a strtod that does not require 53-bit 148165743Sdasrounding precision on systems (such as Intel IA32 systems) that may 149165743Sdassuffer double rounding due to use of extended-precision registers. 150165743SdasFor some conversions this variant of strtod is less efficient than the 151165743Sdasone in strtod.c when the latter is run with 53-bit rounding precision. 152112158Sdas 153165743SdasThe values that the strto* routines return for NaNs are determined by 154165743Sdasgd_qnan.h, which the makefile generates by running the program whose 155165743Sdassource is qnan.c. Note that the rules for distinguishing signaling 156165743Sdasfrom quiet NaNs are system-dependent. For cross-compilation, you need 157165743Sdasto determine arith.h and gd_qnan.h suitably, e.g., using the 158165743Sdasarithmetic of the target machine. 159165743Sdas 160112158SdasC99's hexadecimal floating-point constants are recognized by the 161112158Sdasstrto* routines (but this feature has not yet been heavily tested). 162112158SdasCompiling with NO_HEX_FP #defined disables this feature. 163112158Sdas 164165743SdasWhen compiled with -DINFNAN_CHECK, the strto* routines recognize C99's 165165743SdasNaN and Infinity syntax. Moreover, unless No_Hex_NaN is #defined, the 166165743Sdasstrto* routines also recognize C99's NaN(...) syntax: they accept 167165743Sdas(case insensitively) strings of the form NaN(x), where x is a string 168165743Sdasof hexadecimal digits and spaces; if there is only one string of 169165743Sdashexadecimal digits, it is taken for the fraction bits of the resulting 170165743SdasNaN; if there are two or more strings of hexadecimal digits, each 171165743Sdasstring is assigned to the next available sequence of 32-bit words of 172165743Sdasfractions bits (starting with the most significant), right-aligned in 173165743Sdaseach sequence. 174112158Sdas 175112158SdasFor binary -> decimal conversions, I've provided just one family 176112158Sdasof helper routines: 177112158Sdas 178112158Sdas g_ffmt 179112158Sdas g_dfmt 180112158Sdas g_ddfmt 181112158Sdas g_xfmt 182112158Sdas g_xLfmt 183112158Sdas g_Qfmt 184112158Sdas 185112158Sdaswhich do a "%g" style conversion either to a specified number of decimal 186112158Sdasplaces (if their ndig argument is positive), or to the shortest 187112158Sdasdecimal string that rounds to the given binary floating-point value 188112158Sdas(if ndig <= 0). They write into a buffer supplied as an argument 189112158Sdasand return either a pointer to the end of the string (a null character) 190112158Sdasin the buffer, if the buffer was long enough, or 0. Other forms of 191112158Sdasconversion are easily done with the help of gdtoa(), such as %e or %f 192112158Sdasstyle and conversions with direction of rounding specified (so that, if 193112158Sdasdesired, the decimal value is either >= or <= the binary value). 194112158Sdas 195112158SdasFor an example of more general conversions based on dtoa(), see 196112158Sdasnetlib's "printf.c from ampl/solvers". 197112158Sdas 198112158SdasFor double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic 199112158Sdasof precision max(126, #bits(input)) bits, where #bits(input) is the 200112158Sdasnumber of mantissa bits needed to represent the sum of the two double 201112158Sdasvalues in the input. 202112158Sdas 203112158SdasThe makefile creates a library, gdtoa.a. To use the helper 204112158Sdasroutines, a program only needs to include gdtoa.h. All the 205112158Sdassource files for gdtoa.a include a more extensive gdtoaimp.h; 206112158Sdasamong other things, gdtoaimp.h has #defines that make "internal" 207112158Sdasnames end in _D2A. To make a "system" library, one could modify 208112158Sdasthese #defines to make the names start with __. 209112158Sdas 210112158SdasVarious comments about possible #defines appear in gdtoaimp.h, 211112158Sdasbut for most purposes, arith.h should set suitable #defines. 212112158Sdas 213112158SdasSystems with preemptive scheduling of multiple threads require some 214112158Sdasmanual intervention. On such systems, it's necessary to compile 215112158Sdasdmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined, 216112158Sdasand to provide (or suitably #define) two locks, acquired by 217112158SdasACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1. 218112158Sdas(The second lock, accessed in pow5mult, ensures lazy evaluation of 219112158Sdasonly one copy of high powers of 5; omitting this lock would introduce 220112158Sdasa small probability of wasting memory, but would otherwise be harmless.) 221112158SdasRoutines that call dtoa or gdtoa directly must also invoke freedtoa(s) 222112158Sdasto free the value s returned by dtoa or gdtoa. It's OK to do so whether 223112158Sdasor not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines 224112158Sdaslisted above all do this indirectly (in gfmt_D2A(), which they all call). 225112158Sdas 226112158SdasBy default, there is a private pool of memory of length 2000 bytes 227112158Sdasfor intermediate quantities, and MALLOC (see gdtoaimp.h) is called only 228112158Sdasif the private pool does not suffice. 2000 is large enough that MALLOC 229112158Sdasis called only under very unusual circumstances (decimal -> binary 230112158Sdasconversion of very long strings) for conversions to and from double 231165743Sdasprecision. For systems with preemptively scheduled multiple threads 232112158Sdasor for conversions to extended or quad, it may be appropriate to 233112158Sdas#define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000. 234112158SdasFor extended and quad precisions, -DPRIVATE_MEM=20000 is probably 235112158Sdasplenty even for many digits at the ends of the exponent range. 236112158SdasUse of the private pool avoids some overhead. 237112158Sdas 238112158SdasDirectory test provides some test routines. See its README. 239112158SdasI've also tested this stuff (except double double conversions) 240112158Sdaswith Vern Paxson's testbase program: see 241112158Sdas 242112158Sdas V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal 243112158Sdas Conversion", manuscript, May 1991, 244112158Sdas ftp://ftp.ee.lbl.gov/testbase-report.ps.Z . 245112158Sdas 246112158Sdas(The same ftp directory has source for testbase.) 247112158Sdas 248112158SdasSome system-dependent additions to CFLAGS in the makefile: 249112158Sdas 250112158Sdas HU-UX: -Aa -Ae 251112158Sdas OSF (DEC Unix): -ieee_with_no_inexact 252112158Sdas SunOS 4.1x: -DKR_headers -DBad_float_h 253112158Sdas 254112158SdasIf you want to put this stuff into a shared library and your 255112158Sdasoperating system requires export lists for shared libraries, 256112158Sdasthe following would be an appropriate export list: 257112158Sdas 258112158Sdas dtoa 259112158Sdas freedtoa 260112158Sdas g_Qfmt 261112158Sdas g_ddfmt 262112158Sdas g_dfmt 263112158Sdas g_ffmt 264112158Sdas g_xLfmt 265112158Sdas g_xfmt 266112158Sdas gdtoa 267112158Sdas strtoIQ 268112158Sdas strtoId 269112158Sdas strtoIdd 270112158Sdas strtoIf 271112158Sdas strtoIx 272112158Sdas strtoIxL 273112158Sdas strtod 274112158Sdas strtodI 275112158Sdas strtodg 276112158Sdas strtof 277112158Sdas strtopQ 278112158Sdas strtopd 279112158Sdas strtopdd 280112158Sdas strtopf 281112158Sdas strtopx 282112158Sdas strtopxL 283112158Sdas strtorQ 284112158Sdas strtord 285112158Sdas strtordd 286112158Sdas strtorf 287112158Sdas strtorx 288112158Sdas strtorxL 289112158Sdas 290112158SdasWhen time permits, I (dmg) hope to write in more detail about the 291112158Sdaspresent conversion routines; for now, this README file must suffice. 292112158SdasMeanwhile, if you wish to write helper functions for other kinds of 293112158SdasIEEE-like arithmetic, some explanation of struct FPI and the bits 294112158Sdasarray may be helpful. Both gdtoa and strtodg operate on a bits array 295112158Sdasdescribed by FPI *fpi. The bits array is of type ULong, a 32-bit 296112158Sdasunsigned integer type. Floating-point numbers have fpi->nbits bits, 297112158Sdaswith the least significant 32 bits in bits[0], the next 32 bits in 298112158Sdasbits[1], etc. These numbers are regarded as integers multiplied by 299112158Sdas2^e (i.e., 2 to the power of the exponent e), where e is the second 300112158Sdasargument (be) to gdtoa and is stored in *exp by strtodg. The minimum 301112158Sdasand maximum exponent values fpi->emin and fpi->emax for normalized 302112158Sdasfloating-point numbers reflect this arrangement. For example, the 303112158SdasP754 standard for binary IEEE arithmetic specifies doubles as having 304112158Sdas53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023), 305112158Sdaswith 52 bits (the x's) and the biased exponent b represented explicitly; 306112158Sdasb is an unsigned integer in the range 1 <= b <= 2046 for normalized 307112158Sdasfinite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs. 308112158SdasTo turn an IEEE double into the representation used by strtodg and gdtoa, 309112158Sdaswe multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the 310112158Sdasexponent e = (b-1023) by 52: 311112158Sdas 312112158Sdas fpi->emin = 1 - 1023 - 52 313112158Sdas fpi->emax = 1046 - 1023 - 52 314112158Sdas 315112158SdasIn various wrappers for IEEE double, we actually write -53 + 1 rather 316112158Sdasthan -52, to emphasize that there are 53 bits including one implicit bit. 317112158SdasField fpi->rounding indicates the desired rounding direction, with 318112158Sdaspossible values 319112158Sdas FPI_Round_zero = toward 0, 320112158Sdas FPI_Round_near = unbiased rounding -- the IEEE default, 321112158Sdas FPI_Round_up = toward +Infinity, and 322112158Sdas FPI_Round_down = toward -Infinity 323112158Sdasgiven in gdtoa.h. 324112158Sdas 325112158SdasField fpi->sudden_underflow indicates whether strtodg should return 326112158Sdasdenormals or flush them to zero. Normal floating-point numbers have 327112158Sdasbit fpi->nbits in the bits array on. Denormals have it off, with 328112158Sdasexponent = fpi->emin. Strtodg provides distinct return values for normals 329112158Sdasand denormals; see gdtoa.h. 330112158Sdas 331112415SdasCompiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes 332112415Sdasthe decimal-point character to be taken from the current locale; otherwise 333112415Sdasit is '.'. 334112415Sdas 335182709SdasSource files dtoa.c and strtod.c in this directory are derived from 336182709Sdasnetlib's "dtoa.c from fp" and are meant to function equivalently. 337182709SdasWhen compiled with Honor_FLT_ROUNDS #defined (on systems that provide 338182709SdasFLT_ROUNDS and fegetround() as specified in the C99 standard), they 339182709Sdashonor the current rounding mode. Because FLT_ROUNDS is buggy on some 340182709Sdas(Linux) systems -- not reflecting calls on fesetround(), as the C99 341182709Sdasstandard says it should -- when Honor_FLT_ROUNDS is #defined, the 342182709Sdascurrent rounding mode is obtained from fegetround() rather than from 343182709SdasFLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined. 344182709Sdas 345165743SdasPlease send comments to David M. Gay (dmg at acm dot org, with " at " 346165743Sdaschanged at "@" and " dot " changed to "."). 347