strtod.c revision 174679
1112158Sdas/**************************************************************** 2112158Sdas 3112158SdasThe author of this software is David M. Gay. 4112158Sdas 5112158SdasCopyright (C) 1998-2001 by Lucent Technologies 6112158SdasAll Rights Reserved 7112158Sdas 8112158SdasPermission to use, copy, modify, and distribute this software and 9112158Sdasits documentation for any purpose and without fee is hereby 10112158Sdasgranted, provided that the above copyright notice appear in all 11112158Sdascopies and that both that the copyright notice and this 12112158Sdaspermission notice and warranty disclaimer appear in supporting 13112158Sdasdocumentation, and that the name of Lucent or any of its entities 14112158Sdasnot be used in advertising or publicity pertaining to 15112158Sdasdistribution of the software without specific, written prior 16112158Sdaspermission. 17112158Sdas 18112158SdasLUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19112158SdasINCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20112158SdasIN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21112158SdasSPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22112158SdasWHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23112158SdasIN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24112158SdasARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25112158SdasTHIS SOFTWARE. 26112158Sdas 27112158Sdas****************************************************************/ 28112158Sdas 29165743Sdas/* Please send bug reports to David M. Gay (dmg at acm dot org, 30165743Sdas * with " at " changed at "@" and " dot " changed to "."). */ 31112158Sdas 32174679Sdas/* $FreeBSD: head/contrib/gdtoa/strtod.c 174679 2007-12-16 21:14:33Z das $ */ 33174679Sdas 34112158Sdas#include "gdtoaimp.h" 35165743Sdas#ifndef NO_FENV_H 36165743Sdas#include <fenv.h> 37165743Sdas#endif 38112158Sdas 39112158Sdas#ifdef USE_LOCALE 40112158Sdas#include "locale.h" 41112158Sdas#endif 42112158Sdas 43112158Sdas#ifdef IEEE_Arith 44112158Sdas#ifndef NO_IEEE_Scale 45112158Sdas#define Avoid_Underflow 46112158Sdas#undef tinytens 47112158Sdas/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 48112158Sdas/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 49112158Sdasstatic CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 50112158Sdas 9007199254740992.e-256 51112158Sdas }; 52112158Sdas#endif 53112158Sdas#endif 54112158Sdas 55112158Sdas#ifdef Honor_FLT_ROUNDS 56112158Sdas#define Rounding rounding 57112158Sdas#undef Check_FLT_ROUNDS 58112158Sdas#define Check_FLT_ROUNDS 59112158Sdas#else 60112158Sdas#define Rounding Flt_Rounds 61112158Sdas#endif 62112158Sdas 63112158Sdas double 64112158Sdasstrtod 65112158Sdas#ifdef KR_headers 66112158Sdas (s00, se) CONST char *s00; char **se; 67112158Sdas#else 68112158Sdas (CONST char *s00, char **se) 69112158Sdas#endif 70112158Sdas{ 71112158Sdas#ifdef Avoid_Underflow 72112158Sdas int scale; 73112158Sdas#endif 74165743Sdas int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, 75112158Sdas e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 76112158Sdas CONST char *s, *s0, *s1; 77112158Sdas double aadj, aadj1, adj, rv, rv0; 78112158Sdas Long L; 79112158Sdas ULong y, z; 80112158Sdas Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 81112158Sdas#ifdef SET_INEXACT 82112158Sdas int inexact, oldinexact; 83112158Sdas#endif 84112158Sdas#ifdef Honor_FLT_ROUNDS 85112158Sdas int rounding; 86112158Sdas#endif 87112158Sdas 88165743Sdas sign = nz0 = nz = decpt = 0; 89112158Sdas dval(rv) = 0.; 90112158Sdas for(s = s00;;s++) switch(*s) { 91112158Sdas case '-': 92112158Sdas sign = 1; 93112158Sdas /* no break */ 94112158Sdas case '+': 95112158Sdas if (*++s) 96112158Sdas goto break2; 97112158Sdas /* no break */ 98112158Sdas case 0: 99112158Sdas goto ret0; 100112158Sdas case '\t': 101112158Sdas case '\n': 102112158Sdas case '\v': 103112158Sdas case '\f': 104112158Sdas case '\r': 105112158Sdas case ' ': 106112158Sdas continue; 107112158Sdas default: 108112158Sdas goto break2; 109112158Sdas } 110112158Sdas break2: 111112158Sdas if (*s == '0') { 112112158Sdas#ifndef NO_HEX_FP 113112158Sdas { 114112158Sdas static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 115112158Sdas Long exp; 116112158Sdas ULong bits[2]; 117112158Sdas switch(s[1]) { 118112158Sdas case 'x': 119112158Sdas case 'X': 120165743Sdas { 121165743Sdas#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) 122165743Sdas FPI fpi1 = fpi; 123165743Sdas switch(fegetround()) { 124165743Sdas case FE_TOWARDZERO: fpi1.rounding = 0; break; 125165743Sdas case FE_UPWARD: fpi1.rounding = 2; break; 126165743Sdas case FE_DOWNWARD: fpi1.rounding = 3; 127165743Sdas } 128165743Sdas#else 129165743Sdas#define fpi1 fpi 130165743Sdas#endif 131165743Sdas switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { 132112158Sdas case STRTOG_NoNumber: 133112158Sdas s = s00; 134112158Sdas sign = 0; 135112158Sdas case STRTOG_Zero: 136112158Sdas break; 137112158Sdas default: 138124703Sdas if (bb) { 139124703Sdas copybits(bits, fpi.nbits, bb); 140124703Sdas Bfree(bb); 141124703Sdas } 142112158Sdas ULtod(((U*)&rv)->L, bits, exp, i); 143165743Sdas }} 144112158Sdas goto ret; 145112158Sdas } 146112158Sdas } 147112158Sdas#endif 148112158Sdas nz0 = 1; 149112158Sdas while(*++s == '0') ; 150112158Sdas if (!*s) 151112158Sdas goto ret; 152112158Sdas } 153112158Sdas s0 = s; 154112158Sdas y = z = 0; 155112158Sdas for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 156112158Sdas if (nd < 9) 157112158Sdas y = 10*y + c - '0'; 158112158Sdas else if (nd < 16) 159112158Sdas z = 10*z + c - '0'; 160112158Sdas nd0 = nd; 161112158Sdas#ifdef USE_LOCALE 162112415Sdas if (c == *localeconv()->decimal_point) 163112415Sdas#else 164112415Sdas if (c == '.') 165112158Sdas#endif 166112415Sdas { 167165743Sdas decpt = 1; 168112158Sdas c = *++s; 169112158Sdas if (!nd) { 170112158Sdas for(; c == '0'; c = *++s) 171112158Sdas nz++; 172112158Sdas if (c > '0' && c <= '9') { 173112158Sdas s0 = s; 174112158Sdas nf += nz; 175112158Sdas nz = 0; 176112158Sdas goto have_dig; 177112158Sdas } 178112158Sdas goto dig_done; 179112158Sdas } 180112158Sdas for(; c >= '0' && c <= '9'; c = *++s) { 181112158Sdas have_dig: 182112158Sdas nz++; 183112158Sdas if (c -= '0') { 184112158Sdas nf += nz; 185112158Sdas for(i = 1; i < nz; i++) 186112158Sdas if (nd++ < 9) 187112158Sdas y *= 10; 188112158Sdas else if (nd <= DBL_DIG + 1) 189112158Sdas z *= 10; 190112158Sdas if (nd++ < 9) 191112158Sdas y = 10*y + c; 192112158Sdas else if (nd <= DBL_DIG + 1) 193112158Sdas z = 10*z + c; 194112158Sdas nz = 0; 195112158Sdas } 196112158Sdas } 197112158Sdas } 198112158Sdas dig_done: 199112158Sdas e = 0; 200112158Sdas if (c == 'e' || c == 'E') { 201112158Sdas if (!nd && !nz && !nz0) { 202112158Sdas goto ret0; 203112158Sdas } 204112158Sdas s00 = s; 205112158Sdas esign = 0; 206112158Sdas switch(c = *++s) { 207112158Sdas case '-': 208112158Sdas esign = 1; 209112158Sdas case '+': 210112158Sdas c = *++s; 211112158Sdas } 212112158Sdas if (c >= '0' && c <= '9') { 213112158Sdas while(c == '0') 214112158Sdas c = *++s; 215112158Sdas if (c > '0' && c <= '9') { 216112158Sdas L = c - '0'; 217112158Sdas s1 = s; 218112158Sdas while((c = *++s) >= '0' && c <= '9') 219112158Sdas L = 10*L + c - '0'; 220112158Sdas if (s - s1 > 8 || L > 19999) 221112158Sdas /* Avoid confusion from exponents 222112158Sdas * so large that e might overflow. 223112158Sdas */ 224112158Sdas e = 19999; /* safe for 16 bit ints */ 225112158Sdas else 226112158Sdas e = (int)L; 227112158Sdas if (esign) 228112158Sdas e = -e; 229112158Sdas } 230112158Sdas else 231112158Sdas e = 0; 232112158Sdas } 233112158Sdas else 234112158Sdas s = s00; 235112158Sdas } 236112158Sdas if (!nd) { 237112158Sdas if (!nz && !nz0) { 238112158Sdas#ifdef INFNAN_CHECK 239112158Sdas /* Check for Nan and Infinity */ 240112158Sdas ULong bits[2]; 241112158Sdas static FPI fpinan = /* only 52 explicit bits */ 242112158Sdas { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 243165743Sdas if (!decpt) 244165743Sdas switch(c) { 245112158Sdas case 'i': 246112158Sdas case 'I': 247112158Sdas if (match(&s,"nf")) { 248112158Sdas --s; 249112158Sdas if (!match(&s,"inity")) 250112158Sdas ++s; 251112158Sdas word0(rv) = 0x7ff00000; 252112158Sdas word1(rv) = 0; 253112158Sdas goto ret; 254112158Sdas } 255112158Sdas break; 256112158Sdas case 'n': 257112158Sdas case 'N': 258112158Sdas if (match(&s, "an")) { 259112158Sdas#ifndef No_Hex_NaN 260112158Sdas if (*s == '(' /*)*/ 261112158Sdas && hexnan(&s, &fpinan, bits) 262112158Sdas == STRTOG_NaNbits) { 263174679Sdas word0(rv) = 0x7ff80000 | bits[1]; 264112158Sdas word1(rv) = bits[0]; 265112158Sdas } 266112158Sdas else { 267165743Sdas#endif 268112158Sdas word0(rv) = NAN_WORD0; 269112158Sdas word1(rv) = NAN_WORD1; 270165743Sdas#ifndef No_Hex_NaN 271112158Sdas } 272112158Sdas#endif 273112158Sdas goto ret; 274112158Sdas } 275112158Sdas } 276112158Sdas#endif /* INFNAN_CHECK */ 277112158Sdas ret0: 278112158Sdas s = s00; 279112158Sdas sign = 0; 280112158Sdas } 281112158Sdas goto ret; 282112158Sdas } 283112158Sdas e1 = e -= nf; 284112158Sdas 285112158Sdas /* Now we have nd0 digits, starting at s0, followed by a 286112158Sdas * decimal point, followed by nd-nd0 digits. The number we're 287112158Sdas * after is the integer represented by those digits times 288112158Sdas * 10**e */ 289112158Sdas 290112158Sdas if (!nd0) 291112158Sdas nd0 = nd; 292112158Sdas k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 293112158Sdas dval(rv) = y; 294112158Sdas if (k > 9) { 295112158Sdas#ifdef SET_INEXACT 296112158Sdas if (k > DBL_DIG) 297112158Sdas oldinexact = get_inexact(); 298112158Sdas#endif 299112158Sdas dval(rv) = tens[k - 9] * dval(rv) + z; 300112158Sdas } 301112158Sdas bd0 = 0; 302112158Sdas if (nd <= DBL_DIG 303112158Sdas#ifndef RND_PRODQUOT 304112158Sdas#ifndef Honor_FLT_ROUNDS 305112158Sdas && Flt_Rounds == 1 306112158Sdas#endif 307112158Sdas#endif 308112158Sdas ) { 309112158Sdas if (!e) 310112158Sdas goto ret; 311112158Sdas if (e > 0) { 312112158Sdas if (e <= Ten_pmax) { 313112158Sdas#ifdef VAX 314112158Sdas goto vax_ovfl_check; 315112158Sdas#else 316112158Sdas#ifdef Honor_FLT_ROUNDS 317112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 318112158Sdas if (sign) { 319112158Sdas rv = -rv; 320112158Sdas sign = 0; 321112158Sdas } 322112158Sdas#endif 323112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 324112158Sdas goto ret; 325112158Sdas#endif 326112158Sdas } 327112158Sdas i = DBL_DIG - nd; 328112158Sdas if (e <= Ten_pmax + i) { 329112158Sdas /* A fancier test would sometimes let us do 330112158Sdas * this for larger i values. 331112158Sdas */ 332112158Sdas#ifdef Honor_FLT_ROUNDS 333112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 334112158Sdas if (sign) { 335112158Sdas rv = -rv; 336112158Sdas sign = 0; 337112158Sdas } 338112158Sdas#endif 339112158Sdas e -= i; 340112158Sdas dval(rv) *= tens[i]; 341112158Sdas#ifdef VAX 342112158Sdas /* VAX exponent range is so narrow we must 343112158Sdas * worry about overflow here... 344112158Sdas */ 345112158Sdas vax_ovfl_check: 346112158Sdas word0(rv) -= P*Exp_msk1; 347112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 348112158Sdas if ((word0(rv) & Exp_mask) 349112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 350112158Sdas goto ovfl; 351112158Sdas word0(rv) += P*Exp_msk1; 352112158Sdas#else 353112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 354112158Sdas#endif 355112158Sdas goto ret; 356112158Sdas } 357112158Sdas } 358112158Sdas#ifndef Inaccurate_Divide 359112158Sdas else if (e >= -Ten_pmax) { 360112158Sdas#ifdef Honor_FLT_ROUNDS 361112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 362112158Sdas if (sign) { 363112158Sdas rv = -rv; 364112158Sdas sign = 0; 365112158Sdas } 366112158Sdas#endif 367112158Sdas /* rv = */ rounded_quotient(dval(rv), tens[-e]); 368112158Sdas goto ret; 369112158Sdas } 370112158Sdas#endif 371112158Sdas } 372112158Sdas e1 += nd - k; 373112158Sdas 374112158Sdas#ifdef IEEE_Arith 375112158Sdas#ifdef SET_INEXACT 376112158Sdas inexact = 1; 377112158Sdas if (k <= DBL_DIG) 378112158Sdas oldinexact = get_inexact(); 379112158Sdas#endif 380112158Sdas#ifdef Avoid_Underflow 381112158Sdas scale = 0; 382112158Sdas#endif 383112158Sdas#ifdef Honor_FLT_ROUNDS 384112158Sdas if ((rounding = Flt_Rounds) >= 2) { 385112158Sdas if (sign) 386112158Sdas rounding = rounding == 2 ? 0 : 2; 387112158Sdas else 388112158Sdas if (rounding != 2) 389112158Sdas rounding = 0; 390112158Sdas } 391112158Sdas#endif 392112158Sdas#endif /*IEEE_Arith*/ 393112158Sdas 394112158Sdas /* Get starting approximation = rv * 10**e1 */ 395112158Sdas 396112158Sdas if (e1 > 0) { 397112158Sdas if ( (i = e1 & 15) !=0) 398112158Sdas dval(rv) *= tens[i]; 399112158Sdas if (e1 &= ~15) { 400112158Sdas if (e1 > DBL_MAX_10_EXP) { 401112158Sdas ovfl: 402112158Sdas#ifndef NO_ERRNO 403112158Sdas errno = ERANGE; 404112158Sdas#endif 405112158Sdas /* Can't trust HUGE_VAL */ 406112158Sdas#ifdef IEEE_Arith 407112158Sdas#ifdef Honor_FLT_ROUNDS 408112158Sdas switch(rounding) { 409112158Sdas case 0: /* toward 0 */ 410112158Sdas case 3: /* toward -infinity */ 411112158Sdas word0(rv) = Big0; 412112158Sdas word1(rv) = Big1; 413112158Sdas break; 414112158Sdas default: 415112158Sdas word0(rv) = Exp_mask; 416112158Sdas word1(rv) = 0; 417112158Sdas } 418112158Sdas#else /*Honor_FLT_ROUNDS*/ 419112158Sdas word0(rv) = Exp_mask; 420112158Sdas word1(rv) = 0; 421112158Sdas#endif /*Honor_FLT_ROUNDS*/ 422112158Sdas#ifdef SET_INEXACT 423112158Sdas /* set overflow bit */ 424112158Sdas dval(rv0) = 1e300; 425112158Sdas dval(rv0) *= dval(rv0); 426112158Sdas#endif 427112158Sdas#else /*IEEE_Arith*/ 428112158Sdas word0(rv) = Big0; 429112158Sdas word1(rv) = Big1; 430112158Sdas#endif /*IEEE_Arith*/ 431112158Sdas if (bd0) 432112158Sdas goto retfree; 433112158Sdas goto ret; 434112158Sdas } 435112158Sdas e1 >>= 4; 436112158Sdas for(j = 0; e1 > 1; j++, e1 >>= 1) 437112158Sdas if (e1 & 1) 438112158Sdas dval(rv) *= bigtens[j]; 439112158Sdas /* The last multiplication could overflow. */ 440112158Sdas word0(rv) -= P*Exp_msk1; 441112158Sdas dval(rv) *= bigtens[j]; 442112158Sdas if ((z = word0(rv) & Exp_mask) 443112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 444112158Sdas goto ovfl; 445112158Sdas if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 446112158Sdas /* set to largest number */ 447112158Sdas /* (Can't trust DBL_MAX) */ 448112158Sdas word0(rv) = Big0; 449112158Sdas word1(rv) = Big1; 450112158Sdas } 451112158Sdas else 452112158Sdas word0(rv) += P*Exp_msk1; 453112158Sdas } 454112158Sdas } 455112158Sdas else if (e1 < 0) { 456112158Sdas e1 = -e1; 457112158Sdas if ( (i = e1 & 15) !=0) 458112158Sdas dval(rv) /= tens[i]; 459112158Sdas if (e1 >>= 4) { 460112158Sdas if (e1 >= 1 << n_bigtens) 461112158Sdas goto undfl; 462112158Sdas#ifdef Avoid_Underflow 463112158Sdas if (e1 & Scale_Bit) 464112158Sdas scale = 2*P; 465112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 466112158Sdas if (e1 & 1) 467112158Sdas dval(rv) *= tinytens[j]; 468112158Sdas if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 469112158Sdas >> Exp_shift)) > 0) { 470112158Sdas /* scaled rv is denormal; zap j low bits */ 471112158Sdas if (j >= 32) { 472112158Sdas word1(rv) = 0; 473112158Sdas if (j >= 53) 474112158Sdas word0(rv) = (P+2)*Exp_msk1; 475112158Sdas else 476112158Sdas word0(rv) &= 0xffffffff << j-32; 477112158Sdas } 478112158Sdas else 479112158Sdas word1(rv) &= 0xffffffff << j; 480112158Sdas } 481112158Sdas#else 482112158Sdas for(j = 0; e1 > 1; j++, e1 >>= 1) 483112158Sdas if (e1 & 1) 484112158Sdas dval(rv) *= tinytens[j]; 485112158Sdas /* The last multiplication could underflow. */ 486112158Sdas dval(rv0) = dval(rv); 487112158Sdas dval(rv) *= tinytens[j]; 488112158Sdas if (!dval(rv)) { 489112158Sdas dval(rv) = 2.*dval(rv0); 490112158Sdas dval(rv) *= tinytens[j]; 491112158Sdas#endif 492112158Sdas if (!dval(rv)) { 493112158Sdas undfl: 494112158Sdas dval(rv) = 0.; 495112158Sdas#ifndef NO_ERRNO 496112158Sdas errno = ERANGE; 497112158Sdas#endif 498112158Sdas if (bd0) 499112158Sdas goto retfree; 500112158Sdas goto ret; 501112158Sdas } 502112158Sdas#ifndef Avoid_Underflow 503112158Sdas word0(rv) = Tiny0; 504112158Sdas word1(rv) = Tiny1; 505112158Sdas /* The refinement below will clean 506112158Sdas * this approximation up. 507112158Sdas */ 508112158Sdas } 509112158Sdas#endif 510112158Sdas } 511112158Sdas } 512112158Sdas 513112158Sdas /* Now the hard part -- adjusting rv to the correct value.*/ 514112158Sdas 515112158Sdas /* Put digits into bd: true value = bd * 10^e */ 516112158Sdas 517112158Sdas bd0 = s2b(s0, nd0, nd, y); 518112158Sdas 519112158Sdas for(;;) { 520112158Sdas bd = Balloc(bd0->k); 521112158Sdas Bcopy(bd, bd0); 522112158Sdas bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 523112158Sdas bs = i2b(1); 524112158Sdas 525112158Sdas if (e >= 0) { 526112158Sdas bb2 = bb5 = 0; 527112158Sdas bd2 = bd5 = e; 528112158Sdas } 529112158Sdas else { 530112158Sdas bb2 = bb5 = -e; 531112158Sdas bd2 = bd5 = 0; 532112158Sdas } 533112158Sdas if (bbe >= 0) 534112158Sdas bb2 += bbe; 535112158Sdas else 536112158Sdas bd2 -= bbe; 537112158Sdas bs2 = bb2; 538112158Sdas#ifdef Honor_FLT_ROUNDS 539112158Sdas if (rounding != 1) 540112158Sdas bs2++; 541112158Sdas#endif 542112158Sdas#ifdef Avoid_Underflow 543112158Sdas j = bbe - scale; 544112158Sdas i = j + bbbits - 1; /* logb(rv) */ 545112158Sdas if (i < Emin) /* denormal */ 546112158Sdas j += P - Emin; 547112158Sdas else 548112158Sdas j = P + 1 - bbbits; 549112158Sdas#else /*Avoid_Underflow*/ 550112158Sdas#ifdef Sudden_Underflow 551112158Sdas#ifdef IBM 552112158Sdas j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 553112158Sdas#else 554112158Sdas j = P + 1 - bbbits; 555112158Sdas#endif 556112158Sdas#else /*Sudden_Underflow*/ 557112158Sdas j = bbe; 558112158Sdas i = j + bbbits - 1; /* logb(rv) */ 559112158Sdas if (i < Emin) /* denormal */ 560112158Sdas j += P - Emin; 561112158Sdas else 562112158Sdas j = P + 1 - bbbits; 563112158Sdas#endif /*Sudden_Underflow*/ 564112158Sdas#endif /*Avoid_Underflow*/ 565112158Sdas bb2 += j; 566112158Sdas bd2 += j; 567112158Sdas#ifdef Avoid_Underflow 568112158Sdas bd2 += scale; 569112158Sdas#endif 570112158Sdas i = bb2 < bd2 ? bb2 : bd2; 571112158Sdas if (i > bs2) 572112158Sdas i = bs2; 573112158Sdas if (i > 0) { 574112158Sdas bb2 -= i; 575112158Sdas bd2 -= i; 576112158Sdas bs2 -= i; 577112158Sdas } 578112158Sdas if (bb5 > 0) { 579112158Sdas bs = pow5mult(bs, bb5); 580112158Sdas bb1 = mult(bs, bb); 581112158Sdas Bfree(bb); 582112158Sdas bb = bb1; 583112158Sdas } 584112158Sdas if (bb2 > 0) 585112158Sdas bb = lshift(bb, bb2); 586112158Sdas if (bd5 > 0) 587112158Sdas bd = pow5mult(bd, bd5); 588112158Sdas if (bd2 > 0) 589112158Sdas bd = lshift(bd, bd2); 590112158Sdas if (bs2 > 0) 591112158Sdas bs = lshift(bs, bs2); 592112158Sdas delta = diff(bb, bd); 593112158Sdas dsign = delta->sign; 594112158Sdas delta->sign = 0; 595112158Sdas i = cmp(delta, bs); 596112158Sdas#ifdef Honor_FLT_ROUNDS 597112158Sdas if (rounding != 1) { 598112158Sdas if (i < 0) { 599112158Sdas /* Error is less than an ulp */ 600112158Sdas if (!delta->x[0] && delta->wds <= 1) { 601112158Sdas /* exact */ 602112158Sdas#ifdef SET_INEXACT 603112158Sdas inexact = 0; 604112158Sdas#endif 605112158Sdas break; 606112158Sdas } 607112158Sdas if (rounding) { 608112158Sdas if (dsign) { 609112158Sdas adj = 1.; 610112158Sdas goto apply_adj; 611112158Sdas } 612112158Sdas } 613112158Sdas else if (!dsign) { 614112158Sdas adj = -1.; 615112158Sdas if (!word1(rv) 616112158Sdas && !(word0(rv) & Frac_mask)) { 617112158Sdas y = word0(rv) & Exp_mask; 618112158Sdas#ifdef Avoid_Underflow 619112158Sdas if (!scale || y > 2*P*Exp_msk1) 620112158Sdas#else 621112158Sdas if (y) 622112158Sdas#endif 623112158Sdas { 624112158Sdas delta = lshift(delta,Log2P); 625112158Sdas if (cmp(delta, bs) <= 0) 626112158Sdas adj = -0.5; 627112158Sdas } 628112158Sdas } 629112158Sdas apply_adj: 630112158Sdas#ifdef Avoid_Underflow 631112158Sdas if (scale && (y = word0(rv) & Exp_mask) 632112158Sdas <= 2*P*Exp_msk1) 633112158Sdas word0(adj) += (2*P+1)*Exp_msk1 - y; 634112158Sdas#else 635112158Sdas#ifdef Sudden_Underflow 636112158Sdas if ((word0(rv) & Exp_mask) <= 637112158Sdas P*Exp_msk1) { 638112158Sdas word0(rv) += P*Exp_msk1; 639112158Sdas dval(rv) += adj*ulp(dval(rv)); 640112158Sdas word0(rv) -= P*Exp_msk1; 641112158Sdas } 642112158Sdas else 643112158Sdas#endif /*Sudden_Underflow*/ 644112158Sdas#endif /*Avoid_Underflow*/ 645112158Sdas dval(rv) += adj*ulp(dval(rv)); 646112158Sdas } 647112158Sdas break; 648112158Sdas } 649112158Sdas adj = ratio(delta, bs); 650112158Sdas if (adj < 1.) 651112158Sdas adj = 1.; 652112158Sdas if (adj <= 0x7ffffffe) { 653112158Sdas /* adj = rounding ? ceil(adj) : floor(adj); */ 654112158Sdas y = adj; 655112158Sdas if (y != adj) { 656112158Sdas if (!((rounding>>1) ^ dsign)) 657112158Sdas y++; 658112158Sdas adj = y; 659112158Sdas } 660112158Sdas } 661112158Sdas#ifdef Avoid_Underflow 662112158Sdas if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 663112158Sdas word0(adj) += (2*P+1)*Exp_msk1 - y; 664112158Sdas#else 665112158Sdas#ifdef Sudden_Underflow 666112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 667112158Sdas word0(rv) += P*Exp_msk1; 668112158Sdas adj *= ulp(dval(rv)); 669112158Sdas if (dsign) 670112158Sdas dval(rv) += adj; 671112158Sdas else 672112158Sdas dval(rv) -= adj; 673112158Sdas word0(rv) -= P*Exp_msk1; 674112158Sdas goto cont; 675112158Sdas } 676112158Sdas#endif /*Sudden_Underflow*/ 677112158Sdas#endif /*Avoid_Underflow*/ 678112158Sdas adj *= ulp(dval(rv)); 679112158Sdas if (dsign) 680112158Sdas dval(rv) += adj; 681112158Sdas else 682112158Sdas dval(rv) -= adj; 683112158Sdas goto cont; 684112158Sdas } 685112158Sdas#endif /*Honor_FLT_ROUNDS*/ 686112158Sdas 687112158Sdas if (i < 0) { 688112158Sdas /* Error is less than half an ulp -- check for 689112158Sdas * special case of mantissa a power of two. 690112158Sdas */ 691112158Sdas if (dsign || word1(rv) || word0(rv) & Bndry_mask 692112158Sdas#ifdef IEEE_Arith 693112158Sdas#ifdef Avoid_Underflow 694112158Sdas || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 695112158Sdas#else 696112158Sdas || (word0(rv) & Exp_mask) <= Exp_msk1 697112158Sdas#endif 698112158Sdas#endif 699112158Sdas ) { 700112158Sdas#ifdef SET_INEXACT 701112158Sdas if (!delta->x[0] && delta->wds <= 1) 702112158Sdas inexact = 0; 703112158Sdas#endif 704112158Sdas break; 705112158Sdas } 706112158Sdas if (!delta->x[0] && delta->wds <= 1) { 707112158Sdas /* exact result */ 708112158Sdas#ifdef SET_INEXACT 709112158Sdas inexact = 0; 710112158Sdas#endif 711112158Sdas break; 712112158Sdas } 713112158Sdas delta = lshift(delta,Log2P); 714112158Sdas if (cmp(delta, bs) > 0) 715112158Sdas goto drop_down; 716112158Sdas break; 717112158Sdas } 718112158Sdas if (i == 0) { 719112158Sdas /* exactly half-way between */ 720112158Sdas if (dsign) { 721112158Sdas if ((word0(rv) & Bndry_mask1) == Bndry_mask1 722112158Sdas && word1(rv) == ( 723112158Sdas#ifdef Avoid_Underflow 724112158Sdas (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 725112158Sdas ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 726112158Sdas#endif 727112158Sdas 0xffffffff)) { 728112158Sdas /*boundary case -- increment exponent*/ 729112158Sdas word0(rv) = (word0(rv) & Exp_mask) 730112158Sdas + Exp_msk1 731112158Sdas#ifdef IBM 732112158Sdas | Exp_msk1 >> 4 733112158Sdas#endif 734112158Sdas ; 735112158Sdas word1(rv) = 0; 736112158Sdas#ifdef Avoid_Underflow 737112158Sdas dsign = 0; 738112158Sdas#endif 739112158Sdas break; 740112158Sdas } 741112158Sdas } 742112158Sdas else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 743112158Sdas drop_down: 744112158Sdas /* boundary case -- decrement exponent */ 745112158Sdas#ifdef Sudden_Underflow /*{{*/ 746112158Sdas L = word0(rv) & Exp_mask; 747112158Sdas#ifdef IBM 748112158Sdas if (L < Exp_msk1) 749112158Sdas#else 750112158Sdas#ifdef Avoid_Underflow 751112158Sdas if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 752112158Sdas#else 753112158Sdas if (L <= Exp_msk1) 754112158Sdas#endif /*Avoid_Underflow*/ 755112158Sdas#endif /*IBM*/ 756112158Sdas goto undfl; 757112158Sdas L -= Exp_msk1; 758112158Sdas#else /*Sudden_Underflow}{*/ 759112158Sdas#ifdef Avoid_Underflow 760112158Sdas if (scale) { 761112158Sdas L = word0(rv) & Exp_mask; 762112158Sdas if (L <= (2*P+1)*Exp_msk1) { 763112158Sdas if (L > (P+2)*Exp_msk1) 764112158Sdas /* round even ==> */ 765112158Sdas /* accept rv */ 766112158Sdas break; 767112158Sdas /* rv = smallest denormal */ 768112158Sdas goto undfl; 769112158Sdas } 770112158Sdas } 771112158Sdas#endif /*Avoid_Underflow*/ 772112158Sdas L = (word0(rv) & Exp_mask) - Exp_msk1; 773112158Sdas#endif /*Sudden_Underflow}*/ 774112158Sdas word0(rv) = L | Bndry_mask1; 775112158Sdas word1(rv) = 0xffffffff; 776112158Sdas#ifdef IBM 777112158Sdas goto cont; 778112158Sdas#else 779112158Sdas break; 780112158Sdas#endif 781112158Sdas } 782112158Sdas#ifndef ROUND_BIASED 783112158Sdas if (!(word1(rv) & LSB)) 784112158Sdas break; 785112158Sdas#endif 786112158Sdas if (dsign) 787112158Sdas dval(rv) += ulp(dval(rv)); 788112158Sdas#ifndef ROUND_BIASED 789112158Sdas else { 790112158Sdas dval(rv) -= ulp(dval(rv)); 791112158Sdas#ifndef Sudden_Underflow 792112158Sdas if (!dval(rv)) 793112158Sdas goto undfl; 794112158Sdas#endif 795112158Sdas } 796112158Sdas#ifdef Avoid_Underflow 797112158Sdas dsign = 1 - dsign; 798112158Sdas#endif 799112158Sdas#endif 800112158Sdas break; 801112158Sdas } 802112158Sdas if ((aadj = ratio(delta, bs)) <= 2.) { 803112158Sdas if (dsign) 804112158Sdas aadj = aadj1 = 1.; 805112158Sdas else if (word1(rv) || word0(rv) & Bndry_mask) { 806112158Sdas#ifndef Sudden_Underflow 807112158Sdas if (word1(rv) == Tiny1 && !word0(rv)) 808112158Sdas goto undfl; 809112158Sdas#endif 810112158Sdas aadj = 1.; 811112158Sdas aadj1 = -1.; 812112158Sdas } 813112158Sdas else { 814112158Sdas /* special case -- power of FLT_RADIX to be */ 815112158Sdas /* rounded down... */ 816112158Sdas 817112158Sdas if (aadj < 2./FLT_RADIX) 818112158Sdas aadj = 1./FLT_RADIX; 819112158Sdas else 820112158Sdas aadj *= 0.5; 821112158Sdas aadj1 = -aadj; 822112158Sdas } 823112158Sdas } 824112158Sdas else { 825112158Sdas aadj *= 0.5; 826112158Sdas aadj1 = dsign ? aadj : -aadj; 827112158Sdas#ifdef Check_FLT_ROUNDS 828112158Sdas switch(Rounding) { 829112158Sdas case 2: /* towards +infinity */ 830112158Sdas aadj1 -= 0.5; 831112158Sdas break; 832112158Sdas case 0: /* towards 0 */ 833112158Sdas case 3: /* towards -infinity */ 834112158Sdas aadj1 += 0.5; 835112158Sdas } 836112158Sdas#else 837112158Sdas if (Flt_Rounds == 0) 838112158Sdas aadj1 += 0.5; 839112158Sdas#endif /*Check_FLT_ROUNDS*/ 840112158Sdas } 841112158Sdas y = word0(rv) & Exp_mask; 842112158Sdas 843112158Sdas /* Check for overflow */ 844112158Sdas 845112158Sdas if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 846112158Sdas dval(rv0) = dval(rv); 847112158Sdas word0(rv) -= P*Exp_msk1; 848112158Sdas adj = aadj1 * ulp(dval(rv)); 849112158Sdas dval(rv) += adj; 850112158Sdas if ((word0(rv) & Exp_mask) >= 851112158Sdas Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 852112158Sdas if (word0(rv0) == Big0 && word1(rv0) == Big1) 853112158Sdas goto ovfl; 854112158Sdas word0(rv) = Big0; 855112158Sdas word1(rv) = Big1; 856112158Sdas goto cont; 857112158Sdas } 858112158Sdas else 859112158Sdas word0(rv) += P*Exp_msk1; 860112158Sdas } 861112158Sdas else { 862112158Sdas#ifdef Avoid_Underflow 863112158Sdas if (scale && y <= 2*P*Exp_msk1) { 864112158Sdas if (aadj <= 0x7fffffff) { 865112158Sdas if ((z = aadj) <= 0) 866112158Sdas z = 1; 867112158Sdas aadj = z; 868112158Sdas aadj1 = dsign ? aadj : -aadj; 869112158Sdas } 870112158Sdas word0(aadj1) += (2*P+1)*Exp_msk1 - y; 871112158Sdas } 872112158Sdas adj = aadj1 * ulp(dval(rv)); 873112158Sdas dval(rv) += adj; 874112158Sdas#else 875112158Sdas#ifdef Sudden_Underflow 876112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 877112158Sdas dval(rv0) = dval(rv); 878112158Sdas word0(rv) += P*Exp_msk1; 879112158Sdas adj = aadj1 * ulp(dval(rv)); 880112158Sdas dval(rv) += adj; 881112158Sdas#ifdef IBM 882112158Sdas if ((word0(rv) & Exp_mask) < P*Exp_msk1) 883112158Sdas#else 884112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 885112158Sdas#endif 886112158Sdas { 887112158Sdas if (word0(rv0) == Tiny0 888112158Sdas && word1(rv0) == Tiny1) 889112158Sdas goto undfl; 890112158Sdas word0(rv) = Tiny0; 891112158Sdas word1(rv) = Tiny1; 892112158Sdas goto cont; 893112158Sdas } 894112158Sdas else 895112158Sdas word0(rv) -= P*Exp_msk1; 896112158Sdas } 897112158Sdas else { 898112158Sdas adj = aadj1 * ulp(dval(rv)); 899112158Sdas dval(rv) += adj; 900112158Sdas } 901112158Sdas#else /*Sudden_Underflow*/ 902112158Sdas /* Compute adj so that the IEEE rounding rules will 903112158Sdas * correctly round rv + adj in some half-way cases. 904112158Sdas * If rv * ulp(rv) is denormalized (i.e., 905112158Sdas * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 906112158Sdas * trouble from bits lost to denormalization; 907112158Sdas * example: 1.2e-307 . 908112158Sdas */ 909112158Sdas if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 910112158Sdas aadj1 = (double)(int)(aadj + 0.5); 911112158Sdas if (!dsign) 912112158Sdas aadj1 = -aadj1; 913112158Sdas } 914112158Sdas adj = aadj1 * ulp(dval(rv)); 915112158Sdas dval(rv) += adj; 916112158Sdas#endif /*Sudden_Underflow*/ 917112158Sdas#endif /*Avoid_Underflow*/ 918112158Sdas } 919112158Sdas z = word0(rv) & Exp_mask; 920112158Sdas#ifndef SET_INEXACT 921112158Sdas#ifdef Avoid_Underflow 922112158Sdas if (!scale) 923112158Sdas#endif 924112158Sdas if (y == z) { 925112158Sdas /* Can we stop now? */ 926112158Sdas L = (Long)aadj; 927112158Sdas aadj -= L; 928112158Sdas /* The tolerances below are conservative. */ 929112158Sdas if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 930112158Sdas if (aadj < .4999999 || aadj > .5000001) 931112158Sdas break; 932112158Sdas } 933112158Sdas else if (aadj < .4999999/FLT_RADIX) 934112158Sdas break; 935112158Sdas } 936112158Sdas#endif 937112158Sdas cont: 938112158Sdas Bfree(bb); 939112158Sdas Bfree(bd); 940112158Sdas Bfree(bs); 941112158Sdas Bfree(delta); 942112158Sdas } 943112158Sdas#ifdef SET_INEXACT 944112158Sdas if (inexact) { 945112158Sdas if (!oldinexact) { 946112158Sdas word0(rv0) = Exp_1 + (70 << Exp_shift); 947112158Sdas word1(rv0) = 0; 948112158Sdas dval(rv0) += 1.; 949112158Sdas } 950112158Sdas } 951112158Sdas else if (!oldinexact) 952112158Sdas clear_inexact(); 953112158Sdas#endif 954112158Sdas#ifdef Avoid_Underflow 955112158Sdas if (scale) { 956112158Sdas word0(rv0) = Exp_1 - 2*P*Exp_msk1; 957112158Sdas word1(rv0) = 0; 958112158Sdas dval(rv) *= dval(rv0); 959112158Sdas#ifndef NO_ERRNO 960112158Sdas /* try to avoid the bug of testing an 8087 register value */ 961112158Sdas if (word0(rv) == 0 && word1(rv) == 0) 962112158Sdas errno = ERANGE; 963112158Sdas#endif 964112158Sdas } 965112158Sdas#endif /* Avoid_Underflow */ 966112158Sdas#ifdef SET_INEXACT 967112158Sdas if (inexact && !(word0(rv) & Exp_mask)) { 968112158Sdas /* set underflow bit */ 969112158Sdas dval(rv0) = 1e-300; 970112158Sdas dval(rv0) *= dval(rv0); 971112158Sdas } 972112158Sdas#endif 973112158Sdas retfree: 974112158Sdas Bfree(bb); 975112158Sdas Bfree(bd); 976112158Sdas Bfree(bs); 977112158Sdas Bfree(bd0); 978112158Sdas Bfree(delta); 979112158Sdas ret: 980112158Sdas if (se) 981112158Sdas *se = (char *)s; 982112158Sdas return sign ? -dval(rv) : dval(rv); 983112158Sdas } 984112158Sdas 985