strtodg.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 32112158Sdas#include "gdtoaimp.h" 33112158Sdas 34112158Sdas#ifdef USE_LOCALE 35112158Sdas#include "locale.h" 36112158Sdas#endif 37112158Sdas 38112158Sdas static CONST int 39112158Sdasfivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21, 40112158Sdas 24, 26, 28, 31, 33, 35, 38, 40, 42, 45, 41112158Sdas 47, 49, 52 42112158Sdas#ifdef VAX 43112158Sdas , 54, 56 44112158Sdas#endif 45112158Sdas }; 46112158Sdas 47112158Sdas Bigint * 48112158Sdas#ifdef KR_headers 49112158Sdasincrement(b) Bigint *b; 50112158Sdas#else 51112158Sdasincrement(Bigint *b) 52112158Sdas#endif 53112158Sdas{ 54112158Sdas ULong *x, *xe; 55112158Sdas Bigint *b1; 56112158Sdas#ifdef Pack_16 57112158Sdas ULong carry = 1, y; 58112158Sdas#endif 59112158Sdas 60112158Sdas x = b->x; 61112158Sdas xe = x + b->wds; 62112158Sdas#ifdef Pack_32 63112158Sdas do { 64112158Sdas if (*x < (ULong)0xffffffffL) { 65112158Sdas ++*x; 66112158Sdas return b; 67112158Sdas } 68112158Sdas *x++ = 0; 69112158Sdas } while(x < xe); 70112158Sdas#else 71112158Sdas do { 72112158Sdas y = *x + carry; 73112158Sdas carry = y >> 16; 74112158Sdas *x++ = y & 0xffff; 75112158Sdas if (!carry) 76112158Sdas return b; 77112158Sdas } while(x < xe); 78112158Sdas if (carry) 79112158Sdas#endif 80112158Sdas { 81112158Sdas if (b->wds >= b->maxwds) { 82112158Sdas b1 = Balloc(b->k+1); 83112158Sdas Bcopy(b1,b); 84112158Sdas Bfree(b); 85112158Sdas b = b1; 86112158Sdas } 87112158Sdas b->x[b->wds++] = 1; 88112158Sdas } 89112158Sdas return b; 90112158Sdas } 91112158Sdas 92182709Sdas void 93112158Sdas#ifdef KR_headers 94112158Sdasdecrement(b) Bigint *b; 95112158Sdas#else 96112158Sdasdecrement(Bigint *b) 97112158Sdas#endif 98112158Sdas{ 99112158Sdas ULong *x, *xe; 100112158Sdas#ifdef Pack_16 101112158Sdas ULong borrow = 1, y; 102112158Sdas#endif 103112158Sdas 104112158Sdas x = b->x; 105112158Sdas xe = x + b->wds; 106112158Sdas#ifdef Pack_32 107112158Sdas do { 108112158Sdas if (*x) { 109112158Sdas --*x; 110112158Sdas break; 111112158Sdas } 112112158Sdas *x++ = 0xffffffffL; 113112158Sdas } 114112158Sdas while(x < xe); 115112158Sdas#else 116112158Sdas do { 117112158Sdas y = *x - borrow; 118112158Sdas borrow = (y & 0x10000) >> 16; 119112158Sdas *x++ = y & 0xffff; 120112158Sdas } while(borrow && x < xe); 121112158Sdas#endif 122112158Sdas } 123112158Sdas 124112158Sdas static int 125112158Sdas#ifdef KR_headers 126112158Sdasall_on(b, n) Bigint *b; int n; 127112158Sdas#else 128112158Sdasall_on(Bigint *b, int n) 129112158Sdas#endif 130112158Sdas{ 131112158Sdas ULong *x, *xe; 132112158Sdas 133112158Sdas x = b->x; 134112158Sdas xe = x + (n >> kshift); 135112158Sdas while(x < xe) 136112158Sdas if ((*x++ & ALL_ON) != ALL_ON) 137112158Sdas return 0; 138112158Sdas if (n &= kmask) 139112158Sdas return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON; 140112158Sdas return 1; 141112158Sdas } 142112158Sdas 143112158Sdas Bigint * 144112158Sdas#ifdef KR_headers 145112158Sdasset_ones(b, n) Bigint *b; int n; 146112158Sdas#else 147112158Sdasset_ones(Bigint *b, int n) 148112158Sdas#endif 149112158Sdas{ 150112158Sdas int k; 151112158Sdas ULong *x, *xe; 152112158Sdas 153112158Sdas k = (n + ((1 << kshift) - 1)) >> kshift; 154112158Sdas if (b->k < k) { 155112158Sdas Bfree(b); 156112158Sdas b = Balloc(k); 157112158Sdas } 158112158Sdas k = n >> kshift; 159112158Sdas if (n &= kmask) 160112158Sdas k++; 161112158Sdas b->wds = k; 162112158Sdas x = b->x; 163112158Sdas xe = x + k; 164112158Sdas while(x < xe) 165112158Sdas *x++ = ALL_ON; 166112158Sdas if (n) 167112158Sdas x[-1] >>= ULbits - n; 168112158Sdas return b; 169112158Sdas } 170112158Sdas 171112158Sdas static int 172112158SdasrvOK 173112158Sdas#ifdef KR_headers 174112158Sdas (d, fpi, exp, bits, exact, rd, irv) 175112158Sdas double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv; 176112158Sdas#else 177112158Sdas (double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv) 178112158Sdas#endif 179112158Sdas{ 180112158Sdas Bigint *b; 181112158Sdas ULong carry, inex, lostbits; 182112158Sdas int bdif, e, j, k, k1, nb, rv; 183112158Sdas 184112158Sdas carry = rv = 0; 185112158Sdas b = d2b(d, &e, &bdif); 186112158Sdas bdif -= nb = fpi->nbits; 187112158Sdas e += bdif; 188112158Sdas if (bdif <= 0) { 189112158Sdas if (exact) 190112158Sdas goto trunc; 191112158Sdas goto ret; 192112158Sdas } 193112158Sdas if (P == nb) { 194112158Sdas if ( 195112158Sdas#ifndef IMPRECISE_INEXACT 196112158Sdas exact && 197112158Sdas#endif 198112158Sdas fpi->rounding == 199112158Sdas#ifdef RND_PRODQUOT 200112158Sdas FPI_Round_near 201112158Sdas#else 202112158Sdas Flt_Rounds 203112158Sdas#endif 204112158Sdas ) goto trunc; 205112158Sdas goto ret; 206112158Sdas } 207112158Sdas switch(rd) { 208182709Sdas case 1: /* round down (toward -Infinity) */ 209112158Sdas goto trunc; 210182709Sdas case 2: /* round up (toward +Infinity) */ 211112158Sdas break; 212112158Sdas default: /* round near */ 213112158Sdas k = bdif - 1; 214112158Sdas if (k < 0) 215112158Sdas goto trunc; 216112158Sdas if (!k) { 217112158Sdas if (!exact) 218112158Sdas goto ret; 219112158Sdas if (b->x[0] & 2) 220112158Sdas break; 221112158Sdas goto trunc; 222112158Sdas } 223112158Sdas if (b->x[k>>kshift] & ((ULong)1 << (k & kmask))) 224112158Sdas break; 225112158Sdas goto trunc; 226112158Sdas } 227112158Sdas /* "break" cases: round up 1 bit, then truncate; bdif > 0 */ 228112158Sdas carry = 1; 229112158Sdas trunc: 230112158Sdas inex = lostbits = 0; 231112158Sdas if (bdif > 0) { 232112158Sdas if ( (lostbits = any_on(b, bdif)) !=0) 233112158Sdas inex = STRTOG_Inexlo; 234112158Sdas rshift(b, bdif); 235112158Sdas if (carry) { 236112158Sdas inex = STRTOG_Inexhi; 237112158Sdas b = increment(b); 238112158Sdas if ( (j = nb & kmask) !=0) 239165743Sdas j = ULbits - j; 240112158Sdas if (hi0bits(b->x[b->wds - 1]) != j) { 241112158Sdas if (!lostbits) 242112158Sdas lostbits = b->x[0] & 1; 243112158Sdas rshift(b, 1); 244112158Sdas e++; 245112158Sdas } 246112158Sdas } 247112158Sdas } 248112158Sdas else if (bdif < 0) 249112158Sdas b = lshift(b, -bdif); 250112158Sdas if (e < fpi->emin) { 251112158Sdas k = fpi->emin - e; 252112158Sdas e = fpi->emin; 253112158Sdas if (k > nb || fpi->sudden_underflow) { 254112158Sdas b->wds = inex = 0; 255112158Sdas *irv = STRTOG_Underflow | STRTOG_Inexlo; 256112158Sdas } 257112158Sdas else { 258112158Sdas k1 = k - 1; 259112158Sdas if (k1 > 0 && !lostbits) 260112158Sdas lostbits = any_on(b, k1); 261112158Sdas if (!lostbits && !exact) 262112158Sdas goto ret; 263112158Sdas lostbits |= 264112158Sdas carry = b->x[k1>>kshift] & (1 << (k1 & kmask)); 265112158Sdas rshift(b, k); 266112158Sdas *irv = STRTOG_Denormal; 267112158Sdas if (carry) { 268112158Sdas b = increment(b); 269112158Sdas inex = STRTOG_Inexhi | STRTOG_Underflow; 270112158Sdas } 271112158Sdas else if (lostbits) 272112158Sdas inex = STRTOG_Inexlo | STRTOG_Underflow; 273112158Sdas } 274112158Sdas } 275112158Sdas else if (e > fpi->emax) { 276112158Sdas e = fpi->emax + 1; 277112158Sdas *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 278112158Sdas#ifndef NO_ERRNO 279112158Sdas errno = ERANGE; 280112158Sdas#endif 281112158Sdas b->wds = inex = 0; 282112158Sdas } 283112158Sdas *exp = e; 284112158Sdas copybits(bits, nb, b); 285112158Sdas *irv |= inex; 286112158Sdas rv = 1; 287112158Sdas ret: 288112158Sdas Bfree(b); 289112158Sdas return rv; 290112158Sdas } 291112158Sdas 292112158Sdas static int 293112158Sdas#ifdef KR_headers 294112158Sdasmantbits(d) double d; 295112158Sdas#else 296112158Sdasmantbits(double d) 297112158Sdas#endif 298112158Sdas{ 299112158Sdas ULong L; 300112158Sdas#ifdef VAX 301112158Sdas L = word1(d) << 16 | word1(d) >> 16; 302112158Sdas if (L) 303112158Sdas#else 304112158Sdas if ( (L = word1(d)) !=0) 305112158Sdas#endif 306112158Sdas return P - lo0bits(&L); 307112158Sdas#ifdef VAX 308112158Sdas L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11; 309112158Sdas#else 310112158Sdas L = word0(d) | Exp_msk1; 311112158Sdas#endif 312112158Sdas return P - 32 - lo0bits(&L); 313112158Sdas } 314112158Sdas 315112158Sdas int 316112158Sdasstrtodg 317112158Sdas#ifdef KR_headers 318112158Sdas (s00, se, fpi, exp, bits) 319112158Sdas CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits; 320112158Sdas#else 321112158Sdas (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits) 322112158Sdas#endif 323112158Sdas{ 324112158Sdas int abe, abits, asub; 325165743Sdas int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm; 326165743Sdas int dsign, e, e1, e2, emin, esign, finished, i, inex, irv; 327112158Sdas int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign; 328112158Sdas int sudden_underflow; 329112158Sdas CONST char *s, *s0, *s1; 330112158Sdas double adj, adj0, rv, tol; 331112158Sdas Long L; 332182709Sdas ULong *b, *be, y, z; 333112158Sdas Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; 334112158Sdas 335112158Sdas irv = STRTOG_Zero; 336112158Sdas denorm = sign = nz0 = nz = 0; 337112158Sdas dval(rv) = 0.; 338112158Sdas rvb = 0; 339112158Sdas nbits = fpi->nbits; 340112158Sdas for(s = s00;;s++) switch(*s) { 341112158Sdas case '-': 342112158Sdas sign = 1; 343112158Sdas /* no break */ 344112158Sdas case '+': 345112158Sdas if (*++s) 346112158Sdas goto break2; 347112158Sdas /* no break */ 348112158Sdas case 0: 349112158Sdas sign = 0; 350112158Sdas irv = STRTOG_NoNumber; 351112158Sdas s = s00; 352112158Sdas goto ret; 353112158Sdas case '\t': 354112158Sdas case '\n': 355112158Sdas case '\v': 356112158Sdas case '\f': 357112158Sdas case '\r': 358112158Sdas case ' ': 359112158Sdas continue; 360112158Sdas default: 361112158Sdas goto break2; 362112158Sdas } 363112158Sdas break2: 364112158Sdas if (*s == '0') { 365112158Sdas#ifndef NO_HEX_FP 366112158Sdas switch(s[1]) { 367112158Sdas case 'x': 368112158Sdas case 'X': 369112158Sdas irv = gethex(&s, fpi, exp, &rvb, sign); 370112158Sdas if (irv == STRTOG_NoNumber) { 371112158Sdas s = s00; 372112158Sdas sign = 0; 373112158Sdas } 374112158Sdas goto ret; 375112158Sdas } 376112158Sdas#endif 377112158Sdas nz0 = 1; 378112158Sdas while(*++s == '0') ; 379112158Sdas if (!*s) 380112158Sdas goto ret; 381112158Sdas } 382112158Sdas sudden_underflow = fpi->sudden_underflow; 383112158Sdas s0 = s; 384112158Sdas y = z = 0; 385165743Sdas for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 386112158Sdas if (nd < 9) 387112158Sdas y = 10*y + c - '0'; 388112158Sdas else if (nd < 16) 389112158Sdas z = 10*z + c - '0'; 390112158Sdas nd0 = nd; 391112158Sdas#ifdef USE_LOCALE 392112415Sdas if (c == *localeconv()->decimal_point) 393112415Sdas#else 394112415Sdas if (c == '.') 395112158Sdas#endif 396112415Sdas { 397165743Sdas decpt = 1; 398112158Sdas c = *++s; 399112158Sdas if (!nd) { 400112158Sdas for(; c == '0'; c = *++s) 401112158Sdas nz++; 402112158Sdas if (c > '0' && c <= '9') { 403112158Sdas s0 = s; 404112158Sdas nf += nz; 405112158Sdas nz = 0; 406112158Sdas goto have_dig; 407112158Sdas } 408112158Sdas goto dig_done; 409112158Sdas } 410112158Sdas for(; c >= '0' && c <= '9'; c = *++s) { 411112158Sdas have_dig: 412112158Sdas nz++; 413112158Sdas if (c -= '0') { 414112158Sdas nf += nz; 415112158Sdas for(i = 1; i < nz; i++) 416112158Sdas if (nd++ < 9) 417112158Sdas y *= 10; 418112158Sdas else if (nd <= DBL_DIG + 1) 419112158Sdas z *= 10; 420112158Sdas if (nd++ < 9) 421112158Sdas y = 10*y + c; 422112158Sdas else if (nd <= DBL_DIG + 1) 423112158Sdas z = 10*z + c; 424112158Sdas nz = 0; 425112158Sdas } 426112158Sdas } 427112158Sdas } 428112158Sdas dig_done: 429112158Sdas e = 0; 430112158Sdas if (c == 'e' || c == 'E') { 431112158Sdas if (!nd && !nz && !nz0) { 432112158Sdas irv = STRTOG_NoNumber; 433112158Sdas s = s00; 434112158Sdas goto ret; 435112158Sdas } 436112158Sdas s00 = s; 437112158Sdas esign = 0; 438112158Sdas switch(c = *++s) { 439112158Sdas case '-': 440112158Sdas esign = 1; 441112158Sdas case '+': 442112158Sdas c = *++s; 443112158Sdas } 444112158Sdas if (c >= '0' && c <= '9') { 445112158Sdas while(c == '0') 446112158Sdas c = *++s; 447112158Sdas if (c > '0' && c <= '9') { 448112158Sdas L = c - '0'; 449112158Sdas s1 = s; 450112158Sdas while((c = *++s) >= '0' && c <= '9') 451112158Sdas L = 10*L + c - '0'; 452112158Sdas if (s - s1 > 8 || L > 19999) 453112158Sdas /* Avoid confusion from exponents 454112158Sdas * so large that e might overflow. 455112158Sdas */ 456112158Sdas e = 19999; /* safe for 16 bit ints */ 457112158Sdas else 458112158Sdas e = (int)L; 459112158Sdas if (esign) 460112158Sdas e = -e; 461112158Sdas } 462112158Sdas else 463112158Sdas e = 0; 464112158Sdas } 465112158Sdas else 466112158Sdas s = s00; 467112158Sdas } 468112158Sdas if (!nd) { 469112158Sdas if (!nz && !nz0) { 470112158Sdas#ifdef INFNAN_CHECK 471112158Sdas /* Check for Nan and Infinity */ 472165743Sdas if (!decpt) 473165743Sdas switch(c) { 474112158Sdas case 'i': 475112158Sdas case 'I': 476112158Sdas if (match(&s,"nf")) { 477112158Sdas --s; 478112158Sdas if (!match(&s,"inity")) 479112158Sdas ++s; 480112158Sdas irv = STRTOG_Infinite; 481112158Sdas goto infnanexp; 482112158Sdas } 483112158Sdas break; 484112158Sdas case 'n': 485112158Sdas case 'N': 486112158Sdas if (match(&s, "an")) { 487112158Sdas irv = STRTOG_NaN; 488112158Sdas *exp = fpi->emax + 1; 489112158Sdas#ifndef No_Hex_NaN 490112158Sdas if (*s == '(') /*)*/ 491112158Sdas irv = hexnan(&s, fpi, bits); 492112158Sdas#endif 493112158Sdas goto infnanexp; 494112158Sdas } 495112158Sdas } 496112158Sdas#endif /* INFNAN_CHECK */ 497112158Sdas irv = STRTOG_NoNumber; 498112158Sdas s = s00; 499112158Sdas } 500112158Sdas goto ret; 501112158Sdas } 502112158Sdas 503112158Sdas irv = STRTOG_Normal; 504112158Sdas e1 = e -= nf; 505112158Sdas rd = 0; 506112158Sdas switch(fpi->rounding & 3) { 507112158Sdas case FPI_Round_up: 508112158Sdas rd = 2 - sign; 509112158Sdas break; 510112158Sdas case FPI_Round_zero: 511112158Sdas rd = 1; 512112158Sdas break; 513112158Sdas case FPI_Round_down: 514112158Sdas rd = 1 + sign; 515112158Sdas } 516112158Sdas 517112158Sdas /* Now we have nd0 digits, starting at s0, followed by a 518112158Sdas * decimal point, followed by nd-nd0 digits. The number we're 519112158Sdas * after is the integer represented by those digits times 520112158Sdas * 10**e */ 521112158Sdas 522112158Sdas if (!nd0) 523112158Sdas nd0 = nd; 524112158Sdas k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 525112158Sdas dval(rv) = y; 526112158Sdas if (k > 9) 527112158Sdas dval(rv) = tens[k - 9] * dval(rv) + z; 528112158Sdas bd0 = 0; 529112158Sdas if (nbits <= P && nd <= DBL_DIG) { 530112158Sdas if (!e) { 531112158Sdas if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv)) 532112158Sdas goto ret; 533112158Sdas } 534112158Sdas else if (e > 0) { 535112158Sdas if (e <= Ten_pmax) { 536112158Sdas#ifdef VAX 537112158Sdas goto vax_ovfl_check; 538112158Sdas#else 539112158Sdas i = fivesbits[e] + mantbits(dval(rv)) <= P; 540112158Sdas /* rv = */ rounded_product(dval(rv), tens[e]); 541112158Sdas if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv)) 542112158Sdas goto ret; 543112158Sdas e1 -= e; 544112158Sdas goto rv_notOK; 545112158Sdas#endif 546112158Sdas } 547112158Sdas i = DBL_DIG - nd; 548112158Sdas if (e <= Ten_pmax + i) { 549112158Sdas /* A fancier test would sometimes let us do 550112158Sdas * this for larger i values. 551112158Sdas */ 552112158Sdas e2 = e - i; 553112158Sdas e1 -= i; 554112158Sdas dval(rv) *= tens[i]; 555112158Sdas#ifdef VAX 556112158Sdas /* VAX exponent range is so narrow we must 557112158Sdas * worry about overflow here... 558112158Sdas */ 559112158Sdas vax_ovfl_check: 560112158Sdas dval(adj) = dval(rv); 561112158Sdas word0(adj) -= P*Exp_msk1; 562112158Sdas /* adj = */ rounded_product(dval(adj), tens[e2]); 563112158Sdas if ((word0(adj) & Exp_mask) 564112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 565112158Sdas goto rv_notOK; 566112158Sdas word0(adj) += P*Exp_msk1; 567112158Sdas dval(rv) = dval(adj); 568112158Sdas#else 569112158Sdas /* rv = */ rounded_product(dval(rv), tens[e2]); 570112158Sdas#endif 571112158Sdas if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv)) 572112158Sdas goto ret; 573112158Sdas e1 -= e2; 574112158Sdas } 575112158Sdas } 576112158Sdas#ifndef Inaccurate_Divide 577112158Sdas else if (e >= -Ten_pmax) { 578112158Sdas /* rv = */ rounded_quotient(dval(rv), tens[-e]); 579112158Sdas if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv)) 580112158Sdas goto ret; 581112158Sdas e1 -= e; 582112158Sdas } 583112158Sdas#endif 584112158Sdas } 585112158Sdas rv_notOK: 586112158Sdas e1 += nd - k; 587112158Sdas 588112158Sdas /* Get starting approximation = rv * 10**e1 */ 589112158Sdas 590112158Sdas e2 = 0; 591112158Sdas if (e1 > 0) { 592112158Sdas if ( (i = e1 & 15) !=0) 593112158Sdas dval(rv) *= tens[i]; 594112158Sdas if (e1 &= ~15) { 595112158Sdas e1 >>= 4; 596112158Sdas while(e1 >= (1 << n_bigtens-1)) { 597112158Sdas e2 += ((word0(rv) & Exp_mask) 598112158Sdas >> Exp_shift1) - Bias; 599112158Sdas word0(rv) &= ~Exp_mask; 600112158Sdas word0(rv) |= Bias << Exp_shift1; 601112158Sdas dval(rv) *= bigtens[n_bigtens-1]; 602112158Sdas e1 -= 1 << n_bigtens-1; 603112158Sdas } 604112158Sdas e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; 605112158Sdas word0(rv) &= ~Exp_mask; 606112158Sdas word0(rv) |= Bias << Exp_shift1; 607112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 608112158Sdas if (e1 & 1) 609112158Sdas dval(rv) *= bigtens[j]; 610112158Sdas } 611112158Sdas } 612112158Sdas else if (e1 < 0) { 613112158Sdas e1 = -e1; 614112158Sdas if ( (i = e1 & 15) !=0) 615112158Sdas dval(rv) /= tens[i]; 616112158Sdas if (e1 &= ~15) { 617112158Sdas e1 >>= 4; 618112158Sdas while(e1 >= (1 << n_bigtens-1)) { 619112158Sdas e2 += ((word0(rv) & Exp_mask) 620112158Sdas >> Exp_shift1) - Bias; 621112158Sdas word0(rv) &= ~Exp_mask; 622112158Sdas word0(rv) |= Bias << Exp_shift1; 623112158Sdas dval(rv) *= tinytens[n_bigtens-1]; 624112158Sdas e1 -= 1 << n_bigtens-1; 625112158Sdas } 626112158Sdas e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; 627112158Sdas word0(rv) &= ~Exp_mask; 628112158Sdas word0(rv) |= Bias << Exp_shift1; 629112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 630112158Sdas if (e1 & 1) 631112158Sdas dval(rv) *= tinytens[j]; 632112158Sdas } 633112158Sdas } 634165743Sdas#ifdef IBM 635165743Sdas /* e2 is a correction to the (base 2) exponent of the return 636165743Sdas * value, reflecting adjustments above to avoid overflow in the 637165743Sdas * native arithmetic. For native IBM (base 16) arithmetic, we 638165743Sdas * must multiply e2 by 4 to change from base 16 to 2. 639165743Sdas */ 640165743Sdas e2 <<= 2; 641165743Sdas#endif 642112158Sdas rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */ 643112158Sdas rve += e2; 644112158Sdas if ((j = rvbits - nbits) > 0) { 645112158Sdas rshift(rvb, j); 646112158Sdas rvbits = nbits; 647112158Sdas rve += j; 648112158Sdas } 649112158Sdas bb0 = 0; /* trailing zero bits in rvb */ 650112158Sdas e2 = rve + rvbits - nbits; 651165743Sdas if (e2 > fpi->emax + 1) 652165743Sdas goto huge; 653112158Sdas rve1 = rve + rvbits - nbits; 654112158Sdas if (e2 < (emin = fpi->emin)) { 655112158Sdas denorm = 1; 656112158Sdas j = rve - emin; 657112158Sdas if (j > 0) { 658112158Sdas rvb = lshift(rvb, j); 659112158Sdas rvbits += j; 660112158Sdas } 661112158Sdas else if (j < 0) { 662112158Sdas rvbits += j; 663112158Sdas if (rvbits <= 0) { 664112158Sdas if (rvbits < -1) { 665112158Sdas ufl: 666112158Sdas rvb->wds = 0; 667112158Sdas rvb->x[0] = 0; 668112158Sdas *exp = emin; 669112158Sdas irv = STRTOG_Underflow | STRTOG_Inexlo; 670112158Sdas goto ret; 671112158Sdas } 672112158Sdas rvb->x[0] = rvb->wds = rvbits = 1; 673112158Sdas } 674112158Sdas else 675112158Sdas rshift(rvb, -j); 676112158Sdas } 677112158Sdas rve = rve1 = emin; 678112158Sdas if (sudden_underflow && e2 + 1 < emin) 679112158Sdas goto ufl; 680112158Sdas } 681112158Sdas 682112158Sdas /* Now the hard part -- adjusting rv to the correct value.*/ 683112158Sdas 684112158Sdas /* Put digits into bd: true value = bd * 10^e */ 685112158Sdas 686112158Sdas bd0 = s2b(s0, nd0, nd, y); 687112158Sdas 688112158Sdas for(;;) { 689112158Sdas bd = Balloc(bd0->k); 690112158Sdas Bcopy(bd, bd0); 691112158Sdas bb = Balloc(rvb->k); 692112158Sdas Bcopy(bb, rvb); 693112158Sdas bbbits = rvbits - bb0; 694112158Sdas bbe = rve + bb0; 695112158Sdas bs = i2b(1); 696112158Sdas 697112158Sdas if (e >= 0) { 698112158Sdas bb2 = bb5 = 0; 699112158Sdas bd2 = bd5 = e; 700112158Sdas } 701112158Sdas else { 702112158Sdas bb2 = bb5 = -e; 703112158Sdas bd2 = bd5 = 0; 704112158Sdas } 705112158Sdas if (bbe >= 0) 706112158Sdas bb2 += bbe; 707112158Sdas else 708112158Sdas bd2 -= bbe; 709112158Sdas bs2 = bb2; 710112158Sdas j = nbits + 1 - bbbits; 711112158Sdas i = bbe + bbbits - nbits; 712112158Sdas if (i < emin) /* denormal */ 713112158Sdas j += i - emin; 714112158Sdas bb2 += j; 715112158Sdas bd2 += j; 716112158Sdas i = bb2 < bd2 ? bb2 : bd2; 717112158Sdas if (i > bs2) 718112158Sdas i = bs2; 719112158Sdas if (i > 0) { 720112158Sdas bb2 -= i; 721112158Sdas bd2 -= i; 722112158Sdas bs2 -= i; 723112158Sdas } 724112158Sdas if (bb5 > 0) { 725112158Sdas bs = pow5mult(bs, bb5); 726112158Sdas bb1 = mult(bs, bb); 727112158Sdas Bfree(bb); 728112158Sdas bb = bb1; 729112158Sdas } 730112158Sdas bb2 -= bb0; 731112158Sdas if (bb2 > 0) 732112158Sdas bb = lshift(bb, bb2); 733112158Sdas else if (bb2 < 0) 734112158Sdas rshift(bb, -bb2); 735112158Sdas if (bd5 > 0) 736112158Sdas bd = pow5mult(bd, bd5); 737112158Sdas if (bd2 > 0) 738112158Sdas bd = lshift(bd, bd2); 739112158Sdas if (bs2 > 0) 740112158Sdas bs = lshift(bs, bs2); 741112158Sdas asub = 1; 742112158Sdas inex = STRTOG_Inexhi; 743112158Sdas delta = diff(bb, bd); 744112158Sdas if (delta->wds <= 1 && !delta->x[0]) 745112158Sdas break; 746112158Sdas dsign = delta->sign; 747112158Sdas delta->sign = finished = 0; 748112158Sdas L = 0; 749112158Sdas i = cmp(delta, bs); 750112158Sdas if (rd && i <= 0) { 751112158Sdas irv = STRTOG_Normal; 752112158Sdas if ( (finished = dsign ^ (rd&1)) !=0) { 753112158Sdas if (dsign != 0) { 754112158Sdas irv |= STRTOG_Inexhi; 755112158Sdas goto adj1; 756112158Sdas } 757112158Sdas irv |= STRTOG_Inexlo; 758112158Sdas if (rve1 == emin) 759112158Sdas goto adj1; 760112158Sdas for(i = 0, j = nbits; j >= ULbits; 761112158Sdas i++, j -= ULbits) { 762112158Sdas if (rvb->x[i] & ALL_ON) 763112158Sdas goto adj1; 764112158Sdas } 765112158Sdas if (j > 1 && lo0bits(rvb->x + i) < j - 1) 766112158Sdas goto adj1; 767112158Sdas rve = rve1 - 1; 768112158Sdas rvb = set_ones(rvb, rvbits = nbits); 769112158Sdas break; 770112158Sdas } 771112158Sdas irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; 772112158Sdas break; 773112158Sdas } 774112158Sdas if (i < 0) { 775112158Sdas /* Error is less than half an ulp -- check for 776112158Sdas * special case of mantissa a power of two. 777112158Sdas */ 778112158Sdas irv = dsign 779112158Sdas ? STRTOG_Normal | STRTOG_Inexlo 780112158Sdas : STRTOG_Normal | STRTOG_Inexhi; 781112158Sdas if (dsign || bbbits > 1 || denorm || rve1 == emin) 782112158Sdas break; 783112158Sdas delta = lshift(delta,1); 784112158Sdas if (cmp(delta, bs) > 0) { 785112158Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 786112158Sdas goto drop_down; 787112158Sdas } 788112158Sdas break; 789112158Sdas } 790112158Sdas if (i == 0) { 791112158Sdas /* exactly half-way between */ 792112158Sdas if (dsign) { 793112158Sdas if (denorm && all_on(rvb, rvbits)) { 794112158Sdas /*boundary case -- increment exponent*/ 795112158Sdas rvb->wds = 1; 796112158Sdas rvb->x[0] = 1; 797112158Sdas rve = emin + nbits - (rvbits = 1); 798112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 799112158Sdas denorm = 0; 800112158Sdas break; 801112158Sdas } 802112158Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 803112158Sdas } 804112158Sdas else if (bbbits == 1) { 805112158Sdas irv = STRTOG_Normal; 806112158Sdas drop_down: 807112158Sdas /* boundary case -- decrement exponent */ 808112158Sdas if (rve1 == emin) { 809112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 810112158Sdas if (rvb->wds == 1 && rvb->x[0] == 1) 811112158Sdas sudden_underflow = 1; 812112158Sdas break; 813112158Sdas } 814112158Sdas rve -= nbits; 815112158Sdas rvb = set_ones(rvb, rvbits = nbits); 816112158Sdas break; 817112158Sdas } 818112158Sdas else 819112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 820112158Sdas if (bbbits < nbits && !denorm || !(rvb->x[0] & 1)) 821112158Sdas break; 822112158Sdas if (dsign) { 823112158Sdas rvb = increment(rvb); 824182709Sdas j = kmask & (ULbits - (rvbits & kmask)); 825182709Sdas if (hi0bits(rvb->x[rvb->wds - 1]) != j) 826112158Sdas rvbits++; 827112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 828112158Sdas } 829112158Sdas else { 830112158Sdas if (bbbits == 1) 831112158Sdas goto undfl; 832112158Sdas decrement(rvb); 833112158Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 834112158Sdas } 835112158Sdas break; 836112158Sdas } 837112158Sdas if ((dval(adj) = ratio(delta, bs)) <= 2.) { 838112158Sdas adj1: 839112158Sdas inex = STRTOG_Inexlo; 840112158Sdas if (dsign) { 841112158Sdas asub = 0; 842112158Sdas inex = STRTOG_Inexhi; 843112158Sdas } 844112158Sdas else if (denorm && bbbits <= 1) { 845112158Sdas undfl: 846112158Sdas rvb->wds = 0; 847112158Sdas rve = emin; 848112158Sdas irv = STRTOG_Underflow | STRTOG_Inexlo; 849112158Sdas break; 850112158Sdas } 851112158Sdas adj0 = dval(adj) = 1.; 852112158Sdas } 853112158Sdas else { 854112158Sdas adj0 = dval(adj) *= 0.5; 855112158Sdas if (dsign) { 856112158Sdas asub = 0; 857112158Sdas inex = STRTOG_Inexlo; 858112158Sdas } 859112158Sdas if (dval(adj) < 2147483647.) { 860112158Sdas L = adj0; 861112158Sdas adj0 -= L; 862112158Sdas switch(rd) { 863112158Sdas case 0: 864112158Sdas if (adj0 >= .5) 865112158Sdas goto inc_L; 866112158Sdas break; 867112158Sdas case 1: 868112158Sdas if (asub && adj0 > 0.) 869112158Sdas goto inc_L; 870112158Sdas break; 871112158Sdas case 2: 872112158Sdas if (!asub && adj0 > 0.) { 873112158Sdas inc_L: 874112158Sdas L++; 875112158Sdas inex = STRTOG_Inexact - inex; 876112158Sdas } 877112158Sdas } 878112158Sdas dval(adj) = L; 879112158Sdas } 880112158Sdas } 881112158Sdas y = rve + rvbits; 882112158Sdas 883112158Sdas /* adj *= ulp(dval(rv)); */ 884112158Sdas /* if (asub) rv -= adj; else rv += adj; */ 885112158Sdas 886112158Sdas if (!denorm && rvbits < nbits) { 887112158Sdas rvb = lshift(rvb, j = nbits - rvbits); 888112158Sdas rve -= j; 889112158Sdas rvbits = nbits; 890112158Sdas } 891112158Sdas ab = d2b(dval(adj), &abe, &abits); 892112158Sdas if (abe < 0) 893112158Sdas rshift(ab, -abe); 894112158Sdas else if (abe > 0) 895112158Sdas ab = lshift(ab, abe); 896112158Sdas rvb0 = rvb; 897112158Sdas if (asub) { 898112158Sdas /* rv -= adj; */ 899112158Sdas j = hi0bits(rvb->x[rvb->wds-1]); 900112158Sdas rvb = diff(rvb, ab); 901112158Sdas k = rvb0->wds - 1; 902112158Sdas if (denorm) 903112158Sdas /* do nothing */; 904112158Sdas else if (rvb->wds <= k 905112158Sdas || hi0bits( rvb->x[k]) > 906112158Sdas hi0bits(rvb0->x[k])) { 907112158Sdas /* unlikely; can only have lost 1 high bit */ 908112158Sdas if (rve1 == emin) { 909112158Sdas --rvbits; 910112158Sdas denorm = 1; 911112158Sdas } 912112158Sdas else { 913112158Sdas rvb = lshift(rvb, 1); 914112158Sdas --rve; 915112158Sdas --rve1; 916112158Sdas L = finished = 0; 917112158Sdas } 918112158Sdas } 919112158Sdas } 920112158Sdas else { 921112158Sdas rvb = sum(rvb, ab); 922112158Sdas k = rvb->wds - 1; 923112158Sdas if (k >= rvb0->wds 924112158Sdas || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) { 925112158Sdas if (denorm) { 926112158Sdas if (++rvbits == nbits) 927112158Sdas denorm = 0; 928112158Sdas } 929112158Sdas else { 930112158Sdas rshift(rvb, 1); 931112158Sdas rve++; 932112158Sdas rve1++; 933112158Sdas L = 0; 934112158Sdas } 935112158Sdas } 936112158Sdas } 937112158Sdas Bfree(ab); 938112158Sdas Bfree(rvb0); 939112158Sdas if (finished) 940112158Sdas break; 941112158Sdas 942112158Sdas z = rve + rvbits; 943112158Sdas if (y == z && L) { 944112158Sdas /* Can we stop now? */ 945112158Sdas tol = dval(adj) * 5e-16; /* > max rel error */ 946112158Sdas dval(adj) = adj0 - .5; 947112158Sdas if (dval(adj) < -tol) { 948112158Sdas if (adj0 > tol) { 949112158Sdas irv |= inex; 950112158Sdas break; 951112158Sdas } 952112158Sdas } 953112158Sdas else if (dval(adj) > tol && adj0 < 1. - tol) { 954112158Sdas irv |= inex; 955112158Sdas break; 956112158Sdas } 957112158Sdas } 958112158Sdas bb0 = denorm ? 0 : trailz(rvb); 959112158Sdas Bfree(bb); 960112158Sdas Bfree(bd); 961112158Sdas Bfree(bs); 962112158Sdas Bfree(delta); 963112158Sdas } 964165743Sdas if (!denorm && (j = nbits - rvbits)) { 965165743Sdas if (j > 0) 966165743Sdas rvb = lshift(rvb, j); 967165743Sdas else 968165743Sdas rshift(rvb, -j); 969112158Sdas rve -= j; 970112158Sdas } 971112158Sdas *exp = rve; 972112158Sdas Bfree(bb); 973112158Sdas Bfree(bd); 974112158Sdas Bfree(bs); 975112158Sdas Bfree(bd0); 976112158Sdas Bfree(delta); 977165743Sdas if (rve > fpi->emax) { 978182709Sdas switch(fpi->rounding & 3) { 979182709Sdas case FPI_Round_near: 980182709Sdas goto huge; 981182709Sdas case FPI_Round_up: 982182709Sdas if (!sign) 983182709Sdas goto huge; 984182709Sdas break; 985182709Sdas case FPI_Round_down: 986182709Sdas if (sign) 987182709Sdas goto huge; 988182709Sdas } 989182709Sdas /* Round to largest representable magnitude */ 990182709Sdas Bfree(rvb); 991182709Sdas rvb = 0; 992182709Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 993182709Sdas *exp = fpi->emax; 994182709Sdas b = bits; 995182709Sdas be = b + (fpi->nbits >> 5) + 1; 996182709Sdas while(b < be) 997182709Sdas *b++ = -1; 998182709Sdas if ((j = fpi->nbits & 0x1f)) 999182709Sdas *--be >>= (32 - j); 1000182709Sdas goto ret; 1001165743Sdas huge: 1002165743Sdas rvb->wds = 0; 1003165743Sdas irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 1004165743Sdas#ifndef NO_ERRNO 1005165743Sdas errno = ERANGE; 1006165743Sdas#endif 1007165743Sdas infnanexp: 1008165743Sdas *exp = fpi->emax + 1; 1009165743Sdas } 1010112158Sdas ret: 1011112158Sdas if (denorm) { 1012112158Sdas if (sudden_underflow) { 1013112158Sdas rvb->wds = 0; 1014112158Sdas irv = STRTOG_Underflow | STRTOG_Inexlo; 1015182709Sdas#ifndef NO_ERRNO 1016182709Sdas errno = ERANGE; 1017182709Sdas#endif 1018112158Sdas } 1019112158Sdas else { 1020112158Sdas irv = (irv & ~STRTOG_Retmask) | 1021112158Sdas (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); 1022182709Sdas if (irv & STRTOG_Inexact) { 1023112158Sdas irv |= STRTOG_Underflow; 1024182709Sdas#ifndef NO_ERRNO 1025182709Sdas errno = ERANGE; 1026182709Sdas#endif 1027182709Sdas } 1028112158Sdas } 1029112158Sdas } 1030112158Sdas if (se) 1031112158Sdas *se = (char *)s; 1032112158Sdas if (sign) 1033112158Sdas irv |= STRTOG_Neg; 1034112158Sdas if (rvb) { 1035112158Sdas copybits(bits, nbits, rvb); 1036112158Sdas Bfree(rvb); 1037112158Sdas } 1038112158Sdas return irv; 1039112158Sdas } 1040