| This directory contains source for a library of binary -> decimal |
| and decimal -> binary conversion routines, for single-, double-, |
| and extended-precision IEEE binary floating-point arithmetic, and |
| other IEEE-like binary floating-point, including "double double", |
| as in |
| |
| T. J. Dekker, "A Floating-Point Technique for Extending the |
| Available Precision", Numer. Math. 18 (1971), pp. 224-242 |
| |
| and |
| |
| "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994 |
| |
| The conversion routines use double-precision floating-point arithmetic |
| and, where necessary, high precision integer arithmetic. The routines |
| are generalizations of the strtod and dtoa routines described in |
| |
| David M. Gay, "Correctly Rounded Binary-Decimal and |
| Decimal-Binary Conversions", Numerical Analysis Manuscript |
| No. 90-10, Bell Labs, Murray Hill, 1990; |
| http://cm.bell-labs.com/cm/cs/what/ampl/REFS/rounding.ps.gz |
| |
| (based in part on papers by Clinger and Steele & White: see the |
| references in the above paper). |
| |
| The present conversion routines should be able to use any of IEEE binary, |
| VAX, or IBM-mainframe double-precision arithmetic internally, but I (dmg) |
| have so far only had a chance to test them with IEEE double precision |
| arithmetic. |
| |
| The core conversion routines are strtodg for decimal -> binary conversions |
| and gdtoa for binary -> decimal conversions. These routines operate |
| on arrays of unsigned 32-bit integers of type ULong, a signed 32-bit |
| exponent of type Long, and arithmetic characteristics described in |
| struct FPI; FPI, Long, and ULong are defined in gdtoa.h. File arith.h |
| is supposed to provide #defines that cause gdtoa.h to define its |
| types correctly. File arithchk.c is source for a program that |
| generates a suitable arith.h on all systems where I've been able to |
| test it. |
| |
| The core conversion routines are meant to be called by helper routines |
| that know details of the particular binary arithmetic of interest and |
| convert. The present directory provides helper routines for 5 variants |
| of IEEE binary floating-point arithmetic, each indicated by one or |
| two letters: |
| |
| f IEEE single precision |
| d IEEE double precision |
| x IEEE extended precision, as on Intel 80x87 |
| and software emulations of Motorola 68xxx chips |
| that do not pad the way the 68xxx does, but |
| only store 80 bits |
| xL IEEE extended precision, as on Motorola 68xxx chips |
| Q quad precision, as on Sun Sparc chips |
| dd double double, pairs of IEEE double numbers |
| whose sum is the desired value |
| |
| For decimal -> binary conversions, there are three families of |
| helper routines: one for round-nearest: |
| |
| strtof |
| strtod |
| strtodd |
| strtopd |
| strtopf |
| strtopx |
| strtopxL |
| strtopQ |
| |
| one with rounding direction specified: |
| |
| strtorf |
| strtord |
| strtordd |
| strtorx |
| strtorxL |
| strtorQ |
| |
| and one for computing an interval (at most one bit wide) that contains |
| the decimal number: |
| |
| strtoIf |
| strtoId |
| strtoIdd |
| strtoIx |
| strtoIxL |
| strtoIQ |
| |
| The latter call strtoIg, which makes one call on strtodg and adjusts |
| the result to provide the desired interval. On systems where native |
| arithmetic can easily make one-ulp adjustments on values in the |
| desired floating-point format, it might be more efficient to use the |
| native arithmetic. Routine strtodI is a variant of strtoId that |
| illustrates one way to do this for IEEE binary double-precision |
| arithmetic -- but whether this is more efficient remains to be seen. |
| |
| Functions strtod and strtof have "natural" return types, float and |
| double -- strtod is specified by the C standard, and strtof appears |
| in the stdlib.h of some systems, such as (at least some) Linux systems. |
| The other functions write their results to their final argument(s): |
| to the final two argument for the strtoI... (interval) functions, |
| and to the final argument for the others (strtop... and strtor...). |
| Where possible, these arguments have "natural" return types (double* |
| or float*), to permit at least some type checking. In reality, they |
| are viewed as arrays of ULong (or, for the "x" functions, UShort) |
| values. On systems where long double is the appropriate type, one can |
| pass long double* final argument(s) to these routines. The int value |
| that these routines return is the return value from the call they make |
| on strtodg; see the enum of possible return values in gdtoa.h. |
| |
| Source files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c |
| should use true IEEE double arithmetic (not, e.g., double extended), |
| at least for storing (and viewing the bits of) the variables declared |
| "double" within them. |
| |
| One detail indicated in struct FPI is whether the target binary |
| arithmetic departs from the IEEE standard by flushing denormalized |
| numbers to 0. On systems that do this, the helper routines for |
| conversion to double-double format (when compiled with |
| Sudden_Underflow #defined) penalize the bottom of the exponent |
| range so that they return a nonzero result only when the least |
| significant bit of the less significant member of the pair of |
| double values returned can be expressed as a normalized double |
| value. An alternative would be to drop to 53-bit precision near |
| the bottom of the exponent range. To get correct rounding, this |
| would (in general) require two calls on strtodg (one specifying |
| 126-bit arithmetic, then, if necessary, one specifying 53-bit |
| arithmetic). |
| |
| By default, the core routine strtodg and strtod set errno to ERANGE |
| if the result overflows to +Infinity or underflows to 0. Compile |
| these routines with NO_ERRNO #defined to inhibit errno assignments. |
| |
| Routine strtod is based on netlib's "dtoa.c from fp", and |
| (f = strtod(s,se)) is more efficient for some conversions than, say, |
| strtord(s,se,1,&f). Parts of strtod require true IEEE double |
| arithmetic with the default rounding mode (round-to-nearest) and, on |
| systems with IEEE extended-precision registers, double-precision |
| (53-bit) rounding precision. If the machine uses (the equivalent of) |
| Intel 80x87 arithmetic, the call |
| _control87(PC_53, MCW_PC); |
| does this with many compilers. Whether this or another call is |
| appropriate depends on the compiler; for this to work, it may be |
| necessary to #include "float.h" or another system-dependent header |
| file. |
| |
| The values returned for NaNs may be signaling NaNs on some systems, |
| since the rules for distinguishing signaling from quiet NaNs are |
| system-dependent. You can easily fix this by suitably modifying the |
| ULto* routines in strtor*.c. |
| |
| C99's hexadecimal floating-point constants are recognized by the |
| strto* routines (but this feature has not yet been heavily tested). |
| Compiling with NO_HEX_FP #defined disables this feature. |
| |
| The strto* routines do not (yet) recognize C99's NaN(...) syntax; the |
| strto* routines simply regard '(' as the first unprocessed input |
| character. |
| |
| For binary -> decimal conversions, I've provided just one family |
| of helper routines: |
| |
| g_ffmt |
| g_dfmt |
| g_ddfmt |
| g_xfmt |
| g_xLfmt |
| g_Qfmt |
| |
| which do a "%g" style conversion either to a specified number of decimal |
| places (if their ndig argument is positive), or to the shortest |
| decimal string that rounds to the given binary floating-point value |
| (if ndig <= 0). They write into a buffer supplied as an argument |
| and return either a pointer to the end of the string (a null character) |
| in the buffer, if the buffer was long enough, or 0. Other forms of |
| conversion are easily done with the help of gdtoa(), such as %e or %f |
| style and conversions with direction of rounding specified (so that, if |
| desired, the decimal value is either >= or <= the binary value). |
| |
| For an example of more general conversions based on dtoa(), see |
| netlib's "printf.c from ampl/solvers". |
| |
| For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic |
| of precision max(126, #bits(input)) bits, where #bits(input) is the |
| number of mantissa bits needed to represent the sum of the two double |
| values in the input. |
| |
| The makefile creates a library, gdtoa.a. To use the helper |
| routines, a program only needs to include gdtoa.h. All the |
| source files for gdtoa.a include a more extensive gdtoaimp.h; |
| among other things, gdtoaimp.h has #defines that make "internal" |
| names end in _D2A. To make a "system" library, one could modify |
| these #defines to make the names start with __. |
| |
| Various comments about possible #defines appear in gdtoaimp.h, |
| but for most purposes, arith.h should set suitable #defines. |
| |
| Systems with preemptive scheduling of multiple threads require some |
| manual intervention. On such systems, it's necessary to compile |
| dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined, |
| and to provide (or suitably #define) two locks, acquired by |
| ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1. |
| (The second lock, accessed in pow5mult, ensures lazy evaluation of |
| only one copy of high powers of 5; omitting this lock would introduce |
| a small probability of wasting memory, but would otherwise be harmless.) |
| Routines that call dtoa or gdtoa directly must also invoke freedtoa(s) |
| to free the value s returned by dtoa or gdtoa. It's OK to do so whether |
| or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines |
| listed above all do this indirectly (in gfmt_D2A(), which they all call). |
| |
| By default, there is a private pool of memory of length 2000 bytes |
| for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only |
| if the private pool does not suffice. 2000 is large enough that MALLOC |
| is called only under very unusual circumstances (decimal -> binary |
| conversion of very long strings) for conversions to and from double |
| precision. For systems with preemptivaly scheduled multiple threads |
| or for conversions to extended or quad, it may be appropriate to |
| #define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000. |
| For extended and quad precisions, -DPRIVATE_MEM=20000 is probably |
| plenty even for many digits at the ends of the exponent range. |
| Use of the private pool avoids some overhead. |
| |
| Directory test provides some test routines. See its README. |
| I've also tested this stuff (except double double conversions) |
| with Vern Paxson's testbase program: see |
| |
| V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal |
| Conversion", manuscript, May 1991, |
| ftp://ftp.ee.lbl.gov/testbase-report.ps.Z . |
| |
| (The same ftp directory has source for testbase.) |
| |
| Some system-dependent additions to CFLAGS in the makefile: |
| |
| HU-UX: -Aa -Ae |
| OSF (DEC Unix): -ieee_with_no_inexact |
| SunOS 4.1x: -DKR_headers -DBad_float_h |
| |
| If you want to put this stuff into a shared library and your |
| operating system requires export lists for shared libraries, |
| the following would be an appropriate export list: |
| |
| dtoa |
| freedtoa |
| g_Qfmt |
| g_ddfmt |
| g_dfmt |
| g_ffmt |
| g_xLfmt |
| g_xfmt |
| gdtoa |
| strtoIQ |
| strtoId |
| strtoIdd |
| strtoIf |
| strtoIx |
| strtoIxL |
| strtod |
| strtodI |
| strtodg |
| strtof |
| strtopQ |
| strtopd |
| strtopdd |
| strtopf |
| strtopx |
| strtopxL |
| strtorQ |
| strtord |
| strtordd |
| strtorf |
| strtorx |
| strtorxL |
| |
| When time permits, I (dmg) hope to write in more detail about the |
| present conversion routines; for now, this README file must suffice. |
| Meanwhile, if you wish to write helper functions for other kinds of |
| IEEE-like arithmetic, some explanation of struct FPI and the bits |
| array may be helpful. Both gdtoa and strtodg operate on a bits array |
| described by FPI *fpi. The bits array is of type ULong, a 32-bit |
| unsigned integer type. Floating-point numbers have fpi->nbits bits, |
| with the least significant 32 bits in bits[0], the next 32 bits in |
| bits[1], etc. These numbers are regarded as integers multiplied by |
| 2^e (i.e., 2 to the power of the exponent e), where e is the second |
| argument (be) to gdtoa and is stored in *exp by strtodg. The minimum |
| and maximum exponent values fpi->emin and fpi->emax for normalized |
| floating-point numbers reflect this arrangement. For example, the |
| P754 standard for binary IEEE arithmetic specifies doubles as having |
| 53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023), |
| with 52 bits (the x's) and the biased exponent b represented explicitly; |
| b is an unsigned integer in the range 1 <= b <= 2046 for normalized |
| finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs. |
| To turn an IEEE double into the representation used by strtodg and gdtoa, |
| we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the |
| exponent e = (b-1023) by 52: |
| |
| fpi->emin = 1 - 1023 - 52 |
| fpi->emax = 1046 - 1023 - 52 |
| |
| In various wrappers for IEEE double, we actually write -53 + 1 rather |
| than -52, to emphasize that there are 53 bits including one implicit bit. |
| Field fpi->rounding indicates the desired rounding direction, with |
| possible values |
| FPI_Round_zero = toward 0, |
| FPI_Round_near = unbiased rounding -- the IEEE default, |
| FPI_Round_up = toward +Infinity, and |
| FPI_Round_down = toward -Infinity |
| given in gdtoa.h. |
| |
| Field fpi->sudden_underflow indicates whether strtodg should return |
| denormals or flush them to zero. Normal floating-point numbers have |
| bit fpi->nbits in the bits array on. Denormals have it off, with |
| exponent = fpi->emin. Strtodg provides distinct return values for normals |
| and denormals; see gdtoa.h. |
| |
| Please send comments to |
| |
| David M. Gay |
| Bell Labs, Room 2C-463 |
| 600 Mountain Avenue |
| Murray Hill, NJ 07974-0636, U.S.A. |
| dmg@research.bell-labs.com |