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 59187808Sdashelper routines: one for round-nearest (or the current rounding 60187808Sdasmode on IEEE-arithmetic systems that provide the C99 fegetround() 61187808Sdasfunction, if compiled with -DHonor_FLT_ROUNDS): 62112158Sdas 63112158Sdas strtof 64112158Sdas strtod 65112158Sdas strtodd 66112158Sdas strtopd 67112158Sdas strtopf 68112158Sdas strtopx 69112158Sdas strtopxL 70112158Sdas strtopQ 71112158Sdas 72112158Sdasone with rounding direction specified: 73112158Sdas 74112158Sdas strtorf 75112158Sdas strtord 76112158Sdas strtordd 77112158Sdas strtorx 78112158Sdas strtorxL 79112158Sdas strtorQ 80112158Sdas 81112158Sdasand one for computing an interval (at most one bit wide) that contains 82112158Sdasthe decimal number: 83112158Sdas 84112158Sdas strtoIf 85112158Sdas strtoId 86112158Sdas strtoIdd 87112158Sdas strtoIx 88112158Sdas strtoIxL 89112158Sdas strtoIQ 90112158Sdas 91112158SdasThe latter call strtoIg, which makes one call on strtodg and adjusts 92112158Sdasthe result to provide the desired interval. On systems where native 93112158Sdasarithmetic can easily make one-ulp adjustments on values in the 94112158Sdasdesired floating-point format, it might be more efficient to use the 95112158Sdasnative arithmetic. Routine strtodI is a variant of strtoId that 96112158Sdasillustrates one way to do this for IEEE binary double-precision 97112158Sdasarithmetic -- but whether this is more efficient remains to be seen. 98112158Sdas 99112158SdasFunctions strtod and strtof have "natural" return types, float and 100112158Sdasdouble -- strtod is specified by the C standard, and strtof appears 101112158Sdasin the stdlib.h of some systems, such as (at least some) Linux systems. 102112158SdasThe other functions write their results to their final argument(s): 103112158Sdasto the final two argument for the strtoI... (interval) functions, 104112158Sdasand to the final argument for the others (strtop... and strtor...). 105112158SdasWhere possible, these arguments have "natural" return types (double* 106112158Sdasor float*), to permit at least some type checking. In reality, they 107112158Sdasare viewed as arrays of ULong (or, for the "x" functions, UShort) 108112158Sdasvalues. On systems where long double is the appropriate type, one can 109112158Sdaspass long double* final argument(s) to these routines. The int value 110112158Sdasthat these routines return is the return value from the call they make 111112158Sdason strtodg; see the enum of possible return values in gdtoa.h. 112112158Sdas 113112158SdasSource files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c 114112158Sdasshould use true IEEE double arithmetic (not, e.g., double extended), 115112158Sdasat least for storing (and viewing the bits of) the variables declared 116112158Sdas"double" within them. 117112158Sdas 118112158SdasOne detail indicated in struct FPI is whether the target binary 119112158Sdasarithmetic departs from the IEEE standard by flushing denormalized 120112158Sdasnumbers to 0. On systems that do this, the helper routines for 121112158Sdasconversion to double-double format (when compiled with 122112158SdasSudden_Underflow #defined) penalize the bottom of the exponent 123112158Sdasrange so that they return a nonzero result only when the least 124112158Sdassignificant bit of the less significant member of the pair of 125112158Sdasdouble values returned can be expressed as a normalized double 126112158Sdasvalue. An alternative would be to drop to 53-bit precision near 127112158Sdasthe bottom of the exponent range. To get correct rounding, this 128112158Sdaswould (in general) require two calls on strtodg (one specifying 129112158Sdas126-bit arithmetic, then, if necessary, one specifying 53-bit 130112158Sdasarithmetic). 131112158Sdas 132112158SdasBy default, the core routine strtodg and strtod set errno to ERANGE 133112158Sdasif the result overflows to +Infinity or underflows to 0. Compile 134112158Sdasthese routines with NO_ERRNO #defined to inhibit errno assignments. 135112158Sdas 136112158SdasRoutine strtod is based on netlib's "dtoa.c from fp", and 137112158Sdas(f = strtod(s,se)) is more efficient for some conversions than, say, 138112158Sdasstrtord(s,se,1,&f). Parts of strtod require true IEEE double 139112158Sdasarithmetic with the default rounding mode (round-to-nearest) and, on 140112158Sdassystems with IEEE extended-precision registers, double-precision 141112158Sdas(53-bit) rounding precision. If the machine uses (the equivalent of) 142112158SdasIntel 80x87 arithmetic, the call 143112158Sdas _control87(PC_53, MCW_PC); 144112158Sdasdoes this with many compilers. Whether this or another call is 145112158Sdasappropriate depends on the compiler; for this to work, it may be 146112158Sdasnecessary to #include "float.h" or another system-dependent header 147112158Sdasfile. 148112158Sdas 149165743SdasSource file strtodnrp.c gives a strtod that does not require 53-bit 150165743Sdasrounding precision on systems (such as Intel IA32 systems) that may 151165743Sdassuffer double rounding due to use of extended-precision registers. 152165743SdasFor some conversions this variant of strtod is less efficient than the 153165743Sdasone in strtod.c when the latter is run with 53-bit rounding precision. 154112158Sdas 155165743SdasThe values that the strto* routines return for NaNs are determined by 156165743Sdasgd_qnan.h, which the makefile generates by running the program whose 157165743Sdassource is qnan.c. Note that the rules for distinguishing signaling 158165743Sdasfrom quiet NaNs are system-dependent. For cross-compilation, you need 159165743Sdasto determine arith.h and gd_qnan.h suitably, e.g., using the 160165743Sdasarithmetic of the target machine. 161165743Sdas 162112158SdasC99's hexadecimal floating-point constants are recognized by the 163112158Sdasstrto* routines (but this feature has not yet been heavily tested). 164112158SdasCompiling with NO_HEX_FP #defined disables this feature. 165112158Sdas 166165743SdasWhen compiled with -DINFNAN_CHECK, the strto* routines recognize C99's 167165743SdasNaN and Infinity syntax. Moreover, unless No_Hex_NaN is #defined, the 168165743Sdasstrto* routines also recognize C99's NaN(...) syntax: they accept 169165743Sdas(case insensitively) strings of the form NaN(x), where x is a string 170165743Sdasof hexadecimal digits and spaces; if there is only one string of 171165743Sdashexadecimal digits, it is taken for the fraction bits of the resulting 172165743SdasNaN; if there are two or more strings of hexadecimal digits, each 173165743Sdasstring is assigned to the next available sequence of 32-bit words of 174165743Sdasfractions bits (starting with the most significant), right-aligned in 175165743Sdaseach sequence. 176112158Sdas 177112158SdasFor binary -> decimal conversions, I've provided just one family 178112158Sdasof helper routines: 179112158Sdas 180112158Sdas g_ffmt 181112158Sdas g_dfmt 182112158Sdas g_ddfmt 183112158Sdas g_xfmt 184112158Sdas g_xLfmt 185112158Sdas g_Qfmt 186112158Sdas 187112158Sdaswhich do a "%g" style conversion either to a specified number of decimal 188112158Sdasplaces (if their ndig argument is positive), or to the shortest 189112158Sdasdecimal string that rounds to the given binary floating-point value 190112158Sdas(if ndig <= 0). They write into a buffer supplied as an argument 191112158Sdasand return either a pointer to the end of the string (a null character) 192112158Sdasin the buffer, if the buffer was long enough, or 0. Other forms of 193112158Sdasconversion are easily done with the help of gdtoa(), such as %e or %f 194112158Sdasstyle and conversions with direction of rounding specified (so that, if 195112158Sdasdesired, the decimal value is either >= or <= the binary value). 196187808SdasOn IEEE-arithmetic systems that provide the C99 fegetround() function, 197187808Sdasif compiled with -DHonor_FLT_ROUNDS, these routines honor the current 198187808Sdasrounding mode. 199112158Sdas 200112158SdasFor an example of more general conversions based on dtoa(), see 201112158Sdasnetlib's "printf.c from ampl/solvers". 202112158Sdas 203112158SdasFor double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic 204112158Sdasof precision max(126, #bits(input)) bits, where #bits(input) is the 205112158Sdasnumber of mantissa bits needed to represent the sum of the two double 206112158Sdasvalues in the input. 207112158Sdas 208112158SdasThe makefile creates a library, gdtoa.a. To use the helper 209112158Sdasroutines, a program only needs to include gdtoa.h. All the 210112158Sdassource files for gdtoa.a include a more extensive gdtoaimp.h; 211112158Sdasamong other things, gdtoaimp.h has #defines that make "internal" 212112158Sdasnames end in _D2A. To make a "system" library, one could modify 213112158Sdasthese #defines to make the names start with __. 214112158Sdas 215112158SdasVarious comments about possible #defines appear in gdtoaimp.h, 216112158Sdasbut for most purposes, arith.h should set suitable #defines. 217112158Sdas 218112158SdasSystems with preemptive scheduling of multiple threads require some 219112158Sdasmanual intervention. On such systems, it's necessary to compile 220112158Sdasdmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined, 221112158Sdasand to provide (or suitably #define) two locks, acquired by 222112158SdasACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1. 223112158Sdas(The second lock, accessed in pow5mult, ensures lazy evaluation of 224112158Sdasonly one copy of high powers of 5; omitting this lock would introduce 225112158Sdasa small probability of wasting memory, but would otherwise be harmless.) 226112158SdasRoutines that call dtoa or gdtoa directly must also invoke freedtoa(s) 227112158Sdasto free the value s returned by dtoa or gdtoa. It's OK to do so whether 228112158Sdasor not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines 229112158Sdaslisted above all do this indirectly (in gfmt_D2A(), which they all call). 230112158Sdas 231112158SdasBy default, there is a private pool of memory of length 2000 bytes 232112158Sdasfor intermediate quantities, and MALLOC (see gdtoaimp.h) is called only 233112158Sdasif the private pool does not suffice. 2000 is large enough that MALLOC 234112158Sdasis called only under very unusual circumstances (decimal -> binary 235112158Sdasconversion of very long strings) for conversions to and from double 236165743Sdasprecision. For systems with preemptively scheduled multiple threads 237112158Sdasor for conversions to extended or quad, it may be appropriate to 238112158Sdas#define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000. 239112158SdasFor extended and quad precisions, -DPRIVATE_MEM=20000 is probably 240112158Sdasplenty even for many digits at the ends of the exponent range. 241112158SdasUse of the private pool avoids some overhead. 242112158Sdas 243112158SdasDirectory test provides some test routines. See its README. 244112158SdasI've also tested this stuff (except double double conversions) 245112158Sdaswith Vern Paxson's testbase program: see 246112158Sdas 247112158Sdas V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal 248112158Sdas Conversion", manuscript, May 1991, 249112158Sdas ftp://ftp.ee.lbl.gov/testbase-report.ps.Z . 250112158Sdas 251112158Sdas(The same ftp directory has source for testbase.) 252112158Sdas 253112158SdasSome system-dependent additions to CFLAGS in the makefile: 254112158Sdas 255112158Sdas HU-UX: -Aa -Ae 256112158Sdas OSF (DEC Unix): -ieee_with_no_inexact 257112158Sdas SunOS 4.1x: -DKR_headers -DBad_float_h 258112158Sdas 259112158SdasIf you want to put this stuff into a shared library and your 260112158Sdasoperating system requires export lists for shared libraries, 261112158Sdasthe following would be an appropriate export list: 262112158Sdas 263112158Sdas dtoa 264112158Sdas freedtoa 265112158Sdas g_Qfmt 266112158Sdas g_ddfmt 267112158Sdas g_dfmt 268112158Sdas g_ffmt 269112158Sdas g_xLfmt 270112158Sdas g_xfmt 271112158Sdas gdtoa 272112158Sdas strtoIQ 273112158Sdas strtoId 274112158Sdas strtoIdd 275112158Sdas strtoIf 276112158Sdas strtoIx 277112158Sdas strtoIxL 278112158Sdas strtod 279112158Sdas strtodI 280112158Sdas strtodg 281112158Sdas strtof 282112158Sdas strtopQ 283112158Sdas strtopd 284112158Sdas strtopdd 285112158Sdas strtopf 286112158Sdas strtopx 287112158Sdas strtopxL 288112158Sdas strtorQ 289112158Sdas strtord 290112158Sdas strtordd 291112158Sdas strtorf 292112158Sdas strtorx 293112158Sdas strtorxL 294112158Sdas 295112158SdasWhen time permits, I (dmg) hope to write in more detail about the 296112158Sdaspresent conversion routines; for now, this README file must suffice. 297112158SdasMeanwhile, if you wish to write helper functions for other kinds of 298112158SdasIEEE-like arithmetic, some explanation of struct FPI and the bits 299112158Sdasarray may be helpful. Both gdtoa and strtodg operate on a bits array 300112158Sdasdescribed by FPI *fpi. The bits array is of type ULong, a 32-bit 301112158Sdasunsigned integer type. Floating-point numbers have fpi->nbits bits, 302112158Sdaswith the least significant 32 bits in bits[0], the next 32 bits in 303112158Sdasbits[1], etc. These numbers are regarded as integers multiplied by 304112158Sdas2^e (i.e., 2 to the power of the exponent e), where e is the second 305112158Sdasargument (be) to gdtoa and is stored in *exp by strtodg. The minimum 306112158Sdasand maximum exponent values fpi->emin and fpi->emax for normalized 307112158Sdasfloating-point numbers reflect this arrangement. For example, the 308112158SdasP754 standard for binary IEEE arithmetic specifies doubles as having 309112158Sdas53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023), 310112158Sdaswith 52 bits (the x's) and the biased exponent b represented explicitly; 311112158Sdasb is an unsigned integer in the range 1 <= b <= 2046 for normalized 312112158Sdasfinite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs. 313112158SdasTo turn an IEEE double into the representation used by strtodg and gdtoa, 314112158Sdaswe multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the 315112158Sdasexponent e = (b-1023) by 52: 316112158Sdas 317112158Sdas fpi->emin = 1 - 1023 - 52 318112158Sdas fpi->emax = 1046 - 1023 - 52 319112158Sdas 320112158SdasIn various wrappers for IEEE double, we actually write -53 + 1 rather 321112158Sdasthan -52, to emphasize that there are 53 bits including one implicit bit. 322112158SdasField fpi->rounding indicates the desired rounding direction, with 323112158Sdaspossible values 324112158Sdas FPI_Round_zero = toward 0, 325112158Sdas FPI_Round_near = unbiased rounding -- the IEEE default, 326112158Sdas FPI_Round_up = toward +Infinity, and 327112158Sdas FPI_Round_down = toward -Infinity 328112158Sdasgiven in gdtoa.h. 329112158Sdas 330112158SdasField fpi->sudden_underflow indicates whether strtodg should return 331112158Sdasdenormals or flush them to zero. Normal floating-point numbers have 332112158Sdasbit fpi->nbits in the bits array on. Denormals have it off, with 333112158Sdasexponent = fpi->emin. Strtodg provides distinct return values for normals 334112158Sdasand denormals; see gdtoa.h. 335112158Sdas 336112415SdasCompiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes 337112415Sdasthe decimal-point character to be taken from the current locale; otherwise 338112415Sdasit is '.'. 339112415Sdas 340182709SdasSource files dtoa.c and strtod.c in this directory are derived from 341182709Sdasnetlib's "dtoa.c from fp" and are meant to function equivalently. 342182709SdasWhen compiled with Honor_FLT_ROUNDS #defined (on systems that provide 343182709SdasFLT_ROUNDS and fegetround() as specified in the C99 standard), they 344182709Sdashonor the current rounding mode. Because FLT_ROUNDS is buggy on some 345182709Sdas(Linux) systems -- not reflecting calls on fesetround(), as the C99 346182709Sdasstandard says it should -- when Honor_FLT_ROUNDS is #defined, the 347182709Sdascurrent rounding mode is obtained from fegetround() rather than from 348182709SdasFLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined. 349182709Sdas 350187808SdasCompile with -DUSE_LOCALE to use the current locale; otherwise 351187808Sdasdecimal points are assumed to be '.'. With -DUSE_LOCALE, unless 352187808Sdasyou also compile with -DNO_LOCALE_CACHE, the details about the 353187808Sdascurrent "decimal point" character string are cached and assumed not 354187808Sdasto change during the program's execution. 355187808Sdas 356219557SdasOn machines with a 64-bit long double and perhaps a 113-bit "quad" 357219557Sdastype, you can invoke "make Printf" to add Printf (and variants, such 358219557Sdasas Fprintf) to gdtoa.a. These are analogs, declared in stdio1.h, of 359219557Sdasprintf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long 360219557Sdasdouble and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad 361219557Sdasprecision printing. 362219557Sdas 363165743SdasPlease send comments to David M. Gay (dmg at acm dot org, with " at " 364165743Sdaschanged at "@" and " dot " changed to "."). 365