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 dmg@research.att.com 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 /* 37 * #define IEEE_8087 for IEEE-arithmetic machines where the least 38 * significant byte has the lowest address. 39 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 40 * significant byte has the lowest address. 41 * #define Sudden_Underflow for IEEE-format machines without gradual 42 * underflow (i.e., that flush to zero on underflow). 43 * #define IBM for IBM mainframe-style floating-point arithmetic. 44 * #define VAX for VAX-style floating-point arithmetic. 45 * #define Unsigned_Shifts if >> does treats its left operand as unsigned. 46 * #define No_leftright to omit left-right logic in fast floating-point 47 * computation of dtoa. 48 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. 49 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 50 * that use extended-precision instructions to compute rounded 51 * products and quotients) with IBM. 52 * #define ROUND_BIASED for IEEE-format with biased rounding. 53 * #define Inaccurate_Divide for IEEE-format with correctly rounded 54 * products but inaccurate quotients, e.g., for Intel i860. 55 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision 56 * integer arithmetic. Whether this speeds things up or slows things 57 * down depends on the machine and the number being converted. 58 */ 59 60 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 61 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 62 #endif 63 64 typedef union { 65 double d; 66 unsigned long ul[2]; 67 } Dul; 68 69 #ifdef IEEE_8087 70 #define word0(x) ((x).ul[1]) 71 #define word1(x) ((x).ul[0]) 72 #else 73 #define word0(x) ((x).ul[0]) 74 #define word1(x) ((x).ul[1]) 75 #endif 76 77 #ifdef Unsigned_Shifts 78 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000; 79 #else 80 #define Sign_Extend(a,b) /*no-op*/ 81 #endif 82 83 /* The following definition of Storeinc is appropriate for MIPS processors. 84 * An alternative that might be better on some machines is 85 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 86 */ 87 #if defined(IEEE_8087) + defined(VAX) 88 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 89 ((unsigned short *)a)[0] = (unsigned short)c, a++) 90 #else 91 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 92 ((unsigned short *)a)[1] = (unsigned short)c, a++) 93 #endif 94 95 /* #define P DBL_MANT_DIG */ 96 /* Ten_pmax = floor(P*log(2)/log(5)) */ 97 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 98 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 99 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 100 101 #if defined(IEEE_8087) + defined(IEEE_MC68k) 102 #define Exp_shift 20 103 #define Exp_shift1 20 104 #define Exp_msk1 0x100000 105 #define Exp_msk11 0x100000 106 #define Exp_mask 0x7ff00000 107 #define P 53 108 #define Bias 1023 109 #define IEEE_Arith 110 #define Emin (-1022) 111 #define Exp_1 0x3ff00000 112 #define Exp_11 0x3ff00000 113 #define Ebits 11 114 #define Frac_mask 0xfffff 115 #define Frac_mask1 0xfffff 116 #define Ten_pmax 22 117 #define Bletch 0x10 118 #define Bndry_mask 0xfffff 119 #define Bndry_mask1 0xfffff 120 #define LSB 1 121 #define Sign_bit 0x80000000 122 #define Log2P 1 123 #define Tiny0 0 124 #define Tiny1 1 125 #define Quick_max 14 126 #define Int_max 14 127 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */ 128 #else 129 #undef Sudden_Underflow 130 #define Sudden_Underflow 131 #ifdef IBM 132 #define Exp_shift 24 133 #define Exp_shift1 24 134 #define Exp_msk1 0x1000000 135 #define Exp_msk11 0x1000000 136 #define Exp_mask 0x7f000000 137 #define P 14 138 #define Bias 65 139 #define Exp_1 0x41000000 140 #define Exp_11 0x41000000 141 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 142 #define Frac_mask 0xffffff 143 #define Frac_mask1 0xffffff 144 #define Bletch 4 145 #define Ten_pmax 22 146 #define Bndry_mask 0xefffff 147 #define Bndry_mask1 0xffffff 148 #define LSB 1 149 #define Sign_bit 0x80000000 150 #define Log2P 4 151 #define Tiny0 0x100000 152 #define Tiny1 0 153 #define Quick_max 14 154 #define Int_max 15 155 #else /* VAX */ 156 #define Exp_shift 23 157 #define Exp_shift1 7 158 #define Exp_msk1 0x80 159 #define Exp_msk11 0x800000 160 #define Exp_mask 0x7f80 161 #define P 56 162 #define Bias 129 163 #define Exp_1 0x40800000 164 #define Exp_11 0x4080 165 #define Ebits 8 166 #define Frac_mask 0x7fffff 167 #define Frac_mask1 0xffff007f 168 #define Ten_pmax 24 169 #define Bletch 2 170 #define Bndry_mask 0xffff007f 171 #define Bndry_mask1 0xffff007f 172 #define LSB 0x10000 173 #define Sign_bit 0x8000 174 #define Log2P 1 175 #define Tiny0 0x80 176 #define Tiny1 0 177 #define Quick_max 15 178 #define Int_max 15 179 #endif 180 #endif 181 182 #ifndef IEEE_Arith 183 #define ROUND_BIASED 184 #endif 185 186 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 187 #define Big1 0xffffffff 188 189 #ifndef Just_16 190 /* When Pack_32 is not defined, we store 16 bits per 32-bit long. 191 * This makes some inner loops simpler and sometimes saves work 192 * during multiplications, but it often seems to make things slightly 193 * slower. Hence the default is now to store 32 bits per long. 194 */ 195 #ifndef Pack_32 196 #define Pack_32 197 #endif 198 #endif 199 200 #define Kmax 15 201 202 struct 203 Bigint { 204 struct Bigint *next; 205 int k, maxwds, sign, wds; 206 unsigned long x[1]; 207 }; 208 209 typedef struct Bigint Bigint; 210 211 /* This routines shouldn't be visible externally */ 212 extern Bigint *_Balloc(int); 213 extern void _Bfree(Bigint *); 214 extern Bigint *_multadd(Bigint *, int, int); 215 extern int _hi0bits(unsigned long); 216 extern Bigint *_mult(Bigint *, Bigint *); 217 extern Bigint *_pow5mult(Bigint *, int); 218 extern Bigint *_lshift(Bigint *, int); 219 extern int _cmp(Bigint *, Bigint *); 220 extern Bigint *_diff(Bigint *, Bigint *); 221 extern Bigint *_d2b(double, int *, int *); 222 extern Bigint *_i2b(int); 223 224 extern double _tens[], _bigtens[], _tinytens[]; 225 226 #ifdef IEEE_Arith 227 #define n_bigtens 5 228 #else 229 #ifdef IBM 230 #define n_bigtens 3 231 #else 232 #define n_bigtens 2 233 #endif 234 #endif 235 236 #define Balloc(x) _Balloc(x) 237 #define Bfree(x) _Bfree(x) 238 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 239 y->wds*sizeof(long) + 2*sizeof(int)) 240 #define multadd(x,y,z) _multadd(x,y,z) 241 #define hi0bits(x) _hi0bits(x) 242 #define i2b(x) _i2b(x) 243 #define mult(x,y) _mult(x,y) 244 #define pow5mult(x,y) _pow5mult(x,y) 245 #define lshift(x,y) _lshift(x,y) 246 #define cmp(x,y) _cmp(x,y) 247 #define diff(x,y) _diff(x,y) 248 #define d2b(x,y,z) _d2b(x,y,z) 249 250 #define tens _tens 251 #define bigtens _bigtens 252 #define tinytens _tinytens 253