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