|
|
1.1 root 1: /****************************************************************
2: *
3: * The author of this software (_dtoa, strtod) is David M. Gay.
4: *
5: * Copyright (c) 1991 by AT&T.
6: *
7: * Permission to use, copy, modify, and distribute this software for any
8: * purpose without fee is hereby granted, provided that this entire notice
9: * is included in all copies of any software which is or includes a copy
10: * or modification of this software and in all copies of the supporting
11: * documentation for such software.
12: *
13: * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14: * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
15: * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16: * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17: *
18: ***************************************************************/
19:
20: /* Please send bug reports to
21: David M. Gay
22: AT&T Bell Laboratories, Room 2C-463
23: 600 Mountain Avenue
24: Murray Hill, NJ 07974-2070
25: U.S.A.
26: [email protected] or research!dmg
27: */
28: #include <stdlib.h>
29: #include <string.h>
30: #define _RESEARCH_SOURCE
31: #include <errno.h>
32: #include <float.h>
33: #include <math.h>
34: #define CONST const
35:
36: #ifdef vax
37: #define VAX
38: #else
39: #ifdef mips
40: #define IEEE_MC68k
41: #endif
42: #endif
43:
44: /*
45: * #define IEEE_8087 for IEEE-arithmetic machines where the least
46: * significant byte has the lowest address.
47: * #define IEEE_MC68k for IEEE-arithmetic machines where the most
48: * significant byte has the lowest address.
49: * #define Sudden_Underflow for IEEE-format machines without gradual
50: * underflow (i.e., that flush to zero on underflow).
51: * #define IBM for IBM mainframe-style floating-point arithmetic.
52: * #define VAX for VAX-style floating-point arithmetic.
53: * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
54: * #define No_leftright to omit left-right logic in fast floating-point
55: * computation of dtoa.
56: * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
57: * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
58: * that use extended-precision instructions to compute rounded
59: * products and quotients) with IBM.
60: * #define ROUND_BIASED for IEEE-format with biased rounding.
61: * #define Inaccurate_Divide for IEEE-format with correctly rounded
62: * products but inaccurate quotients, e.g., for Intel i860.
63: * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
64: * integer arithmetic. Whether this speeds things up or slows things
65: * down depends on the machine and the number being converted.
66: */
67:
68: #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
69: Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
70: #endif
71:
72: typedef union {
73: double d;
74: unsigned long ul[2];
75: } Dul;
76:
77: #ifdef IEEE_8087
78: #define word0(x) ((x).ul[1])
79: #define word1(x) ((x).ul[0])
80: #else
81: #define word0(x) ((x).ul[0])
82: #define word1(x) ((x).ul[1])
83: #endif
84:
85: #ifdef Unsigned_Shifts
86: #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
87: #else
88: #define Sign_Extend(a,b) /*no-op*/
89: #endif
90:
91: /* The following definition of Storeinc is appropriate for MIPS processors.
92: * An alternative that might be better on some machines is
93: * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
94: */
95: #if defined(IEEE_8087) + defined(VAX)
96: #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
97: ((unsigned short *)a)[0] = (unsigned short)c, a++)
98: #else
99: #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
100: ((unsigned short *)a)[1] = (unsigned short)c, a++)
101: #endif
102:
103: /* #define P DBL_MANT_DIG */
104: /* Ten_pmax = floor(P*log(2)/log(5)) */
105: /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
106: /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
107: /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
108:
109: #if defined(IEEE_8087) + defined(IEEE_MC68k)
110: #define Exp_shift 20
111: #define Exp_shift1 20
112: #define Exp_msk1 0x100000
113: #define Exp_msk11 0x100000
114: #define Exp_mask 0x7ff00000
115: #define P 53
116: #define Bias 1023
117: #define IEEE_Arith
118: #define Emin (-1022)
119: #define Exp_1 0x3ff00000
120: #define Exp_11 0x3ff00000
121: #define Ebits 11
122: #define Frac_mask 0xfffff
123: #define Frac_mask1 0xfffff
124: #define Ten_pmax 22
125: #define Bletch 0x10
126: #define Bndry_mask 0xfffff
127: #define Bndry_mask1 0xfffff
128: #define LSB 1
129: #define Sign_bit 0x80000000
130: #define Log2P 1
131: #define Tiny0 0
132: #define Tiny1 1
133: #define Quick_max 14
134: #define Int_max 14
135: #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
136: #else
137: #undef Sudden_Underflow
138: #define Sudden_Underflow
139: #ifdef IBM
140: #define Exp_shift 24
141: #define Exp_shift1 24
142: #define Exp_msk1 0x1000000
143: #define Exp_msk11 0x1000000
144: #define Exp_mask 0x7f000000
145: #define P 14
146: #define Bias 65
147: #define Exp_1 0x41000000
148: #define Exp_11 0x41000000
149: #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
150: #define Frac_mask 0xffffff
151: #define Frac_mask1 0xffffff
152: #define Bletch 4
153: #define Ten_pmax 22
154: #define Bndry_mask 0xefffff
155: #define Bndry_mask1 0xffffff
156: #define LSB 1
157: #define Sign_bit 0x80000000
158: #define Log2P 4
159: #define Tiny0 0x100000
160: #define Tiny1 0
161: #define Quick_max 14
162: #define Int_max 15
163: #else /* VAX */
164: #define Exp_shift 23
165: #define Exp_shift1 7
166: #define Exp_msk1 0x80
167: #define Exp_msk11 0x800000
168: #define Exp_mask 0x7f80
169: #define P 56
170: #define Bias 129
171: #define Exp_1 0x40800000
172: #define Exp_11 0x4080
173: #define Ebits 8
174: #define Frac_mask 0x7fffff
175: #define Frac_mask1 0xffff007f
176: #define Ten_pmax 24
177: #define Bletch 2
178: #define Bndry_mask 0xffff007f
179: #define Bndry_mask1 0xffff007f
180: #define LSB 0x10000
181: #define Sign_bit 0x8000
182: #define Log2P 1
183: #define Tiny0 0x80
184: #define Tiny1 0
185: #define Quick_max 15
186: #define Int_max 15
187: #endif
188: #endif
189:
190: #ifndef IEEE_Arith
191: #define ROUND_BIASED
192: #endif
193:
194: #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
195: #define Big1 0xffffffff
196:
197: #ifndef Just_16
198: /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
199: * This makes some inner loops simpler and sometimes saves work
200: * during multiplications, but it often seems to make things slightly
201: * slower. Hence the default is now to store 32 bits per long.
202: */
203: #ifndef Pack_32
204: #define Pack_32
205: #endif
206: #endif
207:
208: #define Kmax 15
209:
210: struct
211: Bigint {
212: struct Bigint *next;
213: int k, maxwds, sign, wds;
214: unsigned long x[1];
215: };
216:
217: typedef struct Bigint Bigint;
218:
219: /* This routines shouldn't be visible externally */
220: extern Bigint *_Balloc(int);
221: extern void _Bfree(Bigint *);
222: extern Bigint *_multadd(Bigint *, int, int);
223: extern int _hi0bits(unsigned long);
224: extern Bigint *_mult(Bigint *, Bigint *);
225: extern Bigint *_pow5mult(Bigint *, int);
226: extern Bigint *_lshift(Bigint *, int);
227: extern int _cmp(Bigint *, Bigint *);
228: extern Bigint *_diff(Bigint *, Bigint *);
229: extern Bigint *_d2b(double, int *, int *);
230: extern Bigint *_i2b(int);
231:
232: extern double _tens[], _bigtens[], _tinytens[];
233:
234: #ifdef IEEE_Arith
235: #define n_bigtens 5
236: #else
237: #ifdef IBM
238: #define n_bigtens 3
239: #else
240: #define n_bigtens 2
241: #endif
242: #endif
243:
244: #define Balloc(x) _Balloc(x)
245: #define Bfree(x) _Bfree(x)
246: #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
247: y->wds*sizeof(long) + 2*sizeof(int))
248: #define multadd(x,y,z) _multadd(x,y,z)
249: #define hi0bits(x) _hi0bits(x)
250: #define i2b(x) _i2b(x)
251: #define mult(x,y) _mult(x,y)
252: #define pow5mult(x,y) _pow5mult(x,y)
253: #define lshift(x,y) _lshift(x,y)
254: #define cmp(x,y) _cmp(x,y)
255: #define diff(x,y) _diff(x,y)
256: #define d2b(x,y,z) _d2b(x,y,z)
257:
258: #define tens _tens
259: #define bigtens _bigtens
260: #define tinytens _tinytens
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.