projects / apple / libc.git / blame
| Commit | Line | Data |
|---|---|---|
| 9385eb3d A |
1 | /**************************************************************** |
| 2 | ||
| 3 | The author of this software is David M. Gay. | |
| 4 | ||
| 5 | Copyright (C) 1998-2001 by Lucent Technologies | |
| 6 | All Rights Reserved | |
| 7 | ||
| 8 | Permission to use, copy, modify, and distribute this software and | |
| 9 | its documentation for any purpose and without fee is hereby | |
| 10 | granted, provided that the above copyright notice appear in all | |
| 11 | copies and that both that the copyright notice and this | |
| 12 | permission notice and warranty disclaimer appear in supporting | |
| 13 | documentation, and that the name of Lucent or any of its entities | |
| 14 | not be used in advertising or publicity pertaining to | |
| 15 | distribution of the software without specific, written prior | |
| 16 | permission. | |
| 17 | ||
| 18 | LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, | |
| 19 | INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. | |
| 20 | IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY | |
| 21 | SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES | |
| 22 | WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER | |
| 23 | IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, | |
| 24 | ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF | |
| 25 | THIS SOFTWARE. | |
| 26 | ||
| 27 | ****************************************************************/ | |
| 28 | ||
| 3d9156a7 A |
29 | /* Please send bug reports to David M. Gay (dmg at acm dot org, |
| 30 | * with " at " changed at "@" and " dot " changed to "."). */ | |
| 9385eb3d A |
31 | |
| 32 | #include "gdtoaimp.h" | |
| eb1cde05 A |
33 | #ifndef NO_FENV_H |
| 34 | #include <fenv.h> | |
| 35 | #endif | |
| 9385eb3d A |
36 | |
| 37 | #ifdef USE_LOCALE | |
| 38 | #include "locale.h" | |
| 39 | #endif | |
| 40 | ||
| 41 | #ifdef IEEE_Arith | |
| 42 | #ifndef NO_IEEE_Scale | |
| 43 | #define Avoid_Underflow | |
| 44 | #undef tinytens | |
| 45 | /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ | |
| 46 | /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ | |
| 47 | static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, | |
| 48 | 9007199254740992.e-256 | |
| 49 | }; | |
| 50 | #endif | |
| 51 | #endif | |
| 52 | ||
| 53 | #ifdef Honor_FLT_ROUNDS | |
| 54 | #define Rounding rounding | |
| 55 | #undef Check_FLT_ROUNDS | |
| 56 | #define Check_FLT_ROUNDS | |
| 57 | #else | |
| 58 | #define Rounding Flt_Rounds | |
| 59 | #endif | |
| 60 | ||
| 61 | double | |
| 62 | strtod | |
| 63 | #ifdef KR_headers | |
| 64 | (s00, se) CONST char *s00; char **se; | |
| 65 | #else | |
| 66 | (CONST char *s00, char **se) | |
| 67 | #endif | |
| 68 | { | |
| 69 | #ifdef Avoid_Underflow | |
| 70 | int scale; | |
| 71 | #endif | |
| 3d9156a7 | 72 | int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, |
| 9385eb3d A |
73 | e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; |
| 74 | CONST char *s, *s0, *s1; | |
| 75 | double aadj, aadj1, adj, rv, rv0; | |
| 76 | Long L; | |
| 77 | ULong y, z; | |
| 78 | Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; | |
| 79 | #ifdef SET_INEXACT | |
| 80 | int inexact, oldinexact; | |
| 81 | #endif | |
| 82 | #ifdef Honor_FLT_ROUNDS | |
| 83 | int rounding; | |
| 84 | #endif | |
| 85 | ||
| 3d9156a7 | 86 | sign = nz0 = nz = decpt = 0; |
| 9385eb3d A |
87 | dval(rv) = 0.; |
| 88 | for(s = s00;;s++) switch(*s) { | |
| 89 | case '-': | |
| 90 | sign = 1; | |
| 91 | /* no break */ | |
| 92 | case '+': | |
| 93 | if (*++s) | |
| 94 | goto break2; | |
| 95 | /* no break */ | |
| 96 | case 0: | |
| 97 | goto ret0; | |
| 98 | case '\t': | |
| 99 | case '\n': | |
| 100 | case '\v': | |
| 101 | case '\f': | |
| 102 | case '\r': | |
| 103 | case ' ': | |
| 104 | continue; | |
| 105 | default: | |
| 106 | goto break2; | |
| 107 | } | |
| 108 | break2: | |
| 109 | if (*s == '0') { | |
| 110 | #ifndef NO_HEX_FP | |
| 111 | { | |
| 112 | static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; | |
| 113 | Long exp; | |
| 114 | ULong bits[2]; | |
| 115 | switch(s[1]) { | |
| 116 | case 'x': | |
| 117 | case 'X': | |
| eb1cde05 A |
118 | { |
| 119 | #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) | |
| 120 | FPI fpi1 = fpi; | |
| 121 | switch(fegetround()) { | |
| 122 | case FE_TOWARDZERO: fpi1.rounding = 0; break; | |
| 123 | case FE_UPWARD: fpi1.rounding = 2; break; | |
| 124 | case FE_DOWNWARD: fpi1.rounding = 3; | |
| 125 | } | |
| 126 | #else | |
| 127 | #define fpi1 fpi | |
| 128 | #endif | |
| 129 | switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { | |
| 9385eb3d A |
130 | case STRTOG_NoNumber: |
| 131 | s = s00; | |
| 132 | sign = 0; | |
| 133 | case STRTOG_Zero: | |
| 134 | break; | |
| 135 | default: | |
| 59e0d9fe A |
136 | if (bb) { |
| 137 | copybits(bits, fpi.nbits, bb); | |
| 138 | Bfree(bb); | |
| 139 | } | |
| 9385eb3d | 140 | ULtod(((U*)&rv)->L, bits, exp, i); |
| eb1cde05 | 141 | }} |
| 9385eb3d A |
142 | goto ret; |
| 143 | } | |
| 144 | } | |
| 145 | #endif | |
| 146 | nz0 = 1; | |
| 147 | while(*++s == '0') ; | |
| 148 | if (!*s) | |
| 149 | goto ret; | |
| 150 | } | |
| 151 | s0 = s; | |
| 152 | y = z = 0; | |
| 153 | for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) | |
| 154 | if (nd < 9) | |
| 155 | y = 10*y + c - '0'; | |
| 156 | else if (nd < 16) | |
| 157 | z = 10*z + c - '0'; | |
| 158 | nd0 = nd; | |
| 159 | #ifdef USE_LOCALE | |
| 160 | if (c == *localeconv()->decimal_point) | |
| 161 | #else | |
| 162 | if (c == '.') | |
| 163 | #endif | |
| 164 | { | |
| 3d9156a7 | 165 | decpt = 1; |
| 9385eb3d A |
166 | c = *++s; |
| 167 | if (!nd) { | |
| 168 | for(; c == '0'; c = *++s) | |
| 169 | nz++; | |
| 170 | if (c > '0' && c <= '9') { | |
| 171 | s0 = s; | |
| 172 | nf += nz; | |
| 173 | nz = 0; | |
| 174 | goto have_dig; | |
| 175 | } | |
| 176 | goto dig_done; | |
| 177 | } | |
| 178 | for(; c >= '0' && c <= '9'; c = *++s) { | |
| 179 | have_dig: | |
| 180 | nz++; | |
| 181 | if (c -= '0') { | |
| 182 | nf += nz; | |
| 183 | for(i = 1; i < nz; i++) | |
| 184 | if (nd++ < 9) | |
| 185 | y *= 10; | |
| 186 | else if (nd <= DBL_DIG + 1) | |
| 187 | z *= 10; | |
| 188 | if (nd++ < 9) | |
| 189 | y = 10*y + c; | |
| 190 | else if (nd <= DBL_DIG + 1) | |
| 191 | z = 10*z + c; | |
| 192 | nz = 0; | |
| 193 | } | |
| 194 | } | |
| 195 | } | |
| 196 | dig_done: | |
| 197 | e = 0; | |
| 198 | if (c == 'e' || c == 'E') { | |
| 199 | if (!nd && !nz && !nz0) { | |
| 200 | goto ret0; | |
| 201 | } | |
| 202 | s00 = s; | |
| 203 | esign = 0; | |
| 204 | switch(c = *++s) { | |
| 205 | case '-': | |
| 206 | esign = 1; | |
| 207 | case '+': | |
| 208 | c = *++s; | |
| 209 | } | |
| 210 | if (c >= '0' && c <= '9') { | |
| 211 | while(c == '0') | |
| 212 | c = *++s; | |
| 213 | if (c > '0' && c <= '9') { | |
| 214 | L = c - '0'; | |
| 215 | s1 = s; | |
| 216 | while((c = *++s) >= '0' && c <= '9') | |
| 217 | L = 10*L + c - '0'; | |
| 218 | if (s - s1 > 8 || L > 19999) | |
| 219 | /* Avoid confusion from exponents | |
| 220 | * so large that e might overflow. | |
| 221 | */ | |
| 222 | e = 19999; /* safe for 16 bit ints */ | |
| 223 | else | |
| 224 | e = (int)L; | |
| 225 | if (esign) | |
| 226 | e = -e; | |
| 227 | } | |
| 228 | else | |
| 229 | e = 0; | |
| 230 | } | |
| 231 | else | |
| 232 | s = s00; | |
| 233 | } | |
| 234 | if (!nd) { | |
| 235 | if (!nz && !nz0) { | |
| 236 | #ifdef INFNAN_CHECK | |
| 237 | /* Check for Nan and Infinity */ | |
| 238 | ULong bits[2]; | |
| 239 | static FPI fpinan = /* only 52 explicit bits */ | |
| 240 | { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; | |
| 3d9156a7 A |
241 | if (!decpt) |
| 242 | switch(c) { | |
| 9385eb3d A |
243 | case 'i': |
| 244 | case 'I': | |
| 245 | if (match(&s,"nf")) { | |
| 246 | --s; | |
| 247 | if (!match(&s,"inity")) | |
| 248 | ++s; | |
| 249 | word0(rv) = 0x7ff00000; | |
| 250 | word1(rv) = 0; | |
| 251 | goto ret; | |
| 252 | } | |
| 253 | break; | |
| 254 | case 'n': | |
| 255 | case 'N': | |
| 256 | if (match(&s, "an")) { | |
| 257 | #ifndef No_Hex_NaN | |
| 258 | if (*s == '(' /*)*/ | |
| 259 | && hexnan(&s, &fpinan, bits) | |
| 260 | == STRTOG_NaNbits) { | |
| 261 | word0(rv) = 0x7ff00000 | bits[1]; | |
| 262 | word1(rv) = bits[0]; | |
| 263 | } | |
| 264 | else { | |
| 3d9156a7 | 265 | #endif |
| 9385eb3d A |
266 | word0(rv) = NAN_WORD0; |
| 267 | word1(rv) = NAN_WORD1; | |
| 3d9156a7 | 268 | #ifndef No_Hex_NaN |
| 9385eb3d A |
269 | } |
| 270 | #endif | |
| 271 | goto ret; | |
| 272 | } | |
| 273 | } | |
| 274 | #endif /* INFNAN_CHECK */ | |
| 275 | ret0: | |
| 276 | s = s00; | |
| 277 | sign = 0; | |
| 278 | } | |
| 279 | goto ret; | |
| 280 | } | |
| 281 | e1 = e -= nf; | |
| 282 | ||
| 283 | /* Now we have nd0 digits, starting at s0, followed by a | |
| 284 | * decimal point, followed by nd-nd0 digits. The number we're | |
| 285 | * after is the integer represented by those digits times | |
| 286 | * 10**e */ | |
| 287 | ||
| 288 | if (!nd0) | |
| 289 | nd0 = nd; | |
| 290 | k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; | |
| 291 | dval(rv) = y; | |
| 292 | if (k > 9) { | |
| 293 | #ifdef SET_INEXACT | |
| 294 | if (k > DBL_DIG) | |
| 295 | oldinexact = get_inexact(); | |
| 296 | #endif | |
| 297 | dval(rv) = tens[k - 9] * dval(rv) + z; | |
| 298 | } | |
| 299 | bd0 = 0; | |
| 300 | if (nd <= DBL_DIG | |
| 301 | #ifndef RND_PRODQUOT | |
| 302 | #ifndef Honor_FLT_ROUNDS | |
| 303 | && Flt_Rounds == 1 | |
| 304 | #endif | |
| 305 | #endif | |
| 306 | ) { | |
| 307 | if (!e) | |
| 308 | goto ret; | |
| 309 | if (e > 0) { | |
| 310 | if (e <= Ten_pmax) { | |
| 311 | #ifdef VAX | |
| 312 | goto vax_ovfl_check; | |
| 313 | #else | |
| 314 | #ifdef Honor_FLT_ROUNDS | |
| 315 | /* round correctly FLT_ROUNDS = 2 or 3 */ | |
| 316 | if (sign) { | |
| 317 | rv = -rv; | |
| 318 | sign = 0; | |
| 319 | } | |
| 320 | #endif | |
| 321 | /* rv = */ rounded_product(dval(rv), tens[e]); | |
| 322 | goto ret; | |
| 323 | #endif | |
| 324 | } | |
| 325 | i = DBL_DIG - nd; | |
| 326 | if (e <= Ten_pmax + i) { | |
| 327 | /* A fancier test would sometimes let us do | |
| 328 | * this for larger i values. | |
| 329 | */ | |
| 330 | #ifdef Honor_FLT_ROUNDS | |
| 331 | /* round correctly FLT_ROUNDS = 2 or 3 */ | |
| 332 | if (sign) { | |
| 333 | rv = -rv; | |
| 334 | sign = 0; | |
| 335 | } | |
| 336 | #endif | |
| 337 | e -= i; | |
| 338 | dval(rv) *= tens[i]; | |
| 339 | #ifdef VAX | |
| 340 | /* VAX exponent range is so narrow we must | |
| 341 | * worry about overflow here... | |
| 342 | */ | |
| 343 | vax_ovfl_check: | |
| 344 | word0(rv) -= P*Exp_msk1; | |
| 345 | /* rv = */ rounded_product(dval(rv), tens[e]); | |
| 346 | if ((word0(rv) & Exp_mask) | |
| 347 | > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) | |
| 348 | goto ovfl; | |
| 349 | word0(rv) += P*Exp_msk1; | |
| 350 | #else | |
| 351 | /* rv = */ rounded_product(dval(rv), tens[e]); | |
| 352 | #endif | |
| 353 | goto ret; | |
| 354 | } | |
| 355 | } | |
| 356 | #ifndef Inaccurate_Divide | |
| 357 | else if (e >= -Ten_pmax) { | |
| 358 | #ifdef Honor_FLT_ROUNDS | |
| 359 | /* round correctly FLT_ROUNDS = 2 or 3 */ | |
| 360 | if (sign) { | |
| 361 | rv = -rv; | |
| 362 | sign = 0; | |
| 363 | } | |
| 364 | #endif | |
| 365 | /* rv = */ rounded_quotient(dval(rv), tens[-e]); | |
| 366 | goto ret; | |
| 367 | } | |
| 368 | #endif | |
| 369 | } | |
| 370 | e1 += nd - k; | |
| 371 | ||
| 372 | #ifdef IEEE_Arith | |
| 373 | #ifdef SET_INEXACT | |
| 374 | inexact = 1; | |
| 375 | if (k <= DBL_DIG) | |
| 376 | oldinexact = get_inexact(); | |
| 377 | #endif | |
| 378 | #ifdef Avoid_Underflow | |
| 379 | scale = 0; | |
| 380 | #endif | |
| 381 | #ifdef Honor_FLT_ROUNDS | |
| 382 | if ((rounding = Flt_Rounds) >= 2) { | |
| 383 | if (sign) | |
| 384 | rounding = rounding == 2 ? 0 : 2; | |
| 385 | else | |
| 386 | if (rounding != 2) | |
| 387 | rounding = 0; | |
| 388 | } | |
| 389 | #endif | |
| 390 | #endif /*IEEE_Arith*/ | |
| 391 | ||
| 392 | /* Get starting approximation = rv * 10**e1 */ | |
| 393 | ||
| 394 | if (e1 > 0) { | |
| 395 | if ( (i = e1 & 15) !=0) | |
| 396 | dval(rv) *= tens[i]; | |
| 397 | if (e1 &= ~15) { | |
| 398 | if (e1 > DBL_MAX_10_EXP) { | |
| 399 | ovfl: | |
| 400 | #ifndef NO_ERRNO | |
| 401 | errno = ERANGE; | |
| 402 | #endif | |
| 403 | /* Can't trust HUGE_VAL */ | |
| 404 | #ifdef IEEE_Arith | |
| 405 | #ifdef Honor_FLT_ROUNDS | |
| 406 | switch(rounding) { | |
| 407 | case 0: /* toward 0 */ | |
| 408 | case 3: /* toward -infinity */ | |
| 409 | word0(rv) = Big0; | |
| 410 | word1(rv) = Big1; | |
| 411 | break; | |
| 412 | default: | |
| 413 | word0(rv) = Exp_mask; | |
| 414 | word1(rv) = 0; | |
| 415 | } | |
| 416 | #else /*Honor_FLT_ROUNDS*/ | |
| 417 | word0(rv) = Exp_mask; | |
| 418 | word1(rv) = 0; | |
| 419 | #endif /*Honor_FLT_ROUNDS*/ | |
| 420 | #ifdef SET_INEXACT | |
| 421 | /* set overflow bit */ | |
| 422 | dval(rv0) = 1e300; | |
| 423 | dval(rv0) *= dval(rv0); | |
| 424 | #endif | |
| 425 | #else /*IEEE_Arith*/ | |
| 426 | word0(rv) = Big0; | |
| 427 | word1(rv) = Big1; | |
| 428 | #endif /*IEEE_Arith*/ | |
| 429 | if (bd0) | |
| 430 | goto retfree; | |
| 431 | goto ret; | |
| 432 | } | |
| 433 | e1 >>= 4; | |
| 434 | for(j = 0; e1 > 1; j++, e1 >>= 1) | |
| 435 | if (e1 & 1) | |
| 436 | dval(rv) *= bigtens[j]; | |
| 437 | /* The last multiplication could overflow. */ | |
| 438 | word0(rv) -= P*Exp_msk1; | |
| 439 | dval(rv) *= bigtens[j]; | |
| 440 | if ((z = word0(rv) & Exp_mask) | |
| 441 | > Exp_msk1*(DBL_MAX_EXP+Bias-P)) | |
| 442 | goto ovfl; | |
| 443 | if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { | |
| 444 | /* set to largest number */ | |
| 445 | /* (Can't trust DBL_MAX) */ | |
| 446 | word0(rv) = Big0; | |
| 447 | word1(rv) = Big1; | |
| 448 | } | |
| 449 | else | |
| 450 | word0(rv) += P*Exp_msk1; | |
| 451 | } | |
| 452 | } | |
| 453 | else if (e1 < 0) { | |
| 454 | e1 = -e1; | |
| 455 | if ( (i = e1 & 15) !=0) | |
| 456 | dval(rv) /= tens[i]; | |
| 457 | if (e1 >>= 4) { | |
| 458 | if (e1 >= 1 << n_bigtens) | |
| 459 | goto undfl; | |
| 460 | #ifdef Avoid_Underflow | |
| 461 | if (e1 & Scale_Bit) | |
| 462 | scale = 2*P; | |
| 463 | for(j = 0; e1 > 0; j++, e1 >>= 1) | |
| 464 | if (e1 & 1) | |
| 465 | dval(rv) *= tinytens[j]; | |
| 466 | if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) | |
| 467 | >> Exp_shift)) > 0) { | |
| 468 | /* scaled rv is denormal; zap j low bits */ | |
| 469 | if (j >= 32) { | |
| 470 | word1(rv) = 0; | |
| 471 | if (j >= 53) | |
| 472 | word0(rv) = (P+2)*Exp_msk1; | |
| 473 | else | |
| 474 | word0(rv) &= 0xffffffff << j-32; | |
| 475 | } | |
| 476 | else | |
| 477 | word1(rv) &= 0xffffffff << j; | |
| 478 | } | |
| 479 | #else | |
| 480 | for(j = 0; e1 > 1; j++, e1 >>= 1) | |
| 481 | if (e1 & 1) | |
| 482 | dval(rv) *= tinytens[j]; | |
| 483 | /* The last multiplication could underflow. */ | |
| 484 | dval(rv0) = dval(rv); | |
| 485 | dval(rv) *= tinytens[j]; | |
| 486 | if (!dval(rv)) { | |
| 487 | dval(rv) = 2.*dval(rv0); | |
| 488 | dval(rv) *= tinytens[j]; | |
| 489 | #endif | |
| 490 | if (!dval(rv)) { | |
| 491 | undfl: | |
| 492 | dval(rv) = 0.; | |
| 493 | #ifndef NO_ERRNO | |
| 494 | errno = ERANGE; | |
| 495 | #endif | |
| 496 | if (bd0) | |
| 497 | goto retfree; | |
| 498 | goto ret; | |
| 499 | } | |
| 500 | #ifndef Avoid_Underflow | |
| 501 | word0(rv) = Tiny0; | |
| 502 | word1(rv) = Tiny1; | |
| 503 | /* The refinement below will clean | |
| 504 | * this approximation up. | |
| 505 | */ | |
| 506 | } | |
| 507 | #endif | |
| 508 | } | |
| 509 | } | |
| 510 | ||
| 511 | /* Now the hard part -- adjusting rv to the correct value.*/ | |
| 512 | ||
| 513 | /* Put digits into bd: true value = bd * 10^e */ | |
| 514 | ||
| 515 | bd0 = s2b(s0, nd0, nd, y); | |
| 516 | ||
| 517 | for(;;) { | |
| 518 | bd = Balloc(bd0->k); | |
| 519 | Bcopy(bd, bd0); | |
| 520 | bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ | |
| 521 | bs = i2b(1); | |
| 522 | ||
| 523 | if (e >= 0) { | |
| 524 | bb2 = bb5 = 0; | |
| 525 | bd2 = bd5 = e; | |
| 526 | } | |
| 527 | else { | |
| 528 | bb2 = bb5 = -e; | |
| 529 | bd2 = bd5 = 0; | |
| 530 | } | |
| 531 | if (bbe >= 0) | |
| 532 | bb2 += bbe; | |
| 533 | else | |
| 534 | bd2 -= bbe; | |
| 535 | bs2 = bb2; | |
| 536 | #ifdef Honor_FLT_ROUNDS | |
| 537 | if (rounding != 1) | |
| 538 | bs2++; | |
| 539 | #endif | |
| 540 | #ifdef Avoid_Underflow | |
| 541 | j = bbe - scale; | |
| 542 | i = j + bbbits - 1; /* logb(rv) */ | |
| 543 | if (i < Emin) /* denormal */ | |
| 544 | j += P - Emin; | |
| 545 | else | |
| 546 | j = P + 1 - bbbits; | |
| 547 | #else /*Avoid_Underflow*/ | |
| 548 | #ifdef Sudden_Underflow | |
| 549 | #ifdef IBM | |
| 550 | j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); | |
| 551 | #else | |
| 552 | j = P + 1 - bbbits; | |
| 553 | #endif | |
| 554 | #else /*Sudden_Underflow*/ | |
| 555 | j = bbe; | |
| 556 | i = j + bbbits - 1; /* logb(rv) */ | |
| 557 | if (i < Emin) /* denormal */ | |
| 558 | j += P - Emin; | |
| 559 | else | |
| 560 | j = P + 1 - bbbits; | |
| 561 | #endif /*Sudden_Underflow*/ | |
| 562 | #endif /*Avoid_Underflow*/ | |
| 563 | bb2 += j; | |
| 564 | bd2 += j; | |
| 565 | #ifdef Avoid_Underflow | |
| 566 | bd2 += scale; | |
| 567 | #endif | |
| 568 | i = bb2 < bd2 ? bb2 : bd2; | |
| 569 | if (i > bs2) | |
| 570 | i = bs2; | |
| 571 | if (i > 0) { | |
| 572 | bb2 -= i; | |
| 573 | bd2 -= i; | |
| 574 | bs2 -= i; | |
| 575 | } | |
| 576 | if (bb5 > 0) { | |
| 577 | bs = pow5mult(bs, bb5); | |
| 578 | bb1 = mult(bs, bb); | |
| 579 | Bfree(bb); | |
| 580 | bb = bb1; | |
| 581 | } | |
| 582 | if (bb2 > 0) | |
| 583 | bb = lshift(bb, bb2); | |
| 584 | if (bd5 > 0) | |
| 585 | bd = pow5mult(bd, bd5); | |
| 586 | if (bd2 > 0) | |
| 587 | bd = lshift(bd, bd2); | |
| 588 | if (bs2 > 0) | |
| 589 | bs = lshift(bs, bs2); | |
| 590 | delta = diff(bb, bd); | |
| 591 | dsign = delta->sign; | |
| 592 | delta->sign = 0; | |
| 593 | i = cmp(delta, bs); | |
| 594 | #ifdef Honor_FLT_ROUNDS | |
| 595 | if (rounding != 1) { | |
| 596 | if (i < 0) { | |
| 597 | /* Error is less than an ulp */ | |
| 598 | if (!delta->x[0] && delta->wds <= 1) { | |
| 599 | /* exact */ | |
| 600 | #ifdef SET_INEXACT | |
| 601 | inexact = 0; | |
| 602 | #endif | |
| 603 | break; | |
| 604 | } | |
| 605 | if (rounding) { | |
| 606 | if (dsign) { | |
| 607 | adj = 1.; | |
| 608 | goto apply_adj; | |
| 609 | } | |
| 610 | } | |
| 611 | else if (!dsign) { | |
| 612 | adj = -1.; | |
| 613 | if (!word1(rv) | |
| 614 | && !(word0(rv) & Frac_mask)) { | |
| 615 | y = word0(rv) & Exp_mask; | |
| 616 | #ifdef Avoid_Underflow | |
| 617 | if (!scale || y > 2*P*Exp_msk1) | |
| 618 | #else | |
| 619 | if (y) | |
| 620 | #endif | |
| 621 | { | |
| 622 | delta = lshift(delta,Log2P); | |
| 623 | if (cmp(delta, bs) <= 0) | |
| 624 | adj = -0.5; | |
| 625 | } | |
| 626 | } | |
| 627 | apply_adj: | |
| 628 | #ifdef Avoid_Underflow | |
| 629 | if (scale && (y = word0(rv) & Exp_mask) | |
| 630 | <= 2*P*Exp_msk1) | |
| 631 | word0(adj) += (2*P+1)*Exp_msk1 - y; | |
| 632 | #else | |
| 633 | #ifdef Sudden_Underflow | |
| 634 | if ((word0(rv) & Exp_mask) <= | |
| 635 | P*Exp_msk1) { | |
| 636 | word0(rv) += P*Exp_msk1; | |
| 637 | dval(rv) += adj*ulp(dval(rv)); | |
| 638 | word0(rv) -= P*Exp_msk1; | |
| 639 | } | |
| 640 | else | |
| 641 | #endif /*Sudden_Underflow*/ | |
| 642 | #endif /*Avoid_Underflow*/ | |
| 643 | dval(rv) += adj*ulp(dval(rv)); | |
| 644 | } | |
| 645 | break; | |
| 646 | } | |
| 647 | adj = ratio(delta, bs); | |
| 648 | if (adj < 1.) | |
| 649 | adj = 1.; | |
| 650 | if (adj <= 0x7ffffffe) { | |
| 651 | /* adj = rounding ? ceil(adj) : floor(adj); */ | |
| 652 | y = adj; | |
| 653 | if (y != adj) { | |
| 654 | if (!((rounding>>1) ^ dsign)) | |
| 655 | y++; | |
| 656 | adj = y; | |
| 657 | } | |
| 658 | } | |
| 659 | #ifdef Avoid_Underflow | |
| 660 | if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) | |
| 661 | word0(adj) += (2*P+1)*Exp_msk1 - y; | |
| 662 | #else | |
| 663 | #ifdef Sudden_Underflow | |
| 664 | if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { | |
| 665 | word0(rv) += P*Exp_msk1; | |
| 666 | adj *= ulp(dval(rv)); | |
| 667 | if (dsign) | |
| 668 | dval(rv) += adj; | |
| 669 | else | |
| 670 | dval(rv) -= adj; | |
| 671 | word0(rv) -= P*Exp_msk1; | |
| 672 | goto cont; | |
| 673 | } | |
| 674 | #endif /*Sudden_Underflow*/ | |
| 675 | #endif /*Avoid_Underflow*/ | |
| 676 | adj *= ulp(dval(rv)); | |
| 677 | if (dsign) | |
| 678 | dval(rv) += adj; | |
| 679 | else | |
| 680 | dval(rv) -= adj; | |
| 681 | goto cont; | |
| 682 | } | |
| 683 | #endif /*Honor_FLT_ROUNDS*/ | |
| 684 | ||
| 685 | if (i < 0) { | |
| 686 | /* Error is less than half an ulp -- check for | |
| 687 | * special case of mantissa a power of two. | |
| 688 | */ | |
| 689 | if (dsign || word1(rv) || word0(rv) & Bndry_mask | |
| 690 | #ifdef IEEE_Arith | |
| 691 | #ifdef Avoid_Underflow | |
| 692 | || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 | |
| 693 | #else | |
| 694 | || (word0(rv) & Exp_mask) <= Exp_msk1 | |
| 695 | #endif | |
| 696 | #endif | |
| 697 | ) { | |
| 698 | #ifdef SET_INEXACT | |
| 699 | if (!delta->x[0] && delta->wds <= 1) | |
| 700 | inexact = 0; | |
| 701 | #endif | |
| 702 | break; | |
| 703 | } | |
| 704 | if (!delta->x[0] && delta->wds <= 1) { | |
| 705 | /* exact result */ | |
| 706 | #ifdef SET_INEXACT | |
| 707 | inexact = 0; | |
| 708 | #endif | |
| 709 | break; | |
| 710 | } | |
| 711 | delta = lshift(delta,Log2P); | |
| 712 | if (cmp(delta, bs) > 0) | |
| 713 | goto drop_down; | |
| 714 | break; | |
| 715 | } | |
| 716 | if (i == 0) { | |
| 717 | /* exactly half-way between */ | |
| 718 | if (dsign) { | |
| 719 | if ((word0(rv) & Bndry_mask1) == Bndry_mask1 | |
| 720 | && word1(rv) == ( | |
| 721 | #ifdef Avoid_Underflow | |
| 722 | (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) | |
| 723 | ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : | |
| 724 | #endif | |
| 725 | 0xffffffff)) { | |
| 726 | /*boundary case -- increment exponent*/ | |
| 727 | word0(rv) = (word0(rv) & Exp_mask) | |
| 728 | + Exp_msk1 | |
| 729 | #ifdef IBM | |
| 730 | | Exp_msk1 >> 4 | |
| 731 | #endif | |
| 732 | ; | |
| 733 | word1(rv) = 0; | |
| 734 | #ifdef Avoid_Underflow | |
| 735 | dsign = 0; | |
| 736 | #endif | |
| 737 | break; | |
| 738 | } | |
| 739 | } | |
| 740 | else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { | |
| 741 | drop_down: | |
| 742 | /* boundary case -- decrement exponent */ | |
| 743 | #ifdef Sudden_Underflow /*{{*/ | |
| 744 | L = word0(rv) & Exp_mask; | |
| 745 | #ifdef IBM | |
| 746 | if (L < Exp_msk1) | |
| 747 | #else | |
| 748 | #ifdef Avoid_Underflow | |
| 749 | if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) | |
| 750 | #else | |
| 751 | if (L <= Exp_msk1) | |
| 752 | #endif /*Avoid_Underflow*/ | |
| 753 | #endif /*IBM*/ | |
| 754 | goto undfl; | |
| 755 | L -= Exp_msk1; | |
| 756 | #else /*Sudden_Underflow}{*/ | |
| 757 | #ifdef Avoid_Underflow | |
| 758 | if (scale) { | |
| 759 | L = word0(rv) & Exp_mask; | |
| 760 | if (L <= (2*P+1)*Exp_msk1) { | |
| 761 | if (L > (P+2)*Exp_msk1) | |
| 762 | /* round even ==> */ | |
| 763 | /* accept rv */ | |
| 764 | break; | |
| 765 | /* rv = smallest denormal */ | |
| 766 | goto undfl; | |
| 767 | } | |
| 768 | } | |
| 769 | #endif /*Avoid_Underflow*/ | |
| 770 | L = (word0(rv) & Exp_mask) - Exp_msk1; | |
| 771 | #endif /*Sudden_Underflow}*/ | |
| 772 | word0(rv) = L | Bndry_mask1; | |
| 773 | word1(rv) = 0xffffffff; | |
| 774 | #ifdef IBM | |
| 775 | goto cont; | |
| 776 | #else | |
| 777 | break; | |
| 778 | #endif | |
| 779 | } | |
| 780 | #ifndef ROUND_BIASED | |
| 781 | if (!(word1(rv) & LSB)) | |
| 782 | break; | |
| 783 | #endif | |
| 784 | if (dsign) | |
| 785 | dval(rv) += ulp(dval(rv)); | |
| 786 | #ifndef ROUND_BIASED | |
| 787 | else { | |
| 788 | dval(rv) -= ulp(dval(rv)); | |
| 789 | #ifndef Sudden_Underflow | |
| 790 | if (!dval(rv)) | |
| 791 | goto undfl; | |
| 792 | #endif | |
| 793 | } | |
| 794 | #ifdef Avoid_Underflow | |
| 795 | dsign = 1 - dsign; | |
| 796 | #endif | |
| 797 | #endif | |
| 798 | break; | |
| 799 | } | |
| 800 | if ((aadj = ratio(delta, bs)) <= 2.) { | |
| 801 | if (dsign) | |
| 802 | aadj = aadj1 = 1.; | |
| 803 | else if (word1(rv) || word0(rv) & Bndry_mask) { | |
| 804 | #ifndef Sudden_Underflow | |
| 805 | if (word1(rv) == Tiny1 && !word0(rv)) | |
| 806 | goto undfl; | |
| 807 | #endif | |
| 808 | aadj = 1.; | |
| 809 | aadj1 = -1.; | |
| 810 | } | |
| 811 | else { | |
| 812 | /* special case -- power of FLT_RADIX to be */ | |
| 813 | /* rounded down... */ | |
| 814 | ||
| 815 | if (aadj < 2./FLT_RADIX) | |
| 816 | aadj = 1./FLT_RADIX; | |
| 817 | else | |
| 818 | aadj *= 0.5; | |
| 819 | aadj1 = -aadj; | |
| 820 | } | |
| 821 | } | |
| 822 | else { | |
| 823 | aadj *= 0.5; | |
| 824 | aadj1 = dsign ? aadj : -aadj; | |
| 825 | #ifdef Check_FLT_ROUNDS | |
| 826 | switch(Rounding) { | |
| 827 | case 2: /* towards +infinity */ | |
| 828 | aadj1 -= 0.5; | |
| 829 | break; | |
| 830 | case 0: /* towards 0 */ | |
| 831 | case 3: /* towards -infinity */ | |
| 832 | aadj1 += 0.5; | |
| 833 | } | |
| 834 | #else | |
| 835 | if (Flt_Rounds == 0) | |
| 836 | aadj1 += 0.5; | |
| 837 | #endif /*Check_FLT_ROUNDS*/ | |
| 838 | } | |
| 839 | y = word0(rv) & Exp_mask; | |
| 840 | ||
| 841 | /* Check for overflow */ | |
| 842 | ||
| 843 | if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { | |
| 844 | dval(rv0) = dval(rv); | |
| 845 | word0(rv) -= P*Exp_msk1; | |
| 846 | adj = aadj1 * ulp(dval(rv)); | |
| 847 | dval(rv) += adj; | |
| 848 | if ((word0(rv) & Exp_mask) >= | |
| 849 | Exp_msk1*(DBL_MAX_EXP+Bias-P)) { | |
| 850 | if (word0(rv0) == Big0 && word1(rv0) == Big1) | |
| 851 | goto ovfl; | |
| 852 | word0(rv) = Big0; | |
| 853 | word1(rv) = Big1; | |
| 854 | goto cont; | |
| 855 | } | |
| 856 | else | |
| 857 | word0(rv) += P*Exp_msk1; | |
| 858 | } | |
| 859 | else { | |
| 860 | #ifdef Avoid_Underflow | |
| 861 | if (scale && y <= 2*P*Exp_msk1) { | |
| 862 | if (aadj <= 0x7fffffff) { | |
| 863 | if ((z = aadj) <= 0) | |
| 864 | z = 1; | |
| 865 | aadj = z; | |
| 866 | aadj1 = dsign ? aadj : -aadj; | |
| 867 | } | |
| 868 | word0(aadj1) += (2*P+1)*Exp_msk1 - y; | |
| 869 | } | |
| 870 | adj = aadj1 * ulp(dval(rv)); | |
| 871 | dval(rv) += adj; | |
| 872 | #else | |
| 873 | #ifdef Sudden_Underflow | |
| 874 | if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { | |
| 875 | dval(rv0) = dval(rv); | |
| 876 | word0(rv) += P*Exp_msk1; | |
| 877 | adj = aadj1 * ulp(dval(rv)); | |
| 878 | dval(rv) += adj; | |
| 879 | #ifdef IBM | |
| 880 | if ((word0(rv) & Exp_mask) < P*Exp_msk1) | |
| 881 | #else | |
| 882 | if ((word0(rv) & Exp_mask) <= P*Exp_msk1) | |
| 883 | #endif | |
| 884 | { | |
| 885 | if (word0(rv0) == Tiny0 | |
| 886 | && word1(rv0) == Tiny1) | |
| 887 | goto undfl; | |
| 888 | word0(rv) = Tiny0; | |
| 889 | word1(rv) = Tiny1; | |
| 890 | goto cont; | |
| 891 | } | |
| 892 | else | |
| 893 | word0(rv) -= P*Exp_msk1; | |
| 894 | } | |
| 895 | else { | |
| 896 | adj = aadj1 * ulp(dval(rv)); | |
| 897 | dval(rv) += adj; | |
| 898 | } | |
| 899 | #else /*Sudden_Underflow*/ | |
| 900 | /* Compute adj so that the IEEE rounding rules will | |
| 901 | * correctly round rv + adj in some half-way cases. | |
| 902 | * If rv * ulp(rv) is denormalized (i.e., | |
| 903 | * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid | |
| 904 | * trouble from bits lost to denormalization; | |
| 905 | * example: 1.2e-307 . | |
| 906 | */ | |
| 907 | if (y <= (P-1)*Exp_msk1 && aadj > 1.) { | |
| 908 | aadj1 = (double)(int)(aadj + 0.5); | |
| 909 | if (!dsign) | |
| 910 | aadj1 = -aadj1; | |
| 911 | } | |
| 912 | adj = aadj1 * ulp(dval(rv)); | |
| 913 | dval(rv) += adj; | |
| 914 | #endif /*Sudden_Underflow*/ | |
| 915 | #endif /*Avoid_Underflow*/ | |
| 916 | } | |
| 917 | z = word0(rv) & Exp_mask; | |
| 918 | #ifndef SET_INEXACT | |
| 919 | #ifdef Avoid_Underflow | |
| 920 | if (!scale) | |
| 921 | #endif | |
| 922 | if (y == z) { | |
| 923 | /* Can we stop now? */ | |
| 924 | L = (Long)aadj; | |
| 925 | aadj -= L; | |
| 926 | /* The tolerances below are conservative. */ | |
| 927 | if (dsign || word1(rv) || word0(rv) & Bndry_mask) { | |
| 928 | if (aadj < .4999999 || aadj > .5000001) | |
| 929 | break; | |
| 930 | } | |
| 931 | else if (aadj < .4999999/FLT_RADIX) | |
| 932 | break; | |
| 933 | } | |
| 934 | #endif | |
| 935 | cont: | |
| 936 | Bfree(bb); | |
| 937 | Bfree(bd); | |
| 938 | Bfree(bs); | |
| 939 | Bfree(delta); | |
| 940 | } | |
| 941 | #ifdef SET_INEXACT | |
| 942 | if (inexact) { | |
| 943 | if (!oldinexact) { | |
| 944 | word0(rv0) = Exp_1 + (70 << Exp_shift); | |
| 945 | word1(rv0) = 0; | |
| 946 | dval(rv0) += 1.; | |
| 947 | } | |
| 948 | } | |
| 949 | else if (!oldinexact) | |
| 950 | clear_inexact(); | |
| 951 | #endif | |
| 952 | #ifdef Avoid_Underflow | |
| 953 | if (scale) { | |
| 954 | word0(rv0) = Exp_1 - 2*P*Exp_msk1; | |
| 955 | word1(rv0) = 0; | |
| 956 | dval(rv) *= dval(rv0); | |
| 957 | #ifndef NO_ERRNO | |
| 958 | /* try to avoid the bug of testing an 8087 register value */ | |
| 959 | if (word0(rv) == 0 && word1(rv) == 0) | |
| 960 | errno = ERANGE; | |
| 961 | #endif | |
| 962 | } | |
| 963 | #endif /* Avoid_Underflow */ | |
| 964 | #ifdef SET_INEXACT | |
| 965 | if (inexact && !(word0(rv) & Exp_mask)) { | |
| 966 | /* set underflow bit */ | |
| 967 | dval(rv0) = 1e-300; | |
| 968 | dval(rv0) *= dval(rv0); | |
| 969 | } | |
| 970 | #endif | |
| 971 | retfree: | |
| 972 | Bfree(bb); | |
| 973 | Bfree(bd); | |
| 974 | Bfree(bs); | |
| 975 | Bfree(bd0); | |
| 976 | Bfree(delta); | |
| 977 | ret: | |
| 978 | if (se) | |
| 979 | *se = (char *)s; | |
| 980 | return sign ? -dval(rv) : dval(rv); | |
| 981 | } | |
| 982 |