1. Summary
  2. Files
  3. Support
  4. Report Spam
  5. Create account
  6. Log in


Revision 2386, 15.0 KB (checked in by sezero, 4 years ago)

set svn:eol-style property.

  • Property svn:mime-type set to text/plain
  • Property svn:eol-style set to native
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
7        T. J. Dekker, "A Floating-Point Technique for Extending the
8        Available Precision", Numer. Math. 18 (1971), pp. 224-242
12        "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994
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
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
23(based in part on papers by Clinger and Steele & White: see the
24references in the above paper).
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
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.
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:
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
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):
63        strtof
64        strtod
65        strtodd
66        strtopd
67        strtopf
68        strtopx
69        strtopxL
70        strtopQ
72one with rounding direction specified:
74        strtorf
75        strtord
76        strtordd
77        strtorx
78        strtorxL
79        strtorQ
81and one for computing an interval (at most one bit wide) that contains
82the decimal number:
84        strtoIf
85        strtoId
86        strtoIdd
87        strtoIx
88        strtoIxL
89        strtoIQ
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.
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.
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.
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
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.
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
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.
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.
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.
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.
177For binary -> decimal conversions, I've provided just one family
178of helper routines:
180        g_ffmt
181        g_dfmt
182        g_ddfmt
183        g_xfmt
184        g_xLfmt
185        g_Qfmt
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.
200For an example of more general conversions based on dtoa(), see
201netlib's "printf.c from ampl/solvers".
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.
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 __.
215Various comments about possible #defines appear in gdtoaimp.h,
216but for most purposes, arith.h should set suitable #defines.
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).
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.
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
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 .
251(The same ftp directory has source for testbase.)
253Some system-dependent additions to CFLAGS in the makefile:
255        HU-UX: -Aa -Ae
256        OSF (DEC Unix): -ieee_with_no_inexact
257        SunOS 4.1x: -DKR_headers -DBad_float_h
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:
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
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:
317        fpi->emin = 1 - 1023 - 52
318        fpi->emax = 1046 - 1023 - 52
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.
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.
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 '.'.
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.
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.
356Please send comments to David M. Gay (dmg at acm dot org, with " at "
357changed at "@" and " dot " changed to ".").
Note: See TracBrowser for help on using the browser.