sys/src/libstdio/dtoa.c

   1 /* derived from /netlib/fp/dtoa.c assuming IEEE, Standard C */
   2 /* kudos to dmg@bell-labs.com, gripes to ehg@bell-labs.com */
   3
   4 /* Let x be the exact mathematical number defined by a decimal
   5  *      string s.  Then atof(s) is the round-nearest-even IEEE
   6  *      floating point value.
   7  * Let y be an IEEE floating point value and let s be the string
   8  *      printed as %.17g.  Then atof(s) is exactly y.
   9  */
  10 #include <u.h>
  11 #include <libc.h>
  12
  13 static Lock _dtoalk[2];
  14 #define ACQUIRE_DTOA_LOCK(n)    lock(&_dtoalk[n])
  15 #define FREE_DTOA_LOCK(n)       unlock(&_dtoalk[n])
  16 #define PRIVATE_mem ((2000+sizeof(double)-1)/sizeof(double))
  17 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
  18 #define FLT_ROUNDS      1
  19 #define DBL_DIG         15
  20 #define DBL_MAX_10_EXP  308
  21 #define DBL_MAX_EXP     1024
  22 #define FLT_RADIX       2
  23 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
  24 #define fpword0(x) ((x).hi)
  25 #define fpword1(x) ((x).lo)
  26
  27 /* Ten_pmax = floor(P*log(2)/log(5)) */
  28 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  29 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  30 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  31
  32 #define Exp_shift  20
  33 #define Exp_shift1 20
  34 #define Exp_msk1    0x100000
  35 #define Exp_msk11   0x100000
  36 #define Exp_mask  0x7ff00000
  37 #define P 53
  38 #define Bias 1023
  39 #define Emin (-1022)
  40 #define Exp_1  0x3ff00000
  41 #define Exp_11 0x3ff00000
  42 #define Ebits 11
  43 #define Frac_mask  0xfffff
  44 #define Frac_mask1 0xfffff
  45 #define Ten_pmax 22
  46 #define Bletch 0x10
  47 #define Bndry_mask  0xfffff
  48 #define Bndry_mask1 0xfffff
  49 #define Sign_bit 0x80000000
  50 #define Log2P 1
  51 #define Tiny0 0
  52 #define Tiny1 1
  53 #define Quick_max 14
  54 #define Int_max 14
  55
  56 #define rounded_product(a,b) a *= b
  57 #define rounded_quotient(a,b) a /= b
  58
  59 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  60 #define Big1 0xffffffff
  61
  62 #define FFFFFFFF 0xffffffffUL
  63
  64 #undef ULint
  65
  66 #define Kmax 15
  67
  68 struct
  69 Bigint {
  70         struct Bigint *next;
  71         int     k, maxwds, sign, wds;
  72         unsigned int x[1];
  73 };
  74
  75 typedef struct Bigint Bigint;
  76
  77 static Bigint *freelist[Kmax+1];
  78
  79 static Bigint *
  80 Balloc(int k)
  81 {
  82         int     x;
  83         Bigint * rv;
  84         unsigned int    len;
  85
  86         assert(k < nelem(freelist));
  87
  88         ACQUIRE_DTOA_LOCK(0);
  89         if (rv = freelist[k]) {
  90                 freelist[k] = rv->next;
  91         } else {
  92                 x = 1 << k;
  93                 len = (sizeof(Bigint) + (x - 1) * sizeof(unsigned int) + sizeof(double) -1)
  94                  / sizeof(double);
  95                 if (pmem_next - private_mem + len <= PRIVATE_mem) {
  96                         rv = (Bigint * )pmem_next;
  97                         pmem_next += len;
  98                 } else
  99                         rv = (Bigint * )malloc(len * sizeof(double));
 100                 rv->k = k;
 101                 rv->maxwds = x;
 102         }
 103         FREE_DTOA_LOCK(0);
 104         rv->sign = rv->wds = 0;
 105         return rv;
 106 }
 107
 108 static void
 109 Bfree(Bigint *v)
 110 {
 111         if (v) {
 112                 ACQUIRE_DTOA_LOCK(0);
 113                 v->next = freelist[v->k];
 114                 freelist[v->k] = v;
 115                 FREE_DTOA_LOCK(0);
 116         }
 117 }
 118
 119 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
 120 y->wds*sizeof(int) + 2*sizeof(int))
 121
 122 static Bigint *
 123 multadd(Bigint *b, int m, int a)        /* multiply by m and add a */
 124 {
 125         int     i, wds;
 126         unsigned int carry, *x, y;
 127         unsigned int xi, z;
 128         Bigint * b1;
 129
 130         wds = b->wds;
 131         x = b->x;
 132         i = 0;
 133         carry = a;
 134         do {
 135                 xi = *x;
 136                 y = (xi & 0xffff) * m + carry;
 137                 z = (xi >> 16) * m + (y >> 16);
 138                 carry = z >> 16;
 139                 *x++ = (z << 16) + (y & 0xffff);
 140         } while (++i < wds);
 141         if (carry) {
 142                 if (wds >= b->maxwds) {
 143                         b1 = Balloc(b->k + 1);
 144                         Bcopy(b1, b);
 145                         Bfree(b);
 146                         b = b1;
 147                 }
 148                 b->x[wds++] = carry;
 149                 b->wds = wds;
 150         }
 151         return b;
 152 }
 153
 154 static Bigint *
 155 s2b(const char *s, int nd0, int nd, unsigned int y9)
 156 {
 157         Bigint * b;
 158         int     i, k;
 159         int x, y;
 160
 161         x = (nd + 8) / 9;
 162         for (k = 0, y = 1; x > y; y <<= 1, k++)
 163                 ;
 164         b = Balloc(k);
 165         b->x[0] = y9;
 166         b->wds = 1;
 167
 168         i = 9;
 169         if (9 < nd0) {
 170                 s += 9;
 171                 do
 172                         b = multadd(b, 10, *s++ - '0');
 173                 while (++i < nd0);
 174                 s++;
 175         } else
 176                 s += 10;
 177         for (; i < nd; i++)
 178                 b = multadd(b, 10, *s++ - '0');
 179         return b;
 180 }
 181
 182 static int
 183 hi0bits(register unsigned int x)
 184 {
 185         register int    k = 0;
 186
 187         if (!(x & 0xffff0000)) {
 188                 k = 16;
 189                 x <<= 16;
 190         }
 191         if (!(x & 0xff000000)) {
 192                 k += 8;
 193                 x <<= 8;
 194         }
 195         if (!(x & 0xf0000000)) {
 196                 k += 4;
 197                 x <<= 4;
 198         }
 199         if (!(x & 0xc0000000)) {
 200                 k += 2;
 201                 x <<= 2;
 202         }
 203         if (!(x & 0x80000000)) {
 204                 k++;
 205                 if (!(x & 0x40000000))
 206                         return 32;
 207         }
 208         return k;
 209 }
 210
 211 static int
 212 lo0bits(unsigned int *y)
 213 {
 214         register int    k;
 215         register unsigned int x = *y;
 216
 217         if (x & 7) {
 218                 if (x & 1)
 219                         return 0;
 220                 if (x & 2) {
 221                         *y = x >> 1;
 222                         return 1;
 223                 }
 224                 *y = x >> 2;
 225                 return 2;
 226         }
 227         k = 0;
 228         if (!(x & 0xffff)) {
 229                 k = 16;
 230                 x >>= 16;
 231         }
 232         if (!(x & 0xff)) {
 233                 k += 8;
 234                 x >>= 8;
 235         }
 236         if (!(x & 0xf)) {
 237                 k += 4;
 238                 x >>= 4;
 239         }
 240         if (!(x & 0x3)) {
 241                 k += 2;
 242                 x >>= 2;
 243         }
 244         if (!(x & 1)) {
 245                 k++;
 246                 x >>= 1;
 247                 if (!x & 1)
 248                         return 32;
 249         }
 250         *y = x;
 251         return k;
 252 }
 253
 254 static Bigint *
 255 i2b(int i)
 256 {
 257         Bigint * b;
 258
 259         b = Balloc(1);
 260         b->x[0] = i;
 261         b->wds = 1;
 262         return b;
 263 }
 264
 265 static Bigint *
 266 mult(Bigint *a, Bigint *b)
 267 {
 268         Bigint * c;
 269         int     k, wa, wb, wc;
 270         unsigned int * x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 271         unsigned int y;
 272         unsigned int carry, z;
 273         unsigned int z2;
 274
 275         if (a->wds < b->wds) {
 276                 c = a;
 277                 a = b;
 278                 b = c;
 279         }
 280         k = a->k;
 281         wa = a->wds;
 282         wb = b->wds;
 283         wc = wa + wb;
 284         if (wc > a->maxwds)
 285                 k++;
 286         c = Balloc(k);
 287         for (x = c->x, xa = x + wc; x < xa; x++)
 288                 *x = 0;
 289         xa = a->x;
 290         xae = xa + wa;
 291         xb = b->x;
 292         xbe = xb + wb;
 293         xc0 = c->x;
 294         for (; xb < xbe; xb++, xc0++) {
 295                 if (y = *xb & 0xffff) {
 296                         x = xa;
 297                         xc = xc0;
 298                         carry = 0;
 299                         do {
 300                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 301                                 carry = z >> 16;
 302                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 303                                 carry = z2 >> 16;
 304                                 Storeinc(xc, z2, z);
 305                         } while (x < xae);
 306                         *xc = carry;
 307                 }
 308                 if (y = *xb >> 16) {
 309                         x = xa;
 310                         xc = xc0;
 311                         carry = 0;
 312                         z2 = *xc;
 313                         do {
 314                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 315                                 carry = z >> 16;
 316                                 Storeinc(xc, z, z2);
 317                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 318                                 carry = z2 >> 16;
 319                         } while (x < xae);
 320                         *xc = z2;
 321                 }
 322         }
 323         for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc)
 324                 ;
 325         c->wds = wc;
 326         return c;
 327 }
 328
 329 static Bigint *p5s;
 330
 331 static Bigint *
 332 pow5mult(Bigint *b, int k)
 333 {
 334         Bigint * b1, *p5, *p51;
 335         int     i;
 336         static int      p05[3] = {
 337                 5, 25, 125      };
 338
 339         if (i = k & 3)
 340                 b = multadd(b, p05[i-1], 0);
 341
 342         if (!(k >>= 2))
 343                 return b;
 344         if (!(p5 = p5s)) {
 345                 /* first time */
 346                 ACQUIRE_DTOA_LOCK(1);
 347                 if (!(p5 = p5s)) {
 348                         p5 = p5s = i2b(625);
 349                         p5->next = 0;
 350                 }
 351                 FREE_DTOA_LOCK(1);
 352         }
 353         for (; ; ) {
 354                 if (k & 1) {
 355                         b1 = mult(b, p5);
 356                         Bfree(b);
 357                         b = b1;
 358                 }
 359                 if (!(k >>= 1))
 360                         break;
 361                 if (!(p51 = p5->next)) {
 362                         ACQUIRE_DTOA_LOCK(1);
 363                         if (!(p51 = p5->next)) {
 364                                 p51 = p5->next = mult(p5, p5);
 365                                 p51->next = 0;
 366                         }
 367                         FREE_DTOA_LOCK(1);
 368                 }
 369                 p5 = p51;
 370         }
 371         return b;
 372 }
 373
 374 static Bigint *
 375 lshift(Bigint *b, int k)
 376 {
 377         int     i, k1, n, n1;
 378         Bigint * b1;
 379         unsigned int * x, *x1, *xe, z;
 380
 381         n = k >> 5;
 382         k1 = b->k;
 383         n1 = n + b->wds + 1;
 384         for (i = b->maxwds; n1 > i; i <<= 1)
 385                 k1++;
 386         b1 = Balloc(k1);
 387         x1 = b1->x;
 388         for (i = 0; i < n; i++)
 389                 *x1++ = 0;
 390         x = b->x;
 391         xe = x + b->wds;
 392         if (k &= 0x1f) {
 393                 k1 = 32 - k;
 394                 z = 0;
 395                 do {
 396                         *x1++ = *x << k | z;
 397                         z = *x++ >> k1;
 398                 } while (x < xe);
 399                 if (*x1 = z)
 400                         ++n1;
 401         } else
 402                 do
 403                         *x1++ = *x++;
 404                 while (x < xe);
 405         b1->wds = n1 - 1;
 406         Bfree(b);
 407         return b1;
 408 }
 409
 410 static int
 411 cmp(Bigint *a, Bigint *b)
 412 {
 413         unsigned int * xa, *xa0, *xb, *xb0;
 414         int     i, j;
 415
 416         i = a->wds;
 417         j = b->wds;
 418         if (i -= j)
 419                 return i;
 420         xa0 = a->x;
 421         xa = xa0 + j;
 422         xb0 = b->x;
 423         xb = xb0 + j;
 424         for (; ; ) {
 425                 if (*--xa != *--xb)
 426                         return * xa < *xb ? -1 : 1;
 427                 if (xa <= xa0)
 428                         break;
 429         }
 430         return 0;
 431 }
 432
 433 static Bigint *
 434 diff(Bigint *a, Bigint *b)
 435 {
 436         Bigint * c;
 437         int     i, wa, wb;
 438         unsigned int * xa, *xae, *xb, *xbe, *xc;
 439         unsigned int borrow, y;
 440         unsigned int z;
 441
 442         i = cmp(a, b);
 443         if (!i) {
 444                 c = Balloc(0);
 445                 c->wds = 1;
 446                 c->x[0] = 0;
 447                 return c;
 448         }
 449         if (i < 0) {
 450                 c = a;
 451                 a = b;
 452                 b = c;
 453                 i = 1;
 454         } else
 455                 i = 0;
 456         c = Balloc(a->k);
 457         c->sign = i;
 458         wa = a->wds;
 459         xa = a->x;
 460         xae = xa + wa;
 461         wb = b->wds;
 462         xb = b->x;
 463         xbe = xb + wb;
 464         xc = c->x;
 465         borrow = 0;
 466         do {
 467                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
 468                 borrow = (y & 0x10000) >> 16;
 469                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
 470                 borrow = (z & 0x10000) >> 16;
 471                 Storeinc(xc, z, y);
 472         } while (xb < xbe);
 473         while (xa < xae) {
 474                 y = (*xa & 0xffff) - borrow;
 475                 borrow = (y & 0x10000) >> 16;
 476                 z = (*xa++ >> 16) - borrow;
 477                 borrow = (z & 0x10000) >> 16;
 478                 Storeinc(xc, z, y);
 479         }
 480         while (!*--xc)
 481                 wa--;
 482         c->wds = wa;
 483         return c;
 484 }
 485
 486 static double
 487 ulp(FPdbleword x)
 488 {
 489         int L;
 490         FPdbleword a;
 491
 492         L = (fpword0(x) & Exp_mask) - (P - 1) * Exp_msk1;
 493         fpword0(a) = L;
 494         fpword1(a) = 0;
 495         return a.x;
 496 }
 497
 498 static FPdbleword
 499 b2d(Bigint *a, int *e)
 500 {
 501         unsigned int * xa, *xa0, w, y, z;
 502         int     k;
 503 #define d0 fpword0(d)
 504 #define d1 fpword1(d)
 505         FPdbleword d;
 506
 507         xa0 = a->x;
 508         xa = xa0 + a->wds;
 509         y = *--xa;
 510         k = hi0bits(y);
 511         *e = 32 - k;
 512         if (k < Ebits) {
 513                 d0 = Exp_1 | y >> Ebits - k;
 514                 w = xa > xa0 ? *--xa : 0;
 515                 d1 = y << (32 - Ebits) + k | w >> Ebits - k;
 516                 goto ret_d;
 517         }
 518         z = xa > xa0 ? *--xa : 0;
 519         if (k -= Ebits) {
 520                 d0 = Exp_1 | y << k | z >> 32 - k;
 521                 y = xa > xa0 ? *--xa : 0;
 522                 d1 = z << k | y >> 32 - k;
 523         } else {
 524                 d0 = Exp_1 | y;
 525                 d1 = z;
 526         }
 527 ret_d:
 528 #undef d0
 529 #undef d1
 530         return d;
 531 }
 532
 533 static Bigint *
 534 d2b(FPdbleword d, int *e, int *bits)
 535 {
 536         Bigint * b;
 537         int     de, k;
 538         unsigned int * x, y, z;
 539 #define d0 d.hi
 540 #define d1 d.lo
 541
 542         b = Balloc(1);
 543         x = b->x;
 544
 545         z = d0 & Frac_mask;
 546         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
 547         de = (int)(d0 >> Exp_shift);
 548         z |= Exp_msk11;
 549         if (y = d1) {
 550                 if (k = lo0bits(&y)) {
 551                         x[0] = y | z << 32 - k;
 552                         z >>= k;
 553                 } else
 554                         x[0] = y;
 555                 b->wds = (x[1] = z) ? 2 : 1;
 556         } else {
 557                 k = lo0bits(&z);
 558                 x[0] = z;
 559                 b->wds = 1;
 560                 k += 32;
 561         }
 562         *e = de - Bias - (P - 1) + k;
 563         *bits = P - k;
 564         return b;
 565 }
 566
 567 #undef d0
 568 #undef d1
 569
 570 static double
 571 ratio(Bigint *a, Bigint *b)
 572 {
 573         FPdbleword da, db;
 574         int     k, ka, kb;
 575
 576         da = b2d(a, &ka);
 577         db = b2d(b, &kb);
 578         k = ka - kb + 32 * (a->wds - b->wds);
 579         if (k > 0)
 580                 fpword0(da) += k * Exp_msk1;
 581         else {
 582                 k = -k;
 583                 fpword0(db) += k * Exp_msk1;
 584         }
 585         return da.x / db.x;
 586 }
 587
 588 static const double
 589 tens[] = {
 590         1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
 591         1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
 592         1e20, 1e21, 1e22
 593 };
 594
 595 static const double
 596 bigtens[] = {
 597         1e16, 1e32, 1e64, 1e128, 1e256 };
 598
 599 static const double tinytens[] = {
 600         1e-16, 1e-32, 1e-64, 1e-128,
 601         9007199254740992.e-256
 602 };
 603
 604 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
 605 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
 606 #define Scale_Bit 0x10
 607 #define n_bigtens 5
 608
 609 #define NAN_WORD0 0x7ff80000
 610
 611 #define NAN_WORD1 0
 612
 613 static int
 614 match(const char **sp, char *t)
 615 {
 616         int     c, d;
 617         const char * s = *sp;
 618
 619         while (d = *t++) {
 620                 if ((c = *++s) >= 'A' && c <= 'Z')
 621                         c += 'a' - 'A';
 622                 if (c != d)
 623                         return 0;
 624         }
 625         *sp = s + 1;
 626         return 1;
 627 }
 628
 629 static void
 630 gethex(FPdbleword *rvp, const char **sp)
 631 {
 632         unsigned int c, x[2];
 633         const char * s;
 634         int     havedig, udx0, xshift;
 635
 636         x[0] = x[1] = 0;
 637         havedig = xshift = 0;
 638         udx0 = 1;
 639         s = *sp;
 640         while (c = *(const unsigned char * )++s) {
 641                 if (c >= '0' && c <= '9')
 642                         c -= '0';
 643                 else if (c >= 'a' && c <= 'f')
 644                         c += 10 - 'a';
 645                 else if (c >= 'A' && c <= 'F')
 646                         c += 10 - 'A';
 647                 else if (c <= ' ') {
 648                         if (udx0 && havedig) {
 649                                 udx0 = 0;
 650                                 xshift = 1;
 651                         }
 652                         continue;
 653                 } else if (/*(*/ c == ')') {
 654                         *sp = s + 1;
 655                         break;
 656                 } else
 657                         return; /* invalid form: don't change *sp */
 658                 havedig = 1;
 659                 if (xshift) {
 660                         xshift = 0;
 661                         x[0] = x[1];
 662                         x[1] = 0;
 663                 }
 664                 if (udx0)
 665                         x[0] = (x[0] << 4) | (x[1] >> 28);
 666                 x[1] = (x[1] << 4) | c;
 667         }
 668         if ((x[0] &= 0xfffff) || x[1]) {
 669                 fpword0(*rvp) = Exp_mask | x[0];
 670                 fpword1(*rvp) = x[1];
 671         }
 672 }
 673
 674 static int
 675 quorem(Bigint *b, Bigint *S)
 676 {
 677         int     n;
 678         unsigned int * bx, *bxe, q, *sx, *sxe;
 679         unsigned int borrow, carry, y, ys;
 680         unsigned int si, z, zs;
 681
 682         n = S->wds;
 683         if (b->wds < n)
 684                 return 0;
 685         sx = S->x;
 686         sxe = sx + --n;
 687         bx = b->x;
 688         bxe = bx + n;
 689         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
 690         if (q) {
 691                 borrow = 0;
 692                 carry = 0;
 693                 do {
 694                         si = *sx++;
 695                         ys = (si & 0xffff) * q + carry;
 696                         zs = (si >> 16) * q + (ys >> 16);
 697                         carry = zs >> 16;
 698                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
 699                         borrow = (y & 0x10000) >> 16;
 700                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
 701                         borrow = (z & 0x10000) >> 16;
 702                         Storeinc(bx, z, y);
 703                 } while (sx <= sxe);
 704                 if (!*bxe) {
 705                         bx = b->x;
 706                         while (--bxe > bx && !*bxe)
 707                                 --n;
 708                         b->wds = n;
 709                 }
 710         }
 711         if (cmp(b, S) >= 0) {
 712                 q++;
 713                 borrow = 0;
 714                 carry = 0;
 715                 bx = b->x;
 716                 sx = S->x;
 717                 do {
 718                         si = *sx++;
 719                         ys = (si & 0xffff) + carry;
 720                         zs = (si >> 16) + (ys >> 16);
 721                         carry = zs >> 16;
 722                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
 723                         borrow = (y & 0x10000) >> 16;
 724                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
 725                         borrow = (z & 0x10000) >> 16;
 726                         Storeinc(bx, z, y);
 727                 } while (sx <= sxe);
 728                 bx = b->x;
 729                 bxe = bx + n;
 730                 if (!*bxe) {
 731                         while (--bxe > bx && !*bxe)
 732                                 --n;
 733                         b->wds = n;
 734                 }
 735         }
 736         return q;
 737 }
 738
 739 static char     *
 740 rv_alloc(int i)
 741 {
 742         int     j, k, *r;
 743
 744         j = sizeof(unsigned int);
 745         for (k = 0;
 746             sizeof(Bigint) - sizeof(unsigned int) - sizeof(int) + j <= i;
 747             j <<= 1)
 748                 k++;
 749         r = (int * )Balloc(k);
 750         *r = k;
 751         return
 752             (char *)(r + 1);
 753 }
 754
 755 static char     *
 756 nrv_alloc(char *s, char **rve, int n)
 757 {
 758         char    *rv, *t;
 759
 760         t = rv = rv_alloc(n);
 761         while (*t = *s++)
 762                 t++;
 763         if (rve)
 764                 *rve = t;
 765         return rv;
 766 }
 767
 768 /* freedtoa(s) must be used to free values s returned by dtoa
 769  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
 770  * but for consistency with earlier versions of dtoa, it is optional
 771  * when MULTIPLE_THREADS is not defined.
 772  */
 773
 774 void
 775 freedtoa(char *s)
 776 {
 777         Bigint * b = (Bigint * )((int *)s - 1);
 778         b->maxwds = 1 << (b->k = *(int * )b);
 779         Bfree(b);
 780 }
 781
 782 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
 783  *
 784  * Inspired by "How to Print Floating-Point Numbers Accurately" by
 785  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
 786  *
 787  * Modifications:
 788  *      1. Rather than iterating, we use a simple numeric overestimate
 789  *         to determine k = floor(log10(d)).  We scale relevant
 790  *         quantities using O(log2(k)) rather than O(k) multiplications.
 791  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
 792  *         try to generate digits strictly left to right.  Instead, we
 793  *         compute with fewer bits and propagate the carry if necessary
 794  *         when rounding the final digit up.  This is often faster.
 795  *      3. Under the assumption that input will be rounded nearest,
 796  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
 797  *         That is, we allow equality in stopping tests when the
 798  *         round-nearest rule will give the same floating-point value
 799  *         as would satisfaction of the stopping test with strict
 800  *         inequality.
 801  *      4. We remove common factors of powers of 2 from relevant
 802  *         quantities.
 803  *      5. When converting floating-point integers less than 1e16,
 804  *         we use floating-point arithmetic rather than resorting
 805  *         to multiple-precision integers.
 806  *      6. When asked to produce fewer than 15 digits, we first try
 807  *         to get by with floating-point arithmetic; we resort to
 808  *         multiple-precision integer arithmetic only if we cannot
 809  *         guarantee that the floating-point calculation has given
 810  *         the correctly rounded result.  For k requested digits and
 811  *         "uniformly" distributed input, the probability is
 812  *         something like 10^(k-15) that we must resort to the int
 813  *         calculation.
 814  */
 815
 816 char    *
 817 dtoa(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
 818 {
 819         /*      Arguments ndigits, decpt, sign are similar to those
 820         of ecvt and fcvt; trailing zeros are suppressed from
 821         the returned string.  If not null, *rve is set to point
 822         to the end of the return value.  If d is +-Infinity or NaN,
 823         then *decpt is set to 9999.
 824
 825         mode:
 826                 0 ==> shortest string that yields d when read in
 827                         and rounded to nearest.
 828                 1 ==> like 0, but with Steele & White stopping rule;
 829                         e.g. with IEEE P754 arithmetic , mode 0 gives
 830                         1e23 whereas mode 1 gives 9.999999999999999e22.
 831                 2 ==> max(1,ndigits) significant digits.  This gives a
 832                         return value similar to that of ecvt, except
 833                         that trailing zeros are suppressed.
 834                 3 ==> through ndigits past the decimal point.  This
 835                         gives a return value similar to that from fcvt,
 836                         except that trailing zeros are suppressed, and
 837                         ndigits can be negative.
 838                 4-9 should give the same return values as 2-3, i.e.,
 839                         4 <= mode <= 9 ==> same return as mode
 840                         2 + (mode & 1).  These modes are mainly for
 841                         debugging; often they run slower but sometimes
 842                         faster than modes 2-3.
 843                 4,5,8,9 ==> left-to-right digit generation.
 844                 6-9 ==> don't try fast floating-point estimate
 845                         (if applicable).
 846
 847                 Values of mode other than 0-9 are treated as mode 0.
 848
 849                 Sufficient space is allocated to the return value
 850                 to hold the suppressed trailing zeros.
 851         */
 852
 853         int     bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
 854         j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
 855         spec_case, try_quick;
 856         int L;
 857         Bigint * b, *b1, *delta, *mlo=nil, *mhi, *S;
 858         double  ds;
 859         FPdbleword d, d2, eps;
 860         char    *s, *s0;
 861
 862         d.x = _d;
 863         if (fpword0(d) & Sign_bit) {
 864                 /* set sign for everything, including 0's and NaNs */
 865                 *sign = 1;
 866                 fpword0(d) &= ~Sign_bit;        /* clear sign bit */
 867         } else
 868                 *sign = 0;
 869
 870         if ((fpword0(d) & Exp_mask) == Exp_mask) {
 871                 /* Infinity or NaN */
 872                 *decpt = 9999;
 873                 if (!fpword1(d) && !(fpword0(d) & 0xfffff))
 874                         return nrv_alloc("Infinity", rve, 8);
 875                 return nrv_alloc("NaN", rve, 3);
 876         }
 877         if (!d.x) {
 878                 *decpt = 1;
 879                 return nrv_alloc("0", rve, 1);
 880         }
 881
 882         b = d2b(d, &be, &bbits);
 883         i = (int)(fpword0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
 884         d2 = d;
 885         fpword0(d2) &= Frac_mask1;
 886         fpword0(d2) |= Exp_11;
 887
 888         /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
 889                  * log10(x)      =  log(x) / log(10)
 890                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 891                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
 892                  *
 893                  * This suggests computing an approximation k to log10(d) by
 894                  *
 895                  * k = (i - Bias)*0.301029995663981
 896                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 897                  *
 898                  * We want k to be too large rather than too small.
 899                  * The error in the first-order Taylor series approximation
 900                  * is in our favor, so we just round up the constant enough
 901                  * to compensate for any error in the multiplication of
 902                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 903                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 904                  * adding 1e-13 to the constant term more than suffices.
 905                  * Hence we adjust the constant term to 0.1760912590558.
 906                  * (We could get a more accurate k by invoking log10,
 907                  *  but this is probably not worthwhile.)
 908                  */
 909
 910         i -= Bias;
 911         ds = (d2.x - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
 912         k = (int)ds;
 913         if (ds < 0. && ds != k)
 914                 k--;    /* want k = floor(ds) */
 915         k_check = 1;
 916         if (k >= 0 && k <= Ten_pmax) {
 917                 if (d.x < tens[k])
 918                         k--;
 919                 k_check = 0;
 920         }
 921         j = bbits - i - 1;
 922         if (j >= 0) {
 923                 b2 = 0;
 924                 s2 = j;
 925         } else {
 926                 b2 = -j;
 927                 s2 = 0;
 928         }
 929         if (k >= 0) {
 930                 b5 = 0;
 931                 s5 = k;
 932                 s2 += k;
 933         } else {
 934                 b2 -= k;
 935                 b5 = -k;
 936                 s5 = 0;
 937         }
 938         assert(k < 100);
 939         if (mode < 0 || mode > 9)
 940                 mode = 0;
 941         try_quick = 1;
 942         if (mode > 5) {
 943                 mode -= 4;
 944                 try_quick = 0;
 945         }
 946         leftright = 1;
 947         switch (mode) {
 948         case 0:
 949         case 1:
 950         default:
 951                 ilim = ilim1 = -1;
 952                 i = 18;
 953                 ndigits = 0;
 954                 break;
 955         case 2:
 956                 leftright = 0;
 957                 /* no break */
 958         case 4:
 959                 if (ndigits <= 0)
 960                         ndigits = 1;
 961                 ilim = ilim1 = i = ndigits;
 962                 break;
 963         case 3:
 964                 leftright = 0;
 965                 /* no break */
 966         case 5:
 967                 i = ndigits + k + 1;
 968                 ilim = i;
 969                 ilim1 = i - 1;
 970                 if (i <= 0)
 971                         i = 1;
 972         }
 973         s = s0 = rv_alloc(i);
 974
 975         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
 976
 977                 /* Try to get by with floating-point arithmetic. */
 978
 979                 i = 0;
 980                 d2 = d;
 981                 k0 = k;
 982                 ilim0 = ilim;
 983                 ieps = 2; /* conservative */
 984                 if (k > 0) {
 985                         ds = tens[k&0xf];
 986                         j = k >> 4;
 987                         if (j & Bletch) {
 988                                 /* prevent overflows */
 989                                 j &= Bletch - 1;
 990                                 d.x /= bigtens[n_bigtens-1];
 991                                 ieps++;
 992                         }
 993                         for (; j; j >>= 1, i++)
 994                                 if (j & 1) {
 995                                         ieps++;
 996                                         ds *= bigtens[i];
 997                                 }
 998                         d.x /= ds;
 999                 } else if (j1 = -k) {
1000                         d.x *= tens[j1 & 0xf];
1001                         for (j = j1 >> 4; j; j >>= 1, i++)
1002                                 if (j & 1) {
1003                                         ieps++;
1004                                         d.x *= bigtens[i];
1005                                 }
1006                 }
1007                 if (k_check && d.x < 1. && ilim > 0) {
1008                         if (ilim1 <= 0)
1009                                 goto fast_failed;
1010                         ilim = ilim1;
1011                         k--;
1012                         d.x *= 10.;
1013                         ieps++;
1014                 }
1015                 eps.x = ieps * d.x + 7.;
1016                 fpword0(eps) -= (P - 1) * Exp_msk1;
1017                 if (ilim == 0) {
1018                         S = mhi = 0;
1019                         d.x -= 5.;
1020                         if (d.x > eps.x)
1021                                 goto one_digit;
1022                         if (d.x < -eps.x)
1023                                 goto no_digits;
1024                         goto fast_failed;
1025                 }
1026                 /* Generate ilim digits, then fix them up. */
1027                 eps.x *= tens[ilim-1];
1028                 for (i = 1; ; i++, d.x *= 10.) {
1029                         L = d.x;
1030                         d.x -= L;
1031                         *s++ = '0' + (int)L;
1032                         if (i == ilim) {
1033                                 if (d.x > 0.5 + eps.x)
1034                                         goto bump_up;
1035                                 else if (d.x < 0.5 - eps.x) {
1036                                         while (*--s == '0')
1037                                                 ;
1038                                         s++;
1039                                         goto ret1;
1040                                 }
1041                                 break;
1042                         }
1043                 }
1044 fast_failed:
1045                 s = s0;
1046                 d.x = d2.x;
1047                 k = k0;
1048                 ilim = ilim0;
1049         }
1050
1051         /* Do we have a "small" integer? */
1052
1053         if (be >= 0 && k <= Int_max) {
1054                 /* Yes. */
1055                 ds = tens[k];
1056                 if (ndigits < 0 && ilim <= 0) {
1057                         S = mhi = 0;
1058                         if (ilim < 0 || d.x <= 5 * ds)
1059                                 goto no_digits;
1060                         goto one_digit;
1061                 }
1062                 for (i = 1; ; i++) {
1063                         L = d.x / ds;
1064                         d.x -= L * ds;
1065                         *s++ = '0' + (int)L;
1066                         if (i == ilim) {
1067                                 d.x += d.x;
1068                                 if (d.x > ds || d.x == ds && L & 1) {
1069 bump_up:
1070                                         while (*--s == '9')
1071                                                 if (s == s0) {
1072                                                         k++;
1073                                                         *s = '0';
1074                                                         break;
1075                                                 }
1076                                         ++ * s++;
1077                                 }
1078                                 break;
1079                         }
1080                         if (!(d.x *= 10.))
1081                                 break;
1082                 }
1083                 goto ret1;
1084         }
1085
1086         m2 = b2;
1087         m5 = b5;
1088         mhi = mlo = 0;
1089         if (leftright) {
1090                 if (mode < 2) {
1091                         i =
1092                             1 + P - bbits;
1093                 } else {
1094                         j = ilim - 1;
1095                         if (m5 >= j)
1096                                 m5 -= j;
1097                         else {
1098                                 s5 += j -= m5;
1099                                 b5 += j;
1100                                 m5 = 0;
1101                         }
1102                         if ((i = ilim) < 0) {
1103                                 m2 -= i;
1104                                 i = 0;
1105                         }
1106                 }
1107                 b2 += i;
1108                 s2 += i;
1109                 mhi = i2b(1);
1110         }
1111         if (m2 > 0 && s2 > 0) {
1112                 i = m2 < s2 ? m2 : s2;
1113                 b2 -= i;
1114                 m2 -= i;
1115                 s2 -= i;
1116         }
1117         if (b5 > 0) {
1118                 if (leftright) {
1119                         if (m5 > 0) {
1120                                 mhi = pow5mult(mhi, m5);
1121                                 b1 = mult(mhi, b);
1122                                 Bfree(b);
1123                                 b = b1;
1124                         }
1125                         if (j = b5 - m5)
1126                                 b = pow5mult(b, j);
1127                 } else
1128                         b = pow5mult(b, b5);
1129         }
1130         S = i2b(1);
1131         if (s5 > 0)
1132                 S = pow5mult(S, s5);
1133
1134         /* Check for special case that d is a normalized power of 2. */
1135
1136         spec_case = 0;
1137         if (mode < 2) {
1138                 if (!fpword1(d) && !(fpword0(d) & Bndry_mask)
1139                     ) {
1140                         /* The special case */
1141                         b2 += Log2P;
1142                         s2 += Log2P;
1143                         spec_case = 1;
1144                 }
1145         }
1146
1147         /* Arrange for convenient computation of quotients:
1148          * shift left if necessary so divisor has 4 leading 0 bits.
1149          *
1150          * Perhaps we should just compute leading 28 bits of S once
1151          * and for all and pass them and a shift to quorem, so it
1152          * can do shifts and ors to compute the numerator for q.
1153          */
1154         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
1155                 i = 32 - i;
1156         if (i > 4) {
1157                 i -= 4;
1158                 b2 += i;
1159                 m2 += i;
1160                 s2 += i;
1161         } else if (i < 4) {
1162                 i += 28;
1163                 b2 += i;
1164                 m2 += i;
1165                 s2 += i;
1166         }
1167         if (b2 > 0)
1168                 b = lshift(b, b2);
1169         if (s2 > 0)
1170                 S = lshift(S, s2);
1171         if (k_check) {
1172                 if (cmp(b, S) < 0) {
1173                         k--;
1174                         b = multadd(b, 10, 0);  /* we botched the k estimate */
1175                         if (leftright)
1176                                 mhi = multadd(mhi, 10, 0);
1177                         ilim = ilim1;
1178                 }
1179         }
1180         if (ilim <= 0 && mode > 2) {
1181                 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) {
1182                         /* no digits, fcvt style */
1183 no_digits:
1184                         k = -1 - ndigits;
1185                         goto ret;
1186                 }
1187 one_digit:
1188                 *s++ = '1';
1189                 k++;
1190                 goto ret;
1191         }
1192         if (leftright) {
1193                 if (m2 > 0)
1194                         mhi = lshift(mhi, m2);
1195
1196                 /* Compute mlo -- check for special case
1197                  * that d is a normalized power of 2.
1198                  */
1199
1200                 mlo = mhi;
1201                 if (spec_case) {
1202                         mhi = Balloc(mhi->k);
1203                         Bcopy(mhi, mlo);
1204                         mhi = lshift(mhi, Log2P);
1205                 }
1206
1207                 for (i = 1; ; i++) {
1208                         dig = quorem(b, S) + '0';
1209                         /* Do we yet have the shortest decimal string
1210                          * that will round to d?
1211                          */
1212                         j = cmp(b, mlo);
1213                         delta = diff(S, mhi);
1214                         j1 = delta->sign ? 1 : cmp(b, delta);
1215                         Bfree(delta);
1216                         if (j1 == 0 && !mode && !(fpword1(d) & 1)) {
1217                                 if (dig == '9')
1218                                         goto round_9_up;
1219                                 if (j > 0)
1220                                         dig++;
1221                                 *s++ = dig;
1222                                 goto ret;
1223                         }
1224                         if (j < 0 || j == 0 && !mode
1225                              && !(fpword1(d) & 1)
1226                             ) {
1227                                 if (j1 > 0) {
1228                                         b = lshift(b, 1);
1229                                         j1 = cmp(b, S);
1230                                         if ((j1 > 0 || j1 == 0 && dig & 1)
1231                                              && dig++ == '9')
1232                                                 goto round_9_up;
1233                                 }
1234                                 *s++ = dig;
1235                                 goto ret;
1236                         }
1237                         if (j1 > 0) {
1238                                 if (dig == '9') { /* possible if i == 1 */
1239 round_9_up:
1240                                         *s++ = '9';
1241                                         goto roundoff;
1242                                 }
1243                                 *s++ = dig + 1;
1244                                 goto ret;
1245                         }
1246                         *s++ = dig;
1247                         if (i == ilim)
1248                                 break;
1249                         b = multadd(b, 10, 0);
1250                         if (mlo == mhi)
1251                                 mlo = mhi = multadd(mhi, 10, 0);
1252                         else {
1253                                 mlo = multadd(mlo, 10, 0);
1254                                 mhi = multadd(mhi, 10, 0);
1255                         }
1256                 }
1257         } else
1258                 for (i = 1; ; i++) {
1259                         *s++ = dig = quorem(b, S) + '0';
1260                         if (i >= ilim)
1261                                 break;
1262                         b = multadd(b, 10, 0);
1263                 }
1264
1265         /* Round off last digit */
1266
1267         b = lshift(b, 1);
1268         j = cmp(b, S);
1269         if (j > 0 || j == 0 && dig & 1) {
1270 roundoff:
1271                 while (*--s == '9')
1272                         if (s == s0) {
1273                                 k++;
1274                                 *s++ = '1';
1275                                 goto ret;
1276                         }
1277                 ++ * s++;
1278         } else {
1279                 while (*--s == '0')
1280                         ;
1281                 s++;
1282         }
1283 ret:
1284         Bfree(S);
1285         if (mhi) {
1286                 if (mlo && mlo != mhi)
1287                         Bfree(mlo);
1288                 Bfree(mhi);
1289         }
1290 ret1:
1291         Bfree(b);
1292         *s = 0;
1293         *decpt = k + 1;
1294         if (rve)
1295                 *rve = s;
1296         return s0;
1297 }
1298