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