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