1 The content below is the README file of the gdtoa distribution, available
2 from http://www.netlib.org/fp/
4 ----------------------------------------------------
6 This directory contains source for a library of binary -> decimal
7 and decimal -> binary conversion routines, for single-, double-,
8 and extended-precision IEEE binary floating-point arithmetic, and
9 other IEEE-like binary floating-point, including "double double",
12 T. J. Dekker, "A Floating-Point Technique for Extending the
13 Available Precision", Numer. Math. 18 (1971), pp. 224-242
17 "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994
19 The conversion routines use double-precision floating-point arithmetic
20 and, where necessary, high precision integer arithmetic. The routines
21 are generalizations of the strtod and dtoa routines described in
23 David M. Gay, "Correctly Rounded Binary-Decimal and
24 Decimal-Binary Conversions", Numerical Analysis Manuscript
25 No. 90-10, Bell Labs, Murray Hill, 1990;
26 http://cm.bell-labs.com/cm/cs/what/ampl/REFS/rounding.ps.gz
28 (based in part on papers by Clinger and Steele & White: see the
29 references in the above paper).
31 The present conversion routines should be able to use any of IEEE binary,
32 VAX, or IBM-mainframe double-precision arithmetic internally, but I (dmg)
33 have so far only had a chance to test them with IEEE double precision
36 The core conversion routines are strtodg for decimal -> binary conversions
37 and gdtoa for binary -> decimal conversions. These routines operate
38 on arrays of unsigned 32-bit integers of type ULong, a signed 32-bit
39 exponent of type Long, and arithmetic characteristics described in
40 struct FPI; FPI, Long, and ULong are defined in gdtoa.h. File arith.h
41 is supposed to provide #defines that cause gdtoa.h to define its
42 types correctly. File arithchk.c is source for a program that
43 generates a suitable arith.h on all systems where I've been able to
46 The core conversion routines are meant to be called by helper routines
47 that know details of the particular binary arithmetic of interest and
48 convert. The present directory provides helper routines for 5 variants
49 of IEEE binary floating-point arithmetic, each indicated by one or
52 f IEEE single precision
53 d IEEE double precision
54 x IEEE extended precision, as on Intel 80x87
55 and software emulations of Motorola 68xxx chips
56 that do not pad the way the 68xxx does, but
58 xL IEEE extended precision, as on Motorola 68xxx chips
59 Q quad precision, as on Sun Sparc chips
60 dd double double, pairs of IEEE double numbers
61 whose sum is the desired value
63 For decimal -> binary conversions, there are three families of
64 helper routines: one for round-nearest (or the current rounding
65 mode on IEEE-arithmetic systems that provide the C99 fegetround()
66 function, if compiled with -DHonor_FLT_ROUNDS):
77 one with rounding direction specified:
86 and one for computing an interval (at most one bit wide) that contains
96 The latter call strtoIg, which makes one call on strtodg and adjusts
97 the result to provide the desired interval. On systems where native
98 arithmetic can easily make one-ulp adjustments on values in the
99 desired floating-point format, it might be more efficient to use the
100 native arithmetic. Routine strtodI is a variant of strtoId that
101 illustrates one way to do this for IEEE binary double-precision
102 arithmetic -- but whether this is more efficient remains to be seen.
104 Functions strtod and strtof have "natural" return types, float and
105 double -- strtod is specified by the C standard, and strtof appears
106 in the stdlib.h of some systems, such as (at least some) Linux systems.
107 The other functions write their results to their final argument(s):
108 to the final two argument for the strtoI... (interval) functions,
109 and to the final argument for the others (strtop... and strtor...).
110 Where possible, these arguments have "natural" return types (double*
111 or float*), to permit at least some type checking. In reality, they
112 are viewed as arrays of ULong (or, for the "x" functions, UShort)
113 values. On systems where long double is the appropriate type, one can
114 pass long double* final argument(s) to these routines. The int value
115 that these routines return is the return value from the call they make
116 on strtodg; see the enum of possible return values in gdtoa.h.
118 Source files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c
119 should use true IEEE double arithmetic (not, e.g., double extended),
120 at least for storing (and viewing the bits of) the variables declared
121 "double" within them.
123 One detail indicated in struct FPI is whether the target binary
124 arithmetic departs from the IEEE standard by flushing denormalized
125 numbers to 0. On systems that do this, the helper routines for
126 conversion to double-double format (when compiled with
127 Sudden_Underflow #defined) penalize the bottom of the exponent
128 range so that they return a nonzero result only when the least
129 significant bit of the less significant member of the pair of
130 double values returned can be expressed as a normalized double
131 value. An alternative would be to drop to 53-bit precision near
132 the bottom of the exponent range. To get correct rounding, this
133 would (in general) require two calls on strtodg (one specifying
134 126-bit arithmetic, then, if necessary, one specifying 53-bit
137 By default, the core routine strtodg and strtod set errno to ERANGE
138 if the result overflows to +Infinity or underflows to 0. Compile
139 these routines with NO_ERRNO #defined to inhibit errno assignments.
141 Routine strtod is based on netlib's "dtoa.c from fp", and
142 (f = strtod(s,se)) is more efficient for some conversions than, say,
143 strtord(s,se,1,&f). Parts of strtod require true IEEE double
144 arithmetic with the default rounding mode (round-to-nearest) and, on
145 systems with IEEE extended-precision registers, double-precision
146 (53-bit) rounding precision. If the machine uses (the equivalent of)
147 Intel 80x87 arithmetic, the call
148 _control87(PC_53, MCW_PC);
149 does this with many compilers. Whether this or another call is
150 appropriate depends on the compiler; for this to work, it may be
151 necessary to #include "float.h" or another system-dependent header
154 Source file strtodnrp.c gives a strtod that does not require 53-bit
155 rounding precision on systems (such as Intel IA32 systems) that may
156 suffer double rounding due to use of extended-precision registers.
157 For some conversions this variant of strtod is less efficient than the
158 one in strtod.c when the latter is run with 53-bit rounding precision.
160 The values that the strto* routines return for NaNs are determined by
161 gd_qnan.h, which the makefile generates by running the program whose
162 source is qnan.c. Note that the rules for distinguishing signaling
163 from quiet NaNs are system-dependent. For cross-compilation, you need
164 to determine arith.h and gd_qnan.h suitably, e.g., using the
165 arithmetic of the target machine.
167 C99's hexadecimal floating-point constants are recognized by the
168 strto* routines (but this feature has not yet been heavily tested).
169 Compiling with NO_HEX_FP #defined disables this feature.
171 When compiled with -DINFNAN_CHECK, the strto* routines recognize C99's
172 NaN and Infinity syntax. Moreover, unless No_Hex_NaN is #defined, the
173 strto* routines also recognize C99's NaN(...) syntax: they accept
174 (case insensitively) strings of the form NaN(x), where x is a string
175 of hexadecimal digits and spaces; if there is only one string of
176 hexadecimal digits, it is taken for the fraction bits of the resulting
177 NaN; if there are two or more strings of hexadecimal digits, each
178 string is assigned to the next available sequence of 32-bit words of
179 fractions bits (starting with the most significant), right-aligned in
182 For binary -> decimal conversions, I've provided just one family
192 which do a "%g" style conversion either to a specified number of decimal
193 places (if their ndig argument is positive), or to the shortest
194 decimal string that rounds to the given binary floating-point value
195 (if ndig <= 0). They write into a buffer supplied as an argument
196 and return either a pointer to the end of the string (a null character)
197 in the buffer, if the buffer was long enough, or 0. Other forms of
198 conversion are easily done with the help of gdtoa(), such as %e or %f
199 style and conversions with direction of rounding specified (so that, if
200 desired, the decimal value is either >= or <= the binary value).
201 On IEEE-arithmetic systems that provide the C99 fegetround() function,
202 if compiled with -DHonor_FLT_ROUNDS, these routines honor the current
205 For an example of more general conversions based on dtoa(), see
206 netlib's "printf.c from ampl/solvers".
208 For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic
209 of precision max(126, #bits(input)) bits, where #bits(input) is the
210 number of mantissa bits needed to represent the sum of the two double
213 The makefile creates a library, gdtoa.a. To use the helper
214 routines, a program only needs to include gdtoa.h. All the
215 source files for gdtoa.a include a more extensive gdtoaimp.h;
216 among other things, gdtoaimp.h has #defines that make "internal"
217 names end in _D2A. To make a "system" library, one could modify
218 these #defines to make the names start with __.
220 Various comments about possible #defines appear in gdtoaimp.h,
221 but for most purposes, arith.h should set suitable #defines.
223 Systems with preemptive scheduling of multiple threads require some
224 manual intervention. On such systems, it's necessary to compile
225 dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined,
226 and to provide (or suitably #define) two locks, acquired by
227 ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1.
228 (The second lock, accessed in pow5mult, ensures lazy evaluation of
229 only one copy of high powers of 5; omitting this lock would introduce
230 a small probability of wasting memory, but would otherwise be harmless.)
231 Routines that call dtoa or gdtoa directly must also invoke freedtoa(s)
232 to free the value s returned by dtoa or gdtoa. It's OK to do so whether
233 or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines
234 listed above all do this indirectly (in gfmt_D2A(), which they all call).
236 By default, there is a private pool of memory of length 2000 bytes
237 for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only
238 if the private pool does not suffice. 2000 is large enough that MALLOC
239 is called only under very unusual circumstances (decimal -> binary
240 conversion of very long strings) for conversions to and from double
241 precision. For systems with preemptively scheduled multiple threads
242 or for conversions to extended or quad, it may be appropriate to
243 #define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000.
244 For extended and quad precisions, -DPRIVATE_MEM=20000 is probably
245 plenty even for many digits at the ends of the exponent range.
246 Use of the private pool avoids some overhead.
248 Directory test provides some test routines. See its README.
249 I've also tested this stuff (except double double conversions)
250 with Vern Paxson's testbase program: see
252 V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal
253 Conversion", manuscript, May 1991,
254 ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
256 (The same ftp directory has source for testbase.)
258 Some system-dependent additions to CFLAGS in the makefile:
261 OSF (DEC Unix): -ieee_with_no_inexact
262 SunOS 4.1x: -DKR_headers -DBad_float_h
264 If you want to put this stuff into a shared library and your
265 operating system requires export lists for shared libraries,
266 the following would be an appropriate export list:
300 When time permits, I (dmg) hope to write in more detail about the
301 present conversion routines; for now, this README file must suffice.
302 Meanwhile, if you wish to write helper functions for other kinds of
303 IEEE-like arithmetic, some explanation of struct FPI and the bits
304 array may be helpful. Both gdtoa and strtodg operate on a bits array
305 described by FPI *fpi. The bits array is of type ULong, a 32-bit
306 unsigned integer type. Floating-point numbers have fpi->nbits bits,
307 with the least significant 32 bits in bits[0], the next 32 bits in
308 bits[1], etc. These numbers are regarded as integers multiplied by
309 2^e (i.e., 2 to the power of the exponent e), where e is the second
310 argument (be) to gdtoa and is stored in *exp by strtodg. The minimum
311 and maximum exponent values fpi->emin and fpi->emax for normalized
312 floating-point numbers reflect this arrangement. For example, the
313 P754 standard for binary IEEE arithmetic specifies doubles as having
314 53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023),
315 with 52 bits (the x's) and the biased exponent b represented explicitly;
316 b is an unsigned integer in the range 1 <= b <= 2046 for normalized
317 finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs.
318 To turn an IEEE double into the representation used by strtodg and gdtoa,
319 we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the
320 exponent e = (b-1023) by 52:
322 fpi->emin = 1 - 1023 - 52
323 fpi->emax = 1046 - 1023 - 52
325 In various wrappers for IEEE double, we actually write -53 + 1 rather
326 than -52, to emphasize that there are 53 bits including one implicit bit.
327 Field fpi->rounding indicates the desired rounding direction, with
329 FPI_Round_zero = toward 0,
330 FPI_Round_near = unbiased rounding -- the IEEE default,
331 FPI_Round_up = toward +Infinity, and
332 FPI_Round_down = toward -Infinity
335 Field fpi->sudden_underflow indicates whether strtodg should return
336 denormals or flush them to zero. Normal floating-point numbers have
337 bit fpi->nbits in the bits array on. Denormals have it off, with
338 exponent = fpi->emin. Strtodg provides distinct return values for normals
339 and denormals; see gdtoa.h.
341 Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
342 the decimal-point character to be taken from the current locale; otherwise
345 Source files dtoa.c and strtod.c in this directory are derived from
346 netlib's "dtoa.c from fp" and are meant to function equivalently.
347 When compiled with Honor_FLT_ROUNDS #defined (on systems that provide
348 FLT_ROUNDS and fegetround() as specified in the C99 standard), they
349 honor the current rounding mode. Because FLT_ROUNDS is buggy on some
350 (Linux) systems -- not reflecting calls on fesetround(), as the C99
351 standard says it should -- when Honor_FLT_ROUNDS is #defined, the
352 current rounding mode is obtained from fegetround() rather than from
353 FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined.
355 Compile with -DUSE_LOCALE to use the current locale; otherwise
356 decimal points are assumed to be '.'. With -DUSE_LOCALE, unless
357 you also compile with -DNO_LOCALE_CACHE, the details about the
358 current "decimal point" character string are cached and assumed not
359 to change during the program's execution.
361 On machines with a 64-bit long double and perhaps a 113-bit "quad"
362 type, you can invoke "make Printf" to add Printf (and variants, such
363 as Fprintf) to gdtoa.a. These are analogs, declared in stdio1.h, of
364 printf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long
365 double and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad
368 Please send comments to David M. Gay (dmg at acm dot org, with " at "
369 changed at "@" and " dot " changed to ".").
OSDN Git Service
