Return to floatconv.c CVS log

Up to [Apple XNU] / GNUtools / libg++ / libio

Annotation of GNUtools/libg++/libio/floatconv.c, revision 1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.