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