13e12c5d1SDavid du Colombier /**************************************************************** 23e12c5d1SDavid du Colombier 3*7dd7cddfSDavid du Colombier The author of this software (_dtoa, strtod) is David M. Gay. 4*7dd7cddfSDavid du Colombier Please send bug reports to 53e12c5d1SDavid du Colombier David M. Gay 6*7dd7cddfSDavid du Colombier Bell Laboratories, Room 2C-463 73e12c5d1SDavid du Colombier 600 Mountain Avenue 83e12c5d1SDavid du Colombier Murray Hill, NJ 07974-2070 93e12c5d1SDavid du Colombier U.S.A. 10*7dd7cddfSDavid du Colombier dmg@research.bell-labs.com 113e12c5d1SDavid du Colombier */ 123e12c5d1SDavid du Colombier #include <stdlib.h> 133e12c5d1SDavid du Colombier #include <string.h> 143e12c5d1SDavid du Colombier #define _RESEARCH_SOURCE 153e12c5d1SDavid du Colombier #include <errno.h> 163e12c5d1SDavid du Colombier #include <float.h> 173e12c5d1SDavid du Colombier #include <math.h> 183e12c5d1SDavid du Colombier #define CONST const 193e12c5d1SDavid du Colombier 203e12c5d1SDavid du Colombier /* 213e12c5d1SDavid du Colombier * #define IEEE_8087 for IEEE-arithmetic machines where the least 223e12c5d1SDavid du Colombier * significant byte has the lowest address. 233e12c5d1SDavid du Colombier * #define IEEE_MC68k for IEEE-arithmetic machines where the most 243e12c5d1SDavid du Colombier * significant byte has the lowest address. 253e12c5d1SDavid du Colombier * #define Sudden_Underflow for IEEE-format machines without gradual 263e12c5d1SDavid du Colombier * underflow (i.e., that flush to zero on underflow). 273e12c5d1SDavid du Colombier * #define IBM for IBM mainframe-style floating-point arithmetic. 283e12c5d1SDavid du Colombier * #define VAX for VAX-style floating-point arithmetic. 293e12c5d1SDavid du Colombier * #define Unsigned_Shifts if >> does treats its left operand as unsigned. 303e12c5d1SDavid du Colombier * #define No_leftright to omit left-right logic in fast floating-point 313e12c5d1SDavid du Colombier * computation of dtoa. 323e12c5d1SDavid du Colombier * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. 333e12c5d1SDavid du Colombier * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 343e12c5d1SDavid du Colombier * that use extended-precision instructions to compute rounded 353e12c5d1SDavid du Colombier * products and quotients) with IBM. 363e12c5d1SDavid du Colombier * #define ROUND_BIASED for IEEE-format with biased rounding. 373e12c5d1SDavid du Colombier * #define Inaccurate_Divide for IEEE-format with correctly rounded 383e12c5d1SDavid du Colombier * products but inaccurate quotients, e.g., for Intel i860. 393e12c5d1SDavid du Colombier * #define Just_16 to store 16 bits per 32-bit long when doing high-precision 403e12c5d1SDavid du Colombier * integer arithmetic. Whether this speeds things up or slows things 413e12c5d1SDavid du Colombier * down depends on the machine and the number being converted. 423e12c5d1SDavid du Colombier */ 433e12c5d1SDavid du Colombier 443e12c5d1SDavid du Colombier #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 453e12c5d1SDavid du Colombier Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 463e12c5d1SDavid du Colombier #endif 473e12c5d1SDavid du Colombier 483e12c5d1SDavid du Colombier typedef union { 493e12c5d1SDavid du Colombier double d; 503e12c5d1SDavid du Colombier unsigned long ul[2]; 513e12c5d1SDavid du Colombier } Dul; 523e12c5d1SDavid du Colombier 533e12c5d1SDavid du Colombier #ifdef IEEE_8087 543e12c5d1SDavid du Colombier #define word0(x) ((x).ul[1]) 553e12c5d1SDavid du Colombier #define word1(x) ((x).ul[0]) 563e12c5d1SDavid du Colombier #else 573e12c5d1SDavid du Colombier #define word0(x) ((x).ul[0]) 583e12c5d1SDavid du Colombier #define word1(x) ((x).ul[1]) 593e12c5d1SDavid du Colombier #endif 603e12c5d1SDavid du Colombier 613e12c5d1SDavid du Colombier #ifdef Unsigned_Shifts 623e12c5d1SDavid du Colombier #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000; 633e12c5d1SDavid du Colombier #else 643e12c5d1SDavid du Colombier #define Sign_Extend(a,b) /*no-op*/ 653e12c5d1SDavid du Colombier #endif 663e12c5d1SDavid du Colombier 673e12c5d1SDavid du Colombier /* The following definition of Storeinc is appropriate for MIPS processors. 683e12c5d1SDavid du Colombier * An alternative that might be better on some machines is 693e12c5d1SDavid du Colombier * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 703e12c5d1SDavid du Colombier */ 713e12c5d1SDavid du Colombier #if defined(IEEE_8087) + defined(VAX) 723e12c5d1SDavid du Colombier #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 733e12c5d1SDavid du Colombier ((unsigned short *)a)[0] = (unsigned short)c, a++) 743e12c5d1SDavid du Colombier #else 753e12c5d1SDavid du Colombier #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 763e12c5d1SDavid du Colombier ((unsigned short *)a)[1] = (unsigned short)c, a++) 773e12c5d1SDavid du Colombier #endif 783e12c5d1SDavid du Colombier 793e12c5d1SDavid du Colombier /* #define P DBL_MANT_DIG */ 803e12c5d1SDavid du Colombier /* Ten_pmax = floor(P*log(2)/log(5)) */ 813e12c5d1SDavid du Colombier /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 823e12c5d1SDavid du Colombier /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 833e12c5d1SDavid du Colombier /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 843e12c5d1SDavid du Colombier 853e12c5d1SDavid du Colombier #if defined(IEEE_8087) + defined(IEEE_MC68k) 863e12c5d1SDavid du Colombier #define Exp_shift 20 873e12c5d1SDavid du Colombier #define Exp_shift1 20 883e12c5d1SDavid du Colombier #define Exp_msk1 0x100000 893e12c5d1SDavid du Colombier #define Exp_msk11 0x100000 903e12c5d1SDavid du Colombier #define Exp_mask 0x7ff00000 913e12c5d1SDavid du Colombier #define P 53 923e12c5d1SDavid du Colombier #define Bias 1023 933e12c5d1SDavid du Colombier #define IEEE_Arith 943e12c5d1SDavid du Colombier #define Emin (-1022) 953e12c5d1SDavid du Colombier #define Exp_1 0x3ff00000 963e12c5d1SDavid du Colombier #define Exp_11 0x3ff00000 973e12c5d1SDavid du Colombier #define Ebits 11 983e12c5d1SDavid du Colombier #define Frac_mask 0xfffff 993e12c5d1SDavid du Colombier #define Frac_mask1 0xfffff 1003e12c5d1SDavid du Colombier #define Ten_pmax 22 1013e12c5d1SDavid du Colombier #define Bletch 0x10 1023e12c5d1SDavid du Colombier #define Bndry_mask 0xfffff 1033e12c5d1SDavid du Colombier #define Bndry_mask1 0xfffff 1043e12c5d1SDavid du Colombier #define LSB 1 1053e12c5d1SDavid du Colombier #define Sign_bit 0x80000000 1063e12c5d1SDavid du Colombier #define Log2P 1 1073e12c5d1SDavid du Colombier #define Tiny0 0 1083e12c5d1SDavid du Colombier #define Tiny1 1 1093e12c5d1SDavid du Colombier #define Quick_max 14 1103e12c5d1SDavid du Colombier #define Int_max 14 1113e12c5d1SDavid du Colombier #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */ 1123e12c5d1SDavid du Colombier #else 1133e12c5d1SDavid du Colombier #undef Sudden_Underflow 1143e12c5d1SDavid du Colombier #define Sudden_Underflow 1153e12c5d1SDavid du Colombier #ifdef IBM 1163e12c5d1SDavid du Colombier #define Exp_shift 24 1173e12c5d1SDavid du Colombier #define Exp_shift1 24 1183e12c5d1SDavid du Colombier #define Exp_msk1 0x1000000 1193e12c5d1SDavid du Colombier #define Exp_msk11 0x1000000 1203e12c5d1SDavid du Colombier #define Exp_mask 0x7f000000 1213e12c5d1SDavid du Colombier #define P 14 1223e12c5d1SDavid du Colombier #define Bias 65 1233e12c5d1SDavid du Colombier #define Exp_1 0x41000000 1243e12c5d1SDavid du Colombier #define Exp_11 0x41000000 1253e12c5d1SDavid du Colombier #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 1263e12c5d1SDavid du Colombier #define Frac_mask 0xffffff 1273e12c5d1SDavid du Colombier #define Frac_mask1 0xffffff 1283e12c5d1SDavid du Colombier #define Bletch 4 1293e12c5d1SDavid du Colombier #define Ten_pmax 22 1303e12c5d1SDavid du Colombier #define Bndry_mask 0xefffff 1313e12c5d1SDavid du Colombier #define Bndry_mask1 0xffffff 1323e12c5d1SDavid du Colombier #define LSB 1 1333e12c5d1SDavid du Colombier #define Sign_bit 0x80000000 1343e12c5d1SDavid du Colombier #define Log2P 4 1353e12c5d1SDavid du Colombier #define Tiny0 0x100000 1363e12c5d1SDavid du Colombier #define Tiny1 0 1373e12c5d1SDavid du Colombier #define Quick_max 14 1383e12c5d1SDavid du Colombier #define Int_max 15 1393e12c5d1SDavid du Colombier #else /* VAX */ 1403e12c5d1SDavid du Colombier #define Exp_shift 23 1413e12c5d1SDavid du Colombier #define Exp_shift1 7 1423e12c5d1SDavid du Colombier #define Exp_msk1 0x80 1433e12c5d1SDavid du Colombier #define Exp_msk11 0x800000 1443e12c5d1SDavid du Colombier #define Exp_mask 0x7f80 1453e12c5d1SDavid du Colombier #define P 56 1463e12c5d1SDavid du Colombier #define Bias 129 1473e12c5d1SDavid du Colombier #define Exp_1 0x40800000 1483e12c5d1SDavid du Colombier #define Exp_11 0x4080 1493e12c5d1SDavid du Colombier #define Ebits 8 1503e12c5d1SDavid du Colombier #define Frac_mask 0x7fffff 1513e12c5d1SDavid du Colombier #define Frac_mask1 0xffff007f 1523e12c5d1SDavid du Colombier #define Ten_pmax 24 1533e12c5d1SDavid du Colombier #define Bletch 2 1543e12c5d1SDavid du Colombier #define Bndry_mask 0xffff007f 1553e12c5d1SDavid du Colombier #define Bndry_mask1 0xffff007f 1563e12c5d1SDavid du Colombier #define LSB 0x10000 1573e12c5d1SDavid du Colombier #define Sign_bit 0x8000 1583e12c5d1SDavid du Colombier #define Log2P 1 1593e12c5d1SDavid du Colombier #define Tiny0 0x80 1603e12c5d1SDavid du Colombier #define Tiny1 0 1613e12c5d1SDavid du Colombier #define Quick_max 15 1623e12c5d1SDavid du Colombier #define Int_max 15 1633e12c5d1SDavid du Colombier #endif 1643e12c5d1SDavid du Colombier #endif 1653e12c5d1SDavid du Colombier 1663e12c5d1SDavid du Colombier #ifndef IEEE_Arith 1673e12c5d1SDavid du Colombier #define ROUND_BIASED 1683e12c5d1SDavid du Colombier #endif 1693e12c5d1SDavid du Colombier 1703e12c5d1SDavid du Colombier #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 1713e12c5d1SDavid du Colombier #define Big1 0xffffffff 1723e12c5d1SDavid du Colombier 1733e12c5d1SDavid du Colombier #ifndef Just_16 1743e12c5d1SDavid du Colombier /* When Pack_32 is not defined, we store 16 bits per 32-bit long. 1753e12c5d1SDavid du Colombier * This makes some inner loops simpler and sometimes saves work 1763e12c5d1SDavid du Colombier * during multiplications, but it often seems to make things slightly 1773e12c5d1SDavid du Colombier * slower. Hence the default is now to store 32 bits per long. 1783e12c5d1SDavid du Colombier */ 1793e12c5d1SDavid du Colombier #ifndef Pack_32 1803e12c5d1SDavid du Colombier #define Pack_32 1813e12c5d1SDavid du Colombier #endif 1823e12c5d1SDavid du Colombier #endif 1833e12c5d1SDavid du Colombier 1843e12c5d1SDavid du Colombier #define Kmax 15 1853e12c5d1SDavid du Colombier 1863e12c5d1SDavid du Colombier struct 1873e12c5d1SDavid du Colombier Bigint { 1883e12c5d1SDavid du Colombier struct Bigint *next; 1893e12c5d1SDavid du Colombier int k, maxwds, sign, wds; 1903e12c5d1SDavid du Colombier unsigned long x[1]; 1913e12c5d1SDavid du Colombier }; 1923e12c5d1SDavid du Colombier 1933e12c5d1SDavid du Colombier typedef struct Bigint Bigint; 1943e12c5d1SDavid du Colombier 1953e12c5d1SDavid du Colombier /* This routines shouldn't be visible externally */ 1963e12c5d1SDavid du Colombier extern Bigint *_Balloc(int); 1973e12c5d1SDavid du Colombier extern void _Bfree(Bigint *); 1983e12c5d1SDavid du Colombier extern Bigint *_multadd(Bigint *, int, int); 1993e12c5d1SDavid du Colombier extern int _hi0bits(unsigned long); 2003e12c5d1SDavid du Colombier extern Bigint *_mult(Bigint *, Bigint *); 2013e12c5d1SDavid du Colombier extern Bigint *_pow5mult(Bigint *, int); 2023e12c5d1SDavid du Colombier extern Bigint *_lshift(Bigint *, int); 2033e12c5d1SDavid du Colombier extern int _cmp(Bigint *, Bigint *); 2043e12c5d1SDavid du Colombier extern Bigint *_diff(Bigint *, Bigint *); 2053e12c5d1SDavid du Colombier extern Bigint *_d2b(double, int *, int *); 2063e12c5d1SDavid du Colombier extern Bigint *_i2b(int); 2073e12c5d1SDavid du Colombier 2083e12c5d1SDavid du Colombier extern double _tens[], _bigtens[], _tinytens[]; 2093e12c5d1SDavid du Colombier 2103e12c5d1SDavid du Colombier #ifdef IEEE_Arith 2113e12c5d1SDavid du Colombier #define n_bigtens 5 2123e12c5d1SDavid du Colombier #else 2133e12c5d1SDavid du Colombier #ifdef IBM 2143e12c5d1SDavid du Colombier #define n_bigtens 3 2153e12c5d1SDavid du Colombier #else 2163e12c5d1SDavid du Colombier #define n_bigtens 2 2173e12c5d1SDavid du Colombier #endif 2183e12c5d1SDavid du Colombier #endif 2193e12c5d1SDavid du Colombier 2203e12c5d1SDavid du Colombier #define Balloc(x) _Balloc(x) 2213e12c5d1SDavid du Colombier #define Bfree(x) _Bfree(x) 2223e12c5d1SDavid du Colombier #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 2233e12c5d1SDavid du Colombier y->wds*sizeof(long) + 2*sizeof(int)) 2243e12c5d1SDavid du Colombier #define multadd(x,y,z) _multadd(x,y,z) 2253e12c5d1SDavid du Colombier #define hi0bits(x) _hi0bits(x) 2263e12c5d1SDavid du Colombier #define i2b(x) _i2b(x) 2273e12c5d1SDavid du Colombier #define mult(x,y) _mult(x,y) 2283e12c5d1SDavid du Colombier #define pow5mult(x,y) _pow5mult(x,y) 2293e12c5d1SDavid du Colombier #define lshift(x,y) _lshift(x,y) 2303e12c5d1SDavid du Colombier #define cmp(x,y) _cmp(x,y) 2313e12c5d1SDavid du Colombier #define diff(x,y) _diff(x,y) 2323e12c5d1SDavid du Colombier #define d2b(x,y,z) _d2b(x,y,z) 2333e12c5d1SDavid du Colombier 2343e12c5d1SDavid du Colombier #define tens _tens 2353e12c5d1SDavid du Colombier #define bigtens _bigtens 2363e12c5d1SDavid du Colombier #define tinytens _tinytens 237