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) 175219557Sdas U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv; 176112158Sdas#else 177219557Sdas (U *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; 185219557Sdas b = d2b(dval(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 294219557Sdasmantbits(d) U *d; 295112158Sdas#else 296219557Sdasmantbits(U *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 316227753Stheravenstrtodg_l 317112158Sdas#ifdef KR_headers 318227753Stheraven (s00, se, fpi, exp, bits, loc) 319227753Stheraven CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits; locale_t loc; 320112158Sdas#else 321227753Stheraven (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits, locale_t loc) 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; 330219557Sdas double adj0, tol; 331112158Sdas Long L; 332219557Sdas U adj, rv; 333182709Sdas ULong *b, *be, y, z; 334112158Sdas Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; 335187808Sdas#ifdef USE_LOCALE /*{{*/ 336187808Sdas#ifdef NO_LOCALE_CACHE 337227753Stheraven char *decimalpoint = localeconv_l(loc)->decimal_point; 338187808Sdas int dplen = strlen(decimalpoint); 339187808Sdas#else 340187808Sdas char *decimalpoint; 341187808Sdas static char *decimalpoint_cache; 342187808Sdas static int dplen; 343187808Sdas if (!(s0 = decimalpoint_cache)) { 344227753Stheraven s0 = localeconv_l(loc)->decimal_point; 345219557Sdas if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { 346187808Sdas strcpy(decimalpoint_cache, s0); 347187808Sdas s0 = decimalpoint_cache; 348187808Sdas } 349187808Sdas dplen = strlen(s0); 350187808Sdas } 351187808Sdas decimalpoint = (char*)s0; 352187808Sdas#endif /*NO_LOCALE_CACHE*/ 353187808Sdas#else /*USE_LOCALE}{*/ 354187808Sdas#define dplen 1 355187808Sdas#endif /*USE_LOCALE}}*/ 356112158Sdas 357112158Sdas irv = STRTOG_Zero; 358112158Sdas denorm = sign = nz0 = nz = 0; 359219557Sdas dval(&rv) = 0.; 360112158Sdas rvb = 0; 361112158Sdas nbits = fpi->nbits; 362112158Sdas for(s = s00;;s++) switch(*s) { 363112158Sdas case '-': 364112158Sdas sign = 1; 365112158Sdas /* no break */ 366112158Sdas case '+': 367112158Sdas if (*++s) 368112158Sdas goto break2; 369112158Sdas /* no break */ 370112158Sdas case 0: 371112158Sdas sign = 0; 372112158Sdas irv = STRTOG_NoNumber; 373112158Sdas s = s00; 374112158Sdas goto ret; 375112158Sdas case '\t': 376112158Sdas case '\n': 377112158Sdas case '\v': 378112158Sdas case '\f': 379112158Sdas case '\r': 380112158Sdas case ' ': 381112158Sdas continue; 382112158Sdas default: 383112158Sdas goto break2; 384112158Sdas } 385112158Sdas break2: 386112158Sdas if (*s == '0') { 387112158Sdas#ifndef NO_HEX_FP 388112158Sdas switch(s[1]) { 389112158Sdas case 'x': 390112158Sdas case 'X': 391112158Sdas irv = gethex(&s, fpi, exp, &rvb, sign); 392112158Sdas if (irv == STRTOG_NoNumber) { 393112158Sdas s = s00; 394112158Sdas sign = 0; 395112158Sdas } 396112158Sdas goto ret; 397112158Sdas } 398112158Sdas#endif 399112158Sdas nz0 = 1; 400112158Sdas while(*++s == '0') ; 401112158Sdas if (!*s) 402112158Sdas goto ret; 403112158Sdas } 404112158Sdas sudden_underflow = fpi->sudden_underflow; 405112158Sdas s0 = s; 406112158Sdas y = z = 0; 407165743Sdas for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 408112158Sdas if (nd < 9) 409112158Sdas y = 10*y + c - '0'; 410112158Sdas else if (nd < 16) 411112158Sdas z = 10*z + c - '0'; 412112158Sdas nd0 = nd; 413112158Sdas#ifdef USE_LOCALE 414187808Sdas if (c == *decimalpoint) { 415187808Sdas for(i = 1; decimalpoint[i]; ++i) 416187808Sdas if (s[i] != decimalpoint[i]) 417187808Sdas goto dig_done; 418187808Sdas s += i; 419187808Sdas c = *s; 420112415Sdas#else 421187808Sdas if (c == '.') { 422187808Sdas c = *++s; 423112158Sdas#endif 424165743Sdas decpt = 1; 425112158Sdas if (!nd) { 426112158Sdas for(; c == '0'; c = *++s) 427112158Sdas nz++; 428112158Sdas if (c > '0' && c <= '9') { 429112158Sdas s0 = s; 430112158Sdas nf += nz; 431112158Sdas nz = 0; 432112158Sdas goto have_dig; 433112158Sdas } 434112158Sdas goto dig_done; 435112158Sdas } 436112158Sdas for(; c >= '0' && c <= '9'; c = *++s) { 437112158Sdas have_dig: 438112158Sdas nz++; 439112158Sdas if (c -= '0') { 440112158Sdas nf += nz; 441112158Sdas for(i = 1; i < nz; i++) 442112158Sdas if (nd++ < 9) 443112158Sdas y *= 10; 444112158Sdas else if (nd <= DBL_DIG + 1) 445112158Sdas z *= 10; 446112158Sdas if (nd++ < 9) 447112158Sdas y = 10*y + c; 448112158Sdas else if (nd <= DBL_DIG + 1) 449112158Sdas z = 10*z + c; 450112158Sdas nz = 0; 451112158Sdas } 452112158Sdas } 453187808Sdas }/*}*/ 454112158Sdas dig_done: 455112158Sdas e = 0; 456112158Sdas if (c == 'e' || c == 'E') { 457112158Sdas if (!nd && !nz && !nz0) { 458112158Sdas irv = STRTOG_NoNumber; 459112158Sdas s = s00; 460112158Sdas goto ret; 461112158Sdas } 462112158Sdas s00 = s; 463112158Sdas esign = 0; 464112158Sdas switch(c = *++s) { 465112158Sdas case '-': 466112158Sdas esign = 1; 467112158Sdas case '+': 468112158Sdas c = *++s; 469112158Sdas } 470112158Sdas if (c >= '0' && c <= '9') { 471112158Sdas while(c == '0') 472112158Sdas c = *++s; 473112158Sdas if (c > '0' && c <= '9') { 474112158Sdas L = c - '0'; 475112158Sdas s1 = s; 476112158Sdas while((c = *++s) >= '0' && c <= '9') 477112158Sdas L = 10*L + c - '0'; 478112158Sdas if (s - s1 > 8 || L > 19999) 479112158Sdas /* Avoid confusion from exponents 480112158Sdas * so large that e might overflow. 481112158Sdas */ 482112158Sdas e = 19999; /* safe for 16 bit ints */ 483112158Sdas else 484112158Sdas e = (int)L; 485112158Sdas if (esign) 486112158Sdas e = -e; 487112158Sdas } 488112158Sdas else 489112158Sdas e = 0; 490112158Sdas } 491112158Sdas else 492112158Sdas s = s00; 493112158Sdas } 494112158Sdas if (!nd) { 495112158Sdas if (!nz && !nz0) { 496112158Sdas#ifdef INFNAN_CHECK 497112158Sdas /* Check for Nan and Infinity */ 498165743Sdas if (!decpt) 499165743Sdas switch(c) { 500112158Sdas case 'i': 501112158Sdas case 'I': 502112158Sdas if (match(&s,"nf")) { 503112158Sdas --s; 504112158Sdas if (!match(&s,"inity")) 505112158Sdas ++s; 506112158Sdas irv = STRTOG_Infinite; 507112158Sdas goto infnanexp; 508112158Sdas } 509112158Sdas break; 510112158Sdas case 'n': 511112158Sdas case 'N': 512112158Sdas if (match(&s, "an")) { 513112158Sdas irv = STRTOG_NaN; 514112158Sdas *exp = fpi->emax + 1; 515112158Sdas#ifndef No_Hex_NaN 516112158Sdas if (*s == '(') /*)*/ 517112158Sdas irv = hexnan(&s, fpi, bits); 518112158Sdas#endif 519112158Sdas goto infnanexp; 520112158Sdas } 521112158Sdas } 522112158Sdas#endif /* INFNAN_CHECK */ 523112158Sdas irv = STRTOG_NoNumber; 524112158Sdas s = s00; 525112158Sdas } 526112158Sdas goto ret; 527112158Sdas } 528112158Sdas 529112158Sdas irv = STRTOG_Normal; 530112158Sdas e1 = e -= nf; 531112158Sdas rd = 0; 532112158Sdas switch(fpi->rounding & 3) { 533112158Sdas case FPI_Round_up: 534112158Sdas rd = 2 - sign; 535112158Sdas break; 536112158Sdas case FPI_Round_zero: 537112158Sdas rd = 1; 538112158Sdas break; 539112158Sdas case FPI_Round_down: 540112158Sdas rd = 1 + sign; 541112158Sdas } 542112158Sdas 543112158Sdas /* Now we have nd0 digits, starting at s0, followed by a 544112158Sdas * decimal point, followed by nd-nd0 digits. The number we're 545112158Sdas * after is the integer represented by those digits times 546112158Sdas * 10**e */ 547112158Sdas 548112158Sdas if (!nd0) 549112158Sdas nd0 = nd; 550112158Sdas k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 551219557Sdas dval(&rv) = y; 552112158Sdas if (k > 9) 553219557Sdas dval(&rv) = tens[k - 9] * dval(&rv) + z; 554112158Sdas bd0 = 0; 555112158Sdas if (nbits <= P && nd <= DBL_DIG) { 556112158Sdas if (!e) { 557219557Sdas if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv)) 558112158Sdas goto ret; 559112158Sdas } 560112158Sdas else if (e > 0) { 561112158Sdas if (e <= Ten_pmax) { 562112158Sdas#ifdef VAX 563112158Sdas goto vax_ovfl_check; 564112158Sdas#else 565219557Sdas i = fivesbits[e] + mantbits(&rv) <= P; 566219557Sdas /* rv = */ rounded_product(dval(&rv), tens[e]); 567219557Sdas if (rvOK(&rv, fpi, exp, bits, i, rd, &irv)) 568112158Sdas goto ret; 569112158Sdas e1 -= e; 570112158Sdas goto rv_notOK; 571112158Sdas#endif 572112158Sdas } 573112158Sdas i = DBL_DIG - nd; 574112158Sdas if (e <= Ten_pmax + i) { 575112158Sdas /* A fancier test would sometimes let us do 576112158Sdas * this for larger i values. 577112158Sdas */ 578112158Sdas e2 = e - i; 579112158Sdas e1 -= i; 580219557Sdas dval(&rv) *= tens[i]; 581112158Sdas#ifdef VAX 582112158Sdas /* VAX exponent range is so narrow we must 583112158Sdas * worry about overflow here... 584112158Sdas */ 585112158Sdas vax_ovfl_check: 586219557Sdas dval(&adj) = dval(&rv); 587219557Sdas word0(&adj) -= P*Exp_msk1; 588219557Sdas /* adj = */ rounded_product(dval(&adj), tens[e2]); 589219557Sdas if ((word0(&adj) & Exp_mask) 590112158Sdas > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 591112158Sdas goto rv_notOK; 592219557Sdas word0(&adj) += P*Exp_msk1; 593219557Sdas dval(&rv) = dval(&adj); 594112158Sdas#else 595219557Sdas /* rv = */ rounded_product(dval(&rv), tens[e2]); 596112158Sdas#endif 597219557Sdas if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv)) 598112158Sdas goto ret; 599112158Sdas e1 -= e2; 600112158Sdas } 601112158Sdas } 602112158Sdas#ifndef Inaccurate_Divide 603112158Sdas else if (e >= -Ten_pmax) { 604219557Sdas /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 605219557Sdas if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv)) 606112158Sdas goto ret; 607112158Sdas e1 -= e; 608112158Sdas } 609112158Sdas#endif 610112158Sdas } 611112158Sdas rv_notOK: 612112158Sdas e1 += nd - k; 613112158Sdas 614112158Sdas /* Get starting approximation = rv * 10**e1 */ 615112158Sdas 616112158Sdas e2 = 0; 617112158Sdas if (e1 > 0) { 618112158Sdas if ( (i = e1 & 15) !=0) 619219557Sdas dval(&rv) *= tens[i]; 620112158Sdas if (e1 &= ~15) { 621112158Sdas e1 >>= 4; 622219557Sdas while(e1 >= (1 << (n_bigtens-1))) { 623219557Sdas e2 += ((word0(&rv) & Exp_mask) 624112158Sdas >> Exp_shift1) - Bias; 625219557Sdas word0(&rv) &= ~Exp_mask; 626219557Sdas word0(&rv) |= Bias << Exp_shift1; 627219557Sdas dval(&rv) *= bigtens[n_bigtens-1]; 628219557Sdas e1 -= 1 << (n_bigtens-1); 629112158Sdas } 630219557Sdas e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias; 631219557Sdas word0(&rv) &= ~Exp_mask; 632219557Sdas word0(&rv) |= Bias << Exp_shift1; 633112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 634112158Sdas if (e1 & 1) 635219557Sdas dval(&rv) *= bigtens[j]; 636112158Sdas } 637112158Sdas } 638112158Sdas else if (e1 < 0) { 639112158Sdas e1 = -e1; 640112158Sdas if ( (i = e1 & 15) !=0) 641219557Sdas dval(&rv) /= tens[i]; 642112158Sdas if (e1 &= ~15) { 643112158Sdas e1 >>= 4; 644219557Sdas while(e1 >= (1 << (n_bigtens-1))) { 645219557Sdas e2 += ((word0(&rv) & Exp_mask) 646112158Sdas >> Exp_shift1) - Bias; 647219557Sdas word0(&rv) &= ~Exp_mask; 648219557Sdas word0(&rv) |= Bias << Exp_shift1; 649219557Sdas dval(&rv) *= tinytens[n_bigtens-1]; 650219557Sdas e1 -= 1 << (n_bigtens-1); 651112158Sdas } 652219557Sdas e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias; 653219557Sdas word0(&rv) &= ~Exp_mask; 654219557Sdas word0(&rv) |= Bias << Exp_shift1; 655112158Sdas for(j = 0; e1 > 0; j++, e1 >>= 1) 656112158Sdas if (e1 & 1) 657219557Sdas dval(&rv) *= tinytens[j]; 658112158Sdas } 659112158Sdas } 660165743Sdas#ifdef IBM 661165743Sdas /* e2 is a correction to the (base 2) exponent of the return 662165743Sdas * value, reflecting adjustments above to avoid overflow in the 663165743Sdas * native arithmetic. For native IBM (base 16) arithmetic, we 664165743Sdas * must multiply e2 by 4 to change from base 16 to 2. 665165743Sdas */ 666165743Sdas e2 <<= 2; 667165743Sdas#endif 668219557Sdas rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */ 669112158Sdas rve += e2; 670112158Sdas if ((j = rvbits - nbits) > 0) { 671112158Sdas rshift(rvb, j); 672112158Sdas rvbits = nbits; 673112158Sdas rve += j; 674112158Sdas } 675112158Sdas bb0 = 0; /* trailing zero bits in rvb */ 676112158Sdas e2 = rve + rvbits - nbits; 677165743Sdas if (e2 > fpi->emax + 1) 678165743Sdas goto huge; 679112158Sdas rve1 = rve + rvbits - nbits; 680112158Sdas if (e2 < (emin = fpi->emin)) { 681112158Sdas denorm = 1; 682112158Sdas j = rve - emin; 683112158Sdas if (j > 0) { 684112158Sdas rvb = lshift(rvb, j); 685112158Sdas rvbits += j; 686112158Sdas } 687112158Sdas else if (j < 0) { 688112158Sdas rvbits += j; 689112158Sdas if (rvbits <= 0) { 690112158Sdas if (rvbits < -1) { 691112158Sdas ufl: 692112158Sdas rvb->wds = 0; 693112158Sdas rvb->x[0] = 0; 694112158Sdas *exp = emin; 695112158Sdas irv = STRTOG_Underflow | STRTOG_Inexlo; 696112158Sdas goto ret; 697112158Sdas } 698112158Sdas rvb->x[0] = rvb->wds = rvbits = 1; 699112158Sdas } 700112158Sdas else 701112158Sdas rshift(rvb, -j); 702112158Sdas } 703112158Sdas rve = rve1 = emin; 704112158Sdas if (sudden_underflow && e2 + 1 < emin) 705112158Sdas goto ufl; 706112158Sdas } 707112158Sdas 708112158Sdas /* Now the hard part -- adjusting rv to the correct value.*/ 709112158Sdas 710112158Sdas /* Put digits into bd: true value = bd * 10^e */ 711112158Sdas 712187808Sdas bd0 = s2b(s0, nd0, nd, y, dplen); 713112158Sdas 714112158Sdas for(;;) { 715112158Sdas bd = Balloc(bd0->k); 716112158Sdas Bcopy(bd, bd0); 717112158Sdas bb = Balloc(rvb->k); 718112158Sdas Bcopy(bb, rvb); 719112158Sdas bbbits = rvbits - bb0; 720112158Sdas bbe = rve + bb0; 721112158Sdas bs = i2b(1); 722112158Sdas 723112158Sdas if (e >= 0) { 724112158Sdas bb2 = bb5 = 0; 725112158Sdas bd2 = bd5 = e; 726112158Sdas } 727112158Sdas else { 728112158Sdas bb2 = bb5 = -e; 729112158Sdas bd2 = bd5 = 0; 730112158Sdas } 731112158Sdas if (bbe >= 0) 732112158Sdas bb2 += bbe; 733112158Sdas else 734112158Sdas bd2 -= bbe; 735112158Sdas bs2 = bb2; 736112158Sdas j = nbits + 1 - bbbits; 737112158Sdas i = bbe + bbbits - nbits; 738112158Sdas if (i < emin) /* denormal */ 739112158Sdas j += i - emin; 740112158Sdas bb2 += j; 741112158Sdas bd2 += j; 742112158Sdas i = bb2 < bd2 ? bb2 : bd2; 743112158Sdas if (i > bs2) 744112158Sdas i = bs2; 745112158Sdas if (i > 0) { 746112158Sdas bb2 -= i; 747112158Sdas bd2 -= i; 748112158Sdas bs2 -= i; 749112158Sdas } 750112158Sdas if (bb5 > 0) { 751112158Sdas bs = pow5mult(bs, bb5); 752112158Sdas bb1 = mult(bs, bb); 753112158Sdas Bfree(bb); 754112158Sdas bb = bb1; 755112158Sdas } 756112158Sdas bb2 -= bb0; 757112158Sdas if (bb2 > 0) 758112158Sdas bb = lshift(bb, bb2); 759112158Sdas else if (bb2 < 0) 760112158Sdas rshift(bb, -bb2); 761112158Sdas if (bd5 > 0) 762112158Sdas bd = pow5mult(bd, bd5); 763112158Sdas if (bd2 > 0) 764112158Sdas bd = lshift(bd, bd2); 765112158Sdas if (bs2 > 0) 766112158Sdas bs = lshift(bs, bs2); 767112158Sdas asub = 1; 768112158Sdas inex = STRTOG_Inexhi; 769112158Sdas delta = diff(bb, bd); 770112158Sdas if (delta->wds <= 1 && !delta->x[0]) 771112158Sdas break; 772112158Sdas dsign = delta->sign; 773112158Sdas delta->sign = finished = 0; 774112158Sdas L = 0; 775112158Sdas i = cmp(delta, bs); 776112158Sdas if (rd && i <= 0) { 777112158Sdas irv = STRTOG_Normal; 778112158Sdas if ( (finished = dsign ^ (rd&1)) !=0) { 779112158Sdas if (dsign != 0) { 780112158Sdas irv |= STRTOG_Inexhi; 781112158Sdas goto adj1; 782112158Sdas } 783112158Sdas irv |= STRTOG_Inexlo; 784112158Sdas if (rve1 == emin) 785112158Sdas goto adj1; 786112158Sdas for(i = 0, j = nbits; j >= ULbits; 787112158Sdas i++, j -= ULbits) { 788112158Sdas if (rvb->x[i] & ALL_ON) 789112158Sdas goto adj1; 790112158Sdas } 791112158Sdas if (j > 1 && lo0bits(rvb->x + i) < j - 1) 792112158Sdas goto adj1; 793112158Sdas rve = rve1 - 1; 794112158Sdas rvb = set_ones(rvb, rvbits = nbits); 795112158Sdas break; 796112158Sdas } 797112158Sdas irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; 798112158Sdas break; 799112158Sdas } 800112158Sdas if (i < 0) { 801112158Sdas /* Error is less than half an ulp -- check for 802112158Sdas * special case of mantissa a power of two. 803112158Sdas */ 804112158Sdas irv = dsign 805112158Sdas ? STRTOG_Normal | STRTOG_Inexlo 806112158Sdas : STRTOG_Normal | STRTOG_Inexhi; 807112158Sdas if (dsign || bbbits > 1 || denorm || rve1 == emin) 808112158Sdas break; 809112158Sdas delta = lshift(delta,1); 810112158Sdas if (cmp(delta, bs) > 0) { 811112158Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 812112158Sdas goto drop_down; 813112158Sdas } 814112158Sdas break; 815112158Sdas } 816112158Sdas if (i == 0) { 817112158Sdas /* exactly half-way between */ 818112158Sdas if (dsign) { 819112158Sdas if (denorm && all_on(rvb, rvbits)) { 820112158Sdas /*boundary case -- increment exponent*/ 821112158Sdas rvb->wds = 1; 822112158Sdas rvb->x[0] = 1; 823112158Sdas rve = emin + nbits - (rvbits = 1); 824112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 825112158Sdas denorm = 0; 826112158Sdas break; 827112158Sdas } 828112158Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 829112158Sdas } 830112158Sdas else if (bbbits == 1) { 831112158Sdas irv = STRTOG_Normal; 832112158Sdas drop_down: 833112158Sdas /* boundary case -- decrement exponent */ 834112158Sdas if (rve1 == emin) { 835112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 836112158Sdas if (rvb->wds == 1 && rvb->x[0] == 1) 837112158Sdas sudden_underflow = 1; 838112158Sdas break; 839112158Sdas } 840112158Sdas rve -= nbits; 841112158Sdas rvb = set_ones(rvb, rvbits = nbits); 842112158Sdas break; 843112158Sdas } 844112158Sdas else 845112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 846219557Sdas if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1)) 847112158Sdas break; 848112158Sdas if (dsign) { 849112158Sdas rvb = increment(rvb); 850182709Sdas j = kmask & (ULbits - (rvbits & kmask)); 851182709Sdas if (hi0bits(rvb->x[rvb->wds - 1]) != j) 852112158Sdas rvbits++; 853112158Sdas irv = STRTOG_Normal | STRTOG_Inexhi; 854112158Sdas } 855112158Sdas else { 856112158Sdas if (bbbits == 1) 857112158Sdas goto undfl; 858112158Sdas decrement(rvb); 859112158Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 860112158Sdas } 861112158Sdas break; 862112158Sdas } 863219557Sdas if ((dval(&adj) = ratio(delta, bs)) <= 2.) { 864112158Sdas adj1: 865112158Sdas inex = STRTOG_Inexlo; 866112158Sdas if (dsign) { 867112158Sdas asub = 0; 868112158Sdas inex = STRTOG_Inexhi; 869112158Sdas } 870112158Sdas else if (denorm && bbbits <= 1) { 871112158Sdas undfl: 872112158Sdas rvb->wds = 0; 873112158Sdas rve = emin; 874112158Sdas irv = STRTOG_Underflow | STRTOG_Inexlo; 875112158Sdas break; 876112158Sdas } 877219557Sdas adj0 = dval(&adj) = 1.; 878112158Sdas } 879112158Sdas else { 880219557Sdas adj0 = dval(&adj) *= 0.5; 881112158Sdas if (dsign) { 882112158Sdas asub = 0; 883112158Sdas inex = STRTOG_Inexlo; 884112158Sdas } 885219557Sdas if (dval(&adj) < 2147483647.) { 886112158Sdas L = adj0; 887112158Sdas adj0 -= L; 888112158Sdas switch(rd) { 889112158Sdas case 0: 890112158Sdas if (adj0 >= .5) 891112158Sdas goto inc_L; 892112158Sdas break; 893112158Sdas case 1: 894112158Sdas if (asub && adj0 > 0.) 895112158Sdas goto inc_L; 896112158Sdas break; 897112158Sdas case 2: 898112158Sdas if (!asub && adj0 > 0.) { 899112158Sdas inc_L: 900112158Sdas L++; 901112158Sdas inex = STRTOG_Inexact - inex; 902112158Sdas } 903112158Sdas } 904219557Sdas dval(&adj) = L; 905112158Sdas } 906112158Sdas } 907112158Sdas y = rve + rvbits; 908112158Sdas 909219557Sdas /* adj *= ulp(dval(&rv)); */ 910112158Sdas /* if (asub) rv -= adj; else rv += adj; */ 911112158Sdas 912112158Sdas if (!denorm && rvbits < nbits) { 913112158Sdas rvb = lshift(rvb, j = nbits - rvbits); 914112158Sdas rve -= j; 915112158Sdas rvbits = nbits; 916112158Sdas } 917219557Sdas ab = d2b(dval(&adj), &abe, &abits); 918112158Sdas if (abe < 0) 919112158Sdas rshift(ab, -abe); 920112158Sdas else if (abe > 0) 921112158Sdas ab = lshift(ab, abe); 922112158Sdas rvb0 = rvb; 923112158Sdas if (asub) { 924112158Sdas /* rv -= adj; */ 925112158Sdas j = hi0bits(rvb->x[rvb->wds-1]); 926112158Sdas rvb = diff(rvb, ab); 927112158Sdas k = rvb0->wds - 1; 928112158Sdas if (denorm) 929112158Sdas /* do nothing */; 930112158Sdas else if (rvb->wds <= k 931112158Sdas || hi0bits( rvb->x[k]) > 932112158Sdas hi0bits(rvb0->x[k])) { 933112158Sdas /* unlikely; can only have lost 1 high bit */ 934112158Sdas if (rve1 == emin) { 935112158Sdas --rvbits; 936112158Sdas denorm = 1; 937112158Sdas } 938112158Sdas else { 939112158Sdas rvb = lshift(rvb, 1); 940112158Sdas --rve; 941112158Sdas --rve1; 942112158Sdas L = finished = 0; 943112158Sdas } 944112158Sdas } 945112158Sdas } 946112158Sdas else { 947112158Sdas rvb = sum(rvb, ab); 948112158Sdas k = rvb->wds - 1; 949112158Sdas if (k >= rvb0->wds 950112158Sdas || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) { 951112158Sdas if (denorm) { 952112158Sdas if (++rvbits == nbits) 953112158Sdas denorm = 0; 954112158Sdas } 955112158Sdas else { 956112158Sdas rshift(rvb, 1); 957112158Sdas rve++; 958112158Sdas rve1++; 959112158Sdas L = 0; 960112158Sdas } 961112158Sdas } 962112158Sdas } 963112158Sdas Bfree(ab); 964112158Sdas Bfree(rvb0); 965112158Sdas if (finished) 966112158Sdas break; 967112158Sdas 968112158Sdas z = rve + rvbits; 969112158Sdas if (y == z && L) { 970112158Sdas /* Can we stop now? */ 971219557Sdas tol = dval(&adj) * 5e-16; /* > max rel error */ 972219557Sdas dval(&adj) = adj0 - .5; 973219557Sdas if (dval(&adj) < -tol) { 974112158Sdas if (adj0 > tol) { 975112158Sdas irv |= inex; 976112158Sdas break; 977112158Sdas } 978112158Sdas } 979219557Sdas else if (dval(&adj) > tol && adj0 < 1. - tol) { 980112158Sdas irv |= inex; 981112158Sdas break; 982112158Sdas } 983112158Sdas } 984112158Sdas bb0 = denorm ? 0 : trailz(rvb); 985112158Sdas Bfree(bb); 986112158Sdas Bfree(bd); 987112158Sdas Bfree(bs); 988112158Sdas Bfree(delta); 989112158Sdas } 990165743Sdas if (!denorm && (j = nbits - rvbits)) { 991165743Sdas if (j > 0) 992165743Sdas rvb = lshift(rvb, j); 993165743Sdas else 994165743Sdas rshift(rvb, -j); 995112158Sdas rve -= j; 996112158Sdas } 997112158Sdas *exp = rve; 998112158Sdas Bfree(bb); 999112158Sdas Bfree(bd); 1000112158Sdas Bfree(bs); 1001112158Sdas Bfree(bd0); 1002112158Sdas Bfree(delta); 1003165743Sdas if (rve > fpi->emax) { 1004182709Sdas switch(fpi->rounding & 3) { 1005182709Sdas case FPI_Round_near: 1006182709Sdas goto huge; 1007182709Sdas case FPI_Round_up: 1008182709Sdas if (!sign) 1009182709Sdas goto huge; 1010182709Sdas break; 1011182709Sdas case FPI_Round_down: 1012182709Sdas if (sign) 1013182709Sdas goto huge; 1014182709Sdas } 1015182709Sdas /* Round to largest representable magnitude */ 1016182709Sdas Bfree(rvb); 1017182709Sdas rvb = 0; 1018182709Sdas irv = STRTOG_Normal | STRTOG_Inexlo; 1019182709Sdas *exp = fpi->emax; 1020182709Sdas b = bits; 1021187808Sdas be = b + ((fpi->nbits + 31) >> 5); 1022182709Sdas while(b < be) 1023182709Sdas *b++ = -1; 1024182709Sdas if ((j = fpi->nbits & 0x1f)) 1025182709Sdas *--be >>= (32 - j); 1026182709Sdas goto ret; 1027165743Sdas huge: 1028165743Sdas rvb->wds = 0; 1029165743Sdas irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 1030165743Sdas#ifndef NO_ERRNO 1031165743Sdas errno = ERANGE; 1032165743Sdas#endif 1033165743Sdas infnanexp: 1034165743Sdas *exp = fpi->emax + 1; 1035165743Sdas } 1036112158Sdas ret: 1037112158Sdas if (denorm) { 1038112158Sdas if (sudden_underflow) { 1039112158Sdas rvb->wds = 0; 1040112158Sdas irv = STRTOG_Underflow | STRTOG_Inexlo; 1041182709Sdas#ifndef NO_ERRNO 1042182709Sdas errno = ERANGE; 1043182709Sdas#endif 1044112158Sdas } 1045112158Sdas else { 1046112158Sdas irv = (irv & ~STRTOG_Retmask) | 1047112158Sdas (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); 1048182709Sdas if (irv & STRTOG_Inexact) { 1049112158Sdas irv |= STRTOG_Underflow; 1050182709Sdas#ifndef NO_ERRNO 1051182709Sdas errno = ERANGE; 1052182709Sdas#endif 1053182709Sdas } 1054112158Sdas } 1055112158Sdas } 1056112158Sdas if (se) 1057112158Sdas *se = (char *)s; 1058112158Sdas if (sign) 1059112158Sdas irv |= STRTOG_Neg; 1060112158Sdas if (rvb) { 1061112158Sdas copybits(bits, nbits, rvb); 1062112158Sdas Bfree(rvb); 1063112158Sdas } 1064112158Sdas return irv; 1065112158Sdas } 1066