strtod.c revision 219557
1235783Skib/**************************************************************** 2235783Skib 3235783SkibThe author of this software is David M. Gay. 4235783Skib 5235783SkibCopyright (C) 1998-2001 by Lucent Technologies 6235783SkibAll Rights Reserved 7235783Skib 8235783SkibPermission to use, copy, modify, and distribute this software and 9235783Skibits documentation for any purpose and without fee is hereby 10235783Skibgranted, provided that the above copyright notice appear in all 11235783Skibcopies and that both that the copyright notice and this 12235783Skibpermission notice and warranty disclaimer appear in supporting 13235783Skibdocumentation, and that the name of Lucent or any of its entities 14235783Skibnot be used in advertising or publicity pertaining to 15235783Skibdistribution of the software without specific, written prior 16235783Skibpermission. 17235783Skib 18235783SkibLUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19235783SkibINCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20235783SkibIN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21235783SkibSPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22235783SkibWHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23235783SkibIN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24235783SkibARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25235783SkibTHIS SOFTWARE. 26235783Skib 27235783Skib****************************************************************/ 28235783Skib 29235783Skib/* Please send bug reports to David M. Gay (dmg at acm dot org, 30235783Skib * with " at " changed at "@" and " dot " changed to "."). */ 31235783Skib 32235783Skib/* $FreeBSD: head/contrib/gdtoa/strtod.c 219557 2011-03-12 07:03:06Z das $ */ 33235783Skib 34235783Skib#include "gdtoaimp.h" 35235783Skib#ifndef NO_FENV_H 36235783Skib#include <fenv.h> 37235783Skib#endif 38235783Skib 39235783Skib#ifdef USE_LOCALE 40235783Skib#include "locale.h" 41235783Skib#endif 42235783Skib 43235783Skib#ifdef IEEE_Arith 44235783Skib#ifndef NO_IEEE_Scale 45235783Skib#define Avoid_Underflow 46235783Skib#undef tinytens 47235783Skib/* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */ 48235783Skib/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 49235783Skibstatic CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 50235783Skib 9007199254740992.*9007199254740992.e-256 51235783Skib }; 52235783Skib#endif 53235783Skib#endif 54235783Skib 55235783Skib#ifdef Honor_FLT_ROUNDS 56235783Skib#undef Check_FLT_ROUNDS 57235783Skib#define Check_FLT_ROUNDS 58235783Skib#else 59235783Skib#define Rounding Flt_Rounds 60235783Skib#endif 61235783Skib 62235783Skib#ifdef Avoid_Underflow /*{*/ 63235783Skib static double 64235783Skibsulp 65235783Skib#ifdef KR_headers 66235783Skib (x, scale) U *x; int scale; 67235783Skib#else 68235783Skib (U *x, int scale) 69235783Skib#endif 70235783Skib{ 71235783Skib U u; 72235783Skib double rv; 73235783Skib int i; 74235783Skib 75235783Skib rv = ulp(x); 76235783Skib if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0) 77235783Skib return rv; /* Is there an example where i <= 0 ? */ 78235783Skib word0(&u) = Exp_1 + (i << Exp_shift); 79235783Skib word1(&u) = 0; 80235783Skib return rv * u.d; 81235783Skib } 82235783Skib#endif /*}*/ 83235783Skib 84235783Skib double 85235783Skibstrtod 86235783Skib#ifdef KR_headers 87235783Skib (s00, se) CONST char *s00; char **se; 88235783Skib#else 89235783Skib (CONST char *s00, char **se) 90235783Skib#endif 91235783Skib{ 92235783Skib#ifdef Avoid_Underflow 93235783Skib int scale; 94235783Skib#endif 95235783Skib int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, 96235783Skib e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 97235783Skib CONST char *s, *s0, *s1; 98235783Skib double aadj; 99235783Skib Long L; 100235783Skib U adj, aadj1, rv, rv0; 101235783Skib ULong y, z; 102235783Skib Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 103235783Skib#ifdef Avoid_Underflow 104235783Skib ULong Lsb, Lsb1; 105235783Skib#endif 106235783Skib#ifdef SET_INEXACT 107235783Skib int inexact, oldinexact; 108235783Skib#endif 109235783Skib#ifdef USE_LOCALE /*{{*/ 110235783Skib#ifdef NO_LOCALE_CACHE 111235783Skib char *decimalpoint = localeconv()->decimal_point; 112235783Skib int dplen = strlen(decimalpoint); 113235783Skib#else 114235783Skib char *decimalpoint; 115235783Skib static char *decimalpoint_cache; 116235783Skib static int dplen; 117235783Skib if (!(s0 = decimalpoint_cache)) { 118235783Skib s0 = localeconv()->decimal_point; 119235783Skib if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { 120235783Skib strcpy(decimalpoint_cache, s0); 121235783Skib s0 = decimalpoint_cache; 122235783Skib } 123235783Skib dplen = strlen(s0); 124235783Skib } 125235783Skib decimalpoint = (char*)s0; 126235783Skib#endif /*NO_LOCALE_CACHE*/ 127235783Skib#else /*USE_LOCALE}{*/ 128235783Skib#define dplen 1 129235783Skib#endif /*USE_LOCALE}}*/ 130235783Skib 131235783Skib#ifdef Honor_FLT_ROUNDS /*{*/ 132235783Skib int Rounding; 133235783Skib#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 134235783Skib Rounding = Flt_Rounds; 135235783Skib#else /*}{*/ 136235783Skib Rounding = 1; 137235783Skib switch(fegetround()) { 138235783Skib case FE_TOWARDZERO: Rounding = 0; break; 139235783Skib case FE_UPWARD: Rounding = 2; break; 140235783Skib case FE_DOWNWARD: Rounding = 3; 141235783Skib } 142235783Skib#endif /*}}*/ 143235783Skib#endif /*}*/ 144235783Skib 145235783Skib sign = nz0 = nz = decpt = 0; 146235783Skib dval(&rv) = 0.; 147235783Skib for(s = s00;;s++) switch(*s) { 148235783Skib case '-': 149235783Skib sign = 1; 150235783Skib /* no break */ 151235783Skib case '+': 152235783Skib if (*++s) 153235783Skib goto break2; 154235783Skib /* no break */ 155235783Skib case 0: 156235783Skib goto ret0; 157235783Skib case '\t': 158235783Skib case '\n': 159235783Skib case '\v': 160235783Skib case '\f': 161235783Skib case '\r': 162235783Skib case ' ': 163235783Skib continue; 164235783Skib default: 165235783Skib goto break2; 166235783Skib } 167235783Skib break2: 168235783Skib if (*s == '0') { 169235783Skib#ifndef NO_HEX_FP /*{*/ 170235783Skib { 171235783Skib static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 172235783Skib Long exp; 173235783Skib ULong bits[2]; 174235783Skib switch(s[1]) { 175235783Skib case 'x': 176235783Skib case 'X': 177235783Skib { 178235783Skib#ifdef Honor_FLT_ROUNDS 179235783Skib FPI fpi1 = fpi; 180235783Skib fpi1.rounding = Rounding; 181235783Skib#else 182235783Skib#define fpi1 fpi 183235783Skib#endif 184235783Skib switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { 185235783Skib case STRTOG_NoNumber: 186235783Skib s = s00; 187235783Skib sign = 0; 188235783Skib case STRTOG_Zero: 189235783Skib break; 190235783Skib default: 191235783Skib if (bb) { 192235783Skib copybits(bits, fpi.nbits, bb); 193235783Skib Bfree(bb); 194235783Skib } 195235783Skib ULtod(((U*)&rv)->L, bits, exp, i); 196235783Skib }} 197235783Skib goto ret; 198235783Skib } 199235783Skib } 200235783Skib#endif /*}*/ 201235783Skib nz0 = 1; 202235783Skib while(*++s == '0') ; 203235783Skib if (!*s) 204235783Skib goto ret; 205235783Skib } 206235783Skib s0 = s; 207235783Skib y = z = 0; 208235783Skib for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 209235783Skib if (nd < 9) 210235783Skib y = 10*y + c - '0'; 211235783Skib else if (nd < 16) 212235783Skib z = 10*z + c - '0'; 213235783Skib nd0 = nd; 214235783Skib#ifdef USE_LOCALE 215235783Skib if (c == *decimalpoint) { 216235783Skib for(i = 1; decimalpoint[i]; ++i) 217235783Skib if (s[i] != decimalpoint[i]) 218235783Skib goto dig_done; 219235783Skib s += i; 220235783Skib c = *s; 221235783Skib#else 222235783Skib if (c == '.') { 223235783Skib c = *++s; 224235783Skib#endif 225235783Skib decpt = 1; 226235783Skib if (!nd) { 227235783Skib for(; c == '0'; c = *++s) 228235783Skib nz++; 229235783Skib if (c > '0' && c <= '9') { 230235783Skib s0 = s; 231235783Skib nf += nz; 232235783Skib nz = 0; 233235783Skib goto have_dig; 234235783Skib } 235235783Skib goto dig_done; 236235783Skib } 237235783Skib for(; c >= '0' && c <= '9'; c = *++s) { 238235783Skib have_dig: 239235783Skib nz++; 240235783Skib if (c -= '0') { 241235783Skib nf += nz; 242235783Skib for(i = 1; i < nz; i++) 243235783Skib if (nd++ < 9) 244235783Skib y *= 10; 245235783Skib else if (nd <= DBL_DIG + 1) 246235783Skib z *= 10; 247235783Skib if (nd++ < 9) 248235783Skib y = 10*y + c; 249235783Skib else if (nd <= DBL_DIG + 1) 250235783Skib z = 10*z + c; 251235783Skib nz = 0; 252235783Skib } 253235783Skib } 254235783Skib }/*}*/ 255235783Skib dig_done: 256235783Skib e = 0; 257235783Skib if (c == 'e' || c == 'E') { 258235783Skib if (!nd && !nz && !nz0) { 259235783Skib goto ret0; 260235783Skib } 261235783Skib s00 = s; 262235783Skib esign = 0; 263235783Skib switch(c = *++s) { 264235783Skib case '-': 265235783Skib esign = 1; 266235783Skib case '+': 267235783Skib c = *++s; 268235783Skib } 269235783Skib if (c >= '0' && c <= '9') { 270235783Skib while(c == '0') 271235783Skib c = *++s; 272235783Skib if (c > '0' && c <= '9') { 273235783Skib L = c - '0'; 274235783Skib s1 = s; 275235783Skib while((c = *++s) >= '0' && c <= '9') 276235783Skib L = 10*L + c - '0'; 277235783Skib if (s - s1 > 8 || L > 19999) 278235783Skib /* Avoid confusion from exponents 279235783Skib * so large that e might overflow. 280235783Skib */ 281235783Skib e = 19999; /* safe for 16 bit ints */ 282235783Skib else 283235783Skib e = (int)L; 284235783Skib if (esign) 285235783Skib e = -e; 286235783Skib } 287235783Skib else 288235783Skib e = 0; 289235783Skib } 290235783Skib else 291235783Skib s = s00; 292235783Skib } 293235783Skib if (!nd) { 294235783Skib if (!nz && !nz0) { 295235783Skib#ifdef INFNAN_CHECK 296235783Skib /* Check for Nan and Infinity */ 297235783Skib ULong bits[2]; 298235783Skib static FPI fpinan = /* only 52 explicit bits */ 299235783Skib { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 300235783Skib if (!decpt) 301235783Skib switch(c) { 302235783Skib case 'i': 303235783Skib case 'I': 304235783Skib if (match(&s,"nf")) { 305235783Skib --s; 306235783Skib if (!match(&s,"inity")) 307235783Skib ++s; 308235783Skib word0(&rv) = 0x7ff00000; 309235783Skib word1(&rv) = 0; 310235783Skib goto ret; 311235783Skib } 312235783Skib break; 313235783Skib case 'n': 314235783Skib case 'N': 315235783Skib if (match(&s, "an")) { 316235783Skib#ifndef No_Hex_NaN 317235783Skib if (*s == '(' /*)*/ 318235783Skib && hexnan(&s, &fpinan, bits) 319235783Skib == STRTOG_NaNbits) { 320235783Skib word0(&rv) = 0x7ff80000 | bits[1]; 321235783Skib word1(&rv) = bits[0]; 322235783Skib } 323235783Skib else { 324235783Skib#endif 325235783Skib word0(&rv) = NAN_WORD0; 326235783Skib word1(&rv) = NAN_WORD1; 327235783Skib#ifndef No_Hex_NaN 328235783Skib } 329235783Skib#endif 330235783Skib goto ret; 331235783Skib } 332235783Skib } 333235783Skib#endif /* INFNAN_CHECK */ 334235783Skib ret0: 335235783Skib s = s00; 336235783Skib sign = 0; 337235783Skib } 338235783Skib goto ret; 339235783Skib } 340235783Skib e1 = e -= nf; 341235783Skib 342235783Skib /* Now we have nd0 digits, starting at s0, followed by a 343235783Skib * decimal point, followed by nd-nd0 digits. The number we're 344235783Skib * after is the integer represented by those digits times 345235783Skib * 10**e */ 346235783Skib 347235783Skib if (!nd0) 348235783Skib nd0 = nd; 349235783Skib k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 350235783Skib dval(&rv) = y; 351235783Skib if (k > 9) { 352235783Skib#ifdef SET_INEXACT 353235783Skib if (k > DBL_DIG) 354235783Skib oldinexact = get_inexact(); 355235783Skib#endif 356235783Skib dval(&rv) = tens[k - 9] * dval(&rv) + z; 357235783Skib } 358235783Skib bd0 = 0; 359235783Skib if (nd <= DBL_DIG 360235783Skib#ifndef RND_PRODQUOT 361235783Skib#ifndef Honor_FLT_ROUNDS 362235783Skib && Flt_Rounds == 1 363235783Skib#endif 364235783Skib#endif 365235783Skib ) { 366235783Skib if (!e) 367235783Skib goto ret; 368235783Skib#ifndef ROUND_BIASED_without_Round_Up 369235783Skib if (e > 0) { 370235783Skib if (e <= Ten_pmax) { 371235783Skib#ifdef VAX 372235783Skib goto vax_ovfl_check; 373235783Skib#else 374235783Skib#ifdef Honor_FLT_ROUNDS 375235783Skib /* round correctly FLT_ROUNDS = 2 or 3 */ 376235783Skib if (sign) { 377235783Skib rv.d = -rv.d; 378235783Skib sign = 0; 379235783Skib } 380235783Skib#endif 381235783Skib /* rv = */ rounded_product(dval(&rv), tens[e]); 382235783Skib goto ret; 383235783Skib#endif 384235783Skib } 385235783Skib i = DBL_DIG - nd; 386235783Skib if (e <= Ten_pmax + i) { 387235783Skib /* A fancier test would sometimes let us do 388235783Skib * this for larger i values. 389235783Skib */ 390235783Skib#ifdef Honor_FLT_ROUNDS 391235783Skib /* round correctly FLT_ROUNDS = 2 or 3 */ 392235783Skib if (sign) { 393235783Skib rv.d = -rv.d; 394235783Skib sign = 0; 395235783Skib } 396235783Skib#endif 397235783Skib e -= i; 398235783Skib dval(&rv) *= tens[i]; 399235783Skib#ifdef VAX 400235783Skib /* VAX exponent range is so narrow we must 401235783Skib * worry about overflow here... 402235783Skib */ 403235783Skib vax_ovfl_check: 404235783Skib word0(&rv) -= P*Exp_msk1; 405235783Skib /* rv = */ rounded_product(dval(&rv), tens[e]); 406235783Skib if ((word0(&rv) & Exp_mask) 407235783Skib > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 408235783Skib goto ovfl; 409235783Skib word0(&rv) += P*Exp_msk1; 410235783Skib#else 411235783Skib /* rv = */ rounded_product(dval(&rv), tens[e]); 412235783Skib#endif 413235783Skib goto ret; 414235783Skib } 415235783Skib } 416235783Skib#ifndef Inaccurate_Divide 417235783Skib else if (e >= -Ten_pmax) { 418235783Skib#ifdef Honor_FLT_ROUNDS 419235783Skib /* round correctly FLT_ROUNDS = 2 or 3 */ 420235783Skib if (sign) { 421235783Skib rv.d = -rv.d; 422235783Skib sign = 0; 423235783Skib } 424235783Skib#endif 425235783Skib /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 426235783Skib goto ret; 427235783Skib } 428235783Skib#endif 429235783Skib#endif /* ROUND_BIASED_without_Round_Up */ 430235783Skib } 431235783Skib e1 += nd - k; 432235783Skib 433235783Skib#ifdef IEEE_Arith 434235783Skib#ifdef SET_INEXACT 435235783Skib inexact = 1; 436235783Skib if (k <= DBL_DIG) 437235783Skib oldinexact = get_inexact(); 438235783Skib#endif 439235783Skib#ifdef Avoid_Underflow 440235783Skib scale = 0; 441235783Skib#endif 442235783Skib#ifdef Honor_FLT_ROUNDS 443235783Skib if (Rounding >= 2) { 444235783Skib if (sign) 445235783Skib Rounding = Rounding == 2 ? 0 : 2; 446235783Skib else 447235783Skib if (Rounding != 2) 448235783Skib Rounding = 0; 449235783Skib } 450235783Skib#endif 451235783Skib#endif /*IEEE_Arith*/ 452235783Skib 453235783Skib /* Get starting approximation = rv * 10**e1 */ 454235783Skib 455235783Skib if (e1 > 0) { 456235783Skib if ( (i = e1 & 15) !=0) 457235783Skib dval(&rv) *= tens[i]; 458235783Skib if (e1 &= ~15) { 459235783Skib if (e1 > DBL_MAX_10_EXP) { 460235783Skib ovfl: 461235783Skib /* Can't trust HUGE_VAL */ 462235783Skib#ifdef IEEE_Arith 463235783Skib#ifdef Honor_FLT_ROUNDS 464235783Skib switch(Rounding) { 465235783Skib case 0: /* toward 0 */ 466235783Skib case 3: /* toward -infinity */ 467235783Skib word0(&rv) = Big0; 468235783Skib word1(&rv) = Big1; 469235783Skib break; 470235783Skib default: 471235783Skib word0(&rv) = Exp_mask; 472235783Skib word1(&rv) = 0; 473235783Skib } 474235783Skib#else /*Honor_FLT_ROUNDS*/ 475235783Skib word0(&rv) = Exp_mask; 476235783Skib word1(&rv) = 0; 477235783Skib#endif /*Honor_FLT_ROUNDS*/ 478235783Skib#ifdef SET_INEXACT 479235783Skib /* set overflow bit */ 480235783Skib dval(&rv0) = 1e300; 481235783Skib dval(&rv0) *= dval(&rv0); 482235783Skib#endif 483235783Skib#else /*IEEE_Arith*/ 484235783Skib word0(&rv) = Big0; 485235783Skib word1(&rv) = Big1; 486235783Skib#endif /*IEEE_Arith*/ 487235783Skib range_err: 488235783Skib if (bd0) { 489235783Skib Bfree(bb); 490235783Skib Bfree(bd); 491235783Skib Bfree(bs); 492235783Skib Bfree(bd0); 493235783Skib Bfree(delta); 494235783Skib } 495235783Skib#ifndef NO_ERRNO 496235783Skib errno = ERANGE; 497235783Skib#endif 498235783Skib goto ret; 499235783Skib } 500235783Skib e1 >>= 4; 501235783Skib for(j = 0; e1 > 1; j++, e1 >>= 1) 502235783Skib if (e1 & 1) 503235783Skib dval(&rv) *= bigtens[j]; 504235783Skib /* The last multiplication could overflow. */ 505235783Skib word0(&rv) -= P*Exp_msk1; 506235783Skib dval(&rv) *= bigtens[j]; 507235783Skib if ((z = word0(&rv) & Exp_mask) 508235783Skib > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 509235783Skib goto ovfl; 510235783Skib if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 511235783Skib /* set to largest number */ 512235783Skib /* (Can't trust DBL_MAX) */ 513235783Skib word0(&rv) = Big0; 514235783Skib word1(&rv) = Big1; 515235783Skib } 516235783Skib else 517235783Skib word0(&rv) += P*Exp_msk1; 518235783Skib } 519235783Skib } 520235783Skib else if (e1 < 0) { 521235783Skib e1 = -e1; 522235783Skib if ( (i = e1 & 15) !=0) 523235783Skib dval(&rv) /= tens[i]; 524235783Skib if (e1 >>= 4) { 525235783Skib if (e1 >= 1 << n_bigtens) 526235783Skib goto undfl; 527235783Skib#ifdef Avoid_Underflow 528235783Skib if (e1 & Scale_Bit) 529235783Skib scale = 2*P; 530235783Skib for(j = 0; e1 > 0; j++, e1 >>= 1) 531235783Skib if (e1 & 1) 532235783Skib dval(&rv) *= tinytens[j]; 533235783Skib if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 534235783Skib >> Exp_shift)) > 0) { 535235783Skib /* scaled rv is denormal; zap j low bits */ 536235783Skib if (j >= 32) { 537235783Skib word1(&rv) = 0; 538235783Skib if (j >= 53) 539235783Skib word0(&rv) = (P+2)*Exp_msk1; 540235783Skib else 541235783Skib word0(&rv) &= 0xffffffff << (j-32); 542235783Skib } 543235783Skib else 544235783Skib word1(&rv) &= 0xffffffff << j; 545235783Skib } 546235783Skib#else 547235783Skib for(j = 0; e1 > 1; j++, e1 >>= 1) 548235783Skib if (e1 & 1) 549235783Skib dval(&rv) *= tinytens[j]; 550235783Skib /* The last multiplication could underflow. */ 551235783Skib dval(&rv0) = dval(&rv); 552235783Skib dval(&rv) *= tinytens[j]; 553235783Skib if (!dval(&rv)) { 554235783Skib dval(&rv) = 2.*dval(&rv0); 555235783Skib dval(&rv) *= tinytens[j]; 556235783Skib#endif 557235783Skib if (!dval(&rv)) { 558235783Skib undfl: 559235783Skib dval(&rv) = 0.; 560235783Skib goto range_err; 561235783Skib } 562235783Skib#ifndef Avoid_Underflow 563235783Skib word0(&rv) = Tiny0; 564235783Skib word1(&rv) = Tiny1; 565235783Skib /* The refinement below will clean 566235783Skib * this approximation up. 567235783Skib */ 568235783Skib } 569235783Skib#endif 570235783Skib } 571235783Skib } 572235783Skib 573235783Skib /* Now the hard part -- adjusting rv to the correct value.*/ 574235783Skib 575235783Skib /* Put digits into bd: true value = bd * 10^e */ 576235783Skib 577235783Skib bd0 = s2b(s0, nd0, nd, y, dplen); 578235783Skib 579235783Skib for(;;) { 580235783Skib bd = Balloc(bd0->k); 581235783Skib Bcopy(bd, bd0); 582235783Skib bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 583235783Skib bs = i2b(1); 584235783Skib 585235783Skib if (e >= 0) { 586235783Skib bb2 = bb5 = 0; 587235783Skib bd2 = bd5 = e; 588235783Skib } 589235783Skib else { 590235783Skib bb2 = bb5 = -e; 591235783Skib bd2 = bd5 = 0; 592235783Skib } 593235783Skib if (bbe >= 0) 594235783Skib bb2 += bbe; 595235783Skib else 596235783Skib bd2 -= bbe; 597235783Skib bs2 = bb2; 598235783Skib#ifdef Honor_FLT_ROUNDS 599235783Skib if (Rounding != 1) 600235783Skib bs2++; 601235783Skib#endif 602235783Skib#ifdef Avoid_Underflow 603235783Skib Lsb = LSB; 604235783Skib Lsb1 = 0; 605235783Skib j = bbe - scale; 606235783Skib i = j + bbbits - 1; /* logb(rv) */ 607235783Skib j = P + 1 - bbbits; 608235783Skib if (i < Emin) { /* denormal */ 609235783Skib i = Emin - i; 610235783Skib j -= i; 611235783Skib if (i < 32) 612235783Skib Lsb <<= i; 613235783Skib else 614235783Skib Lsb1 = Lsb << (i-32); 615235783Skib } 616235783Skib#else /*Avoid_Underflow*/ 617235783Skib#ifdef Sudden_Underflow 618235783Skib#ifdef IBM 619235783Skib j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 620235783Skib#else 621 j = P + 1 - bbbits; 622#endif 623#else /*Sudden_Underflow*/ 624 j = bbe; 625 i = j + bbbits - 1; /* logb(&rv) */ 626 if (i < Emin) /* denormal */ 627 j += P - Emin; 628 else 629 j = P + 1 - bbbits; 630#endif /*Sudden_Underflow*/ 631#endif /*Avoid_Underflow*/ 632 bb2 += j; 633 bd2 += j; 634#ifdef Avoid_Underflow 635 bd2 += scale; 636#endif 637 i = bb2 < bd2 ? bb2 : bd2; 638 if (i > bs2) 639 i = bs2; 640 if (i > 0) { 641 bb2 -= i; 642 bd2 -= i; 643 bs2 -= i; 644 } 645 if (bb5 > 0) { 646 bs = pow5mult(bs, bb5); 647 bb1 = mult(bs, bb); 648 Bfree(bb); 649 bb = bb1; 650 } 651 if (bb2 > 0) 652 bb = lshift(bb, bb2); 653 if (bd5 > 0) 654 bd = pow5mult(bd, bd5); 655 if (bd2 > 0) 656 bd = lshift(bd, bd2); 657 if (bs2 > 0) 658 bs = lshift(bs, bs2); 659 delta = diff(bb, bd); 660 dsign = delta->sign; 661 delta->sign = 0; 662 i = cmp(delta, bs); 663#ifdef Honor_FLT_ROUNDS 664 if (Rounding != 1) { 665 if (i < 0) { 666 /* Error is less than an ulp */ 667 if (!delta->x[0] && delta->wds <= 1) { 668 /* exact */ 669#ifdef SET_INEXACT 670 inexact = 0; 671#endif 672 break; 673 } 674 if (Rounding) { 675 if (dsign) { 676 dval(&adj) = 1.; 677 goto apply_adj; 678 } 679 } 680 else if (!dsign) { 681 dval(&adj) = -1.; 682 if (!word1(&rv) 683 && !(word0(&rv) & Frac_mask)) { 684 y = word0(&rv) & Exp_mask; 685#ifdef Avoid_Underflow 686 if (!scale || y > 2*P*Exp_msk1) 687#else 688 if (y) 689#endif 690 { 691 delta = lshift(delta,Log2P); 692 if (cmp(delta, bs) <= 0) 693 dval(&adj) = -0.5; 694 } 695 } 696 apply_adj: 697#ifdef Avoid_Underflow 698 if (scale && (y = word0(&rv) & Exp_mask) 699 <= 2*P*Exp_msk1) 700 word0(&adj) += (2*P+1)*Exp_msk1 - y; 701#else 702#ifdef Sudden_Underflow 703 if ((word0(&rv) & Exp_mask) <= 704 P*Exp_msk1) { 705 word0(&rv) += P*Exp_msk1; 706 dval(&rv) += adj*ulp(&rv); 707 word0(&rv) -= P*Exp_msk1; 708 } 709 else 710#endif /*Sudden_Underflow*/ 711#endif /*Avoid_Underflow*/ 712 dval(&rv) += adj.d*ulp(&rv); 713 } 714 break; 715 } 716 dval(&adj) = ratio(delta, bs); 717 if (adj.d < 1.) 718 dval(&adj) = 1.; 719 if (adj.d <= 0x7ffffffe) { 720 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */ 721 y = adj.d; 722 if (y != adj.d) { 723 if (!((Rounding>>1) ^ dsign)) 724 y++; 725 dval(&adj) = y; 726 } 727 } 728#ifdef Avoid_Underflow 729 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 730 word0(&adj) += (2*P+1)*Exp_msk1 - y; 731#else 732#ifdef Sudden_Underflow 733 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 734 word0(&rv) += P*Exp_msk1; 735 dval(&adj) *= ulp(&rv); 736 if (dsign) 737 dval(&rv) += adj; 738 else 739 dval(&rv) -= adj; 740 word0(&rv) -= P*Exp_msk1; 741 goto cont; 742 } 743#endif /*Sudden_Underflow*/ 744#endif /*Avoid_Underflow*/ 745 dval(&adj) *= ulp(&rv); 746 if (dsign) { 747 if (word0(&rv) == Big0 && word1(&rv) == Big1) 748 goto ovfl; 749 dval(&rv) += adj.d; 750 } 751 else 752 dval(&rv) -= adj.d; 753 goto cont; 754 } 755#endif /*Honor_FLT_ROUNDS*/ 756 757 if (i < 0) { 758 /* Error is less than half an ulp -- check for 759 * special case of mantissa a power of two. 760 */ 761 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask 762#ifdef IEEE_Arith 763#ifdef Avoid_Underflow 764 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 765#else 766 || (word0(&rv) & Exp_mask) <= Exp_msk1 767#endif 768#endif 769 ) { 770#ifdef SET_INEXACT 771 if (!delta->x[0] && delta->wds <= 1) 772 inexact = 0; 773#endif 774 break; 775 } 776 if (!delta->x[0] && delta->wds <= 1) { 777 /* exact result */ 778#ifdef SET_INEXACT 779 inexact = 0; 780#endif 781 break; 782 } 783 delta = lshift(delta,Log2P); 784 if (cmp(delta, bs) > 0) 785 goto drop_down; 786 break; 787 } 788 if (i == 0) { 789 /* exactly half-way between */ 790 if (dsign) { 791 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 792 && word1(&rv) == ( 793#ifdef Avoid_Underflow 794 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 795 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 796#endif 797 0xffffffff)) { 798 /*boundary case -- increment exponent*/ 799 if (word0(&rv) == Big0 && word1(&rv) == Big1) 800 goto ovfl; 801 word0(&rv) = (word0(&rv) & Exp_mask) 802 + Exp_msk1 803#ifdef IBM 804 | Exp_msk1 >> 4 805#endif 806 ; 807 word1(&rv) = 0; 808#ifdef Avoid_Underflow 809 dsign = 0; 810#endif 811 break; 812 } 813 } 814 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 815 drop_down: 816 /* boundary case -- decrement exponent */ 817#ifdef Sudden_Underflow /*{{*/ 818 L = word0(&rv) & Exp_mask; 819#ifdef IBM 820 if (L < Exp_msk1) 821#else 822#ifdef Avoid_Underflow 823 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 824#else 825 if (L <= Exp_msk1) 826#endif /*Avoid_Underflow*/ 827#endif /*IBM*/ 828 goto undfl; 829 L -= Exp_msk1; 830#else /*Sudden_Underflow}{*/ 831#ifdef Avoid_Underflow 832 if (scale) { 833 L = word0(&rv) & Exp_mask; 834 if (L <= (2*P+1)*Exp_msk1) { 835 if (L > (P+2)*Exp_msk1) 836 /* round even ==> */ 837 /* accept rv */ 838 break; 839 /* rv = smallest denormal */ 840 goto undfl; 841 } 842 } 843#endif /*Avoid_Underflow*/ 844 L = (word0(&rv) & Exp_mask) - Exp_msk1; 845#endif /*Sudden_Underflow}}*/ 846 word0(&rv) = L | Bndry_mask1; 847 word1(&rv) = 0xffffffff; 848#ifdef IBM 849 goto cont; 850#else 851 break; 852#endif 853 } 854#ifndef ROUND_BIASED 855#ifdef Avoid_Underflow 856 if (Lsb1) { 857 if (!(word0(&rv) & Lsb1)) 858 break; 859 } 860 else if (!(word1(&rv) & Lsb)) 861 break; 862#else 863 if (!(word1(&rv) & LSB)) 864 break; 865#endif 866#endif 867 if (dsign) 868#ifdef Avoid_Underflow 869 dval(&rv) += sulp(&rv, scale); 870#else 871 dval(&rv) += ulp(&rv); 872#endif 873#ifndef ROUND_BIASED 874 else { 875#ifdef Avoid_Underflow 876 dval(&rv) -= sulp(&rv, scale); 877#else 878 dval(&rv) -= ulp(&rv); 879#endif 880#ifndef Sudden_Underflow 881 if (!dval(&rv)) 882 goto undfl; 883#endif 884 } 885#ifdef Avoid_Underflow 886 dsign = 1 - dsign; 887#endif 888#endif 889 break; 890 } 891 if ((aadj = ratio(delta, bs)) <= 2.) { 892 if (dsign) 893 aadj = dval(&aadj1) = 1.; 894 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 895#ifndef Sudden_Underflow 896 if (word1(&rv) == Tiny1 && !word0(&rv)) 897 goto undfl; 898#endif 899 aadj = 1.; 900 dval(&aadj1) = -1.; 901 } 902 else { 903 /* special case -- power of FLT_RADIX to be */ 904 /* rounded down... */ 905 906 if (aadj < 2./FLT_RADIX) 907 aadj = 1./FLT_RADIX; 908 else 909 aadj *= 0.5; 910 dval(&aadj1) = -aadj; 911 } 912 } 913 else { 914 aadj *= 0.5; 915 dval(&aadj1) = dsign ? aadj : -aadj; 916#ifdef Check_FLT_ROUNDS 917 switch(Rounding) { 918 case 2: /* towards +infinity */ 919 dval(&aadj1) -= 0.5; 920 break; 921 case 0: /* towards 0 */ 922 case 3: /* towards -infinity */ 923 dval(&aadj1) += 0.5; 924 } 925#else 926 if (Flt_Rounds == 0) 927 dval(&aadj1) += 0.5; 928#endif /*Check_FLT_ROUNDS*/ 929 } 930 y = word0(&rv) & Exp_mask; 931 932 /* Check for overflow */ 933 934 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 935 dval(&rv0) = dval(&rv); 936 word0(&rv) -= P*Exp_msk1; 937 dval(&adj) = dval(&aadj1) * ulp(&rv); 938 dval(&rv) += dval(&adj); 939 if ((word0(&rv) & Exp_mask) >= 940 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 941 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 942 goto ovfl; 943 word0(&rv) = Big0; 944 word1(&rv) = Big1; 945 goto cont; 946 } 947 else 948 word0(&rv) += P*Exp_msk1; 949 } 950 else { 951#ifdef Avoid_Underflow 952 if (scale && y <= 2*P*Exp_msk1) { 953 if (aadj <= 0x7fffffff) { 954 if ((z = aadj) <= 0) 955 z = 1; 956 aadj = z; 957 dval(&aadj1) = dsign ? aadj : -aadj; 958 } 959 word0(&aadj1) += (2*P+1)*Exp_msk1 - y; 960 } 961 dval(&adj) = dval(&aadj1) * ulp(&rv); 962 dval(&rv) += dval(&adj); 963#else 964#ifdef Sudden_Underflow 965 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 966 dval(&rv0) = dval(&rv); 967 word0(&rv) += P*Exp_msk1; 968 dval(&adj) = dval(&aadj1) * ulp(&rv); 969 dval(&rv) += adj; 970#ifdef IBM 971 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) 972#else 973 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) 974#endif 975 { 976 if (word0(&rv0) == Tiny0 977 && word1(&rv0) == Tiny1) 978 goto undfl; 979 word0(&rv) = Tiny0; 980 word1(&rv) = Tiny1; 981 goto cont; 982 } 983 else 984 word0(&rv) -= P*Exp_msk1; 985 } 986 else { 987 dval(&adj) = dval(&aadj1) * ulp(&rv); 988 dval(&rv) += adj; 989 } 990#else /*Sudden_Underflow*/ 991 /* Compute dval(&adj) so that the IEEE rounding rules will 992 * correctly round rv + dval(&adj) in some half-way cases. 993 * If rv * ulp(&rv) is denormalized (i.e., 994 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 995 * trouble from bits lost to denormalization; 996 * example: 1.2e-307 . 997 */ 998 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 999 dval(&aadj1) = (double)(int)(aadj + 0.5); 1000 if (!dsign) 1001 dval(&aadj1) = -dval(&aadj1); 1002 } 1003 dval(&adj) = dval(&aadj1) * ulp(&rv); 1004 dval(&rv) += adj; 1005#endif /*Sudden_Underflow*/ 1006#endif /*Avoid_Underflow*/ 1007 } 1008 z = word0(&rv) & Exp_mask; 1009#ifndef SET_INEXACT 1010#ifdef Avoid_Underflow 1011 if (!scale) 1012#endif 1013 if (y == z) { 1014 /* Can we stop now? */ 1015 L = (Long)aadj; 1016 aadj -= L; 1017 /* The tolerances below are conservative. */ 1018 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 1019 if (aadj < .4999999 || aadj > .5000001) 1020 break; 1021 } 1022 else if (aadj < .4999999/FLT_RADIX) 1023 break; 1024 } 1025#endif 1026 cont: 1027 Bfree(bb); 1028 Bfree(bd); 1029 Bfree(bs); 1030 Bfree(delta); 1031 } 1032 Bfree(bb); 1033 Bfree(bd); 1034 Bfree(bs); 1035 Bfree(bd0); 1036 Bfree(delta); 1037#ifdef SET_INEXACT 1038 if (inexact) { 1039 if (!oldinexact) { 1040 word0(&rv0) = Exp_1 + (70 << Exp_shift); 1041 word1(&rv0) = 0; 1042 dval(&rv0) += 1.; 1043 } 1044 } 1045 else if (!oldinexact) 1046 clear_inexact(); 1047#endif 1048#ifdef Avoid_Underflow 1049 if (scale) { 1050 word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 1051 word1(&rv0) = 0; 1052 dval(&rv) *= dval(&rv0); 1053#ifndef NO_ERRNO 1054 /* try to avoid the bug of testing an 8087 register value */ 1055#ifdef IEEE_Arith 1056 if (!(word0(&rv) & Exp_mask)) 1057#else 1058 if (word0(&rv) == 0 && word1(&rv) == 0) 1059#endif 1060 errno = ERANGE; 1061#endif 1062 } 1063#endif /* Avoid_Underflow */ 1064#ifdef SET_INEXACT 1065 if (inexact && !(word0(&rv) & Exp_mask)) { 1066 /* set underflow bit */ 1067 dval(&rv0) = 1e-300; 1068 dval(&rv0) *= dval(&rv0); 1069 } 1070#endif 1071 ret: 1072 if (se) 1073 *se = (char *)s; 1074 return sign ? -dval(&rv) : dval(&rv); 1075 } 1076 1077