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