NameDateSize

..27-Jun-202090

arithchk.c06-Nov-20124 KiB

changes06-Nov-201230.3 KiB

dmisc.c06-Nov-20124.6 KiB

dtoa.c06-Nov-201217.1 KiB

g__fmt.c06-Nov-20123.3 KiB

g_ddfmt.c06-Nov-20124.1 KiB

g_dfmt.c06-Nov-20122.6 KiB

g_ffmt.c06-Nov-20122.5 KiB

g_Qfmt.c06-Nov-20122.9 KiB

g_xfmt.c06-Nov-20122.8 KiB

g_xLfmt.c06-Nov-20122.7 KiB

gdtoa.c06-Nov-201216.9 KiB

gdtoa.h06-Nov-20124.8 KiB

gdtoa_fltrnds.h06-Nov-2012421

gdtoaimp.h06-Nov-201222.7 KiB

gethex.c06-Nov-20126.9 KiB

gmisc.c06-Nov-20122 KiB

hd_init.c06-Nov-20121.8 KiB

hexnan.c06-Nov-20123.4 KiB

makefile16-Dec-20123.1 KiB

misc.c06-Nov-201214.1 KiB

printf.c06-Nov-2012158

printf.c006-Nov-201229 KiB

qnan.c06-Nov-20123.3 KiB

README06-Nov-201215.4 KiB

smisc.c06-Nov-20123.6 KiB

stdio1.h06-Nov-20122.8 KiB

strtod.c06-Nov-201222.4 KiB

strtodg.c06-Nov-201220.7 KiB

strtodI.c06-Nov-20123.9 KiB

strtodnrp.c06-Nov-20122.5 KiB

strtof.c06-Nov-20122.4 KiB

strtoId.c06-Nov-20121.9 KiB

strtoIdd.c06-Nov-20122.1 KiB

strtoIf.c06-Nov-20121.8 KiB

strtoIg.c06-Nov-20123.5 KiB

strtoIQ.c06-Nov-20121.9 KiB

strtoIx.c06-Nov-20121.9 KiB

strtoIxL.c06-Nov-20121.9 KiB

strtopd.c06-Nov-20121.7 KiB

strtopdd.c06-Nov-20124.5 KiB

strtopf.c06-Nov-20122.1 KiB

strtopQ.c06-Nov-20122.6 KiB

strtopx.c06-Nov-20122.6 KiB

strtopxL.c06-Nov-20122.4 KiB

strtord.c06-Nov-20122.5 KiB

strtordd.c06-Nov-20124.9 KiB

strtorf.c06-Nov-20122.3 KiB

strtorQ.c06-Nov-20123 KiB

strtorx.c06-Nov-20123.1 KiB

strtorxL.c06-Nov-20122.7 KiB

sum.c06-Nov-20122.4 KiB

ulp.c06-Nov-20121.8 KiB

README

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
155The values that the strto* routines return for NaNs are determined by
156gd_qnan.h, which the makefile generates by running the program whose
157source is qnan.c.  Note that the rules for distinguishing signaling
158from quiet NaNs are system-dependent.  For cross-compilation, you need
159to determine arith.h and gd_qnan.h suitably, e.g., using the
160arithmetic of the target machine.
161
162C99's hexadecimal floating-point constants are recognized by the
163strto* routines (but this feature has not yet been heavily tested).
164Compiling with NO_HEX_FP #defined disables this feature.
165
166When compiled with -DINFNAN_CHECK, the strto* routines recognize C99's
167NaN and Infinity syntax.  Moreover, unless No_Hex_NaN is #defined, the
168strto* routines also recognize C99's NaN(...) syntax: they accept
169(case insensitively) strings of the form NaN(x), where x is a string
170of hexadecimal digits and spaces; if there is only one string of
171hexadecimal digits, it is taken for the fraction bits of the resulting
172NaN; if there are two or more strings of hexadecimal digits, each
173string is assigned to the next available sequence of 32-bit words of
174fractions bits (starting with the most significant), right-aligned in
175each sequence.
176
177For binary -> decimal conversions, I've provided just one family
178of helper routines:
179
180	g_ffmt
181	g_dfmt
182	g_ddfmt
183	g_xfmt
184	g_xLfmt
185	g_Qfmt
186
187which do a "%g" style conversion either to a specified number of decimal
188places (if their ndig argument is positive), or to the shortest
189decimal string that rounds to the given binary floating-point value
190(if ndig <= 0).  They write into a buffer supplied as an argument
191and return either a pointer to the end of the string (a null character)
192in the buffer, if the buffer was long enough, or 0.  Other forms of
193conversion are easily done with the help of gdtoa(), such as %e or %f
194style and conversions with direction of rounding specified (so that, if
195desired, the decimal value is either >= or <= the binary value).
196On IEEE-arithmetic systems that provide the C99 fegetround() function,
197if compiled with -DHonor_FLT_ROUNDS, these routines honor the current
198rounding mode.
199
200For an example of more general conversions based on dtoa(), see
201netlib's "printf.c from ampl/solvers".
202
203For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic
204of precision max(126, #bits(input)) bits, where #bits(input) is the
205number of mantissa bits needed to represent the sum of the two double
206values in the input.
207
208The makefile creates a library, gdtoa.a.  To use the helper
209routines, a program only needs to include gdtoa.h.  All the
210source files for gdtoa.a include a more extensive gdtoaimp.h;
211among other things, gdtoaimp.h has #defines that make "internal"
212names end in _D2A.  To make a "system" library, one could modify
213these #defines to make the names start with __.
214
215Various comments about possible #defines appear in gdtoaimp.h,
216but for most purposes, arith.h should set suitable #defines.
217
218Systems with preemptive scheduling of multiple threads require some
219manual intervention.  On such systems, it's necessary to compile
220dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined,
221and to provide (or suitably #define) two locks, acquired by
222ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1.
223(The second lock, accessed in pow5mult, ensures lazy evaluation of
224only one copy of high powers of 5; omitting this lock would introduce
225a small probability of wasting memory, but would otherwise be harmless.)
226Routines that call dtoa or gdtoa directly must also invoke freedtoa(s)
227to free the value s returned by dtoa or gdtoa.  It's OK to do so whether
228or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines
229listed above all do this indirectly (in gfmt_D2A(), which they all call).
230
231By default, there is a private pool of memory of length 2000 bytes
232for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only
233if the private pool does not suffice.   2000 is large enough that MALLOC
234is called only under very unusual circumstances (decimal -> binary
235conversion of very long strings) for conversions to and from double
236precision.  For systems with preemptively scheduled multiple threads
237or for conversions to extended or quad, it may be appropriate to
238#define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000.
239For extended and quad precisions, -DPRIVATE_MEM=20000 is probably
240plenty even for many digits at the ends of the exponent range.
241Use of the private pool avoids some overhead.
242
243Directory test provides some test routines.  See its README.
244I've also tested this stuff (except double double conversions)
245with Vern Paxson's testbase program: see
246
247	V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal
248	Conversion", manuscript, May 1991,
249	ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
250
251(The same ftp directory has source for testbase.)
252
253Some system-dependent additions to CFLAGS in the makefile:
254
255	HU-UX: -Aa -Ae
256	OSF (DEC Unix): -ieee_with_no_inexact
257	SunOS 4.1x: -DKR_headers -DBad_float_h
258
259If you want to put this stuff into a shared library and your
260operating system requires export lists for shared libraries,
261the following would be an appropriate export list:
262
263	dtoa
264	freedtoa
265	g_Qfmt
266	g_ddfmt
267	g_dfmt
268	g_ffmt
269	g_xLfmt
270	g_xfmt
271	gdtoa
272	strtoIQ
273	strtoId
274	strtoIdd
275	strtoIf
276	strtoIx
277	strtoIxL
278	strtod
279	strtodI
280	strtodg
281	strtof
282	strtopQ
283	strtopd
284	strtopdd
285	strtopf
286	strtopx
287	strtopxL
288	strtorQ
289	strtord
290	strtordd
291	strtorf
292	strtorx
293	strtorxL
294
295When time permits, I (dmg) hope to write in more detail about the
296present conversion routines; for now, this README file must suffice.
297Meanwhile, if you wish to write helper functions for other kinds of
298IEEE-like arithmetic, some explanation of struct FPI and the bits
299array may be helpful.  Both gdtoa and strtodg operate on a bits array
300described by FPI *fpi.  The bits array is of type ULong, a 32-bit
301unsigned integer type.  Floating-point numbers have fpi->nbits bits,
302with the least significant 32 bits in bits[0], the next 32 bits in
303bits[1], etc.  These numbers are regarded as integers multiplied by
3042^e (i.e., 2 to the power of the exponent e), where e is the second
305argument (be) to gdtoa and is stored in *exp by strtodg.  The minimum
306and maximum exponent values fpi->emin and fpi->emax for normalized
307floating-point numbers reflect this arrangement.  For example, the
308P754 standard for binary IEEE arithmetic specifies doubles as having
30953 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023),
310with 52 bits (the x's) and the biased exponent b represented explicitly;
311b is an unsigned integer in the range 1 <= b <= 2046 for normalized
312finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs.
313To turn an IEEE double into the representation used by strtodg and gdtoa,
314we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the
315exponent e = (b-1023) by 52:
316
317	fpi->emin = 1 - 1023 - 52
318	fpi->emax = 1046 - 1023 - 52
319
320In various wrappers for IEEE double, we actually write -53 + 1 rather
321than -52, to emphasize that there are 53 bits including one implicit bit.
322Field fpi->rounding indicates the desired rounding direction, with
323possible values
324	FPI_Round_zero = toward 0,
325	FPI_Round_near = unbiased rounding -- the IEEE default,
326	FPI_Round_up = toward +Infinity, and
327	FPI_Round_down = toward -Infinity
328given in gdtoa.h.
329
330Field fpi->sudden_underflow indicates whether strtodg should return
331denormals or flush them to zero.  Normal floating-point numbers have
332bit fpi->nbits in the bits array on.  Denormals have it off, with
333exponent = fpi->emin.  Strtodg provides distinct return values for normals
334and denormals; see gdtoa.h.
335
336Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
337the decimal-point character to be taken from the current locale; otherwise
338it is '.'.
339
340Source files dtoa.c and strtod.c in this directory are derived from
341netlib's "dtoa.c from fp" and are meant to function equivalently.
342When compiled with Honor_FLT_ROUNDS #defined (on systems that provide
343FLT_ROUNDS and fegetround() as specified in the C99 standard), they
344honor the current rounding mode.  Because FLT_ROUNDS is buggy on some
345(Linux) systems -- not reflecting calls on fesetround(), as the C99
346standard says it should -- when Honor_FLT_ROUNDS is #defined, the
347current rounding mode is obtained from fegetround() rather than from
348FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined.
349
350Compile with -DUSE_LOCALE to use the current locale; otherwise
351decimal points are assumed to be '.'.  With -DUSE_LOCALE, unless
352you also compile with -DNO_LOCALE_CACHE, the details about the
353current "decimal point" character string are cached and assumed not
354to change during the program's execution.
355
356On machines with a 64-bit long double and perhaps a 113-bit "quad"
357type, you can invoke "make Printf" to add Printf (and variants, such
358as Fprintf) to gdtoa.a.  These are analogs, declared in stdio1.h, of
359printf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long
360double and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad
361precision printing.
362
363Please send comments to	David M. Gay (dmg at acm dot org, with " at "
364changed at "@" and " dot " changed to ".").
365