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