strtod.c revision 182709
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 182709 2008-09-03 07:23:57Z 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 47182709Sdas/* The factor of 2^106 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, 50182709Sdas 9007199254740992.*9007199254740992.e-256 51112158Sdas }; 52112158Sdas#endif 53112158Sdas#endif 54112158Sdas 55112158Sdas#ifdef Honor_FLT_ROUNDS 56112158Sdas#undef Check_FLT_ROUNDS 57112158Sdas#define Check_FLT_ROUNDS 58112158Sdas#else 59112158Sdas#define Rounding Flt_Rounds 60112158Sdas#endif 61112158Sdas 62112158Sdas double 63112158Sdasstrtod 64112158Sdas#ifdef KR_headers 65112158Sdas (s00, se) CONST char *s00; char **se; 66112158Sdas#else 67112158Sdas (CONST char *s00, char **se) 68112158Sdas#endif 69112158Sdas{ 70112158Sdas#ifdef Avoid_Underflow 71112158Sdas int scale; 72112158Sdas#endif 73165743Sdas int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, 74112158Sdas e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 75112158Sdas CONST char *s, *s0, *s1; 76112158Sdas double aadj, aadj1, adj, rv, rv0; 77112158Sdas Long L; 78112158Sdas ULong y, z; 79112158Sdas Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 80112158Sdas#ifdef SET_INEXACT 81112158Sdas int inexact, oldinexact; 82112158Sdas#endif 83182709Sdas#ifdef Honor_FLT_ROUNDS /*{*/ 84182709Sdas int Rounding; 85182709Sdas#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 86182709Sdas Rounding = Flt_Rounds; 87182709Sdas#else /*}{*/ 88182709Sdas Rounding = 1; 89182709Sdas switch(fegetround()) { 90182709Sdas case FE_TOWARDZERO: Rounding = 0; break; 91182709Sdas case FE_UPWARD: Rounding = 2; break; 92182709Sdas case FE_DOWNWARD: Rounding = 3; 93182709Sdas } 94182709Sdas#endif /*}}*/ 95182709Sdas#endif /*}*/ 96112158Sdas 97165743Sdas sign = nz0 = nz = decpt = 0; 98112158Sdas dval(rv) = 0.; 99112158Sdas for(s = s00;;s++) switch(*s) { 100112158Sdas case '-': 101112158Sdas sign = 1; 102112158Sdas /* no break */ 103112158Sdas case '+': 104112158Sdas if (*++s) 105112158Sdas goto break2; 106112158Sdas /* no break */ 107112158Sdas case 0: 108112158Sdas goto ret0; 109112158Sdas case '\t': 110112158Sdas case '\n': 111112158Sdas case '\v': 112112158Sdas case '\f': 113112158Sdas case '\r': 114112158Sdas case ' ': 115112158Sdas continue; 116112158Sdas default: 117112158Sdas goto break2; 118112158Sdas } 119112158Sdas break2: 120112158Sdas if (*s == '0') { 121182709Sdas#ifndef NO_HEX_FP /*{{*/ 122112158Sdas { 123112158Sdas static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 124112158Sdas Long exp; 125112158Sdas ULong bits[2]; 126112158Sdas switch(s[1]) { 127112158Sdas case 'x': 128112158Sdas case 'X': 129165743Sdas { 130182709Sdas#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/ 131165743Sdas FPI fpi1 = fpi; 132182709Sdas#ifdef Honor_FLT_ROUNDS /*{{*/ 133182709Sdas fpi1.rounding = Rounding; 134182709Sdas#else /*}{*/ 135165743Sdas switch(fegetround()) { 136165743Sdas case FE_TOWARDZERO: fpi1.rounding = 0; break; 137165743Sdas case FE_UPWARD: fpi1.rounding = 2; break; 138165743Sdas case FE_DOWNWARD: fpi1.rounding = 3; 139165743Sdas } 140182709Sdas#endif /*}}*/ 141182709Sdas#else /*}{*/ 142165743Sdas#define fpi1 fpi 143182709Sdas#endif /*}}*/ 144165743Sdas switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { 145112158Sdas case STRTOG_NoNumber: 146112158Sdas s = s00; 147112158Sdas sign = 0; 148112158Sdas case STRTOG_Zero: 149112158Sdas break; 150112158Sdas default: 151124703Sdas if (bb) { 152124703Sdas copybits(bits, fpi.nbits, bb); 153124703Sdas Bfree(bb); 154124703Sdas } 155112158Sdas ULtod(((U*)&rv)->L, bits, exp, i); 156165743Sdas }} 157112158Sdas goto ret; 158112158Sdas } 159112158Sdas } 160112158Sdas#endif 161112158Sdas nz0 = 1; 162112158Sdas while(*++s == '0') ; 163112158Sdas if (!*s) 164112158Sdas goto ret; 165112158Sdas } 166112158Sdas s0 = s; 167112158Sdas y = z = 0; 168112158Sdas for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 169112158Sdas if (nd < 9) 170112158Sdas y = 10*y + c - '0'; 171112158Sdas else if (nd < 16) 172112158Sdas z = 10*z + c - '0'; 173112158Sdas nd0 = nd; 174112158Sdas#ifdef USE_LOCALE 175112415Sdas if (c == *localeconv()->decimal_point) 176112415Sdas#else 177112415Sdas if (c == '.') 178112158Sdas#endif 179112415Sdas { 180165743Sdas decpt = 1; 181112158Sdas c = *++s; 182112158Sdas if (!nd) { 183112158Sdas for(; c == '0'; c = *++s) 184112158Sdas nz++; 185112158Sdas if (c > '0' && c <= '9') { 186112158Sdas s0 = s; 187112158Sdas nf += nz; 188112158Sdas nz = 0; 189112158Sdas goto have_dig; 190112158Sdas } 191112158Sdas goto dig_done; 192112158Sdas } 193112158Sdas for(; c >= '0' && c <= '9'; c = *++s) { 194112158Sdas have_dig: 195112158Sdas nz++; 196112158Sdas if (c -= '0') { 197112158Sdas nf += nz; 198112158Sdas for(i = 1; i < nz; i++) 199112158Sdas if (nd++ < 9) 200112158Sdas y *= 10; 201112158Sdas else if (nd <= DBL_DIG + 1) 202112158Sdas z *= 10; 203112158Sdas if (nd++ < 9) 204112158Sdas y = 10*y + c; 205112158Sdas else if (nd <= DBL_DIG + 1) 206112158Sdas z = 10*z + c; 207112158Sdas nz = 0; 208112158Sdas } 209112158Sdas } 210112158Sdas } 211112158Sdas dig_done: 212112158Sdas e = 0; 213112158Sdas if (c == 'e' || c == 'E') { 214112158Sdas if (!nd && !nz && !nz0) { 215112158Sdas goto ret0; 216112158Sdas } 217112158Sdas s00 = s; 218112158Sdas esign = 0; 219112158Sdas switch(c = *++s) { 220112158Sdas case '-': 221112158Sdas esign = 1; 222112158Sdas case '+': 223112158Sdas c = *++s; 224112158Sdas } 225112158Sdas if (c >= '0' && c <= '9') { 226112158Sdas while(c == '0') 227112158Sdas c = *++s; 228112158Sdas if (c > '0' && c <= '9') { 229112158Sdas L = c - '0'; 230112158Sdas s1 = s; 231112158Sdas while((c = *++s) >= '0' && c <= '9') 232112158Sdas L = 10*L + c - '0'; 233112158Sdas if (s - s1 > 8 || L > 19999) 234112158Sdas /* Avoid confusion from exponents 235112158Sdas * so large that e might overflow. 236112158Sdas */ 237112158Sdas e = 19999; /* safe for 16 bit ints */ 238112158Sdas else 239112158Sdas e = (int)L; 240112158Sdas if (esign) 241112158Sdas e = -e; 242112158Sdas } 243112158Sdas else 244112158Sdas e = 0; 245112158Sdas } 246112158Sdas else 247112158Sdas s = s00; 248112158Sdas } 249112158Sdas if (!nd) { 250112158Sdas if (!nz && !nz0) { 251112158Sdas#ifdef INFNAN_CHECK 252112158Sdas /* Check for Nan and Infinity */ 253112158Sdas ULong bits[2]; 254112158Sdas static FPI fpinan = /* only 52 explicit bits */ 255112158Sdas { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 256165743Sdas if (!decpt) 257165743Sdas switch(c) { 258112158Sdas case 'i': 259112158Sdas case 'I': 260112158Sdas if (match(&s,"nf")) { 261112158Sdas --s; 262112158Sdas if (!match(&s,"inity")) 263112158Sdas ++s; 264112158Sdas word0(rv) = 0x7ff00000; 265112158Sdas word1(rv) = 0; 266112158Sdas goto ret; 267112158Sdas } 268112158Sdas break; 269112158Sdas case 'n': 270112158Sdas case 'N': 271112158Sdas if (match(&s, "an")) { 272112158Sdas#ifndef No_Hex_NaN 273112158Sdas if (*s == '(' /*)*/ 274112158Sdas && hexnan(&s, &fpinan, bits) 275112158Sdas == STRTOG_NaNbits) { 276174679Sdas word0(rv) = 0x7ff80000 | bits[1]; 277112158Sdas word1(rv) = bits[0]; 278112158Sdas } 279112158Sdas else { 280165743Sdas#endif 281112158Sdas word0(rv) = NAN_WORD0; 282112158Sdas word1(rv) = NAN_WORD1; 283165743Sdas#ifndef No_Hex_NaN 284112158Sdas } 285112158Sdas#endif 286112158Sdas goto ret; 287112158Sdas } 288112158Sdas } 289112158Sdas#endif /* INFNAN_CHECK */ 290112158Sdas ret0: 291112158Sdas s = s00; 292112158Sdas sign = 0; 293112158Sdas } 294112158Sdas goto ret; 295112158Sdas } 296112158Sdas e1 = e -= nf; 297112158Sdas 298112158Sdas /* Now we have nd0 digits, starting at s0, followed by a 299112158Sdas * decimal point, followed by nd-nd0 digits. The number we're 300112158Sdas * after is the integer represented by those digits times 301112158Sdas * 10**e */ 302112158Sdas 303112158Sdas if (!nd0) 304112158Sdas nd0 = nd; 305112158Sdas k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 306112158Sdas dval(rv) = y; 307112158Sdas if (k > 9) { 308112158Sdas#ifdef SET_INEXACT 309112158Sdas if (k > DBL_DIG) 310112158Sdas oldinexact = get_inexact(); 311112158Sdas#endif 312112158Sdas dval(rv) = tens[k - 9] * dval(rv) + z; 313112158Sdas } 314112158Sdas bd0 = 0; 315112158Sdas if (nd <= DBL_DIG 316112158Sdas#ifndef RND_PRODQUOT 317112158Sdas#ifndef Honor_FLT_ROUNDS 318112158Sdas && Flt_Rounds == 1 319112158Sdas#endif 320112158Sdas#endif 321112158Sdas ) { 322112158Sdas if (!e) 323112158Sdas goto ret; 324112158Sdas if (e > 0) { 325112158Sdas if (e <= Ten_pmax) { 326112158Sdas#ifdef VAX 327112158Sdas goto vax_ovfl_check; 328112158Sdas#else 329112158Sdas#ifdef Honor_FLT_ROUNDS 330112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 331112158Sdas if (sign) { 332112158Sdas rv = -rv; 333112158Sdas sign = 0; 334112158Sdas } 335112158Sdas#endif 336112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 337112158Sdas goto ret; 338112158Sdas#endif 339112158Sdas } 340112158Sdas i = DBL_DIG - nd; 341112158Sdas if (e <= Ten_pmax + i) { 342112158Sdas /* A fancier test would sometimes let us do 343112158Sdas * this for larger i values. 344112158Sdas */ 345112158Sdas#ifdef Honor_FLT_ROUNDS 346112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 347112158Sdas if (sign) { 348112158Sdas rv = -rv; 349112158Sdas sign = 0; 350112158Sdas } 351112158Sdas#endif 352112158Sdas e -= i; 353112158Sdas dval(rv) *= tens[i]; 354112158Sdas#ifdef VAX 355112158Sdas /* VAX exponent range is so narrow we must 356112158Sdas * worry about overflow here... 357112158Sdas */ 358112158Sdas vax_ovfl_check: 359112158Sdas word0(rv) -= P*Exp_msk1; 360112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 361112158Sdas if ((word0(rv) & Exp_mask) 362112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 363112158Sdas goto ovfl; 364112158Sdas word0(rv) += P*Exp_msk1; 365112158Sdas#else 366112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 367112158Sdas#endif 368112158Sdas goto ret; 369112158Sdas } 370112158Sdas } 371112158Sdas#ifndef Inaccurate_Divide 372112158Sdas else if (e >= -Ten_pmax) { 373112158Sdas#ifdef Honor_FLT_ROUNDS 374112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 375112158Sdas if (sign) { 376112158Sdas rv = -rv; 377112158Sdas sign = 0; 378112158Sdas } 379112158Sdas#endif 380112158Sdas /* rv = */ rounded_quotient(dval(rv), tens[-e]); 381112158Sdas goto ret; 382112158Sdas } 383112158Sdas#endif 384112158Sdas } 385112158Sdas e1 += nd - k; 386112158Sdas 387112158Sdas#ifdef IEEE_Arith 388112158Sdas#ifdef SET_INEXACT 389112158Sdas inexact = 1; 390112158Sdas if (k <= DBL_DIG) 391112158Sdas oldinexact = get_inexact(); 392112158Sdas#endif 393112158Sdas#ifdef Avoid_Underflow 394112158Sdas scale = 0; 395112158Sdas#endif 396112158Sdas#ifdef Honor_FLT_ROUNDS 397182709Sdas if (Rounding >= 2) { 398112158Sdas if (sign) 399182709Sdas Rounding = Rounding == 2 ? 0 : 2; 400112158Sdas else 401182709Sdas if (Rounding != 2) 402182709Sdas Rounding = 0; 403112158Sdas } 404112158Sdas#endif 405112158Sdas#endif /*IEEE_Arith*/ 406112158Sdas 407112158Sdas /* Get starting approximation = rv * 10**e1 */ 408112158Sdas 409112158Sdas if (e1 > 0) { 410112158Sdas if ( (i = e1 & 15) !=0) 411112158Sdas dval(rv) *= tens[i]; 412112158Sdas if (e1 &= ~15) { 413112158Sdas if (e1 > DBL_MAX_10_EXP) { 414112158Sdas ovfl: 415112158Sdas#ifndef NO_ERRNO 416112158Sdas errno = ERANGE; 417112158Sdas#endif 418112158Sdas /* Can't trust HUGE_VAL */ 419112158Sdas#ifdef IEEE_Arith 420112158Sdas#ifdef Honor_FLT_ROUNDS 421182709Sdas switch(Rounding) { 422112158Sdas case 0: /* toward 0 */ 423112158Sdas case 3: /* toward -infinity */ 424112158Sdas word0(rv) = Big0; 425112158Sdas word1(rv) = Big1; 426112158Sdas break; 427112158Sdas default: 428112158Sdas word0(rv) = Exp_mask; 429112158Sdas word1(rv) = 0; 430112158Sdas } 431112158Sdas#else /*Honor_FLT_ROUNDS*/ 432112158Sdas word0(rv) = Exp_mask; 433112158Sdas word1(rv) = 0; 434112158Sdas#endif /*Honor_FLT_ROUNDS*/ 435112158Sdas#ifdef SET_INEXACT 436112158Sdas /* set overflow bit */ 437112158Sdas dval(rv0) = 1e300; 438112158Sdas dval(rv0) *= dval(rv0); 439112158Sdas#endif 440112158Sdas#else /*IEEE_Arith*/ 441112158Sdas word0(rv) = Big0; 442112158Sdas word1(rv) = Big1; 443112158Sdas#endif /*IEEE_Arith*/ 444112158Sdas if (bd0) 445112158Sdas goto retfree; 446112158Sdas goto ret; 447112158Sdas } 448112158Sdas e1 >>= 4; 449112158Sdas for(j = 0; e1 > 1; j++, e1 >>= 1) 450112158Sdas if (e1 & 1) 451112158Sdas dval(rv) *= bigtens[j]; 452112158Sdas /* The last multiplication could overflow. */ 453112158Sdas word0(rv) -= P*Exp_msk1; 454112158Sdas dval(rv) *= bigtens[j]; 455112158Sdas if ((z = word0(rv) & Exp_mask) 456112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 457112158Sdas goto ovfl; 458112158Sdas if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 459112158Sdas /* set to largest number */ 460112158Sdas /* (Can't trust DBL_MAX) */ 461112158Sdas word0(rv) = Big0; 462112158Sdas word1(rv) = Big1; 463112158Sdas } 464112158Sdas else 465112158Sdas word0(rv) += P*Exp_msk1; 466112158Sdas } 467112158Sdas } 468112158Sdas else if (e1 < 0) { 469112158Sdas e1 = -e1; 470112158Sdas if ( (i = e1 & 15) !=0) 471112158Sdas dval(rv) /= tens[i]; 472112158Sdas if (e1 >>= 4) { 473112158Sdas if (e1 >= 1 << n_bigtens) 474112158Sdas goto undfl; 475112158Sdas#ifdef Avoid_Underflow 476112158Sdas if (e1 & Scale_Bit) 477112158Sdas scale = 2*P; 478112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 479112158Sdas if (e1 & 1) 480112158Sdas dval(rv) *= tinytens[j]; 481112158Sdas if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 482112158Sdas >> Exp_shift)) > 0) { 483112158Sdas /* scaled rv is denormal; zap j low bits */ 484112158Sdas if (j >= 32) { 485112158Sdas word1(rv) = 0; 486112158Sdas if (j >= 53) 487112158Sdas word0(rv) = (P+2)*Exp_msk1; 488112158Sdas else 489112158Sdas word0(rv) &= 0xffffffff << j-32; 490112158Sdas } 491112158Sdas else 492112158Sdas word1(rv) &= 0xffffffff << j; 493112158Sdas } 494112158Sdas#else 495112158Sdas for(j = 0; e1 > 1; j++, e1 >>= 1) 496112158Sdas if (e1 & 1) 497112158Sdas dval(rv) *= tinytens[j]; 498112158Sdas /* The last multiplication could underflow. */ 499112158Sdas dval(rv0) = dval(rv); 500112158Sdas dval(rv) *= tinytens[j]; 501112158Sdas if (!dval(rv)) { 502112158Sdas dval(rv) = 2.*dval(rv0); 503112158Sdas dval(rv) *= tinytens[j]; 504112158Sdas#endif 505112158Sdas if (!dval(rv)) { 506112158Sdas undfl: 507112158Sdas dval(rv) = 0.; 508112158Sdas#ifndef NO_ERRNO 509112158Sdas errno = ERANGE; 510112158Sdas#endif 511112158Sdas if (bd0) 512112158Sdas goto retfree; 513112158Sdas goto ret; 514112158Sdas } 515112158Sdas#ifndef Avoid_Underflow 516112158Sdas word0(rv) = Tiny0; 517112158Sdas word1(rv) = Tiny1; 518112158Sdas /* The refinement below will clean 519112158Sdas * this approximation up. 520112158Sdas */ 521112158Sdas } 522112158Sdas#endif 523112158Sdas } 524112158Sdas } 525112158Sdas 526112158Sdas /* Now the hard part -- adjusting rv to the correct value.*/ 527112158Sdas 528112158Sdas /* Put digits into bd: true value = bd * 10^e */ 529112158Sdas 530112158Sdas bd0 = s2b(s0, nd0, nd, y); 531112158Sdas 532112158Sdas for(;;) { 533112158Sdas bd = Balloc(bd0->k); 534112158Sdas Bcopy(bd, bd0); 535112158Sdas bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 536112158Sdas bs = i2b(1); 537112158Sdas 538112158Sdas if (e >= 0) { 539112158Sdas bb2 = bb5 = 0; 540112158Sdas bd2 = bd5 = e; 541112158Sdas } 542112158Sdas else { 543112158Sdas bb2 = bb5 = -e; 544112158Sdas bd2 = bd5 = 0; 545112158Sdas } 546112158Sdas if (bbe >= 0) 547112158Sdas bb2 += bbe; 548112158Sdas else 549112158Sdas bd2 -= bbe; 550112158Sdas bs2 = bb2; 551112158Sdas#ifdef Honor_FLT_ROUNDS 552182709Sdas if (Rounding != 1) 553112158Sdas bs2++; 554112158Sdas#endif 555112158Sdas#ifdef Avoid_Underflow 556112158Sdas j = bbe - scale; 557112158Sdas i = j + bbbits - 1; /* logb(rv) */ 558112158Sdas if (i < Emin) /* denormal */ 559112158Sdas j += P - Emin; 560112158Sdas else 561112158Sdas j = P + 1 - bbbits; 562112158Sdas#else /*Avoid_Underflow*/ 563112158Sdas#ifdef Sudden_Underflow 564112158Sdas#ifdef IBM 565112158Sdas j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 566112158Sdas#else 567112158Sdas j = P + 1 - bbbits; 568112158Sdas#endif 569112158Sdas#else /*Sudden_Underflow*/ 570112158Sdas j = bbe; 571112158Sdas i = j + bbbits - 1; /* logb(rv) */ 572112158Sdas if (i < Emin) /* denormal */ 573112158Sdas j += P - Emin; 574112158Sdas else 575112158Sdas j = P + 1 - bbbits; 576112158Sdas#endif /*Sudden_Underflow*/ 577112158Sdas#endif /*Avoid_Underflow*/ 578112158Sdas bb2 += j; 579112158Sdas bd2 += j; 580112158Sdas#ifdef Avoid_Underflow 581112158Sdas bd2 += scale; 582112158Sdas#endif 583112158Sdas i = bb2 < bd2 ? bb2 : bd2; 584112158Sdas if (i > bs2) 585112158Sdas i = bs2; 586112158Sdas if (i > 0) { 587112158Sdas bb2 -= i; 588112158Sdas bd2 -= i; 589112158Sdas bs2 -= i; 590112158Sdas } 591112158Sdas if (bb5 > 0) { 592112158Sdas bs = pow5mult(bs, bb5); 593112158Sdas bb1 = mult(bs, bb); 594112158Sdas Bfree(bb); 595112158Sdas bb = bb1; 596112158Sdas } 597112158Sdas if (bb2 > 0) 598112158Sdas bb = lshift(bb, bb2); 599112158Sdas if (bd5 > 0) 600112158Sdas bd = pow5mult(bd, bd5); 601112158Sdas if (bd2 > 0) 602112158Sdas bd = lshift(bd, bd2); 603112158Sdas if (bs2 > 0) 604112158Sdas bs = lshift(bs, bs2); 605112158Sdas delta = diff(bb, bd); 606112158Sdas dsign = delta->sign; 607112158Sdas delta->sign = 0; 608112158Sdas i = cmp(delta, bs); 609112158Sdas#ifdef Honor_FLT_ROUNDS 610182709Sdas if (Rounding != 1) { 611112158Sdas if (i < 0) { 612112158Sdas /* Error is less than an ulp */ 613112158Sdas if (!delta->x[0] && delta->wds <= 1) { 614112158Sdas /* exact */ 615112158Sdas#ifdef SET_INEXACT 616112158Sdas inexact = 0; 617112158Sdas#endif 618112158Sdas break; 619112158Sdas } 620182709Sdas if (Rounding) { 621112158Sdas if (dsign) { 622112158Sdas adj = 1.; 623112158Sdas goto apply_adj; 624112158Sdas } 625112158Sdas } 626112158Sdas else if (!dsign) { 627112158Sdas adj = -1.; 628112158Sdas if (!word1(rv) 629112158Sdas && !(word0(rv) & Frac_mask)) { 630112158Sdas y = word0(rv) & Exp_mask; 631112158Sdas#ifdef Avoid_Underflow 632112158Sdas if (!scale || y > 2*P*Exp_msk1) 633112158Sdas#else 634112158Sdas if (y) 635112158Sdas#endif 636112158Sdas { 637112158Sdas delta = lshift(delta,Log2P); 638112158Sdas if (cmp(delta, bs) <= 0) 639112158Sdas adj = -0.5; 640112158Sdas } 641112158Sdas } 642112158Sdas apply_adj: 643112158Sdas#ifdef Avoid_Underflow 644112158Sdas if (scale && (y = word0(rv) & Exp_mask) 645112158Sdas <= 2*P*Exp_msk1) 646112158Sdas word0(adj) += (2*P+1)*Exp_msk1 - y; 647112158Sdas#else 648112158Sdas#ifdef Sudden_Underflow 649112158Sdas if ((word0(rv) & Exp_mask) <= 650112158Sdas P*Exp_msk1) { 651112158Sdas word0(rv) += P*Exp_msk1; 652112158Sdas dval(rv) += adj*ulp(dval(rv)); 653112158Sdas word0(rv) -= P*Exp_msk1; 654112158Sdas } 655112158Sdas else 656112158Sdas#endif /*Sudden_Underflow*/ 657112158Sdas#endif /*Avoid_Underflow*/ 658112158Sdas dval(rv) += adj*ulp(dval(rv)); 659112158Sdas } 660112158Sdas break; 661112158Sdas } 662112158Sdas adj = ratio(delta, bs); 663112158Sdas if (adj < 1.) 664112158Sdas adj = 1.; 665112158Sdas if (adj <= 0x7ffffffe) { 666182709Sdas /* adj = Rounding ? ceil(adj) : floor(adj); */ 667112158Sdas y = adj; 668112158Sdas if (y != adj) { 669182709Sdas if (!((Rounding>>1) ^ dsign)) 670112158Sdas y++; 671112158Sdas adj = y; 672112158Sdas } 673112158Sdas } 674112158Sdas#ifdef Avoid_Underflow 675112158Sdas if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 676112158Sdas word0(adj) += (2*P+1)*Exp_msk1 - y; 677112158Sdas#else 678112158Sdas#ifdef Sudden_Underflow 679112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 680112158Sdas word0(rv) += P*Exp_msk1; 681112158Sdas adj *= ulp(dval(rv)); 682112158Sdas if (dsign) 683112158Sdas dval(rv) += adj; 684112158Sdas else 685112158Sdas dval(rv) -= adj; 686112158Sdas word0(rv) -= P*Exp_msk1; 687112158Sdas goto cont; 688112158Sdas } 689112158Sdas#endif /*Sudden_Underflow*/ 690112158Sdas#endif /*Avoid_Underflow*/ 691112158Sdas adj *= ulp(dval(rv)); 692182709Sdas if (dsign) { 693182709Sdas if (word0(rv) == Big0 && word1(rv) == Big1) 694182709Sdas goto ovfl; 695112158Sdas dval(rv) += adj; 696182709Sdas } 697112158Sdas else 698112158Sdas dval(rv) -= adj; 699112158Sdas goto cont; 700112158Sdas } 701112158Sdas#endif /*Honor_FLT_ROUNDS*/ 702112158Sdas 703112158Sdas if (i < 0) { 704112158Sdas /* Error is less than half an ulp -- check for 705112158Sdas * special case of mantissa a power of two. 706112158Sdas */ 707112158Sdas if (dsign || word1(rv) || word0(rv) & Bndry_mask 708112158Sdas#ifdef IEEE_Arith 709112158Sdas#ifdef Avoid_Underflow 710112158Sdas || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 711112158Sdas#else 712112158Sdas || (word0(rv) & Exp_mask) <= Exp_msk1 713112158Sdas#endif 714112158Sdas#endif 715112158Sdas ) { 716112158Sdas#ifdef SET_INEXACT 717112158Sdas if (!delta->x[0] && delta->wds <= 1) 718112158Sdas inexact = 0; 719112158Sdas#endif 720112158Sdas break; 721112158Sdas } 722112158Sdas if (!delta->x[0] && delta->wds <= 1) { 723112158Sdas /* exact result */ 724112158Sdas#ifdef SET_INEXACT 725112158Sdas inexact = 0; 726112158Sdas#endif 727112158Sdas break; 728112158Sdas } 729112158Sdas delta = lshift(delta,Log2P); 730112158Sdas if (cmp(delta, bs) > 0) 731112158Sdas goto drop_down; 732112158Sdas break; 733112158Sdas } 734112158Sdas if (i == 0) { 735112158Sdas /* exactly half-way between */ 736112158Sdas if (dsign) { 737112158Sdas if ((word0(rv) & Bndry_mask1) == Bndry_mask1 738112158Sdas && word1(rv) == ( 739112158Sdas#ifdef Avoid_Underflow 740112158Sdas (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 741112158Sdas ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 742112158Sdas#endif 743112158Sdas 0xffffffff)) { 744112158Sdas /*boundary case -- increment exponent*/ 745112158Sdas word0(rv) = (word0(rv) & Exp_mask) 746112158Sdas + Exp_msk1 747112158Sdas#ifdef IBM 748112158Sdas | Exp_msk1 >> 4 749112158Sdas#endif 750112158Sdas ; 751112158Sdas word1(rv) = 0; 752112158Sdas#ifdef Avoid_Underflow 753112158Sdas dsign = 0; 754112158Sdas#endif 755112158Sdas break; 756112158Sdas } 757112158Sdas } 758112158Sdas else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 759112158Sdas drop_down: 760112158Sdas /* boundary case -- decrement exponent */ 761112158Sdas#ifdef Sudden_Underflow /*{{*/ 762112158Sdas L = word0(rv) & Exp_mask; 763112158Sdas#ifdef IBM 764112158Sdas if (L < Exp_msk1) 765112158Sdas#else 766112158Sdas#ifdef Avoid_Underflow 767112158Sdas if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 768112158Sdas#else 769112158Sdas if (L <= Exp_msk1) 770112158Sdas#endif /*Avoid_Underflow*/ 771112158Sdas#endif /*IBM*/ 772112158Sdas goto undfl; 773112158Sdas L -= Exp_msk1; 774112158Sdas#else /*Sudden_Underflow}{*/ 775112158Sdas#ifdef Avoid_Underflow 776112158Sdas if (scale) { 777112158Sdas L = word0(rv) & Exp_mask; 778112158Sdas if (L <= (2*P+1)*Exp_msk1) { 779112158Sdas if (L > (P+2)*Exp_msk1) 780112158Sdas /* round even ==> */ 781112158Sdas /* accept rv */ 782112158Sdas break; 783112158Sdas /* rv = smallest denormal */ 784112158Sdas goto undfl; 785112158Sdas } 786112158Sdas } 787112158Sdas#endif /*Avoid_Underflow*/ 788112158Sdas L = (word0(rv) & Exp_mask) - Exp_msk1; 789182709Sdas#endif /*Sudden_Underflow}}*/ 790112158Sdas word0(rv) = L | Bndry_mask1; 791112158Sdas word1(rv) = 0xffffffff; 792112158Sdas#ifdef IBM 793112158Sdas goto cont; 794112158Sdas#else 795112158Sdas break; 796112158Sdas#endif 797112158Sdas } 798112158Sdas#ifndef ROUND_BIASED 799112158Sdas if (!(word1(rv) & LSB)) 800112158Sdas break; 801112158Sdas#endif 802112158Sdas if (dsign) 803112158Sdas dval(rv) += ulp(dval(rv)); 804112158Sdas#ifndef ROUND_BIASED 805112158Sdas else { 806112158Sdas dval(rv) -= ulp(dval(rv)); 807112158Sdas#ifndef Sudden_Underflow 808112158Sdas if (!dval(rv)) 809112158Sdas goto undfl; 810112158Sdas#endif 811112158Sdas } 812112158Sdas#ifdef Avoid_Underflow 813112158Sdas dsign = 1 - dsign; 814112158Sdas#endif 815112158Sdas#endif 816112158Sdas break; 817112158Sdas } 818112158Sdas if ((aadj = ratio(delta, bs)) <= 2.) { 819112158Sdas if (dsign) 820112158Sdas aadj = aadj1 = 1.; 821112158Sdas else if (word1(rv) || word0(rv) & Bndry_mask) { 822112158Sdas#ifndef Sudden_Underflow 823112158Sdas if (word1(rv) == Tiny1 && !word0(rv)) 824112158Sdas goto undfl; 825112158Sdas#endif 826112158Sdas aadj = 1.; 827112158Sdas aadj1 = -1.; 828112158Sdas } 829112158Sdas else { 830112158Sdas /* special case -- power of FLT_RADIX to be */ 831112158Sdas /* rounded down... */ 832112158Sdas 833112158Sdas if (aadj < 2./FLT_RADIX) 834112158Sdas aadj = 1./FLT_RADIX; 835112158Sdas else 836112158Sdas aadj *= 0.5; 837112158Sdas aadj1 = -aadj; 838112158Sdas } 839112158Sdas } 840112158Sdas else { 841112158Sdas aadj *= 0.5; 842112158Sdas aadj1 = dsign ? aadj : -aadj; 843112158Sdas#ifdef Check_FLT_ROUNDS 844112158Sdas switch(Rounding) { 845112158Sdas case 2: /* towards +infinity */ 846112158Sdas aadj1 -= 0.5; 847112158Sdas break; 848112158Sdas case 0: /* towards 0 */ 849112158Sdas case 3: /* towards -infinity */ 850112158Sdas aadj1 += 0.5; 851112158Sdas } 852112158Sdas#else 853112158Sdas if (Flt_Rounds == 0) 854112158Sdas aadj1 += 0.5; 855112158Sdas#endif /*Check_FLT_ROUNDS*/ 856112158Sdas } 857112158Sdas y = word0(rv) & Exp_mask; 858112158Sdas 859112158Sdas /* Check for overflow */ 860112158Sdas 861112158Sdas if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 862112158Sdas dval(rv0) = dval(rv); 863112158Sdas word0(rv) -= P*Exp_msk1; 864112158Sdas adj = aadj1 * ulp(dval(rv)); 865112158Sdas dval(rv) += adj; 866112158Sdas if ((word0(rv) & Exp_mask) >= 867112158Sdas Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 868112158Sdas if (word0(rv0) == Big0 && word1(rv0) == Big1) 869112158Sdas goto ovfl; 870112158Sdas word0(rv) = Big0; 871112158Sdas word1(rv) = Big1; 872112158Sdas goto cont; 873112158Sdas } 874112158Sdas else 875112158Sdas word0(rv) += P*Exp_msk1; 876112158Sdas } 877112158Sdas else { 878112158Sdas#ifdef Avoid_Underflow 879112158Sdas if (scale && y <= 2*P*Exp_msk1) { 880112158Sdas if (aadj <= 0x7fffffff) { 881112158Sdas if ((z = aadj) <= 0) 882112158Sdas z = 1; 883112158Sdas aadj = z; 884112158Sdas aadj1 = dsign ? aadj : -aadj; 885112158Sdas } 886112158Sdas word0(aadj1) += (2*P+1)*Exp_msk1 - y; 887112158Sdas } 888112158Sdas adj = aadj1 * ulp(dval(rv)); 889112158Sdas dval(rv) += adj; 890112158Sdas#else 891112158Sdas#ifdef Sudden_Underflow 892112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 893112158Sdas dval(rv0) = dval(rv); 894112158Sdas word0(rv) += P*Exp_msk1; 895112158Sdas adj = aadj1 * ulp(dval(rv)); 896112158Sdas dval(rv) += adj; 897112158Sdas#ifdef IBM 898112158Sdas if ((word0(rv) & Exp_mask) < P*Exp_msk1) 899112158Sdas#else 900112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 901112158Sdas#endif 902112158Sdas { 903112158Sdas if (word0(rv0) == Tiny0 904112158Sdas && word1(rv0) == Tiny1) 905112158Sdas goto undfl; 906112158Sdas word0(rv) = Tiny0; 907112158Sdas word1(rv) = Tiny1; 908112158Sdas goto cont; 909112158Sdas } 910112158Sdas else 911112158Sdas word0(rv) -= P*Exp_msk1; 912112158Sdas } 913112158Sdas else { 914112158Sdas adj = aadj1 * ulp(dval(rv)); 915112158Sdas dval(rv) += adj; 916112158Sdas } 917112158Sdas#else /*Sudden_Underflow*/ 918112158Sdas /* Compute adj so that the IEEE rounding rules will 919112158Sdas * correctly round rv + adj in some half-way cases. 920112158Sdas * If rv * ulp(rv) is denormalized (i.e., 921112158Sdas * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 922112158Sdas * trouble from bits lost to denormalization; 923112158Sdas * example: 1.2e-307 . 924112158Sdas */ 925112158Sdas if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 926112158Sdas aadj1 = (double)(int)(aadj + 0.5); 927112158Sdas if (!dsign) 928112158Sdas aadj1 = -aadj1; 929112158Sdas } 930112158Sdas adj = aadj1 * ulp(dval(rv)); 931112158Sdas dval(rv) += adj; 932112158Sdas#endif /*Sudden_Underflow*/ 933112158Sdas#endif /*Avoid_Underflow*/ 934112158Sdas } 935112158Sdas z = word0(rv) & Exp_mask; 936112158Sdas#ifndef SET_INEXACT 937112158Sdas#ifdef Avoid_Underflow 938112158Sdas if (!scale) 939112158Sdas#endif 940112158Sdas if (y == z) { 941112158Sdas /* Can we stop now? */ 942112158Sdas L = (Long)aadj; 943112158Sdas aadj -= L; 944112158Sdas /* The tolerances below are conservative. */ 945112158Sdas if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 946112158Sdas if (aadj < .4999999 || aadj > .5000001) 947112158Sdas break; 948112158Sdas } 949112158Sdas else if (aadj < .4999999/FLT_RADIX) 950112158Sdas break; 951112158Sdas } 952112158Sdas#endif 953112158Sdas cont: 954112158Sdas Bfree(bb); 955112158Sdas Bfree(bd); 956112158Sdas Bfree(bs); 957112158Sdas Bfree(delta); 958112158Sdas } 959112158Sdas#ifdef SET_INEXACT 960112158Sdas if (inexact) { 961112158Sdas if (!oldinexact) { 962112158Sdas word0(rv0) = Exp_1 + (70 << Exp_shift); 963112158Sdas word1(rv0) = 0; 964112158Sdas dval(rv0) += 1.; 965112158Sdas } 966112158Sdas } 967112158Sdas else if (!oldinexact) 968112158Sdas clear_inexact(); 969112158Sdas#endif 970112158Sdas#ifdef Avoid_Underflow 971112158Sdas if (scale) { 972112158Sdas word0(rv0) = Exp_1 - 2*P*Exp_msk1; 973112158Sdas word1(rv0) = 0; 974112158Sdas dval(rv) *= dval(rv0); 975112158Sdas#ifndef NO_ERRNO 976112158Sdas /* try to avoid the bug of testing an 8087 register value */ 977112158Sdas if (word0(rv) == 0 && word1(rv) == 0) 978112158Sdas errno = ERANGE; 979112158Sdas#endif 980112158Sdas } 981112158Sdas#endif /* Avoid_Underflow */ 982112158Sdas#ifdef SET_INEXACT 983112158Sdas if (inexact && !(word0(rv) & Exp_mask)) { 984112158Sdas /* set underflow bit */ 985112158Sdas dval(rv0) = 1e-300; 986112158Sdas dval(rv0) *= dval(rv0); 987112158Sdas } 988112158Sdas#endif 989112158Sdas retfree: 990112158Sdas Bfree(bb); 991112158Sdas Bfree(bd); 992112158Sdas Bfree(bs); 993112158Sdas Bfree(bd0); 994112158Sdas Bfree(delta); 995112158Sdas ret: 996112158Sdas if (se) 997112158Sdas *se = (char *)s; 998112158Sdas return sign ? -dval(rv) : dval(rv); 999112158Sdas } 1000112158Sdas 1001