strtod.c revision 187808
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 187808 2009-01-28 04:36:34Z 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 83187808Sdas#ifdef USE_LOCALE /*{{*/ 84187808Sdas#ifdef NO_LOCALE_CACHE 85187808Sdas char *decimalpoint = localeconv()->decimal_point; 86187808Sdas int dplen = strlen(decimalpoint); 87187808Sdas#else 88187808Sdas char *decimalpoint; 89187808Sdas static char *decimalpoint_cache; 90187808Sdas static int dplen; 91187808Sdas if (!(s0 = decimalpoint_cache)) { 92187808Sdas s0 = localeconv()->decimal_point; 93187808Sdas if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) { 94187808Sdas strcpy(decimalpoint_cache, s0); 95187808Sdas s0 = decimalpoint_cache; 96187808Sdas } 97187808Sdas dplen = strlen(s0); 98187808Sdas } 99187808Sdas decimalpoint = (char*)s0; 100187808Sdas#endif /*NO_LOCALE_CACHE*/ 101187808Sdas#else /*USE_LOCALE}{*/ 102187808Sdas#define dplen 1 103187808Sdas#endif /*USE_LOCALE}}*/ 104187808Sdas 105182709Sdas#ifdef Honor_FLT_ROUNDS /*{*/ 106182709Sdas int Rounding; 107182709Sdas#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 108182709Sdas Rounding = Flt_Rounds; 109182709Sdas#else /*}{*/ 110182709Sdas Rounding = 1; 111182709Sdas switch(fegetround()) { 112182709Sdas case FE_TOWARDZERO: Rounding = 0; break; 113182709Sdas case FE_UPWARD: Rounding = 2; break; 114182709Sdas case FE_DOWNWARD: Rounding = 3; 115182709Sdas } 116182709Sdas#endif /*}}*/ 117182709Sdas#endif /*}*/ 118112158Sdas 119165743Sdas sign = nz0 = nz = decpt = 0; 120112158Sdas dval(rv) = 0.; 121112158Sdas for(s = s00;;s++) switch(*s) { 122112158Sdas case '-': 123112158Sdas sign = 1; 124112158Sdas /* no break */ 125112158Sdas case '+': 126112158Sdas if (*++s) 127112158Sdas goto break2; 128112158Sdas /* no break */ 129112158Sdas case 0: 130112158Sdas goto ret0; 131112158Sdas case '\t': 132112158Sdas case '\n': 133112158Sdas case '\v': 134112158Sdas case '\f': 135112158Sdas case '\r': 136112158Sdas case ' ': 137112158Sdas continue; 138112158Sdas default: 139112158Sdas goto break2; 140112158Sdas } 141112158Sdas break2: 142112158Sdas if (*s == '0') { 143187808Sdas#ifndef NO_HEX_FP /*{*/ 144112158Sdas { 145112158Sdas static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 146112158Sdas Long exp; 147112158Sdas ULong bits[2]; 148112158Sdas switch(s[1]) { 149112158Sdas case 'x': 150112158Sdas case 'X': 151165743Sdas { 152182709Sdas#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/ 153165743Sdas FPI fpi1 = fpi; 154182709Sdas#ifdef Honor_FLT_ROUNDS /*{{*/ 155182709Sdas fpi1.rounding = Rounding; 156182709Sdas#else /*}{*/ 157165743Sdas switch(fegetround()) { 158165743Sdas case FE_TOWARDZERO: fpi1.rounding = 0; break; 159165743Sdas case FE_UPWARD: fpi1.rounding = 2; break; 160165743Sdas case FE_DOWNWARD: fpi1.rounding = 3; 161165743Sdas } 162182709Sdas#endif /*}}*/ 163182709Sdas#else /*}{*/ 164165743Sdas#define fpi1 fpi 165182709Sdas#endif /*}}*/ 166165743Sdas switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { 167112158Sdas case STRTOG_NoNumber: 168112158Sdas s = s00; 169112158Sdas sign = 0; 170112158Sdas case STRTOG_Zero: 171112158Sdas break; 172112158Sdas default: 173124703Sdas if (bb) { 174124703Sdas copybits(bits, fpi.nbits, bb); 175124703Sdas Bfree(bb); 176124703Sdas } 177112158Sdas ULtod(((U*)&rv)->L, bits, exp, i); 178165743Sdas }} 179112158Sdas goto ret; 180112158Sdas } 181112158Sdas } 182187808Sdas#endif /*}*/ 183112158Sdas nz0 = 1; 184112158Sdas while(*++s == '0') ; 185112158Sdas if (!*s) 186112158Sdas goto ret; 187112158Sdas } 188112158Sdas s0 = s; 189112158Sdas y = z = 0; 190112158Sdas for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 191112158Sdas if (nd < 9) 192112158Sdas y = 10*y + c - '0'; 193112158Sdas else if (nd < 16) 194112158Sdas z = 10*z + c - '0'; 195112158Sdas nd0 = nd; 196112158Sdas#ifdef USE_LOCALE 197187808Sdas if (c == *decimalpoint) { 198187808Sdas for(i = 1; decimalpoint[i]; ++i) 199187808Sdas if (s[i] != decimalpoint[i]) 200187808Sdas goto dig_done; 201187808Sdas s += i; 202187808Sdas c = *s; 203112415Sdas#else 204187808Sdas if (c == '.') { 205187808Sdas c = *++s; 206112158Sdas#endif 207165743Sdas decpt = 1; 208112158Sdas if (!nd) { 209112158Sdas for(; c == '0'; c = *++s) 210112158Sdas nz++; 211112158Sdas if (c > '0' && c <= '9') { 212112158Sdas s0 = s; 213112158Sdas nf += nz; 214112158Sdas nz = 0; 215112158Sdas goto have_dig; 216112158Sdas } 217112158Sdas goto dig_done; 218112158Sdas } 219112158Sdas for(; c >= '0' && c <= '9'; c = *++s) { 220112158Sdas have_dig: 221112158Sdas nz++; 222112158Sdas if (c -= '0') { 223112158Sdas nf += nz; 224112158Sdas for(i = 1; i < nz; i++) 225112158Sdas if (nd++ < 9) 226112158Sdas y *= 10; 227112158Sdas else if (nd <= DBL_DIG + 1) 228112158Sdas z *= 10; 229112158Sdas if (nd++ < 9) 230112158Sdas y = 10*y + c; 231112158Sdas else if (nd <= DBL_DIG + 1) 232112158Sdas z = 10*z + c; 233112158Sdas nz = 0; 234112158Sdas } 235112158Sdas } 236187808Sdas }/*}*/ 237112158Sdas dig_done: 238112158Sdas e = 0; 239112158Sdas if (c == 'e' || c == 'E') { 240112158Sdas if (!nd && !nz && !nz0) { 241112158Sdas goto ret0; 242112158Sdas } 243112158Sdas s00 = s; 244112158Sdas esign = 0; 245112158Sdas switch(c = *++s) { 246112158Sdas case '-': 247112158Sdas esign = 1; 248112158Sdas case '+': 249112158Sdas c = *++s; 250112158Sdas } 251112158Sdas if (c >= '0' && c <= '9') { 252112158Sdas while(c == '0') 253112158Sdas c = *++s; 254112158Sdas if (c > '0' && c <= '9') { 255112158Sdas L = c - '0'; 256112158Sdas s1 = s; 257112158Sdas while((c = *++s) >= '0' && c <= '9') 258112158Sdas L = 10*L + c - '0'; 259112158Sdas if (s - s1 > 8 || L > 19999) 260112158Sdas /* Avoid confusion from exponents 261112158Sdas * so large that e might overflow. 262112158Sdas */ 263112158Sdas e = 19999; /* safe for 16 bit ints */ 264112158Sdas else 265112158Sdas e = (int)L; 266112158Sdas if (esign) 267112158Sdas e = -e; 268112158Sdas } 269112158Sdas else 270112158Sdas e = 0; 271112158Sdas } 272112158Sdas else 273112158Sdas s = s00; 274112158Sdas } 275112158Sdas if (!nd) { 276112158Sdas if (!nz && !nz0) { 277112158Sdas#ifdef INFNAN_CHECK 278112158Sdas /* Check for Nan and Infinity */ 279112158Sdas ULong bits[2]; 280112158Sdas static FPI fpinan = /* only 52 explicit bits */ 281112158Sdas { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 282165743Sdas if (!decpt) 283165743Sdas switch(c) { 284112158Sdas case 'i': 285112158Sdas case 'I': 286112158Sdas if (match(&s,"nf")) { 287112158Sdas --s; 288112158Sdas if (!match(&s,"inity")) 289112158Sdas ++s; 290112158Sdas word0(rv) = 0x7ff00000; 291112158Sdas word1(rv) = 0; 292112158Sdas goto ret; 293112158Sdas } 294112158Sdas break; 295112158Sdas case 'n': 296112158Sdas case 'N': 297112158Sdas if (match(&s, "an")) { 298112158Sdas#ifndef No_Hex_NaN 299112158Sdas if (*s == '(' /*)*/ 300112158Sdas && hexnan(&s, &fpinan, bits) 301112158Sdas == STRTOG_NaNbits) { 302174679Sdas word0(rv) = 0x7ff80000 | bits[1]; 303112158Sdas word1(rv) = bits[0]; 304112158Sdas } 305112158Sdas else { 306165743Sdas#endif 307112158Sdas word0(rv) = NAN_WORD0; 308112158Sdas word1(rv) = NAN_WORD1; 309165743Sdas#ifndef No_Hex_NaN 310112158Sdas } 311112158Sdas#endif 312112158Sdas goto ret; 313112158Sdas } 314112158Sdas } 315112158Sdas#endif /* INFNAN_CHECK */ 316112158Sdas ret0: 317112158Sdas s = s00; 318112158Sdas sign = 0; 319112158Sdas } 320112158Sdas goto ret; 321112158Sdas } 322112158Sdas e1 = e -= nf; 323112158Sdas 324112158Sdas /* Now we have nd0 digits, starting at s0, followed by a 325112158Sdas * decimal point, followed by nd-nd0 digits. The number we're 326112158Sdas * after is the integer represented by those digits times 327112158Sdas * 10**e */ 328112158Sdas 329112158Sdas if (!nd0) 330112158Sdas nd0 = nd; 331112158Sdas k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 332112158Sdas dval(rv) = y; 333112158Sdas if (k > 9) { 334112158Sdas#ifdef SET_INEXACT 335112158Sdas if (k > DBL_DIG) 336112158Sdas oldinexact = get_inexact(); 337112158Sdas#endif 338112158Sdas dval(rv) = tens[k - 9] * dval(rv) + z; 339112158Sdas } 340112158Sdas bd0 = 0; 341112158Sdas if (nd <= DBL_DIG 342112158Sdas#ifndef RND_PRODQUOT 343112158Sdas#ifndef Honor_FLT_ROUNDS 344112158Sdas && Flt_Rounds == 1 345112158Sdas#endif 346112158Sdas#endif 347112158Sdas ) { 348112158Sdas if (!e) 349112158Sdas goto ret; 350112158Sdas if (e > 0) { 351112158Sdas if (e <= Ten_pmax) { 352112158Sdas#ifdef VAX 353112158Sdas goto vax_ovfl_check; 354112158Sdas#else 355112158Sdas#ifdef Honor_FLT_ROUNDS 356112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 357112158Sdas if (sign) { 358112158Sdas rv = -rv; 359112158Sdas sign = 0; 360112158Sdas } 361112158Sdas#endif 362112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 363112158Sdas goto ret; 364112158Sdas#endif 365112158Sdas } 366112158Sdas i = DBL_DIG - nd; 367112158Sdas if (e <= Ten_pmax + i) { 368112158Sdas /* A fancier test would sometimes let us do 369112158Sdas * this for larger i values. 370112158Sdas */ 371112158Sdas#ifdef Honor_FLT_ROUNDS 372112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 373112158Sdas if (sign) { 374112158Sdas rv = -rv; 375112158Sdas sign = 0; 376112158Sdas } 377112158Sdas#endif 378112158Sdas e -= i; 379112158Sdas dval(rv) *= tens[i]; 380112158Sdas#ifdef VAX 381112158Sdas /* VAX exponent range is so narrow we must 382112158Sdas * worry about overflow here... 383112158Sdas */ 384112158Sdas vax_ovfl_check: 385112158Sdas word0(rv) -= P*Exp_msk1; 386112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 387112158Sdas if ((word0(rv) & Exp_mask) 388112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 389112158Sdas goto ovfl; 390112158Sdas word0(rv) += P*Exp_msk1; 391112158Sdas#else 392112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 393112158Sdas#endif 394112158Sdas goto ret; 395112158Sdas } 396112158Sdas } 397112158Sdas#ifndef Inaccurate_Divide 398112158Sdas else if (e >= -Ten_pmax) { 399112158Sdas#ifdef Honor_FLT_ROUNDS 400112158Sdas /* round correctly FLT_ROUNDS = 2 or 3 */ 401112158Sdas if (sign) { 402112158Sdas rv = -rv; 403112158Sdas sign = 0; 404112158Sdas } 405112158Sdas#endif 406112158Sdas /* rv = */ rounded_quotient(dval(rv), tens[-e]); 407112158Sdas goto ret; 408112158Sdas } 409112158Sdas#endif 410112158Sdas } 411112158Sdas e1 += nd - k; 412112158Sdas 413112158Sdas#ifdef IEEE_Arith 414112158Sdas#ifdef SET_INEXACT 415112158Sdas inexact = 1; 416112158Sdas if (k <= DBL_DIG) 417112158Sdas oldinexact = get_inexact(); 418112158Sdas#endif 419112158Sdas#ifdef Avoid_Underflow 420112158Sdas scale = 0; 421112158Sdas#endif 422112158Sdas#ifdef Honor_FLT_ROUNDS 423182709Sdas if (Rounding >= 2) { 424112158Sdas if (sign) 425182709Sdas Rounding = Rounding == 2 ? 0 : 2; 426112158Sdas else 427182709Sdas if (Rounding != 2) 428182709Sdas Rounding = 0; 429112158Sdas } 430112158Sdas#endif 431112158Sdas#endif /*IEEE_Arith*/ 432112158Sdas 433112158Sdas /* Get starting approximation = rv * 10**e1 */ 434112158Sdas 435112158Sdas if (e1 > 0) { 436112158Sdas if ( (i = e1 & 15) !=0) 437112158Sdas dval(rv) *= tens[i]; 438112158Sdas if (e1 &= ~15) { 439112158Sdas if (e1 > DBL_MAX_10_EXP) { 440112158Sdas ovfl: 441112158Sdas#ifndef NO_ERRNO 442112158Sdas errno = ERANGE; 443112158Sdas#endif 444112158Sdas /* Can't trust HUGE_VAL */ 445112158Sdas#ifdef IEEE_Arith 446112158Sdas#ifdef Honor_FLT_ROUNDS 447182709Sdas switch(Rounding) { 448112158Sdas case 0: /* toward 0 */ 449112158Sdas case 3: /* toward -infinity */ 450112158Sdas word0(rv) = Big0; 451112158Sdas word1(rv) = Big1; 452112158Sdas break; 453112158Sdas default: 454112158Sdas word0(rv) = Exp_mask; 455112158Sdas word1(rv) = 0; 456112158Sdas } 457112158Sdas#else /*Honor_FLT_ROUNDS*/ 458112158Sdas word0(rv) = Exp_mask; 459112158Sdas word1(rv) = 0; 460112158Sdas#endif /*Honor_FLT_ROUNDS*/ 461112158Sdas#ifdef SET_INEXACT 462112158Sdas /* set overflow bit */ 463112158Sdas dval(rv0) = 1e300; 464112158Sdas dval(rv0) *= dval(rv0); 465112158Sdas#endif 466112158Sdas#else /*IEEE_Arith*/ 467112158Sdas word0(rv) = Big0; 468112158Sdas word1(rv) = Big1; 469112158Sdas#endif /*IEEE_Arith*/ 470112158Sdas if (bd0) 471112158Sdas goto retfree; 472112158Sdas goto ret; 473112158Sdas } 474112158Sdas e1 >>= 4; 475112158Sdas for(j = 0; e1 > 1; j++, e1 >>= 1) 476112158Sdas if (e1 & 1) 477112158Sdas dval(rv) *= bigtens[j]; 478112158Sdas /* The last multiplication could overflow. */ 479112158Sdas word0(rv) -= P*Exp_msk1; 480112158Sdas dval(rv) *= bigtens[j]; 481112158Sdas if ((z = word0(rv) & Exp_mask) 482112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 483112158Sdas goto ovfl; 484112158Sdas if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 485112158Sdas /* set to largest number */ 486112158Sdas /* (Can't trust DBL_MAX) */ 487112158Sdas word0(rv) = Big0; 488112158Sdas word1(rv) = Big1; 489112158Sdas } 490112158Sdas else 491112158Sdas word0(rv) += P*Exp_msk1; 492112158Sdas } 493112158Sdas } 494112158Sdas else if (e1 < 0) { 495112158Sdas e1 = -e1; 496112158Sdas if ( (i = e1 & 15) !=0) 497112158Sdas dval(rv) /= tens[i]; 498112158Sdas if (e1 >>= 4) { 499112158Sdas if (e1 >= 1 << n_bigtens) 500112158Sdas goto undfl; 501112158Sdas#ifdef Avoid_Underflow 502112158Sdas if (e1 & Scale_Bit) 503112158Sdas scale = 2*P; 504112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 505112158Sdas if (e1 & 1) 506112158Sdas dval(rv) *= tinytens[j]; 507112158Sdas if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask) 508112158Sdas >> Exp_shift)) > 0) { 509112158Sdas /* scaled rv is denormal; zap j low bits */ 510112158Sdas if (j >= 32) { 511112158Sdas word1(rv) = 0; 512112158Sdas if (j >= 53) 513112158Sdas word0(rv) = (P+2)*Exp_msk1; 514112158Sdas else 515112158Sdas word0(rv) &= 0xffffffff << j-32; 516112158Sdas } 517112158Sdas else 518112158Sdas word1(rv) &= 0xffffffff << j; 519112158Sdas } 520112158Sdas#else 521112158Sdas for(j = 0; e1 > 1; j++, e1 >>= 1) 522112158Sdas if (e1 & 1) 523112158Sdas dval(rv) *= tinytens[j]; 524112158Sdas /* The last multiplication could underflow. */ 525112158Sdas dval(rv0) = dval(rv); 526112158Sdas dval(rv) *= tinytens[j]; 527112158Sdas if (!dval(rv)) { 528112158Sdas dval(rv) = 2.*dval(rv0); 529112158Sdas dval(rv) *= tinytens[j]; 530112158Sdas#endif 531112158Sdas if (!dval(rv)) { 532112158Sdas undfl: 533112158Sdas dval(rv) = 0.; 534112158Sdas#ifndef NO_ERRNO 535112158Sdas errno = ERANGE; 536112158Sdas#endif 537112158Sdas if (bd0) 538112158Sdas goto retfree; 539112158Sdas goto ret; 540112158Sdas } 541112158Sdas#ifndef Avoid_Underflow 542112158Sdas word0(rv) = Tiny0; 543112158Sdas word1(rv) = Tiny1; 544112158Sdas /* The refinement below will clean 545112158Sdas * this approximation up. 546112158Sdas */ 547112158Sdas } 548112158Sdas#endif 549112158Sdas } 550112158Sdas } 551112158Sdas 552112158Sdas /* Now the hard part -- adjusting rv to the correct value.*/ 553112158Sdas 554112158Sdas /* Put digits into bd: true value = bd * 10^e */ 555112158Sdas 556187808Sdas bd0 = s2b(s0, nd0, nd, y, dplen); 557112158Sdas 558112158Sdas for(;;) { 559112158Sdas bd = Balloc(bd0->k); 560112158Sdas Bcopy(bd, bd0); 561112158Sdas bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 562112158Sdas bs = i2b(1); 563112158Sdas 564112158Sdas if (e >= 0) { 565112158Sdas bb2 = bb5 = 0; 566112158Sdas bd2 = bd5 = e; 567112158Sdas } 568112158Sdas else { 569112158Sdas bb2 = bb5 = -e; 570112158Sdas bd2 = bd5 = 0; 571112158Sdas } 572112158Sdas if (bbe >= 0) 573112158Sdas bb2 += bbe; 574112158Sdas else 575112158Sdas bd2 -= bbe; 576112158Sdas bs2 = bb2; 577112158Sdas#ifdef Honor_FLT_ROUNDS 578182709Sdas if (Rounding != 1) 579112158Sdas bs2++; 580112158Sdas#endif 581112158Sdas#ifdef Avoid_Underflow 582112158Sdas j = bbe - scale; 583112158Sdas i = j + bbbits - 1; /* logb(rv) */ 584112158Sdas if (i < Emin) /* denormal */ 585112158Sdas j += P - Emin; 586112158Sdas else 587112158Sdas j = P + 1 - bbbits; 588112158Sdas#else /*Avoid_Underflow*/ 589112158Sdas#ifdef Sudden_Underflow 590112158Sdas#ifdef IBM 591112158Sdas j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 592112158Sdas#else 593112158Sdas j = P + 1 - bbbits; 594112158Sdas#endif 595112158Sdas#else /*Sudden_Underflow*/ 596112158Sdas j = bbe; 597112158Sdas i = j + bbbits - 1; /* logb(rv) */ 598112158Sdas if (i < Emin) /* denormal */ 599112158Sdas j += P - Emin; 600112158Sdas else 601112158Sdas j = P + 1 - bbbits; 602112158Sdas#endif /*Sudden_Underflow*/ 603112158Sdas#endif /*Avoid_Underflow*/ 604112158Sdas bb2 += j; 605112158Sdas bd2 += j; 606112158Sdas#ifdef Avoid_Underflow 607112158Sdas bd2 += scale; 608112158Sdas#endif 609112158Sdas i = bb2 < bd2 ? bb2 : bd2; 610112158Sdas if (i > bs2) 611112158Sdas i = bs2; 612112158Sdas if (i > 0) { 613112158Sdas bb2 -= i; 614112158Sdas bd2 -= i; 615112158Sdas bs2 -= i; 616112158Sdas } 617112158Sdas if (bb5 > 0) { 618112158Sdas bs = pow5mult(bs, bb5); 619112158Sdas bb1 = mult(bs, bb); 620112158Sdas Bfree(bb); 621112158Sdas bb = bb1; 622112158Sdas } 623112158Sdas if (bb2 > 0) 624112158Sdas bb = lshift(bb, bb2); 625112158Sdas if (bd5 > 0) 626112158Sdas bd = pow5mult(bd, bd5); 627112158Sdas if (bd2 > 0) 628112158Sdas bd = lshift(bd, bd2); 629112158Sdas if (bs2 > 0) 630112158Sdas bs = lshift(bs, bs2); 631112158Sdas delta = diff(bb, bd); 632112158Sdas dsign = delta->sign; 633112158Sdas delta->sign = 0; 634112158Sdas i = cmp(delta, bs); 635112158Sdas#ifdef Honor_FLT_ROUNDS 636182709Sdas if (Rounding != 1) { 637112158Sdas if (i < 0) { 638112158Sdas /* Error is less than an ulp */ 639112158Sdas if (!delta->x[0] && delta->wds <= 1) { 640112158Sdas /* exact */ 641112158Sdas#ifdef SET_INEXACT 642112158Sdas inexact = 0; 643112158Sdas#endif 644112158Sdas break; 645112158Sdas } 646182709Sdas if (Rounding) { 647112158Sdas if (dsign) { 648112158Sdas adj = 1.; 649112158Sdas goto apply_adj; 650112158Sdas } 651112158Sdas } 652112158Sdas else if (!dsign) { 653112158Sdas adj = -1.; 654112158Sdas if (!word1(rv) 655112158Sdas && !(word0(rv) & Frac_mask)) { 656112158Sdas y = word0(rv) & Exp_mask; 657112158Sdas#ifdef Avoid_Underflow 658112158Sdas if (!scale || y > 2*P*Exp_msk1) 659112158Sdas#else 660112158Sdas if (y) 661112158Sdas#endif 662112158Sdas { 663112158Sdas delta = lshift(delta,Log2P); 664112158Sdas if (cmp(delta, bs) <= 0) 665112158Sdas adj = -0.5; 666112158Sdas } 667112158Sdas } 668112158Sdas apply_adj: 669112158Sdas#ifdef Avoid_Underflow 670112158Sdas if (scale && (y = word0(rv) & Exp_mask) 671112158Sdas <= 2*P*Exp_msk1) 672112158Sdas word0(adj) += (2*P+1)*Exp_msk1 - y; 673112158Sdas#else 674112158Sdas#ifdef Sudden_Underflow 675112158Sdas if ((word0(rv) & Exp_mask) <= 676112158Sdas P*Exp_msk1) { 677112158Sdas word0(rv) += P*Exp_msk1; 678112158Sdas dval(rv) += adj*ulp(dval(rv)); 679112158Sdas word0(rv) -= P*Exp_msk1; 680112158Sdas } 681112158Sdas else 682112158Sdas#endif /*Sudden_Underflow*/ 683112158Sdas#endif /*Avoid_Underflow*/ 684112158Sdas dval(rv) += adj*ulp(dval(rv)); 685112158Sdas } 686112158Sdas break; 687112158Sdas } 688112158Sdas adj = ratio(delta, bs); 689112158Sdas if (adj < 1.) 690112158Sdas adj = 1.; 691112158Sdas if (adj <= 0x7ffffffe) { 692182709Sdas /* adj = Rounding ? ceil(adj) : floor(adj); */ 693112158Sdas y = adj; 694112158Sdas if (y != adj) { 695182709Sdas if (!((Rounding>>1) ^ dsign)) 696112158Sdas y++; 697112158Sdas adj = y; 698112158Sdas } 699112158Sdas } 700112158Sdas#ifdef Avoid_Underflow 701112158Sdas if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 702112158Sdas word0(adj) += (2*P+1)*Exp_msk1 - y; 703112158Sdas#else 704112158Sdas#ifdef Sudden_Underflow 705112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 706112158Sdas word0(rv) += P*Exp_msk1; 707112158Sdas adj *= ulp(dval(rv)); 708112158Sdas if (dsign) 709112158Sdas dval(rv) += adj; 710112158Sdas else 711112158Sdas dval(rv) -= adj; 712112158Sdas word0(rv) -= P*Exp_msk1; 713112158Sdas goto cont; 714112158Sdas } 715112158Sdas#endif /*Sudden_Underflow*/ 716112158Sdas#endif /*Avoid_Underflow*/ 717112158Sdas adj *= ulp(dval(rv)); 718182709Sdas if (dsign) { 719182709Sdas if (word0(rv) == Big0 && word1(rv) == Big1) 720182709Sdas goto ovfl; 721112158Sdas dval(rv) += adj; 722182709Sdas } 723112158Sdas else 724112158Sdas dval(rv) -= adj; 725112158Sdas goto cont; 726112158Sdas } 727112158Sdas#endif /*Honor_FLT_ROUNDS*/ 728112158Sdas 729112158Sdas if (i < 0) { 730112158Sdas /* Error is less than half an ulp -- check for 731112158Sdas * special case of mantissa a power of two. 732112158Sdas */ 733112158Sdas if (dsign || word1(rv) || word0(rv) & Bndry_mask 734112158Sdas#ifdef IEEE_Arith 735112158Sdas#ifdef Avoid_Underflow 736112158Sdas || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 737112158Sdas#else 738112158Sdas || (word0(rv) & Exp_mask) <= Exp_msk1 739112158Sdas#endif 740112158Sdas#endif 741112158Sdas ) { 742112158Sdas#ifdef SET_INEXACT 743112158Sdas if (!delta->x[0] && delta->wds <= 1) 744112158Sdas inexact = 0; 745112158Sdas#endif 746112158Sdas break; 747112158Sdas } 748112158Sdas if (!delta->x[0] && delta->wds <= 1) { 749112158Sdas /* exact result */ 750112158Sdas#ifdef SET_INEXACT 751112158Sdas inexact = 0; 752112158Sdas#endif 753112158Sdas break; 754112158Sdas } 755112158Sdas delta = lshift(delta,Log2P); 756112158Sdas if (cmp(delta, bs) > 0) 757112158Sdas goto drop_down; 758112158Sdas break; 759112158Sdas } 760112158Sdas if (i == 0) { 761112158Sdas /* exactly half-way between */ 762112158Sdas if (dsign) { 763112158Sdas if ((word0(rv) & Bndry_mask1) == Bndry_mask1 764112158Sdas && word1(rv) == ( 765112158Sdas#ifdef Avoid_Underflow 766112158Sdas (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 767112158Sdas ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 768112158Sdas#endif 769112158Sdas 0xffffffff)) { 770112158Sdas /*boundary case -- increment exponent*/ 771112158Sdas word0(rv) = (word0(rv) & Exp_mask) 772112158Sdas + Exp_msk1 773112158Sdas#ifdef IBM 774112158Sdas | Exp_msk1 >> 4 775112158Sdas#endif 776112158Sdas ; 777112158Sdas word1(rv) = 0; 778112158Sdas#ifdef Avoid_Underflow 779112158Sdas dsign = 0; 780112158Sdas#endif 781112158Sdas break; 782112158Sdas } 783112158Sdas } 784112158Sdas else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 785112158Sdas drop_down: 786112158Sdas /* boundary case -- decrement exponent */ 787112158Sdas#ifdef Sudden_Underflow /*{{*/ 788112158Sdas L = word0(rv) & Exp_mask; 789112158Sdas#ifdef IBM 790112158Sdas if (L < Exp_msk1) 791112158Sdas#else 792112158Sdas#ifdef Avoid_Underflow 793112158Sdas if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 794112158Sdas#else 795112158Sdas if (L <= Exp_msk1) 796112158Sdas#endif /*Avoid_Underflow*/ 797112158Sdas#endif /*IBM*/ 798112158Sdas goto undfl; 799112158Sdas L -= Exp_msk1; 800112158Sdas#else /*Sudden_Underflow}{*/ 801112158Sdas#ifdef Avoid_Underflow 802112158Sdas if (scale) { 803112158Sdas L = word0(rv) & Exp_mask; 804112158Sdas if (L <= (2*P+1)*Exp_msk1) { 805112158Sdas if (L > (P+2)*Exp_msk1) 806112158Sdas /* round even ==> */ 807112158Sdas /* accept rv */ 808112158Sdas break; 809112158Sdas /* rv = smallest denormal */ 810112158Sdas goto undfl; 811112158Sdas } 812112158Sdas } 813112158Sdas#endif /*Avoid_Underflow*/ 814112158Sdas L = (word0(rv) & Exp_mask) - Exp_msk1; 815182709Sdas#endif /*Sudden_Underflow}}*/ 816112158Sdas word0(rv) = L | Bndry_mask1; 817112158Sdas word1(rv) = 0xffffffff; 818112158Sdas#ifdef IBM 819112158Sdas goto cont; 820112158Sdas#else 821112158Sdas break; 822112158Sdas#endif 823112158Sdas } 824112158Sdas#ifndef ROUND_BIASED 825112158Sdas if (!(word1(rv) & LSB)) 826112158Sdas break; 827112158Sdas#endif 828112158Sdas if (dsign) 829112158Sdas dval(rv) += ulp(dval(rv)); 830112158Sdas#ifndef ROUND_BIASED 831112158Sdas else { 832112158Sdas dval(rv) -= ulp(dval(rv)); 833112158Sdas#ifndef Sudden_Underflow 834112158Sdas if (!dval(rv)) 835112158Sdas goto undfl; 836112158Sdas#endif 837112158Sdas } 838112158Sdas#ifdef Avoid_Underflow 839112158Sdas dsign = 1 - dsign; 840112158Sdas#endif 841112158Sdas#endif 842112158Sdas break; 843112158Sdas } 844112158Sdas if ((aadj = ratio(delta, bs)) <= 2.) { 845112158Sdas if (dsign) 846112158Sdas aadj = aadj1 = 1.; 847112158Sdas else if (word1(rv) || word0(rv) & Bndry_mask) { 848112158Sdas#ifndef Sudden_Underflow 849112158Sdas if (word1(rv) == Tiny1 && !word0(rv)) 850112158Sdas goto undfl; 851112158Sdas#endif 852112158Sdas aadj = 1.; 853112158Sdas aadj1 = -1.; 854112158Sdas } 855112158Sdas else { 856112158Sdas /* special case -- power of FLT_RADIX to be */ 857112158Sdas /* rounded down... */ 858112158Sdas 859112158Sdas if (aadj < 2./FLT_RADIX) 860112158Sdas aadj = 1./FLT_RADIX; 861112158Sdas else 862112158Sdas aadj *= 0.5; 863112158Sdas aadj1 = -aadj; 864112158Sdas } 865112158Sdas } 866112158Sdas else { 867112158Sdas aadj *= 0.5; 868112158Sdas aadj1 = dsign ? aadj : -aadj; 869112158Sdas#ifdef Check_FLT_ROUNDS 870112158Sdas switch(Rounding) { 871112158Sdas case 2: /* towards +infinity */ 872112158Sdas aadj1 -= 0.5; 873112158Sdas break; 874112158Sdas case 0: /* towards 0 */ 875112158Sdas case 3: /* towards -infinity */ 876112158Sdas aadj1 += 0.5; 877112158Sdas } 878112158Sdas#else 879112158Sdas if (Flt_Rounds == 0) 880112158Sdas aadj1 += 0.5; 881112158Sdas#endif /*Check_FLT_ROUNDS*/ 882112158Sdas } 883112158Sdas y = word0(rv) & Exp_mask; 884112158Sdas 885112158Sdas /* Check for overflow */ 886112158Sdas 887112158Sdas if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 888112158Sdas dval(rv0) = dval(rv); 889112158Sdas word0(rv) -= P*Exp_msk1; 890112158Sdas adj = aadj1 * ulp(dval(rv)); 891112158Sdas dval(rv) += adj; 892112158Sdas if ((word0(rv) & Exp_mask) >= 893112158Sdas Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 894112158Sdas if (word0(rv0) == Big0 && word1(rv0) == Big1) 895112158Sdas goto ovfl; 896112158Sdas word0(rv) = Big0; 897112158Sdas word1(rv) = Big1; 898112158Sdas goto cont; 899112158Sdas } 900112158Sdas else 901112158Sdas word0(rv) += P*Exp_msk1; 902112158Sdas } 903112158Sdas else { 904112158Sdas#ifdef Avoid_Underflow 905112158Sdas if (scale && y <= 2*P*Exp_msk1) { 906112158Sdas if (aadj <= 0x7fffffff) { 907112158Sdas if ((z = aadj) <= 0) 908112158Sdas z = 1; 909112158Sdas aadj = z; 910112158Sdas aadj1 = dsign ? aadj : -aadj; 911112158Sdas } 912112158Sdas word0(aadj1) += (2*P+1)*Exp_msk1 - y; 913112158Sdas } 914112158Sdas adj = aadj1 * ulp(dval(rv)); 915112158Sdas dval(rv) += adj; 916112158Sdas#else 917112158Sdas#ifdef Sudden_Underflow 918112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 919112158Sdas dval(rv0) = dval(rv); 920112158Sdas word0(rv) += P*Exp_msk1; 921112158Sdas adj = aadj1 * ulp(dval(rv)); 922112158Sdas dval(rv) += adj; 923112158Sdas#ifdef IBM 924112158Sdas if ((word0(rv) & Exp_mask) < P*Exp_msk1) 925112158Sdas#else 926112158Sdas if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 927112158Sdas#endif 928112158Sdas { 929112158Sdas if (word0(rv0) == Tiny0 930112158Sdas && word1(rv0) == Tiny1) 931112158Sdas goto undfl; 932112158Sdas word0(rv) = Tiny0; 933112158Sdas word1(rv) = Tiny1; 934112158Sdas goto cont; 935112158Sdas } 936112158Sdas else 937112158Sdas word0(rv) -= P*Exp_msk1; 938112158Sdas } 939112158Sdas else { 940112158Sdas adj = aadj1 * ulp(dval(rv)); 941112158Sdas dval(rv) += adj; 942112158Sdas } 943112158Sdas#else /*Sudden_Underflow*/ 944112158Sdas /* Compute adj so that the IEEE rounding rules will 945112158Sdas * correctly round rv + adj in some half-way cases. 946112158Sdas * If rv * ulp(rv) is denormalized (i.e., 947112158Sdas * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 948112158Sdas * trouble from bits lost to denormalization; 949112158Sdas * example: 1.2e-307 . 950112158Sdas */ 951112158Sdas if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 952112158Sdas aadj1 = (double)(int)(aadj + 0.5); 953112158Sdas if (!dsign) 954112158Sdas aadj1 = -aadj1; 955112158Sdas } 956112158Sdas adj = aadj1 * ulp(dval(rv)); 957112158Sdas dval(rv) += adj; 958112158Sdas#endif /*Sudden_Underflow*/ 959112158Sdas#endif /*Avoid_Underflow*/ 960112158Sdas } 961112158Sdas z = word0(rv) & Exp_mask; 962112158Sdas#ifndef SET_INEXACT 963112158Sdas#ifdef Avoid_Underflow 964112158Sdas if (!scale) 965112158Sdas#endif 966112158Sdas if (y == z) { 967112158Sdas /* Can we stop now? */ 968112158Sdas L = (Long)aadj; 969112158Sdas aadj -= L; 970112158Sdas /* The tolerances below are conservative. */ 971112158Sdas if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 972112158Sdas if (aadj < .4999999 || aadj > .5000001) 973112158Sdas break; 974112158Sdas } 975112158Sdas else if (aadj < .4999999/FLT_RADIX) 976112158Sdas break; 977112158Sdas } 978112158Sdas#endif 979112158Sdas cont: 980112158Sdas Bfree(bb); 981112158Sdas Bfree(bd); 982112158Sdas Bfree(bs); 983112158Sdas Bfree(delta); 984112158Sdas } 985112158Sdas#ifdef SET_INEXACT 986112158Sdas if (inexact) { 987112158Sdas if (!oldinexact) { 988112158Sdas word0(rv0) = Exp_1 + (70 << Exp_shift); 989112158Sdas word1(rv0) = 0; 990112158Sdas dval(rv0) += 1.; 991112158Sdas } 992112158Sdas } 993112158Sdas else if (!oldinexact) 994112158Sdas clear_inexact(); 995112158Sdas#endif 996112158Sdas#ifdef Avoid_Underflow 997112158Sdas if (scale) { 998112158Sdas word0(rv0) = Exp_1 - 2*P*Exp_msk1; 999112158Sdas word1(rv0) = 0; 1000112158Sdas dval(rv) *= dval(rv0); 1001112158Sdas#ifndef NO_ERRNO 1002112158Sdas /* try to avoid the bug of testing an 8087 register value */ 1003187808Sdas#ifdef IEEE_Arith 1004187808Sdas if (!(word0(rv) & Exp_mask)) 1005187808Sdas#else 1006112158Sdas if (word0(rv) == 0 && word1(rv) == 0) 1007187808Sdas#endif 1008112158Sdas errno = ERANGE; 1009112158Sdas#endif 1010112158Sdas } 1011112158Sdas#endif /* Avoid_Underflow */ 1012112158Sdas#ifdef SET_INEXACT 1013112158Sdas if (inexact && !(word0(rv) & Exp_mask)) { 1014112158Sdas /* set underflow bit */ 1015112158Sdas dval(rv0) = 1e-300; 1016112158Sdas dval(rv0) *= dval(rv0); 1017112158Sdas } 1018112158Sdas#endif 1019112158Sdas retfree: 1020112158Sdas Bfree(bb); 1021112158Sdas Bfree(bd); 1022112158Sdas Bfree(bs); 1023112158Sdas Bfree(bd0); 1024112158Sdas Bfree(delta); 1025112158Sdas ret: 1026112158Sdas if (se) 1027112158Sdas *se = (char *)s; 1028112158Sdas return sign ? -dval(rv) : dval(rv); 1029112158Sdas } 1030112158Sdas 1031