1 /**************************************************************** 2 * 3 * The author of this software 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 29 #include <machine/ieeefp.h> 30 #include <math.h> 31 #include <float.h> 32 #include <errno.h> 33 #include <assert.h> 34 #include <sys/types.h> 35 #include "../locale/setlocale.h" 36 37 #ifdef __IEEE_LITTLE_ENDIAN 38 #define IEEE_8087 39 #endif 40 41 #ifdef __IEEE_BIG_ENDIAN 42 #define IEEE_MC68k 43 #endif 44 45 #ifdef __Z8000__ 46 #define Just_16 47 #endif 48 49 #ifdef DEBUG 50 #include "stdio.h" 51 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 52 #endif 53 54 #ifdef Unsigned_Shifts 55 #define Sign_Extend(a,b) if (b < 0) a |= (__uint32_t)0xffff0000; 56 #else 57 #define Sign_Extend(a,b) /*no-op*/ 58 #endif 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 /* If we are going to examine or modify specific bits in a double using 65 the word0 and/or word1 macros, then we must wrap the double inside 66 a union. This is necessary to avoid undefined behavior according to 67 the ANSI C spec. */ 68 union double_union 69 { 70 double d; 71 __uint32_t i[2]; 72 }; 73 74 #ifdef IEEE_8087 75 #define word0(x) (x.i[1]) 76 #define word1(x) (x.i[0]) 77 #else 78 #define word0(x) (x.i[0]) 79 #define word1(x) (x.i[1]) 80 #endif 81 82 /* The following is taken from gdtoaimp.h for use with new strtod, but 83 adjusted to avoid invalid type-punning. */ 84 typedef __int32_t Long; 85 86 /* Unfortunately, because __ULong might be a different type than 87 __uint32_t, we can't re-use union double_union as-is without 88 further edits in strtod.c. */ 89 typedef union { double d; __ULong i[2]; } U; 90 91 #define dword0(x) word0(x) 92 #define dword1(x) word1(x) 93 #define dval(x) (x.d) 94 95 #undef SI 96 #ifdef Sudden_Underflow 97 #define SI 1 98 #else 99 #define SI 0 100 #endif 101 102 #define Storeinc(a,b,c) (*(a)++ = ((b) << 16) | ((c) & 0xffff)) 103 104 /* #define P DBL_MANT_DIG */ 105 /* Ten_pmax = floor(P*log(2)/log(5)) */ 106 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 107 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 108 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 109 110 #if defined(IEEE_8087) + defined(IEEE_MC68k) 111 #if defined (_DOUBLE_IS_32BITS) 112 #define Exp_shift 23 113 #define Exp_shift1 23 114 #define Exp_msk1 ((__uint32_t)0x00800000L) 115 #define Exp_msk11 ((__uint32_t)0x00800000L) 116 #define Exp_mask ((__uint32_t)0x7f800000L) 117 #define P 24 118 #define Bias 127 119 #define NO_HEX_FP /* not supported in this case */ 120 #define IEEE_Arith 121 #define Emin (-126) 122 #define Exp_1 ((__uint32_t)0x3f800000L) 123 #define Exp_11 ((__uint32_t)0x3f800000L) 124 #define Ebits 8 125 #define Frac_mask ((__uint32_t)0x007fffffL) 126 #define Frac_mask1 ((__uint32_t)0x007fffffL) 127 #define Ten_pmax 10 128 #define Sign_bit ((__uint32_t)0x80000000L) 129 #define Ten_pmax 10 130 #define Bletch 2 131 #define Bndry_mask ((__uint32_t)0x007fffffL) 132 #define Bndry_mask1 ((__uint32_t)0x007fffffL) 133 #define LSB 1 134 #define Sign_bit ((__uint32_t)0x80000000L) 135 #define Log2P 1 136 #define Tiny0 0 137 #define Tiny1 1 138 #define Quick_max 5 139 #define Int_max 6 140 #define Infinite(x) (word0(x) == ((__uint32_t)0x7f800000L)) 141 #undef word0 142 #undef word1 143 #undef dword0 144 #undef dword1 145 146 #define word0(x) (x.i[0]) 147 #define word1(x) 0 148 #define dword0(x) word0(x) 149 #define dword1(x) 0 150 #else 151 152 #define Exp_shift 20 153 #define Exp_shift1 20 154 #define Exp_msk1 ((__uint32_t)0x100000L) 155 #define Exp_msk11 ((__uint32_t)0x100000L) 156 #define Exp_mask ((__uint32_t)0x7ff00000L) 157 #define P 53 158 #define Bias 1023 159 #define IEEE_Arith 160 #define Emin (-1022) 161 #define Exp_1 ((__uint32_t)0x3ff00000L) 162 #define Exp_11 ((__uint32_t)0x3ff00000L) 163 #define Ebits 11 164 #define Frac_mask ((__uint32_t)0xfffffL) 165 #define Frac_mask1 ((__uint32_t)0xfffffL) 166 #define Ten_pmax 22 167 #define Bletch 0x10 168 #define Bndry_mask ((__uint32_t)0xfffffL) 169 #define Bndry_mask1 ((__uint32_t)0xfffffL) 170 #define LSB 1 171 #define Sign_bit ((__uint32_t)0x80000000L) 172 #define Log2P 1 173 #define Tiny0 0 174 #define Tiny1 1 175 #define Quick_max 14 176 #define Int_max 14 177 #define Infinite(x) (word0(x) == ((__uint32_t)0x7ff00000L)) /* sufficient test for here */ 178 179 #endif /* !_DOUBLE_IS_32BITS */ 180 181 #ifndef Flt_Rounds 182 #ifdef FLT_ROUNDS 183 #define Flt_Rounds FLT_ROUNDS 184 #else 185 #define Flt_Rounds 1 186 #endif 187 #endif /*Flt_Rounds*/ 188 189 #else /* !IEEE_8087 && !IEEE_MC68k */ 190 #undef Sudden_Underflow 191 #define Sudden_Underflow 192 #ifdef IBM 193 #define Flt_Rounds 0 194 #define Exp_shift 24 195 #define Exp_shift1 24 196 #define Exp_msk1 ((__uint32_t)0x1000000L) 197 #define Exp_msk11 ((__uint32_t)0x1000000L) 198 #define Exp_mask ((__uint32_t)0x7f000000L) 199 #define P 14 200 #define Bias 65 201 #define Exp_1 ((__uint32_t)0x41000000L) 202 #define Exp_11 ((__uint32_t)0x41000000L) 203 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 204 #define Frac_mask ((__uint32_t)0xffffffL) 205 #define Frac_mask1 ((__uint32_t)0xffffffL) 206 #define Bletch 4 207 #define Ten_pmax 22 208 #define Bndry_mask ((__uint32_t)0xefffffL) 209 #define Bndry_mask1 ((__uint32_t)0xffffffL) 210 #define LSB 1 211 #define Sign_bit ((__uint32_t)0x80000000L) 212 #define Log2P 4 213 #define Tiny0 ((__uint32_t)0x100000L) 214 #define Tiny1 0 215 #define Quick_max 14 216 #define Int_max 15 217 #else /* VAX */ 218 #define Flt_Rounds 1 219 #define Exp_shift 23 220 #define Exp_shift1 7 221 #define Exp_msk1 0x80 222 #define Exp_msk11 ((__uint32_t)0x800000L) 223 #define Exp_mask ((__uint32_t)0x7f80L) 224 #define P 56 225 #define Bias 129 226 #define Exp_1 ((__uint32_t)0x40800000L) 227 #define Exp_11 ((__uint32_t)0x4080L) 228 #define Ebits 8 229 #define Frac_mask ((__uint32_t)0x7fffffL) 230 #define Frac_mask1 ((__uint32_t)0xffff007fL) 231 #define Ten_pmax 24 232 #define Bletch 2 233 #define Bndry_mask ((__uint32_t)0xffff007fL) 234 #define Bndry_mask1 ((__uint32_t)0xffff007fL) 235 #define LSB ((__uint32_t)0x10000L) 236 #define Sign_bit ((__uint32_t)0x8000L) 237 #define Log2P 1 238 #define Tiny0 0x80 239 #define Tiny1 0 240 #define Quick_max 15 241 #define Int_max 15 242 #endif 243 #endif 244 245 #ifndef IEEE_Arith 246 #define ROUND_BIASED 247 #else 248 #define Scale_Bit 0x10 249 #if defined(_DOUBLE_IS_32BITS) && defined(__v800) 250 #define n_bigtens 2 251 #else 252 #define n_bigtens 5 253 #endif 254 #endif 255 256 #ifdef IBM 257 #define n_bigtens 3 258 #endif 259 260 #ifdef VAX 261 #define n_bigtens 2 262 #endif 263 264 #ifndef __NO_INFNAN_CHECK 265 #define INFNAN_CHECK 266 #endif 267 268 #ifdef RND_PRODQUOT 269 #define rounded_product(a,b) a = rnd_prod(a, b) 270 #define rounded_quotient(a,b) a = rnd_quot(a, b) 271 #ifdef KR_headers 272 extern double rnd_prod(), rnd_quot(); 273 #else 274 extern double rnd_prod(double, double), rnd_quot(double, double); 275 #endif 276 #else 277 #define rounded_product(a,b) a *= b 278 #define rounded_quotient(a,b) a /= b 279 #endif 280 281 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 282 #define Big1 ((__uint32_t)0xffffffffL) 283 284 #ifndef Just_16 285 /* When Pack_32 is not defined, we store 16 bits per 32-bit long. 286 * This makes some inner loops simpler and sometimes saves work 287 * during multiplications, but it often seems to make things slightly 288 * slower. Hence the default is now to store 32 bits per long. 289 */ 290 291 #ifndef Pack_32 292 #define Pack_32 293 #endif 294 #else /* Just_16 */ 295 #ifndef Pack_16 296 #define Pack_16 297 #endif 298 #endif /* Just_16 */ 299 300 #ifdef Pack_32 301 #define ULbits 32 302 #define kshift 5 303 #define kmask 31 304 #define ALL_ON 0xffffffff 305 #else 306 #define ULbits 16 307 #define kshift 4 308 #define kmask 15 309 #define ALL_ON 0xffff 310 #endif 311 312 #ifdef __cplusplus 313 extern "C" double strtod(const char *s00, char **se); 314 extern "C" char *dtoa(double d, int mode, int ndigits, 315 int *decpt, int *sign, char **rve); 316 #endif 317 318 319 typedef struct _Bigint _Bigint; 320 321 #define Balloc _Balloc 322 #define Bfree _Bfree 323 #define multadd __multadd 324 #define s2b __s2b 325 #define lo0bits __lo0bits 326 #define hi0bits __hi0bits 327 #define i2b __i2b 328 #define mult __multiply 329 #define pow5mult __pow5mult 330 #define lshift __lshift 331 #define match __match 332 #define cmp __mcmp 333 #define diff __mdiff 334 #define ulp __ulp 335 #define b2d __b2d 336 #define d2b __d2b 337 #define ratio __ratio 338 #define any_on __any_on 339 #define gethex __gethex 340 #define copybits __copybits 341 #define hexnan __hexnan 342 343 #if !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) 344 #define __get_hexdig(x) __hexdig[x] /* NOTE: must evaluate arg only once */ 345 #else /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ 346 #define __get_hexdig(x) __hexdig_fun(x) 347 #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ 348 349 #define tens __mprec_tens 350 #define bigtens __mprec_bigtens 351 #define tinytens __mprec_tinytens 352 353 char *__alloc_dtoa_result(int len); 354 355 int __mprec_register_exit(void); 356 357 struct _Bigint 358 { 359 int _k, _maxwds, _sign, _wds; 360 uint32_t _x[]; 361 }; 362 363 struct FPI; 364 double ulp (double x); 365 double b2d (_Bigint *a , int *e); 366 _Bigint * Balloc (int k); 367 void Bfree (_Bigint *v); 368 _Bigint * multadd (_Bigint *, int, int); 369 _Bigint * s2b (const char*, int, int, __ULong); 370 _Bigint * i2b (int); 371 _Bigint * mult (_Bigint *, _Bigint *); 372 _Bigint * pow5mult (_Bigint *, int k); 373 int hi0bits (__ULong); 374 int lo0bits (__ULong *); 375 _Bigint * d2b (double d, int *e, int *bits); 376 _Bigint * lshift (_Bigint *b, int k); 377 int match (const char**, char*); 378 _Bigint * diff (_Bigint *a, _Bigint *b); 379 int cmp (_Bigint *a, _Bigint *b); 380 int gethex (const char **sp, const struct FPI *fpi, Long *exp, _Bigint **bp, int sign, locale_t loc); 381 double ratio (_Bigint *a, _Bigint *b); 382 __ULong any_on (_Bigint *b, int k); 383 void copybits (__ULong *c, int n, _Bigint *b); 384 double _strtod_l (const char *__restrict s00, 385 char **__restrict se, locale_t loc); 386 #if defined (_HAVE_LONG_DOUBLE) && !defined (_LDBL_EQ_DBL) || 1 387 int _strtorx_l (const char *, char **, int, 388 void *, locale_t); 389 int _strtodg_l (const char *s00, char **se, 390 struct FPI *fpi, Long *exp, __ULong *bits, 391 locale_t); 392 #endif /* _HAVE_LONG_DOUBLE && !_LDBL_EQ_DBL */ 393 394 #if defined(PREFER_SIZE_OVER_SPEED) || defined(__OPTIMIZE_SIZE__) || defined(_SMALL_HEXDIG) 395 unsigned char __hexdig_fun (unsigned char); 396 #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ 397 #ifdef INFNAN_CHECK 398 int hexnan (const char **sp, const struct FPI *fpi, __ULong *x0); 399 #endif 400 401 #define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int)) 402 403 extern const double tinytens[]; 404 extern const double bigtens[]; 405 extern const double tens[]; 406 #if !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) 407 extern const unsigned char __hexdig[]; 408 #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */ 409 410 411 double _mprec_log10 (int); 412
