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