1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322

This directory contains source for a library of binary > decimal
and decimal > binary conversion routines, for single, double,
and extendedprecision IEEE binary floatingpoint arithmetic, and
other IEEElike binary floatingpoint, including "double double",
as in
T. J. Dekker, "A FloatingPoint Technique for Extending the
Available Precision", Numer. Math. 18 (1971), pp. 224242
and
"Inside Macintosh: PowerPC Numerics", AddisonWesley, 1994
The conversion routines use doubleprecision floatingpoint 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 BinaryDecimal and
DecimalBinary Conversions", Numerical Analysis Manuscript
No. 9010, Bell Labs, Murray Hill, 1990;
http://cm.belllabs.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 IBMmainframe doubleprecision 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 32bit integers of type ULong, a signed 32bit
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 floatingpoint 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 roundnearest:
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 oneulp adjustments on values in the
desired floatingpoint 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 doubleprecision
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 doubledouble 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 53bit precision near
the bottom of the exponent range. To get correct rounding, this
would (in general) require two calls on strtodg (one specifying
126bit arithmetic, then, if necessary, one specifying 53bit
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 (roundtonearest) and, on
systems with IEEE extendedprecision registers, doubleprecision
(53bit) 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 systemdependent header
file.
The values returned for NaNs may be signaling NaNs on some systems,
since the rules for distinguishing signaling from quiet NaNs are
systemdependent. You can easily fix this by suitably modifying the
ULto* routines in strtor*.c.
C99's hexadecimal floatingpoint 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 floatingpoint 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 doubledouble > decimal, g_ddfmt() assumes IEEElike 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 BinaryDecimal
Conversion", manuscript, May 1991,
ftp://ftp.ee.lbl.gov/testbasereport.ps.Z .
(The same ftp directory has source for testbase.)
Some systemdependent additions to CFLAGS in the makefile:
HUUX: 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
IEEElike 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 32bit
unsigned integer type. Floatingpoint 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
floatingpoint 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^(b1023),
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 = (b1023) 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 floatingpoint 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 2C463
600 Mountain Avenue
Murray Hill, NJ 079740636, U.S.A.
dmg@research.belllabs.com
