contrib/libio/floatconv.c

   1 /*
   2 Copyright (C) 1993, 1994 Free Software Foundation
   3
   4 This file is part of the GNU IO Library.  This library is free
   5 software; you can redistribute it and/or modify it under the
   6 terms of the GNU General Public License as published by the
   7 Free Software Foundation; either version 2, or (at your option)
   8 any later version.
   9
  10 This library is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this library; see the file COPYING.  If not, write to the Free
  17 Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  18
  19 As a special exception, if you link this library with files
  20 compiled with a GNU compiler to produce an executable, this does not cause
  21 the resulting executable to be covered by the GNU General Public License.
  22 This exception does not however invalidate any other reasons why
  23 the executable file might be covered by the GNU General Public License. */
  24
  25 #include <libioP.h>
  26 #ifdef _IO_USE_DTOA
  27 /****************************************************************
  28  *
  29  * The author of this software is David M. Gay.
  30  *
  31  * Copyright (c) 1991 by AT&T.
  32  *
  33  * Permission to use, copy, modify, and distribute this software for any
  34  * purpose without fee is hereby granted, provided that this entire notice
  35  * is included in all copies of any software which is or includes a copy
  36  * or modification of this software and in all copies of the supporting
  37  * documentation for such software.
  38  *
  39  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  40  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  41  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  42  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  43  *
  44  ***************************************************************/
  45
  46 /* Some cleaning up by Per Bothner, bothner@cygnus.com, 1992, 1993.
  47    Re-written to not need static variables
  48    (except result, result_k, HIWORD, LOWORD). */
  49
  50 /* Note that the checking of _DOUBLE_IS_32BITS is for use with the
  51    cross targets that employ the newlib ieeefp.h header.  -- brendan */
  52
  53 /* Please send bug reports to
  54         David M. Gay
  55         AT&T Bell Laboratories, Room 2C-463
  56         600 Mountain Avenue
  57         Murray Hill, NJ 07974-2070
  58         U.S.A.
  59         dmg@research.att.com or research!dmg
  60  */
  61
  62 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  63  *
  64  * This strtod returns a nearest machine number to the input decimal
  65  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  66  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  67  * biased rounding (add half and chop).
  68  *
  69  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  70  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  71  *
  72  * Modifications:
  73  *
  74  *      1. We only require IEEE, IBM, or VAX double-precision
  75  *              arithmetic (not IEEE double-extended).
  76  *      2. We get by with floating-point arithmetic in a case that
  77  *              Clinger missed -- when we're computing d * 10^n
  78  *              for a small integer d and the integer n is not too
  79  *              much larger than 22 (the maximum integer k for which
  80  *              we can represent 10^k exactly), we may be able to
  81  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  82  *      3. Rather than a bit-at-a-time adjustment of the binary
  83  *              result in the hard case, we use floating-point
  84  *              arithmetic to determine the adjustment to within
  85  *              one bit; only in really hard cases do we need to
  86  *              compute a second residual.
  87  *      4. Because of 3., we don't need a large table of powers of 10
  88  *              for ten-to-e (just some small tables, e.g. of 10^k
  89  *              for 0 <= k <= 22).
  90  */
  91
  92 /*
  93  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  94  *      significant byte has the lowest address.
  95  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  96  *      significant byte has the lowest address.
  97  * #define Sudden_Underflow for IEEE-format machines without gradual
  98  *      underflow (i.e., that flush to zero on underflow).
  99  * #define IBM for IBM mainframe-style floating-point arithmetic.
 100  * #define VAX for VAX-style floating-point arithmetic.
 101  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
 102  * #define No_leftright to omit left-right logic in fast floating-point
 103  *      computation of dtoa.
 104  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
 105  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
 106  *      that use extended-precision instructions to compute rounded
 107  *      products and quotients) with IBM.
 108  * #define ROUND_BIASED for IEEE-format with biased rounding.
 109  * #define Inaccurate_Divide for IEEE-format with correctly rounded
 110  *      products but inaccurate quotients, e.g., for Intel i860.
 111  * #define KR_headers for old-style C function headers.
 112  */
 113
 114 #ifdef DEBUG
 115 #include <stdio.h>
 116 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
 117 #endif
 118
 119 #ifdef __STDC__
 120 #include <stdlib.h>
 121 #include <string.h>
 122 #include <float.h>
 123 #define CONST const
 124 #else
 125 #define CONST
 126 #define KR_headers
 127
 128 /* In this case, we assume IEEE floats. */
 129 #define FLT_ROUNDS 1
 130 #define FLT_RADIX 2
 131 #define DBL_MANT_DIG 53
 132 #define DBL_DIG 15
 133 #define DBL_MAX_10_EXP 308
 134 #define DBL_MAX_EXP 1024
 135 #endif
 136
 137 #include <errno.h>
 138 #ifndef __MATH_H__
 139 #include <math.h>
 140 #endif
 141
 142 #ifdef Unsigned_Shifts
 143 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
 144 #else
 145 #define Sign_Extend(a,b) /*no-op*/
 146 #endif
 147
 148 #if defined(__i386__) || defined(__i860__) || defined(clipper)
 149 #define IEEE_8087
 150 #endif
 151 #if defined(MIPSEL) || defined(__alpha__)
 152 #define IEEE_8087
 153 #endif
 154 #if defined(__sparc__) || defined(sparc) || defined(MIPSEB)
 155 #define IEEE_MC68k
 156 #endif
 157
 158 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
 159
 160 #ifndef _DOUBLE_IS_32BITS
 161 #if FLT_RADIX==16
 162 #define IBM
 163 #else
 164 #if DBL_MANT_DIG==56
 165 #define VAX
 166 #else
 167 #if DBL_MANT_DIG==53 && DBL_MAX_10_EXP==308
 168 #define IEEE_Unknown
 169 #else
 170 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
 171 #endif
 172 #endif
 173 #endif
 174 #endif /* !_DOUBLE_IS_32BITS */
 175 #endif
 176
 177 typedef _G_uint32_t unsigned32;
 178
 179 union doubleword {
 180   double d;
 181   unsigned32 u[2];
 182 };
 183
 184 #ifdef IEEE_8087
 185 #define HIWORD 1
 186 #define LOWORD 0
 187 #define TEST_ENDIANNESS  /* nothing */
 188 #else
 189 #if defined(IEEE_MC68k)
 190 #define HIWORD 0
 191 #define LOWORD 1
 192 #define TEST_ENDIANNESS  /* nothing */
 193 #else
 194 static int HIWORD = -1, LOWORD;
 195 static void test_endianness()
 196 {
 197     union doubleword dw;
 198     dw.d = 10;
 199     if (dw.u[0] != 0) /* big-endian */
 200         HIWORD=0, LOWORD=1;
 201     else
 202         HIWORD=1, LOWORD=0;
 203 }
 204 #define TEST_ENDIANNESS  if (HIWORD<0) test_endianness();
 205 #endif
 206 #endif
 207
 208 #if 0
 209 union doubleword _temp;
 210 #endif
 211 #if defined(__GNUC__) && !defined(_DOUBLE_IS_32BITS)
 212 #define word0(x) ({ union doubleword _du; _du.d = (x); _du.u[HIWORD]; })
 213 #define word1(x) ({ union doubleword _du; _du.d = (x); _du.u[LOWORD]; })
 214 #define setword0(D,W) \
 215   ({ union doubleword _du; _du.d = (D); _du.u[HIWORD]=(W); (D)=_du.d; })
 216 #define setword1(D,W) \
 217   ({ union doubleword _du; _du.d = (D); _du.u[LOWORD]=(W); (D)=_du.d; })
 218 #define setwords(D,W0,W1) ({ union doubleword _du; \
 219   _du.u[HIWORD]=(W0); _du.u[LOWORD]=(W1); (D)=_du.d; })
 220 #define addword0(D,W) \
 221   ({ union doubleword _du; _du.d = (D); _du.u[HIWORD]+=(W); (D)=_du.d; })
 222 #else
 223 #define word0(x) ((unsigned32 *)&x)[HIWORD]
 224 #ifndef _DOUBLE_IS_32BITS
 225 #define word1(x) ((unsigned32 *)&x)[LOWORD]
 226 #else
 227 #define word1(x) 0
 228 #endif
 229 #define setword0(D,W) word0(D) = (W)
 230 #ifndef _DOUBLE_IS_32BITS
 231 #define setword1(D,W) word1(D) = (W)
 232 #define setwords(D,W0,W1) (setword0(D,W0),setword1(D,W1))
 233 #else
 234 #define setword1(D,W)
 235 #define setwords(D,W0,W1) (setword0(D,W0))
 236 #endif
 237 #define addword0(D,X) (word0(D) += (X))
 238 #endif
 239
 240 /* The following definition of Storeinc is appropriate for MIPS processors. */
 241 #if defined(IEEE_8087) + defined(VAX)
 242 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
 243 ((unsigned short *)a)[0] = (unsigned short)c, a++)
 244 #else
 245 #if defined(IEEE_MC68k)
 246 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
 247 ((unsigned short *)a)[1] = (unsigned short)c, a++)
 248 #else
 249 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
 250 #endif
 251 #endif
 252
 253 /* #define P DBL_MANT_DIG */
 254 /* Ten_pmax = floor(P*log(2)/log(5)) */
 255 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
 256 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
 257 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
 258
 259 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_Unknown)
 260 #define Exp_shift  20
 261 #define Exp_shift1 20
 262 #define Exp_msk1    0x100000
 263 #define Exp_msk11   0x100000
 264 #define Exp_mask  0x7ff00000
 265 #define P 53
 266 #define Bias 1023
 267 #define IEEE_Arith
 268 #define Emin (-1022)
 269 #define Exp_1  0x3ff00000
 270 #define Exp_11 0x3ff00000
 271 #define Ebits 11
 272 #define Frac_mask  0xfffff
 273 #define Frac_mask1 0xfffff
 274 #define Ten_pmax 22
 275 #define Bletch 0x10
 276 #define Bndry_mask  0xfffff
 277 #define Bndry_mask1 0xfffff
 278 #define LSB 1
 279 #define Sign_bit 0x80000000
 280 #define Log2P 1
 281 #define Tiny0 0
 282 #define Tiny1 1
 283 #define Quick_max 14
 284 #define Int_max 14
 285 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
 286 #else
 287 #undef  Sudden_Underflow
 288 #define Sudden_Underflow
 289 #ifdef IBM
 290 #define Exp_shift  24
 291 #define Exp_shift1 24
 292 #define Exp_msk1   0x1000000
 293 #define Exp_msk11  0x1000000
 294 #define Exp_mask  0x7f000000
 295 #define P 14
 296 #define Bias 65
 297 #define Exp_1  0x41000000
 298 #define Exp_11 0x41000000
 299 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
 300 #define Frac_mask  0xffffff
 301 #define Frac_mask1 0xffffff
 302 #define Bletch 4
 303 #define Ten_pmax 22
 304 #define Bndry_mask  0xefffff
 305 #define Bndry_mask1 0xffffff
 306 #define LSB 1
 307 #define Sign_bit 0x80000000
 308 #define Log2P 4
 309 #define Tiny0 0x100000
 310 #define Tiny1 0
 311 #define Quick_max 14
 312 #define Int_max 15
 313 #else /* VAX */
 314 #define Exp_shift  23
 315 #define Exp_shift1 7
 316 #define Exp_msk1    0x80
 317 #define Exp_msk11   0x800000
 318 #define Exp_mask  0x7f80
 319 #define P 56
 320 #define Bias 129
 321 #define Exp_1  0x40800000
 322 #define Exp_11 0x4080
 323 #define Ebits 8
 324 #define Frac_mask  0x7fffff
 325 #define Frac_mask1 0xffff007f
 326 #define Ten_pmax 24
 327 #define Bletch 2
 328 #define Bndry_mask  0xffff007f
 329 #define Bndry_mask1 0xffff007f
 330 #define LSB 0x10000
 331 #define Sign_bit 0x8000
 332 #define Log2P 1
 333 #define Tiny0 0x80
 334 #define Tiny1 0
 335 #define Quick_max 15
 336 #define Int_max 15
 337 #endif
 338 #endif
 339
 340 #ifndef IEEE_Arith
 341 #define ROUND_BIASED
 342 #endif
 343
 344 #ifdef RND_PRODQUOT
 345 #define rounded_product(a,b) a = rnd_prod(a, b)
 346 #define rounded_quotient(a,b) a = rnd_quot(a, b)
 347 extern double rnd_prod(double, double), rnd_quot(double, double);
 348 #else
 349 #define rounded_product(a,b) a *= b
 350 #define rounded_quotient(a,b) a /= b
 351 #endif
 352
 353 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
 354 #define Big1 0xffffffff
 355
 356 #define Kmax 15
 357
 358 /* (1<<BIGINT_MINIMUM_K) is the minimum number of words to allocate
 359    in a Bigint.  dtoa usually manages with 1<<2, and has not been
 360    known to need more than 1<<3.  */
 361
 362 #define BIGINT_MINIMUM_K 3
 363
 364 struct Bigint {
 365   struct Bigint *next;
 366   int k;                /* Parameter given to Balloc(k) */
 367   int maxwds;           /* Allocated space: equals 1<<k. */
 368   short on_stack;       /* 1 if stack-allocated. */
 369   short sign;           /* 0 if value is positive or zero; 1 if negative. */
 370   int wds;              /* Current length. */
 371   unsigned32 x[1<<BIGINT_MINIMUM_K]; /* Actually: x[maxwds] */
 372 };
 373
 374 #define BIGINT_HEADER_SIZE \
 375   (sizeof(Bigint) - (1<<BIGINT_MINIMUM_K) * sizeof(unsigned32))
 376
 377 typedef struct Bigint Bigint;
 378
 379 /* Initialize a stack-allocated Bigint. */
 380
 381 static Bigint *
 382 Binit
 383 #ifdef KR_headers
 384         (v) Bigint *v;
 385 #else
 386         (Bigint *v)
 387 #endif
 388 {
 389   v->on_stack = 1;
 390   v->k = BIGINT_MINIMUM_K;
 391   v->maxwds = 1 << BIGINT_MINIMUM_K;
 392   v->sign = v->wds = 0;
 393   return v;
 394 }
 395
 396 /* Allocate a Bigint with '1<<k' big digits. */
 397
 398 static Bigint *
 399 Balloc
 400 #ifdef KR_headers
 401         (k) int k;
 402 #else
 403         (int k)
 404 #endif
 405 {
 406   int x;
 407   Bigint *rv;
 408
 409   if (k < BIGINT_MINIMUM_K)
 410     k = BIGINT_MINIMUM_K;
 411
 412   x = 1 << k;
 413   rv = (Bigint *)
 414     malloc(BIGINT_HEADER_SIZE + x * sizeof(unsigned32));
 415   rv->k = k;
 416   rv->maxwds = x;
 417   rv->sign = rv->wds = 0;
 418   rv->on_stack = 0;
 419   return rv;
 420 }
 421
 422 static void
 423 Bfree
 424 #ifdef KR_headers
 425         (v) Bigint *v;
 426 #else
 427         (Bigint *v)
 428 #endif
 429 {
 430   if (v && !v->on_stack)
 431     free (v);
 432 }
 433
 434 static void
 435 Bcopy
 436 #ifdef KR_headers
 437         (x, y) Bigint *x, *y;
 438 #else
 439         (Bigint *x, Bigint *y)
 440 #endif
 441 {
 442   register unsigned32 *xp, *yp;
 443   register int i = y->wds;
 444   x->sign = y->sign;
 445   x->wds = i;
 446   for (xp = x->x, yp = y->x; --i >= 0; )
 447     *xp++ = *yp++;
 448 }
 449
 450 /* Make sure b has room for at least 1<<k big digits. */
 451
 452 static Bigint *
 453 Brealloc
 454 #ifdef KR_headers
 455         (b, k) Bigint *b; int k;
 456 #else
 457         (Bigint * b, int k)
 458 #endif
 459 {
 460   if (b == NULL)
 461     return Balloc(k);
 462   if (b->k >= k)
 463     return b;
 464   else
 465     {
 466       Bigint *rv = Balloc (k);
 467       Bcopy(rv, b);
 468       Bfree(b);
 469       return rv;
 470     }
 471 }
 472
 473 /* Return b*m+a.  b is modified.
 474    Assumption:  0xFFFF*m+a fits in 32 bits. */
 475
 476 static Bigint *
 477 multadd
 478 #ifdef KR_headers
 479         (b, m, a) Bigint *b; int m, a;
 480 #else
 481         (Bigint *b, int m, int a)
 482 #endif
 483 {
 484         int i, wds;
 485         unsigned32 *x, y;
 486         unsigned32 xi, z;
 487
 488         wds = b->wds;
 489         x = b->x;
 490         i = 0;
 491         do {
 492                 xi = *x;
 493                 y = (xi & 0xffff) * m + a;
 494                 z = (xi >> 16) * m + (y >> 16);
 495                 a = (int)(z >> 16);
 496                 *x++ = (z << 16) + (y & 0xffff);
 497                 }
 498                 while(++i < wds);
 499         if (a) {
 500                 if (wds >= b->maxwds)
 501                         b = Brealloc(b, b->k+1);
 502                 b->x[wds++] = a;
 503                 b->wds = wds;
 504                 }
 505         return b;
 506         }
 507
 508 static Bigint *
 509 s2b
 510 #ifdef KR_headers
 511         (result, s, nd0, nd, y9)
 512         Bigint *result; CONST char *s; int nd0, nd; unsigned32 y9;
 513 #else
 514         (Bigint *result, CONST char *s, int nd0, int nd, unsigned32 y9)
 515 #endif
 516 {
 517   int i, k;
 518   _G_int32_t x, y;
 519
 520   x = (nd + 8) / 9;
 521   for(k = 0, y = 1; x > y; y <<= 1, k++) ;
 522   result = Brealloc(result, k);
 523   result->x[0] = y9;
 524   result->wds = 1;
 525
 526   i = 9;
 527   if (9 < nd0)
 528     {
 529       s += 9;
 530       do
 531         result = multadd(result, 10, *s++ - '0');
 532       while (++i < nd0);
 533       s++;
 534     }
 535   else
 536     s += 10;
 537   for(; i < nd; i++)
 538     result = multadd(result, 10, *s++ - '0');
 539   return result;
 540 }
 541
 542 static int
 543 hi0bits
 544 #ifdef KR_headers
 545         (x) register unsigned32 x;
 546 #else
 547         (register unsigned32 x)
 548 #endif
 549 {
 550         register int k = 0;
 551
 552         if (!(x & 0xffff0000)) {
 553                 k = 16;
 554                 x <<= 16;
 555                 }
 556         if (!(x & 0xff000000)) {
 557                 k += 8;
 558                 x <<= 8;
 559                 }
 560         if (!(x & 0xf0000000)) {
 561                 k += 4;
 562                 x <<= 4;
 563                 }
 564         if (!(x & 0xc0000000)) {
 565                 k += 2;
 566                 x <<= 2;
 567                 }
 568         if (!(x & 0x80000000)) {
 569                 k++;
 570                 if (!(x & 0x40000000))
 571                         return 32;
 572                 }
 573         return k;
 574         }
 575
 576 static int
 577 lo0bits
 578 #ifdef KR_headers
 579         (y) unsigned32 *y;
 580 #else
 581         (unsigned32 *y)
 582 #endif
 583 {
 584         register int k;
 585         register unsigned32 x = *y;
 586
 587         if (x & 7) {
 588                 if (x & 1)
 589                         return 0;
 590                 if (x & 2) {
 591                         *y = x >> 1;
 592                         return 1;
 593                         }
 594                 *y = x >> 2;
 595                 return 2;
 596                 }
 597         k = 0;
 598         if (!(x & 0xffff)) {
 599                 k = 16;
 600                 x >>= 16;
 601                 }
 602         if (!(x & 0xff)) {
 603                 k += 8;
 604                 x >>= 8;
 605                 }
 606         if (!(x & 0xf)) {
 607                 k += 4;
 608                 x >>= 4;
 609                 }
 610         if (!(x & 0x3)) {
 611                 k += 2;
 612                 x >>= 2;
 613                 }
 614         if (!(x & 1)) {
 615                 k++;
 616                 x >>= 1;
 617                 if (!x & 1)
 618                         return 32;
 619                 }
 620         *y = x;
 621         return k;
 622         }
 623
 624 static Bigint *
 625 i2b
 626 #ifdef KR_headers
 627         (result, i) Bigint *result; int i;
 628 #else
 629         (Bigint* result, int i)
 630 #endif
 631 {
 632   result = Brealloc(result, 1);
 633   result->x[0] = i;
 634   result->wds = 1;
 635   return result;
 636 }
 637
 638 /* Do: c = a * b. */
 639
 640 static Bigint *
 641 mult
 642 #ifdef KR_headers
 643         (c, a, b) Bigint *a, *b, *c;
 644 #else
 645         (Bigint *c, Bigint *a, Bigint *b)
 646 #endif
 647 {
 648         int k, wa, wb, wc;
 649         unsigned32 carry, y, z;
 650         unsigned32 *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 651         unsigned32 z2;
 652         if (a->wds < b->wds) {
 653                 Bigint *tmp = a;
 654                 a = b;
 655                 b = tmp;
 656                 }
 657         k = a->k;
 658         wa = a->wds;
 659         wb = b->wds;
 660         wc = wa + wb;
 661         if (wc > a->maxwds)
 662                 k++;
 663         c = Brealloc(c, k);
 664         for(x = c->x, xa = x + wc; x < xa; x++)
 665                 *x = 0;
 666         xa = a->x;
 667         xae = xa + wa;
 668         xb = b->x;
 669         xbe = xb + wb;
 670         xc0 = c->x;
 671         for(; xb < xbe; xb++, xc0++) {
 672                 if ((y = *xb & 0xffff)) {
 673                         x = xa;
 674                         xc = xc0;
 675                         carry = 0;
 676                         do {
 677                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 678                                 carry = z >> 16;
 679                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 680                                 carry = z2 >> 16;
 681                                 Storeinc(xc, z2, z);
 682                                 }
 683                                 while(x < xae);
 684                         *xc = carry;
 685                         }
 686                 if ((y = *xb >> 16)) {
 687                         x = xa;
 688                         xc = xc0;
 689                         carry = 0;
 690                         z2 = *xc;
 691                         do {
 692                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 693                                 carry = z >> 16;
 694                                 Storeinc(xc, z, z2);
 695                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 696                                 carry = z2 >> 16;
 697                                 }
 698                                 while(x < xae);
 699                         *xc = z2;
 700                         }
 701                 }
 702         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
 703         c->wds = wc;
 704         return c;
 705         }
 706
 707 /* Returns b*(5**k).  b is modified. */
 708 /* Re-written by Per Bothner to not need a static list. */
 709
 710 static Bigint *
 711 pow5mult
 712 #ifdef KR_headers
 713         (b, k) Bigint *b; int k;
 714 #else
 715         (Bigint *b, int k)
 716 #endif
 717 {
 718   static int p05[6] = { 5, 25, 125, 625, 3125, 15625 };
 719
 720   for (; k > 6; k -= 6)
 721     b = multadd(b, 15625, 0); /* b *= 5**6 */
 722   if (k == 0)
 723     return b;
 724   else
 725     return multadd(b, p05[k-1], 0);
 726 }
 727
 728 /* Re-written by Per Bothner so shift can be in place. */
 729
 730 static Bigint *
 731 lshift
 732 #ifdef KR_headers
 733         (b, k) Bigint *b; int k;
 734 #else
 735         (Bigint *b, int k)
 736 #endif
 737 {
 738   int i;
 739   unsigned32 *x, *x1, *xe;
 740   int old_wds = b->wds;
 741   int n = k >> 5;
 742   int k1 = b->k;
 743   int n1 = n + old_wds + 1;
 744
 745   if (k == 0)
 746     return b;
 747
 748   for(i = b->maxwds; n1 > i; i <<= 1)
 749     k1++;
 750   b = Brealloc(b, k1);
 751
 752   xe = b->x; /* Source limit */
 753   x = xe + old_wds; /* Source pointer */
 754   x1 = x + n; /* Destination pointer */
 755   if (k &= 0x1f) {
 756     int k1 = 32 - k;
 757     unsigned32 z = *--x;
 758     if ((*x1 = (z >> k1)) != 0) {
 759       ++n1;
 760     }
 761     while (x > xe) {
 762       unsigned32 w = *--x;
 763       *--x1 = (z << k) | (w >> k1);
 764       z = w;
 765     }
 766     *--x1 = z << k;
 767   }
 768   else
 769     do {
 770       *--x1 = *--x;
 771     } while(x > xe);
 772   while (x1 > xe)
 773     *--x1 = 0;
 774   b->wds = n1 - 1;
 775   return b;
 776 }
 777
 778 static int
 779 cmp
 780 #ifdef KR_headers
 781         (a, b) Bigint *a, *b;
 782 #else
 783         (Bigint *a, Bigint *b)
 784 #endif
 785 {
 786         unsigned32 *xa, *xa0, *xb, *xb0;
 787         int i, j;
 788
 789         i = a->wds;
 790         j = b->wds;
 791 #ifdef DEBUG
 792         if (i > 1 && !a->x[i-1])
 793                 Bug("cmp called with a->x[a->wds-1] == 0");
 794         if (j > 1 && !b->x[j-1])
 795                 Bug("cmp called with b->x[b->wds-1] == 0");
 796 #endif
 797         if (i -= j)
 798                 return i;
 799         xa0 = a->x;
 800         xa = xa0 + j;
 801         xb0 = b->x;
 802         xb = xb0 + j;
 803         for(;;) {
 804                 if (*--xa != *--xb)
 805                         return *xa < *xb ? -1 : 1;
 806                 if (xa <= xa0)
 807                         break;
 808                 }
 809         return 0;
 810         }
 811
 812 /* Do: c = a-b. */
 813
 814 static Bigint *
 815 diff
 816 #ifdef KR_headers
 817         (c, a, b) Bigint *c, *a, *b;
 818 #else
 819         (Bigint *c, Bigint *a, Bigint *b)
 820 #endif
 821 {
 822         int i, wa, wb;
 823         _G_int32_t borrow, y; /* We need signed shifts here. */
 824         unsigned32 *xa, *xae, *xb, *xbe, *xc;
 825         _G_int32_t z;
 826
 827         i = cmp(a,b);
 828         if (!i) {
 829                 c = Brealloc(c, 0);
 830                 c->wds = 1;
 831                 c->x[0] = 0;
 832                 return c;
 833                 }
 834         if (i < 0) {
 835                 Bigint *tmp = a;
 836                 a = b;
 837                 b = tmp;
 838                 i = 1;
 839                 }
 840         else
 841                 i = 0;
 842         c = Brealloc(c, a->k);
 843         c->sign = i;
 844         wa = a->wds;
 845         xa = a->x;
 846         xae = xa + wa;
 847         wb = b->wds;
 848         xb = b->x;
 849         xbe = xb + wb;
 850         xc = c->x;
 851         borrow = 0;
 852         do {
 853                 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
 854                 borrow = y >> 16;
 855                 Sign_Extend(borrow, y);
 856                 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
 857                 borrow = z >> 16;
 858                 Sign_Extend(borrow, z);
 859                 Storeinc(xc, z, y);
 860                 }
 861                 while(xb < xbe);
 862         while(xa < xae) {
 863                 y = (*xa & 0xffff) + borrow;
 864                 borrow = y >> 16;
 865                 Sign_Extend(borrow, y);
 866                 z = (*xa++ >> 16) + borrow;
 867                 borrow = z >> 16;
 868                 Sign_Extend(borrow, z);
 869                 Storeinc(xc, z, y);
 870                 }
 871         while(!*--xc)
 872                 wa--;
 873         c->wds = wa;
 874         return c;
 875         }
 876
 877 static double
 878 ulp
 879 #ifdef KR_headers
 880         (x) double x;
 881 #else
 882         (double x)
 883 #endif
 884 {
 885         register _G_int32_t L;
 886         double a;
 887
 888         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
 889 #ifndef Sudden_Underflow
 890         if (L > 0) {
 891 #endif
 892 #ifdef IBM
 893                 L |= Exp_msk1 >> 4;
 894 #endif
 895                 setwords(a, L, 0);
 896 #ifndef Sudden_Underflow
 897                 }
 898         else {
 899                 L = -L >> Exp_shift;
 900                 if (L < Exp_shift)
 901                         setwords(a, 0x80000 >> L, 0);
 902                 else {
 903                         L -= Exp_shift;
 904                         setwords(a, 0, L >= 31 ? 1 : 1 << (31 - L));
 905                         }
 906                 }
 907 #endif
 908         return a;
 909         }
 910
 911 static double
 912 b2d
 913 #ifdef KR_headers
 914         (a, e) Bigint *a; int *e;
 915 #else
 916         (Bigint *a, int *e)
 917 #endif
 918 {
 919         unsigned32 *xa, *xa0, w, y, z;
 920         int k;
 921         double d;
 922         unsigned32 d0, d1;
 923
 924         xa0 = a->x;
 925         xa = xa0 + a->wds;
 926         y = *--xa;
 927 #ifdef DEBUG
 928         if (!y) Bug("zero y in b2d");
 929 #endif
 930         k = hi0bits(y);
 931         *e = 32 - k;
 932         if (k < Ebits) {
 933                 d0 = Exp_1 | y >> (Ebits - k);
 934                 w = xa > xa0 ? *--xa : 0;
 935 #ifndef _DOUBLE_IS_32BITS
 936                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
 937 #endif
 938                 goto ret_d;
 939                 }
 940         z = xa > xa0 ? *--xa : 0;
 941         if (k -= Ebits) {
 942                 d0 = Exp_1 | y << k | z >> (32 - k);
 943                 y = xa > xa0 ? *--xa : 0;
 944 #ifndef _DOUBLE_IS_32BITS
 945                 d1 = z << k | y >> (32 - k);
 946 #endif
 947                 }
 948         else {
 949                 d0 = Exp_1 | y;
 950 #ifndef _DOUBLE_IS_32BITS
 951                 d1 = z;
 952 #endif
 953                 }
 954  ret_d:
 955 #ifdef VAX
 956         setwords(d, d0 >> 16 | d0 << 16, d1 >> 16 | d1 << 16);
 957 #else
 958         setwords (d, d0, d1);
 959 #endif
 960         return d;
 961         }
 962
 963 static Bigint *
 964 d2b
 965 #ifdef KR_headers
 966         (result, d, e, bits) Bigint *result; double d; _G_int32_t *e, *bits;
 967 #else
 968         (Bigint *result, double d, _G_int32_t *e, _G_int32_t *bits)
 969 #endif
 970 {
 971         int de, i, k;
 972         unsigned32 *x, y, z;
 973         unsigned32 d0, d1;
 974 #ifdef VAX
 975         d0 = word0(d) >> 16 | word0(d) << 16;
 976         d1 = word1(d) >> 16 | word1(d) << 16;
 977 #else
 978         d0 = word0(d);
 979         d1 = word1(d);
 980 #endif
 981
 982         result = Brealloc(result, 1);
 983         x = result->x;
 984
 985         z = d0 & Frac_mask;
 986         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
 987
 988         de = (int)(d0 >> Exp_shift);  /* The exponent part of d. */
 989
 990         /* Put back the suppressed high-order bit, if normalized. */
 991 #ifndef IBM
 992 #ifndef Sudden_Underflow
 993         if (de)
 994 #endif
 995           z |= Exp_msk11;
 996 #endif
 997
 998 #ifndef _DOUBLE_IS_32BITS
 999         if ((y = d1)) {
1000                 if ((k = lo0bits(&y))) {
1001                         x[0] = y | z << (32 - k);
1002                         z >>= k;
1003                         }
1004                 else
1005                         x[0] = y;
1006                 i = result->wds = (x[1] = z) ? 2 : 1;
1007                 }
1008         else {
1009 #endif /* !_DOUBLE_IS_32BITS */
1010 #ifdef DEBUG
1011                 if (!z)
1012                         Bug("Zero passed to d2b");
1013 #endif
1014                 k = lo0bits(&z);
1015                 x[0] = z;
1016                 i = result->wds = 1;
1017 #ifndef _DOUBLE_IS_32BITS
1018                 k += 32;
1019                 }
1020 #endif
1021 #ifndef Sudden_Underflow
1022         if (de) {
1023 #endif
1024 #ifdef IBM
1025                 *e = (de - Bias - (P-1) << 2) + k;
1026                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1027 #else
1028                 *e = de - Bias - (P-1) + k;
1029                 *bits = P - k;
1030 #endif
1031 #ifndef Sudden_Underflow
1032                 }
1033         else {
1034                 *e = de - Bias - (P-1) + 1 + k;
1035                 *bits = 32*i - hi0bits(x[i-1]);
1036                 }
1037 #endif
1038         return result;
1039         }
1040
1041 static double
1042 ratio
1043 #ifdef KR_headers
1044         (a, b) Bigint *a, *b;
1045 #else
1046         (Bigint *a, Bigint *b)
1047 #endif
1048 {
1049         double da, db;
1050         int k, ka, kb;
1051
1052         da = b2d(a, &ka);
1053         db = b2d(b, &kb);
1054         k = ka - kb + 32*(a->wds - b->wds);
1055 #ifdef IBM
1056         if (k > 0) {
1057                 addword0(da, (k >> 2)*Exp_msk1);
1058                 if (k &= 3)
1059                         da *= 1 << k;
1060                 }
1061         else {
1062                 k = -k;
1063                 addword0(db,(k >> 2)*Exp_msk1);
1064                 if (k &= 3)
1065                         db *= 1 << k;
1066                 }
1067 #else
1068         if (k > 0)
1069                 addword0(da, k*Exp_msk1);
1070         else {
1071                 k = -k;
1072                 addword0(db, k*Exp_msk1);
1073                 }
1074 #endif
1075         return da / db;
1076         }
1077
1078 static CONST double
1079 tens[] = {
1080                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1081                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1082                 1e20, 1e21, 1e22
1083 #ifdef VAX
1084                 , 1e23, 1e24
1085 #endif
1086                 };
1087
1088 #ifdef IEEE_Arith
1089 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1090 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1091 #define n_bigtens 5
1092 #else
1093 #ifdef IBM
1094 static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1095 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1096 #define n_bigtens 3
1097 #else
1098 /* Also used for the case when !_DOUBLE_IS_32BITS.  */
1099 static CONST double bigtens[] = { 1e16, 1e32 };
1100 static CONST double tinytens[] = { 1e-16, 1e-32 };
1101 #define n_bigtens 2
1102 #endif
1103 #endif
1104
1105  double
1106 _IO_strtod
1107 #ifdef KR_headers
1108         (s00, se) CONST char *s00; char **se;
1109 #else
1110         (CONST char *s00, char **se)
1111 #endif
1112 {
1113         _G_int32_t bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1114                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1115         CONST char *s, *s0, *s1;
1116         double aadj, aadj1, adj, rv, rv0;
1117         _G_int32_t L;
1118         unsigned32 y, z;
1119         Bigint _bb, _b_avail, _bd, _bd0, _bs, _delta;
1120         Bigint *bb = Binit(&_bb);
1121         Bigint *bd = Binit(&_bd);
1122         Bigint *bd0 = Binit(&_bd0);
1123         Bigint *bs = Binit(&_bs);
1124         Bigint *b_avail = Binit(&_b_avail);
1125         Bigint *delta = Binit(&_delta);
1126
1127         TEST_ENDIANNESS;
1128         sign = nz0 = nz = 0;
1129         rv = 0.;
1130         (void)&rv;              /* Force rv into the stack */
1131         for(s = s00;;s++) switch(*s) {
1132                 case '-':
1133                         sign = 1;
1134                         /* no break */
1135                 case '+':
1136                         if (*++s)
1137                                 goto break2;
1138                         /* no break */
1139                 case 0:
1140                         /* "+" and "-" should be reported as an error? */
1141                         sign = 0;
1142                         s = s00;
1143                         goto ret;
1144                 case '\t':
1145                 case '\n':
1146                 case '\v':
1147                 case '\f':
1148                 case '\r':
1149                 case ' ':
1150                         continue;
1151                 default:
1152                         goto break2;
1153                 }
1154  break2:
1155         if (*s == '0') {
1156                 nz0 = 1;
1157                 while(*++s == '0') ;
1158                 if (!*s)
1159                         goto ret;
1160                 }
1161         s0 = s;
1162         y = z = 0;
1163         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1164                 if (nd < 9)
1165                         y = 10*y + c - '0';
1166                 else if (nd < 16)
1167                         z = 10*z + c - '0';
1168         nd0 = nd;
1169         if (c == '.') {
1170                 c = *++s;
1171                 if (!nd) {
1172                         for(; c == '0'; c = *++s)
1173                                 nz++;
1174                         if (c > '0' && c <= '9') {
1175                                 s0 = s;
1176                                 nf += nz;
1177                                 nz = 0;
1178                                 goto have_dig;
1179                                 }
1180                         goto dig_done;
1181                         }
1182                 for(; c >= '0' && c <= '9'; c = *++s) {
1183  have_dig:
1184                         nz++;
1185                         if (c -= '0') {
1186                                 nf += nz;
1187                                 for(i = 1; i < nz; i++)
1188                                         if (nd++ < 9)
1189                                                 y *= 10;
1190                                         else if (nd <= DBL_DIG + 1)
1191                                                 z *= 10;
1192                                 if (nd++ < 9)
1193                                         y = 10*y + c;
1194                                 else if (nd <= DBL_DIG + 1)
1195                                         z = 10*z + c;
1196                                 nz = 0;
1197                                 }
1198                         }
1199                 }
1200  dig_done:
1201         e = 0;
1202         if (c == 'e' || c == 'E') {
1203                 if (!nd && !nz && !nz0) {
1204                         s = s00;
1205                         goto ret;
1206                         }
1207                 s00 = s;
1208                 esign = 0;
1209                 switch(c = *++s) {
1210                         case '-':
1211                                 esign = 1;
1212                         case '+':
1213                                 c = *++s;
1214                         }
1215                 if (c >= '0' && c <= '9') {
1216                         while(c == '0')
1217                                 c = *++s;
1218                         if (c > '0' && c <= '9') {
1219                                 e = c - '0';
1220                                 s1 = s;
1221                                 while((c = *++s) >= '0' && c <= '9')
1222                                         e = 10*e + c - '0';
1223                                 if (s - s1 > 8)
1224                                         /* Avoid confusion from exponents
1225                                          * so large that e might overflow.
1226                                          */
1227                                         e = 9999999;
1228                                 if (esign)
1229                                         e = -e;
1230                                 }
1231                         else
1232                                 e = 0;
1233                         }
1234                 else
1235                         s = s00;
1236                 }
1237         if (!nd) {
1238                 if (!nz && !nz0)
1239                         s = s00;
1240                 goto ret;
1241                 }
1242         e1 = e -= nf;
1243
1244         /* Now we have nd0 digits, starting at s0, followed by a
1245          * decimal point, followed by nd-nd0 digits.  The number we're
1246          * after is the integer represented by those digits times
1247          * 10**e */
1248
1249         if (!nd0)
1250                 nd0 = nd;
1251         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1252         rv = y;
1253         if (k > 9)
1254                 rv = tens[k - 9] * rv + z;
1255         if (nd <= DBL_DIG
1256 #ifndef RND_PRODQUOT
1257                 && FLT_ROUNDS == 1
1258 #endif
1259                         ) {
1260                 if (!e)
1261                         goto ret;
1262                 if (e > 0) {
1263                         if (e <= Ten_pmax) {
1264 #ifdef VAX
1265                                 goto vax_ovfl_check;
1266 #else
1267                                 /* rv = */ rounded_product(rv, tens[e]);
1268                                 goto ret;
1269 #endif
1270                                 }
1271                         i = DBL_DIG - nd;
1272                         if (e <= Ten_pmax + i) {
1273                                 /* A fancier test would sometimes let us do
1274                                  * this for larger i values.
1275                                  */
1276                                 e -= i;
1277                                 rv *= tens[i];
1278 #ifdef VAX
1279                                 /* VAX exponent range is so narrow we must
1280                                  * worry about overflow here...
1281                                  */
1282  vax_ovfl_check:
1283                                 addword0(rv, - P*Exp_msk1);
1284                                 /* rv = */ rounded_product(rv, tens[e]);
1285                                 if ((word0(rv) & Exp_mask)
1286                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1287                                         goto ovfl;
1288                                 addword0(rv, P*Exp_msk1);
1289 #else
1290                                 /* rv = */ rounded_product(rv, tens[e]);
1291 #endif
1292                                 goto ret;
1293                                 }
1294                         }
1295 #ifndef Inaccurate_Divide
1296                 else if (e >= -Ten_pmax) {
1297                         /* rv = */ rounded_quotient(rv, tens[-e]);
1298                         goto ret;
1299                         }
1300 #endif
1301                 }
1302         e1 += nd - k;
1303
1304         /* Get starting approximation = rv * 10**e1 */
1305
1306         if (e1 > 0) {
1307                 if ((i = e1 & 15))
1308                         rv *= tens[i];
1309                 if (e1 &= ~15) {
1310                         if (e1 > DBL_MAX_10_EXP) {
1311  ovfl:
1312                                 errno = ERANGE;
1313 #if defined(sun) && !defined(__svr4__)
1314 /* SunOS defines HUGE_VAL as __infinity(), which is in libm. */
1315 #undef HUGE_VAL
1316 #endif
1317 #ifndef HUGE_VAL
1318 #define HUGE_VAL        1.7976931348623157E+308
1319 #endif
1320                                 rv = HUGE_VAL;
1321                                 goto ret;
1322                                 }
1323                         if (e1 >>= 4) {
1324                                 for(j = 0; e1 > 1; j++, e1 >>= 1)
1325                                         if (e1 & 1)
1326                                                 rv *= bigtens[j];
1327                         /* The last multiplication could overflow. */
1328                                 addword0(rv, -P*Exp_msk1);
1329                                 rv *= bigtens[j];
1330                                 if ((z = word0(rv) & Exp_mask)
1331                                  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1332                                         goto ovfl;
1333                                 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1334                                         /* set to largest number */
1335                                         /* (Can't trust DBL_MAX) */
1336                                         setwords(rv, Big0, Big1);
1337                                         }
1338                                 else
1339                                         addword0(rv, P*Exp_msk1);
1340                                 }
1341
1342                         }
1343                 }
1344         else if (e1 < 0) {
1345                 e1 = -e1;
1346                 if ((i = e1 & 15))
1347                         rv /= tens[i];
1348                 if (e1 &= ~15) {
1349                         e1 >>= 4;
1350                         for(j = 0; e1 > 1; j++, e1 >>= 1)
1351                                 if (e1 & 1)
1352                                         rv *= tinytens[j];
1353                         /* The last multiplication could underflow. */
1354                         rv0 = rv;
1355                         rv *= tinytens[j];
1356                         if (!rv) {
1357                                 rv = 2.*rv0;
1358                                 rv *= tinytens[j];
1359                                 if (!rv) {
1360  undfl:
1361                                         rv = 0.;
1362                                         errno = ERANGE;
1363                                         goto ret;
1364                                         }
1365                                 setwords(rv, Tiny0, Tiny1);
1366                                 /* The refinement below will clean
1367                                  * this approximation up.
1368                                  */
1369                                 }
1370                         }
1371                 }
1372
1373         /* Now the hard part -- adjusting rv to the correct value.*/
1374
1375         /* Put digits into bd: true value = bd * 10^e */
1376
1377         bd0 = s2b(bd0, s0, nd0, nd, y);
1378         bd = Brealloc(bd, bd0->k);
1379
1380         for(;;) {
1381                 Bcopy(bd, bd0);
1382                 bb = d2b(bb, rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
1383                 bs = i2b(bs, 1);
1384
1385                 if (e >= 0) {
1386                         bb2 = bb5 = 0;
1387                         bd2 = bd5 = e;
1388                         }
1389                 else {
1390                         bb2 = bb5 = -e;
1391                         bd2 = bd5 = 0;
1392                         }
1393                 if (bbe >= 0)
1394                         bb2 += bbe;
1395                 else
1396                         bd2 -= bbe;
1397                 bs2 = bb2;
1398 #ifdef Sudden_Underflow
1399 #ifdef IBM
1400                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1401 #else
1402                 j = P + 1 - bbbits;
1403 #endif
1404 #else
1405                 i = bbe + bbbits - 1;   /* logb(rv) */
1406                 if (i < Emin)   /* denormal */
1407                         j = bbe + (P-Emin);
1408                 else
1409                         j = P + 1 - bbbits;
1410 #endif
1411                 bb2 += j;
1412                 bd2 += j;
1413                 i = bb2 < bd2 ? bb2 : bd2;
1414                 if (i > bs2)
1415                         i = bs2;
1416                 if (i > 0) {
1417                         bb2 -= i;
1418                         bd2 -= i;
1419                         bs2 -= i;
1420                         }
1421                 if (bb5 > 0) {
1422                         Bigint *b_tmp;
1423                         bs = pow5mult(bs, bb5);
1424                         b_tmp = mult(b_avail, bs, bb);
1425                         b_avail = bb;
1426                         bb = b_tmp;
1427                         }
1428                 if (bb2 > 0)
1429                         bb = lshift(bb, bb2);
1430                 if (bd5 > 0)
1431                         bd = pow5mult(bd, bd5);
1432                 if (bd2 > 0)
1433                         bd = lshift(bd, bd2);
1434                 if (bs2 > 0)
1435                         bs = lshift(bs, bs2);
1436                 delta = diff(delta, bb, bd);
1437                 dsign = delta->sign;
1438                 delta->sign = 0;
1439                 i = cmp(delta, bs);
1440                 if (i < 0) {
1441                         /* Error is less than half an ulp -- check for
1442                          * special case of mantissa a power of two.
1443                          */
1444                         if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1445                                 break;
1446                         delta = lshift(delta,Log2P);
1447                         if (cmp(delta, bs) > 0)
1448                                 goto drop_down;
1449                         break;
1450                         }
1451                 if (i == 0) {
1452                         /* exactly half-way between */
1453                         if (dsign) {
1454                                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1455                                  &&  word1(rv) == 0xffffffff) {
1456                                         /*boundary case -- increment exponent*/
1457                                         setword0(rv, (word0(rv) & Exp_mask)
1458                                                  + Exp_msk1);
1459 #ifdef IBM
1460                                         setword0 (rv,
1461                                                   word0(rv) | (Exp_msk1 >> 4));
1462 #endif
1463                                         setword1(rv, 0);
1464                                         break;
1465                                         }
1466                                 }
1467                         else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1468  drop_down:
1469                                 /* boundary case -- decrement exponent */
1470 #ifdef Sudden_Underflow
1471                                 L = word0(rv) & Exp_mask;
1472 #ifdef IBM
1473                                 if (L <  Exp_msk1)
1474 #else
1475                                 if (L <= Exp_msk1)
1476 #endif
1477                                         goto undfl;
1478                                 L -= Exp_msk1;
1479 #else
1480                                 L = (word0(rv) & Exp_mask) - Exp_msk1;
1481 #endif
1482                                 setwords(rv, L | Bndry_mask1, 0xffffffff);
1483 #ifdef IBM
1484                                 continue;
1485 #else
1486                                 break;
1487 #endif
1488                                 }
1489 #ifndef ROUND_BIASED
1490                         if (!(word1(rv) & LSB))
1491                                 break;
1492 #endif
1493                         if (dsign)
1494                                 rv += ulp(rv);
1495 #ifndef ROUND_BIASED
1496                         else {
1497                                 rv -= ulp(rv);
1498 #ifndef Sudden_Underflow
1499                                 if (!rv)
1500                                         goto undfl;
1501 #endif
1502                                 }
1503 #endif
1504                         break;
1505                         }
1506                 if ((aadj = ratio(delta, bs)) <= 2.) {
1507                         if (dsign)
1508                                 aadj = aadj1 = 1.;
1509                         else if (word1(rv) || word0(rv) & Bndry_mask) {
1510 #ifndef Sudden_Underflow
1511                                 if (word1(rv) == Tiny1 && !word0(rv))
1512                                         goto undfl;
1513 #endif
1514                                 aadj = 1.;
1515                                 aadj1 = -1.;
1516                                 }
1517                         else {
1518                                 /* special case -- power of FLT_RADIX to be */
1519                                 /* rounded down... */
1520
1521                                 if (aadj < 2./FLT_RADIX)
1522                                         aadj = 1./FLT_RADIX;
1523                                 else
1524                                         aadj *= 0.5;
1525                                 aadj1 = -aadj;
1526                                 }
1527                         }
1528                 else {
1529                         aadj *= 0.5;
1530                         aadj1 = dsign ? aadj : -aadj;
1531 #ifdef Check_FLT_ROUNDS
1532                         switch(FLT_ROUNDS) {
1533                                 case 2: /* towards +infinity */
1534                                         aadj1 -= 0.5;
1535                                         break;
1536                                 case 0: /* towards 0 */
1537                                 case 3: /* towards -infinity */
1538                                         aadj1 += 0.5;
1539                                 }
1540 #else
1541                         if (FLT_ROUNDS == 0)
1542                                 aadj1 += 0.5;
1543 #endif
1544                         }
1545                 y = word0(rv) & Exp_mask;
1546
1547                 /* Check for overflow */
1548
1549                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1550                         rv0 = rv;
1551                         addword0(rv, - P*Exp_msk1);
1552                         adj = aadj1 * ulp(rv);
1553                         rv += adj;
1554                         if ((word0(rv) & Exp_mask) >=
1555                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1556                                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1557                                         goto ovfl;
1558                                 setwords(rv, Big0, Big1);
1559                                 continue;
1560                                 }
1561                         else
1562                                 addword0(rv, P*Exp_msk1);
1563                         }
1564                 else {
1565 #ifdef Sudden_Underflow
1566                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1567                                 rv0 = rv;
1568                                 addword0(rv, P*Exp_msk1);
1569                                 adj = aadj1 * ulp(rv);
1570                                 rv += adj;
1571 #ifdef IBM
1572                                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
1573 #else
1574                                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1575 #endif
1576                                         {
1577                                         if (word0(rv0) == Tiny0
1578                                          && word1(rv0) == Tiny1)
1579                                                 goto undfl;
1580                                         setwords(rv, Tiny0, Tiny1);
1581                                         continue;
1582                                         }
1583                                 else
1584                                         addword0(rv, -P*Exp_msk1);
1585                                 }
1586                         else {
1587                                 adj = aadj1 * ulp(rv);
1588                                 rv += adj;
1589                                 }
1590 #else
1591                         /* Compute adj so that the IEEE rounding rules will
1592                          * correctly round rv + adj in some half-way cases.
1593                          * If rv * ulp(rv) is denormalized (i.e.,
1594                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1595                          * trouble from bits lost to denormalization;
1596                          * example: 1.2e-307 .
1597                          */
1598                         if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1599                                 aadj1 = (double)(int)(aadj + 0.5);
1600                                 if (!dsign)
1601                                         aadj1 = -aadj1;
1602                                 }
1603                         adj = aadj1 * ulp(rv);
1604                         rv += adj;
1605 #endif
1606                         }
1607                 z = word0(rv) & Exp_mask;
1608                 if (y == z) {
1609                         /* Can we stop now? */
1610                         L = (_G_int32_t)aadj;
1611                         aadj -= L;
1612                         /* The tolerances below are conservative. */
1613                         if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1614                                 if (aadj < .4999999 || aadj > .5000001)
1615                                         break;
1616                                 }
1617                         else if (aadj < .4999999/FLT_RADIX)
1618                                 break;
1619                         }
1620                 }
1621         Bfree(bb);
1622         Bfree(bd);
1623         Bfree(bs);
1624         Bfree(bd0);
1625         Bfree(delta);
1626         Bfree(b_avail);
1627  ret:
1628         if (se)
1629                 *se = (char *)s;
1630         return sign ? -rv : rv;
1631         }
1632
1633 static int
1634 quorem
1635 #ifdef KR_headers
1636         (b, S) Bigint *b, *S;
1637 #else
1638         (Bigint *b, Bigint *S)
1639 #endif
1640 {
1641         int n;
1642         _G_int32_t borrow, y;
1643         unsigned32 carry, q, ys;
1644         unsigned32 *bx, *bxe, *sx, *sxe;
1645         _G_int32_t z;
1646         unsigned32 si, zs;
1647
1648         n = S->wds;
1649 #ifdef DEBUG
1650         /*debug*/ if (b->wds > n)
1651         /*debug*/       Bug("oversize b in quorem");
1652 #endif
1653         if (b->wds < n)
1654                 return 0;
1655         sx = S->x;
1656         sxe = sx + --n;
1657         bx = b->x;
1658         bxe = bx + n;
1659         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
1660 #ifdef DEBUG
1661         /*debug*/ if (q > 9)
1662         /*debug*/       Bug("oversized quotient in quorem");
1663 #endif
1664         if (q) {
1665                 borrow = 0;
1666                 carry = 0;
1667                 do {
1668                         si = *sx++;
1669                         ys = (si & 0xffff) * q + carry;
1670                         zs = (si >> 16) * q + (ys >> 16);
1671                         carry = zs >> 16;
1672                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1673                         borrow = y >> 16;
1674                         Sign_Extend(borrow, y);
1675                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
1676                         borrow = z >> 16;
1677                         Sign_Extend(borrow, z);
1678                         Storeinc(bx, z, y);
1679                         }
1680                         while(sx <= sxe);
1681                 if (!*bxe) {
1682                         bx = b->x;
1683                         while(--bxe > bx && !*bxe)
1684                                 --n;
1685                         b->wds = n;
1686                         }
1687                 }
1688         if (cmp(b, S) >= 0) {
1689                 q++;
1690                 borrow = 0;
1691                 carry = 0;
1692                 bx = b->x;
1693                 sx = S->x;
1694                 do {
1695                         si = *sx++;
1696                         ys = (si & 0xffff) + carry;
1697                         zs = (si >> 16) + (ys >> 16);
1698                         carry = zs >> 16;
1699                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1700                         borrow = y >> 16;
1701                         Sign_Extend(borrow, y);
1702                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
1703                         borrow = z >> 16;
1704                         Sign_Extend(borrow, z);
1705                         Storeinc(bx, z, y);
1706                         }
1707                         while(sx <= sxe);
1708                 bx = b->x;
1709                 bxe = bx + n;
1710                 if (!*bxe) {
1711                         while(--bxe > bx && !*bxe)
1712                                 --n;
1713                         b->wds = n;
1714                         }
1715                 }
1716         return q;
1717         }
1718
1719 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1720  *
1721  * Inspired by "How to Print Floating-Point Numbers Accurately" by
1722  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1723  *
1724  * Modifications:
1725  *      1. Rather than iterating, we use a simple numeric overestimate
1726  *         to determine k = floor(log10(d)).  We scale relevant
1727  *         quantities using O(log2(k)) rather than O(k) multiplications.
1728  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1729  *         try to generate digits strictly left to right.  Instead, we
1730  *         compute with fewer bits and propagate the carry if necessary
1731  *         when rounding the final digit up.  This is often faster.
1732  *      3. Under the assumption that input will be rounded nearest,
1733  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1734  *         That is, we allow equality in stopping tests when the
1735  *         round-nearest rule will give the same floating-point value
1736  *         as would satisfaction of the stopping test with strict
1737  *         inequality.
1738  *      4. We remove common factors of powers of 2 from relevant
1739  *         quantities.
1740  *      5. When converting floating-point integers less than 1e16,
1741  *         we use floating-point arithmetic rather than resorting
1742  *         to multiple-precision integers.
1743  *      6. When asked to produce fewer than 15 digits, we first try
1744  *         to get by with floating-point arithmetic; we resort to
1745  *         multiple-precision integer arithmetic only if we cannot
1746  *         guarantee that the floating-point calculation has given
1747  *         the correctly rounded result.  For k requested digits and
1748  *         "uniformly" distributed input, the probability is
1749  *         something like 10^(k-15) that we must resort to the long
1750  *         calculation.
1751  */
1752
1753  char *
1754 _IO_dtoa
1755 #ifdef KR_headers
1756         (d, mode, ndigits, decpt, sign, rve)
1757         double d; int mode, ndigits, *decpt, *sign; char **rve;
1758 #else
1759         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1760 #endif
1761 {
1762  /*     Arguments ndigits, decpt, sign are similar to those
1763         of ecvt and fcvt; trailing zeros are suppressed from
1764         the returned string.  If not null, *rve is set to point
1765         to the end of the return value.  If d is +-Infinity or NaN,
1766         then *decpt is set to 9999.
1767
1768         mode:
1769                 0 ==> shortest string that yields d when read in
1770                         and rounded to nearest.
1771                 1 ==> like 0, but with Steele & White stopping rule;
1772                         e.g. with IEEE P754 arithmetic , mode 0 gives
1773                         1e23 whereas mode 1 gives 9.999999999999999e22.
1774                 2 ==> max(1,ndigits) significant digits.  This gives a
1775                         return value similar to that of ecvt, except
1776                         that trailing zeros are suppressed.
1777                 3 ==> through ndigits past the decimal point.  This
1778                         gives a return value similar to that from fcvt,
1779                         except that trailing zeros are suppressed, and
1780                         ndigits can be negative.
1781                 4-9 should give the same return values as 2-3, i.e.,
1782                         4 <= mode <= 9 ==> same return as mode
1783                         2 + (mode & 1).  These modes are mainly for
1784                         debugging; often they run slower but sometimes
1785                         faster than modes 2-3.
1786                 4,5,8,9 ==> left-to-right digit generation.
1787                 6-9 ==> don't try fast floating-point estimate
1788                         (if applicable).
1789
1790                 Values of mode other than 0-9 are treated as mode 0.
1791
1792                 Sufficient space is allocated to the return value
1793                 to hold the suppressed trailing zeros.
1794         */
1795
1796         _G_int32_t bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1797                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1798                 spec_case, try_quick;
1799         _G_int32_t L;
1800 #ifndef Sudden_Underflow
1801         int denorm;
1802 #endif
1803         Bigint _b_avail, _b, _mhi, _mlo, _S;
1804         Bigint *b_avail = Binit(&_b_avail);
1805         Bigint *b = Binit(&_b);
1806         Bigint *S = Binit(&_S);
1807         /* mhi and mlo are only set and used if leftright. */
1808         Bigint *mhi = NULL, *mlo = NULL;
1809         double d2, ds, eps;
1810         char *s, *s0;
1811         static Bigint *result = NULL;
1812         static int result_k;
1813
1814         TEST_ENDIANNESS;
1815         if (result) {
1816                 /* result is contains a string, so its fields (interpreted
1817                    as a Bigint have been trashed.  Restore them.
1818                    This is a really ugly interface - result should
1819                    not be static, since that is not thread-safe.  FIXME. */
1820                 result->k = result_k;
1821                 result->maxwds = 1 << result_k;
1822                 result->on_stack = 0;
1823                 }
1824
1825         if (word0(d) & Sign_bit) {
1826                 /* set sign for everything, including 0's and NaNs */
1827                 *sign = 1;
1828                 setword0(d, word0(d) & ~Sign_bit);  /* clear sign bit */
1829                 }
1830         else
1831                 *sign = 0;
1832
1833 #if defined(IEEE_Arith) + defined(VAX)
1834 #ifdef IEEE_Arith
1835         if ((word0(d) & Exp_mask) == Exp_mask)
1836 #else
1837         if (word0(d)  == 0x8000)
1838 #endif
1839                 {
1840                 /* Infinity or NaN */
1841                 *decpt = 9999;
1842 #ifdef IEEE_Arith
1843                 if (!word1(d) && !(word0(d) & 0xfffff))
1844                   {
1845                     s = "Infinity";
1846                     if (rve)
1847                       *rve = s + 8;
1848                   }
1849                 else
1850 #endif
1851                   {
1852                     s = "NaN";
1853                     if (rve)
1854                       *rve = s +3;
1855                   }
1856                 return s;
1857                 }
1858 #endif
1859 #ifdef IBM
1860         d += 0; /* normalize */
1861 #endif
1862         if (!d) {
1863                 *decpt = 1;
1864                 s = "0";
1865                 if (rve)
1866                         *rve = s + 1;
1867                 return s;
1868                 }
1869
1870         b = d2b(b, d, &be, &bbits);
1871         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1872 #ifndef Sudden_Underflow
1873         if (i) {
1874 #endif
1875                 d2 = d;
1876                 setword0(d2, (word0(d2) & Frac_mask1) | Exp_11);
1877 #ifdef IBM
1878                 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
1879                         d2 /= 1 << j;
1880 #endif
1881
1882                 i -= Bias;
1883 #ifdef IBM
1884                 i <<= 2;
1885                 i += j;
1886 #endif
1887 #ifndef Sudden_Underflow
1888                 denorm = 0;
1889                 }
1890         else {
1891                 /* d is denormalized */
1892                 unsigned32 x;
1893
1894                 i = bbits + be + (Bias + (P-1) - 1);
1895                 x = i > 32  ? word0(d) << (64 - i) | word1(d) >> (i - 32)
1896                             : word1(d) << (32 - i);
1897                 d2 = x;
1898                 addword0(d2, - 31*Exp_msk1); /* adjust exponent */
1899                 i -= (Bias + (P-1) - 1) + 1;
1900                 denorm = 1;
1901                 }
1902 #endif
1903
1904         /* Now i is the unbiased base-2 exponent. */
1905
1906         /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
1907          * log10(x)      =  log(x) / log(10)
1908          *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1909          * log10(d) = i*log(2)/log(10) + log10(d2)
1910          *
1911          * This suggests computing an approximation k to log10(d) by
1912          *
1913          * k = i*0.301029995663981
1914          *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1915          *
1916          * We want k to be too large rather than too small.
1917          * The error in the first-order Taylor series approximation
1918          * is in our favor, so we just round up the constant enough
1919          * to compensate for any error in the multiplication of
1920          * (i) by 0.301029995663981; since |i| <= 1077,
1921          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1922          * adding 1e-13 to the constant term more than suffices.
1923          * Hence we adjust the constant term to 0.1760912590558.
1924          * (We could get a more accurate k by invoking log10,
1925          *  but this is probably not worthwhile.)
1926          */
1927
1928         ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1929         k = (int)ds;
1930         if (ds < 0. && ds != k)
1931                 k--;    /* want k = floor(ds) */
1932         k_check = 1;
1933         if (k >= 0 && k <= Ten_pmax) {
1934                 if (d < tens[k])
1935                         k--;
1936                 k_check = 0;
1937                 }
1938         j = bbits - i - 1;
1939         if (j >= 0) {
1940                 b2 = 0;
1941                 s2 = j;
1942                 }
1943         else {
1944                 b2 = -j;
1945                 s2 = 0;
1946                 }
1947         if (k >= 0) {
1948                 b5 = 0;
1949                 s5 = k;
1950                 s2 += k;
1951                 }
1952         else {
1953                 b2 -= k;
1954                 b5 = -k;
1955                 s5 = 0;
1956                 }
1957         if (mode < 0 || mode > 9)
1958                 mode = 0;
1959         try_quick = 1;
1960         if (mode > 5) {
1961                 mode -= 4;
1962                 try_quick = 0;
1963                 }
1964         leftright = 1;
1965         switch(mode) {
1966                 case 0:
1967                 case 1:
1968                         ilim = ilim1 = -1;
1969                         i = 18;
1970                         ndigits = 0;
1971                         break;
1972                 case 2:
1973                         leftright = 0;
1974                         /* no break */
1975                 case 4:
1976                         if (ndigits <= 0)
1977                                 ndigits = 1;
1978                         ilim = ilim1 = i = ndigits;
1979                         break;
1980                 case 3:
1981                         leftright = 0;
1982                         /* no break */
1983                 case 5:
1984                         i = ndigits + k + 1;
1985                         ilim = i;
1986                         ilim1 = i - 1;
1987                         if (i <= 0)
1988                                 i = 1;
1989                 }
1990         /* i is now an upper bound of the number of digits to generate. */
1991         j = sizeof(unsigned32) * (1<<BIGINT_MINIMUM_K);
1992         /* The test is <= so as to allow room for the final '\0'. */
1993         for(result_k = BIGINT_MINIMUM_K; BIGINT_HEADER_SIZE + j <= i;
1994                 j <<= 1) result_k++;
1995         if (!result || result_k > result->k)
1996         {
1997           Bfree (result);
1998           result = Balloc(result_k);
1999         }
2000         s = s0 = (char *)result;
2001
2002         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2003
2004                 /* Try to get by with floating-point arithmetic. */
2005
2006                 i = 0;
2007                 d2 = d;
2008                 k0 = k;
2009                 ilim0 = ilim;
2010                 ieps = 2; /* conservative */
2011                 if (k > 0) {
2012                         ds = tens[k&0xf];
2013                         j = k >> 4;
2014                         if (j & Bletch) {
2015                                 /* prevent overflows */
2016                                 j &= Bletch - 1;
2017                                 d /= bigtens[n_bigtens-1];
2018                                 ieps++;
2019                                 }
2020                         for(; j; j >>= 1, i++)
2021                                 if (j & 1) {
2022                                         ieps++;
2023                                         ds *= bigtens[i];
2024                                         }
2025                         d /= ds;
2026                         }
2027                 else if ((j1 = -k)) {
2028                         d *= tens[j1 & 0xf];
2029                         for(j = j1 >> 4; j; j >>= 1, i++)
2030                                 if (j & 1) {
2031                                         ieps++;
2032                                         d *= bigtens[i];
2033                                         }
2034                         }
2035                 if (k_check && d < 1. && ilim > 0) {
2036                         if (ilim1 <= 0)
2037                                 goto fast_failed;
2038                         ilim = ilim1;
2039                         k--;
2040                         d *= 10.;
2041                         ieps++;
2042                         }
2043                 eps = ieps*d + 7.;
2044                 addword0(eps, - (P-1)*Exp_msk1);
2045                 if (ilim == 0) {
2046                         d -= 5.;
2047                         if (d > eps)
2048                                 goto one_digit;
2049                         if (d < -eps)
2050                                 goto no_digits;
2051                         goto fast_failed;
2052                         }
2053 #ifndef No_leftright
2054                 if (leftright) {
2055                         /* Use Steele & White method of only
2056                          * generating digits needed.
2057                          */
2058                         eps = 0.5/tens[ilim-1] - eps;
2059                         for(i = 0;;) {
2060                                 L = (_G_int32_t)d;
2061                                 d -= L;
2062                                 *s++ = '0' + (int)L;
2063                                 if (d < eps)
2064                                         goto ret1;
2065                                 if (1. - d < eps)
2066                                         goto bump_up;
2067                                 if (++i >= ilim)
2068                                         break;
2069                                 eps *= 10.;
2070                                 d *= 10.;
2071                                 }
2072                         }
2073                 else {
2074 #endif
2075                         /* Generate ilim digits, then fix them up. */
2076                         eps *= tens[ilim-1];
2077                         for(i = 1;; i++, d *= 10.) {
2078                                 L = (_G_int32_t)d;
2079                                 d -= L;
2080                                 *s++ = '0' + (int)L;
2081                                 if (i == ilim) {
2082                                         if (d > 0.5 + eps)
2083                                                 goto bump_up;
2084                                         else if (d < 0.5 - eps) {
2085                                                 while(*--s == '0');
2086                                                 s++;
2087                                                 goto ret1;
2088                                                 }
2089                                         break;
2090                                         }
2091                                 }
2092 #ifndef No_leftright
2093                         }
2094 #endif
2095  fast_failed:
2096                 s = s0;
2097                 d = d2;
2098                 k = k0;
2099                 ilim = ilim0;
2100                 }
2101
2102         /* Do we have a "small" integer? */
2103
2104         if (be >= 0 && k <= Int_max) {
2105                 /* Yes. */
2106                 ds = tens[k];
2107                 if (ndigits < 0 && ilim <= 0) {
2108                         if (ilim < 0 || d <= 5*ds)
2109                                 goto no_digits;
2110                         goto one_digit;
2111                         }
2112                 for(i = 1;; i++) {
2113                         L = (_G_int32_t)(d / ds);
2114                         d -= L*ds;
2115 #ifdef Check_FLT_ROUNDS
2116                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2117                         if (d < 0) {
2118                                 L--;
2119                                 d += ds;
2120                                 }
2121 #endif
2122                         *s++ = '0' + (int)L;
2123                         if (i == ilim) {
2124                                 d += d;
2125                                 if (d > ds || (d == ds && L & 1)) {
2126  bump_up:
2127                                         while(*--s == '9')
2128                                                 if (s == s0) {
2129                                                         k++;
2130                                                         *s = '0';
2131                                                         break;
2132                                                         }
2133                                         ++*s++;
2134                                         }
2135                                 break;
2136                                 }
2137                         if (!(d *= 10.))
2138                                 break;
2139                         }
2140                 goto ret1;
2141                 }
2142
2143         m2 = b2;
2144         m5 = b5;
2145         if (leftright) {
2146                 if (mode < 2) {
2147                         i =
2148 #ifndef Sudden_Underflow
2149                                 denorm ? be + (Bias + (P-1) - 1 + 1) :
2150 #endif
2151 #ifdef IBM
2152                                 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2153 #else
2154                                 1 + P - bbits;
2155 #endif
2156                         }
2157                 else {
2158                         j = ilim - 1;
2159                         if (m5 >= j)
2160                                 m5 -= j;
2161                         else {
2162                                 s5 += j -= m5;
2163                                 b5 += j;
2164                                 m5 = 0;
2165                                 }
2166                         if ((i = ilim) < 0) {
2167                                 m2 -= i;
2168                                 i = 0;
2169                                 }
2170                         }
2171                 b2 += i;
2172                 s2 += i;
2173                 mhi = i2b(Binit(&_mhi), 1);
2174                 }
2175         if (m2 > 0 && s2 > 0) {
2176                 i = m2 < s2 ? m2 : s2;
2177                 b2 -= i;
2178                 m2 -= i;
2179                 s2 -= i;
2180                 }
2181         if (b5 > 0) {
2182                 if (leftright) {
2183                         if (m5 > 0) {
2184                                 Bigint *b_tmp;
2185                                 mhi = pow5mult(mhi, m5);
2186                                 b_tmp = mult(b_avail, mhi, b);
2187                                 b_avail = b;
2188                                 b = b_tmp;
2189                                 }
2190                         if ((j = b5 - m5))
2191                                 b = pow5mult(b, j);
2192                         }
2193                 else
2194                         b = pow5mult(b, b5);
2195                 }
2196         S = i2b(S, 1);
2197         if (s5 > 0)
2198                 S = pow5mult(S, s5);
2199
2200         /* Check for special case that d is a normalized power of 2. */
2201
2202         if (mode < 2) {
2203                 if (!word1(d) && !(word0(d) & Bndry_mask)
2204 #ifndef Sudden_Underflow
2205                  && word0(d) & Exp_mask
2206 #endif
2207                                 ) {
2208                         /* The special case */
2209                         b2 += Log2P;
2210                         s2 += Log2P;
2211                         spec_case = 1;
2212                         }
2213                 else
2214                         spec_case = 0;
2215                 }
2216
2217         /* Arrange for convenient computation of quotients:
2218          * shift left if necessary so divisor has 4 leading 0 bits.
2219          *
2220          * Perhaps we should just compute leading 28 bits of S once
2221          * and for all and pass them and a shift to quorem, so it
2222          * can do shifts and ors to compute the numerator for q.
2223          */
2224         if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
2225                 i = 32 - i;
2226         if (i > 4) {
2227                 i -= 4;
2228                 b2 += i;
2229                 m2 += i;
2230                 s2 += i;
2231                 }
2232         else if (i < 4) {
2233                 i += 28;
2234                 b2 += i;
2235                 m2 += i;
2236                 s2 += i;
2237                 }
2238         if (b2 > 0)
2239                 b = lshift(b, b2);
2240         if (s2 > 0)
2241                 S = lshift(S, s2);
2242         if (k_check) {
2243                 if (cmp(b,S) < 0) {
2244                         k--;
2245                         b = multadd(b, 10, 0);  /* we botched the k estimate */
2246                         if (leftright)
2247                                 mhi = multadd(mhi, 10, 0);
2248                         ilim = ilim1;
2249                         }
2250                 }
2251         if (ilim <= 0 && mode > 2) {
2252                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2253                         /* no digits, fcvt style */
2254  no_digits:
2255                         k = -1 - ndigits;
2256                         goto ret;
2257                         }
2258  one_digit:
2259                 *s++ = '1';
2260                 k++;
2261                 goto ret;
2262                 }
2263         if (leftright) {
2264                 if (m2 > 0)
2265                         mhi = lshift(mhi, m2);
2266
2267                 /* Compute mlo -- check for special case
2268                  * that d is a normalized power of 2.
2269                  */
2270
2271                 if (spec_case) {
2272                         mlo = Brealloc(Binit(&_mlo), mhi->k);
2273                         Bcopy(mlo, mhi);
2274                         mhi = lshift(mhi, Log2P);
2275                         }
2276                 else
2277                         mlo = mhi;
2278
2279                 for(i = 1;;i++) {
2280                         dig = quorem(b,S) + '0';
2281                         /* Do we yet have the shortest decimal string
2282                          * that will round to d?
2283                          */
2284                         j = cmp(b, mlo);
2285                         b_avail = diff(b_avail, S, mhi); /* b_avail = S - mi */
2286                         j1 = b_avail->sign ? 1 : cmp(b, b_avail);
2287 #ifndef ROUND_BIASED
2288                         if (j1 == 0 && !mode && !(word1(d) & 1)) {
2289                                 if (dig == '9')
2290                                         goto round_9_up;
2291                                 if (j > 0)
2292                                         dig++;
2293                                 *s++ = dig;
2294                                 goto ret;
2295                                 }
2296 #endif
2297                         if (j < 0 || (j == 0 && !mode
2298 #ifndef ROUND_BIASED
2299                                                         && !(word1(d) & 1)
2300 #endif
2301                                         )) {
2302                                 if (j1 > 0) {
2303                                         b = lshift(b, 1);
2304                                         j1 = cmp(b, S);
2305                                         if ((j1 > 0 || (j1 == 0 && dig & 1))
2306                                         && dig++ == '9')
2307                                                 goto round_9_up;
2308                                         }
2309                                 *s++ = dig;
2310                                 goto ret;
2311                                 }
2312                         if (j1 > 0) {
2313                                 if (dig == '9') { /* possible if i == 1 */
2314  round_9_up:
2315                                         *s++ = '9';
2316                                         goto roundoff;
2317                                         }
2318                                 *s++ = dig + 1;
2319                                 goto ret;
2320                                 }
2321                         *s++ = dig;
2322                         if (i == ilim)
2323                                 break;
2324                         b = multadd(b, 10, 0);
2325                         if (mlo == mhi)
2326                                 mlo = mhi = multadd(mhi, 10, 0);
2327                         else {
2328                                 mlo = multadd(mlo, 10, 0);
2329                                 mhi = multadd(mhi, 10, 0);
2330                                 }
2331                         }
2332                 }
2333         else
2334                 for(i = 1;; i++) {
2335                         *s++ = dig = quorem(b,S) + '0';
2336                         if (i >= ilim)
2337                                 break;
2338                         b = multadd(b, 10, 0);
2339                         }
2340
2341         /* Round off last digit */
2342
2343         b = lshift(b, 1);
2344         j = cmp(b, S);
2345         if (j > 0 || (j == 0 && dig & 1)) {
2346  roundoff:
2347                 while(*--s == '9')
2348                         if (s == s0) {
2349                                 k++;
2350                                 *s++ = '1';
2351                                 goto ret;
2352                                 }
2353                 ++*s++;
2354                 }
2355         else {
2356                 while(*--s == '0');
2357                 s++;
2358                 }
2359  ret:
2360         Bfree(b_avail);
2361         Bfree(S);
2362         if (mhi) {
2363                 if (mlo && mlo != mhi)
2364                         Bfree(mlo);
2365                 Bfree(mhi);
2366                 }
2367  ret1:
2368         Bfree(b);
2369         *s = 0;
2370         *decpt = k + 1;
2371         if (rve)
2372                 *rve = s;
2373         return s0;
2374         }
2375 #endif /* _IO_USE_DTOA */