dtoa.cpp

00001 /****************************************************************
00002  *
00003  * The author of this software is David M. Gay.
00004  *
00005  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
00006  *
00007  * Permission to use, copy, modify, and distribute this software for any
00008  * purpose without fee is hereby granted, provided that this entire notice
00009  * is included in all copies of any software which is or includes a copy
00010  * or modification of this software and in all copies of the supporting
00011  * documentation for such software.
00012  *
00013  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00014  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
00015  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00016  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00017  *
00018  ***************************************************************/
00019 
00020 /* Please send bug reports to
00021     David M. Gay
00022     Bell Laboratories, Room 2C-463
00023     600 Mountain Avenue
00024     Murray Hill, NJ 07974-0636
00025     U.S.A.
00026     dmg@bell-labs.com
00027  */
00028 
00029 /* On a machine with IEEE extended-precision registers, it is
00030  * necessary to specify double-precision (53-bit) rounding precision
00031  * before invoking strtod or dtoa.  If the machine uses (the equivalent
00032  * of) Intel 80x87 arithmetic, the call
00033  *  _control87(PC_53, MCW_PC);
00034  * does this with many compilers.  Whether this or another call is
00035  * appropriate depends on the compiler; for this to work, it may be
00036  * necessary to #include "float.h" or another system-dependent header
00037  * file.
00038  */
00039 
00040 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
00041  *
00042  * This strtod returns a nearest machine number to the input decimal
00043  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
00044  * broken by the IEEE round-even rule.  Otherwise ties are broken by
00045  * biased rounding (add half and chop).
00046  *
00047  * Inspired loosely by William D. Clinger's paper "How to Read Floating
00048  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
00049  *
00050  * Modifications:
00051  *
00052  *  1. We only require IEEE, IBM, or VAX double-precision
00053  *      arithmetic (not IEEE double-extended).
00054  *  2. We get by with floating-point arithmetic in a case that
00055  *      Clinger missed -- when we're computing d * 10^n
00056  *      for a small integer d and the integer n is not too
00057  *      much larger than 22 (the maximum integer k for which
00058  *      we can represent 10^k exactly), we may be able to
00059  *      compute (d*10^k) * 10^(e-k) with just one roundoff.
00060  *  3. Rather than a bit-at-a-time adjustment of the binary
00061  *      result in the hard case, we use floating-point
00062  *      arithmetic to determine the adjustment to within
00063  *      one bit; only in really hard cases do we need to
00064  *      compute a second residual.
00065  *  4. Because of 3., we don't need a large table of powers of 10
00066  *      for ten-to-e (just some small tables, e.g. of 10^k
00067  *      for 0 <= k <= 22).
00068  */
00069 
00070 /*
00071  * #define IEEE_8087 for IEEE-arithmetic machines where the least
00072  *  significant byte has the lowest address.
00073  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
00074  *  significant byte has the lowest address.
00075  * #define Long int on machines with 32-bit ints and 64-bit longs.
00076  * #define IBM for IBM mainframe-style floating-point arithmetic.
00077  * #define VAX for VAX-style floating-point arithmetic (D_floating).
00078  * #define No_leftright to omit left-right logic in fast floating-point
00079  *  computation of dtoa.
00080  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
00081  *  and strtod and dtoa should round accordingly.
00082  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
00083  *  and Honor_FLT_ROUNDS is not #defined.
00084  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
00085  *  that use extended-precision instructions to compute rounded
00086  *  products and quotients) with IBM.
00087  * #define ROUND_BIASED for IEEE-format with biased rounding.
00088  * #define Inaccurate_Divide for IEEE-format with correctly rounded
00089  *  products but inaccurate quotients, e.g., for Intel i860.
00090  * #define NO_LONG_LONG on machines that do not have a "long long"
00091  *  integer type (of >= 64 bits).  On such machines, you can
00092  *  #define Just_16 to store 16 bits per 32-bit Long when doing
00093  *  high-precision integer arithmetic.  Whether this speeds things
00094  *  up or slows things down depends on the machine and the number
00095  *  being converted.  If long long is available and the name is
00096  *  something other than "long long", #define Llong to be the name,
00097  *  and if "unsigned Llong" does not work as an unsigned version of
00098  *  Llong, #define #ULLong to be the corresponding unsigned type.
00099  * #define KR_headers for old-style C function headers.
00100  * #define Bad_float_h if your system lacks a float.h or if it does not
00101  *  define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
00102  *  FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
00103  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
00104  *  if memory is available and otherwise does something you deem
00105  *  appropriate.  If MALLOC is undefined, malloc will be invoked
00106  *  directly -- and assumed always to succeed.
00107  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
00108  *  memory allocations from a private pool of memory when possible.
00109  *  When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
00110  *  unless #defined to be a different length.  This default length
00111  *  suffices to get rid of MALLOC calls except for unusual cases,
00112  *  such as decimal-to-binary conversion of a very long string of
00113  *  digits.  The longest string dtoa can return is about 751 bytes
00114  *  long.  For conversions by strtod of strings of 800 digits and
00115  *  all dtoa conversions in single-threaded executions with 8-byte
00116  *  pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
00117  *  pointers, PRIVATE_MEM >= 7112 appears adequate.
00118  * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
00119  *  Infinity and NaN (case insensitively).  On some systems (e.g.,
00120  *  some HP systems), it may be necessary to #define NAN_WORD0
00121  *  appropriately -- to the most significant word of a quiet NaN.
00122  *  (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
00123  *  When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
00124  *  strtod also accepts (case insensitively) strings of the form
00125  *  NaN(x), where x is a string of hexadecimal digits and spaces;
00126  *  if there is only one string of hexadecimal digits, it is taken
00127  *  for the 52 fraction bits of the resulting NaN; if there are two
00128  *  or more strings of hex digits, the first is for the high 20 bits,
00129  *  the second and subsequent for the low 32 bits, with intervening
00130  *  white space ignored; but if this results in none of the 52
00131  *  fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
00132  *  and NAN_WORD1 are used instead.
00133  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
00134  *  multiple threads.  In this case, you must provide (or suitably
00135  *  #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
00136  *  by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
00137  *  in pow5mult, ensures lazy evaluation of only one copy of high
00138  *  powers of 5; omitting this lock would introduce a small
00139  *  probability of wasting memory, but would otherwise be harmless.)
00140  *  You must also invoke freedtoa(s) to free the value s returned by
00141  *  dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
00142  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
00143  *  avoids underflows on inputs whose result does not underflow.
00144  *  If you #define NO_IEEE_Scale on a machine that uses IEEE-format
00145  *  floating-point numbers and flushes underflows to zero rather
00146  *  than implementing gradual underflow, then you must also #define
00147  *  Sudden_Underflow.
00148  * #define YES_ALIAS to permit aliasing certain double values with
00149  *  arrays of ULongs.  This leads to slightly better code with
00150  *  some compilers and was always used prior to 19990916, but it
00151  *  is not strictly legal and can cause trouble with aggressively
00152  *  optimizing compilers (e.g., gcc 2.95.1 under -O2).
00153  * #define USE_LOCALE to use the current locale's decimal_point value.
00154  * #define SET_INEXACT if IEEE arithmetic is being used and extra
00155  *  computation should be done to set the inexact flag when the
00156  *  result is inexact and avoid setting inexact when the result
00157  *  is exact.  In this case, dtoa.c must be compiled in
00158  *  an environment, perhaps provided by #include "dtoa.c" in a
00159  *  suitable wrapper, that defines two functions,
00160  *      int get_inexact(void);
00161  *      void clear_inexact(void);
00162  *  such that get_inexact() returns a nonzero value if the
00163  *  inexact bit is already set, and clear_inexact() sets the
00164  *  inexact bit to 0.  When SET_INEXACT is #defined, strtod
00165  *  also does extra computations to set the underflow and overflow
00166  *  flags when appropriate (i.e., when the result is tiny and
00167  *  inexact or when it is a numeric value rounded to +-infinity).
00168  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
00169  *  the result overflows to +-Infinity or underflows to 0.
00170  */
00171 
00172 // Put this before anything else that may import a definition of CONST. CONST from grammar.cpp conflicts with this.
00173 #ifdef KDE_USE_FINAL
00174 #undef CONST
00175 #endif
00176 
00177 #include "dtoa.h"
00178 #include <config.h>
00179 
00180 #include "global.h"
00181 
00182 #ifdef WORDS_BIGENDIAN
00183 #define IEEE_MC68k
00184 #else
00185 #define IEEE_8087
00186 #endif
00187 #define INFNAN_CHECK
00188 
00189 
00190 
00191 #ifndef Long
00192 #define Long int
00193 #endif
00194 #ifndef ULong
00195 typedef unsigned Long ULong;
00196 #endif
00197 
00198 #ifdef DEBUG
00199 #include <stdio.h>
00200 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
00201 #endif
00202 
00203 #include <stdlib.h>
00204 #include <string.h>
00205 
00206 #ifdef USE_LOCALE
00207 #include <locale.h>
00208 #endif
00209 
00210 #ifdef MALLOC
00211 extern void *MALLOC(size_t);
00212 #else
00213 #define MALLOC malloc
00214 #endif
00215 
00216 #ifndef Omit_Private_Memory
00217 #ifndef PRIVATE_MEM
00218 #define PRIVATE_MEM 2304
00219 #endif
00220 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
00221 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
00222 #endif
00223 
00224 #undef IEEE_Arith
00225 #undef Avoid_Underflow
00226 #ifdef IEEE_MC68k
00227 #define IEEE_Arith
00228 #endif
00229 #ifdef IEEE_8087
00230 #define IEEE_Arith
00231 #endif
00232 
00233 #include <errno.h>
00234 
00235 #ifdef Bad_float_h
00236 
00237 #ifdef IEEE_Arith
00238 #define DBL_DIG 15
00239 #define DBL_MAX_10_EXP 308
00240 #define DBL_MAX_EXP 1024
00241 #define FLT_RADIX 2
00242 #endif /*IEEE_Arith*/
00243 
00244 #ifdef IBM
00245 #define DBL_DIG 16
00246 #define DBL_MAX_10_EXP 75
00247 #define DBL_MAX_EXP 63
00248 #define FLT_RADIX 16
00249 #define DBL_MAX 7.2370055773322621e+75
00250 #endif
00251 
00252 #ifdef VAX
00253 #define DBL_DIG 16
00254 #define DBL_MAX_10_EXP 38
00255 #define DBL_MAX_EXP 127
00256 #define FLT_RADIX 2
00257 #define DBL_MAX 1.7014118346046923e+38
00258 #endif
00259 
00260 #ifndef LONG_MAX
00261 #define LONG_MAX 2147483647
00262 #endif
00263 
00264 #else /* ifndef Bad_float_h */
00265 #include <float.h>
00266 #endif /* Bad_float_h */
00267 
00268 #ifndef __MATH_H__
00269 #include <math.h>
00270 #endif
00271 
00272 #define strtod kjs_strtod
00273 #define dtoa kjs_dtoa
00274 #define freedtoa kjs_freedtoa
00275 
00276 #ifdef __cplusplus
00277 extern "C" {
00278 #endif
00279 
00280 #ifndef CONST
00281 #define CONST const
00282 #endif
00283 
00284 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
00285 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
00286 #endif
00287 
00288 typedef union { double d; ULong L[2]; } U;
00289 
00290 #define dval(x) (x).d
00291 #ifdef IEEE_8087
00292 #define word0(x) (x).L[1]
00293 #define word1(x) (x).L[0]
00294 #else
00295 #define word0(x) (x).L[0]
00296 #define word1(x) (x).L[1]
00297 #endif
00298 
00299 /* The following definition of Storeinc is appropriate for MIPS processors.
00300  * An alternative that might be better on some machines is
00301  */
00302 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
00303 
00304 /* #define P DBL_MANT_DIG */
00305 /* Ten_pmax = floor(P*log(2)/log(5)) */
00306 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
00307 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
00308 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
00309 
00310 #ifdef IEEE_Arith
00311 #define Exp_shift  20
00312 #define Exp_shift1 20
00313 #define Exp_msk1    0x100000
00314 #define Exp_msk11   0x100000
00315 #define Exp_mask  0x7ff00000
00316 #define P 53
00317 #define Bias 1023
00318 #define Emin (-1022)
00319 #define Exp_1  0x3ff00000
00320 #define Exp_11 0x3ff00000
00321 #define Ebits 11
00322 #define Frac_mask  0xfffff
00323 #define Frac_mask1 0xfffff
00324 #define Ten_pmax 22
00325 #define Bletch 0x10
00326 #define Bndry_mask  0xfffff
00327 #define Bndry_mask1 0xfffff
00328 #define LSB 1
00329 #define Sign_bit 0x80000000
00330 #define Log2P 1
00331 #define Tiny0 0
00332 #define Tiny1 1
00333 #define Quick_max 14
00334 #define Int_max 14
00335 #ifndef NO_IEEE_Scale
00336 #define Avoid_Underflow
00337 #ifdef Flush_Denorm /* debugging option */
00338 #undef Sudden_Underflow
00339 #endif
00340 #endif
00341 
00342 #ifndef Flt_Rounds
00343 #ifdef FLT_ROUNDS
00344 #define Flt_Rounds FLT_ROUNDS
00345 #else
00346 #define Flt_Rounds 1
00347 #endif
00348 #endif /*Flt_Rounds*/
00349 
00350 #ifdef Honor_FLT_ROUNDS
00351 #define Rounding rounding
00352 #undef Check_FLT_ROUNDS
00353 #define Check_FLT_ROUNDS
00354 #else
00355 #define Rounding Flt_Rounds
00356 #endif
00357 
00358 #else /* ifndef IEEE_Arith */
00359 #undef Check_FLT_ROUNDS
00360 #undef Honor_FLT_ROUNDS
00361 #undef SET_INEXACT
00362 #undef  Sudden_Underflow
00363 #define Sudden_Underflow
00364 #ifdef IBM
00365 #undef Flt_Rounds
00366 #define Flt_Rounds 0
00367 #define Exp_shift  24
00368 #define Exp_shift1 24
00369 #define Exp_msk1   0x1000000
00370 #define Exp_msk11  0x1000000
00371 #define Exp_mask  0x7f000000
00372 #define P 14
00373 #define Bias 65
00374 #define Exp_1  0x41000000
00375 #define Exp_11 0x41000000
00376 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
00377 #define Frac_mask  0xffffff
00378 #define Frac_mask1 0xffffff
00379 #define Bletch 4
00380 #define Ten_pmax 22
00381 #define Bndry_mask  0xefffff
00382 #define Bndry_mask1 0xffffff
00383 #define LSB 1
00384 #define Sign_bit 0x80000000
00385 #define Log2P 4
00386 #define Tiny0 0x100000
00387 #define Tiny1 0
00388 #define Quick_max 14
00389 #define Int_max 15
00390 #else /* VAX */
00391 #undef Flt_Rounds
00392 #define Flt_Rounds 1
00393 #define Exp_shift  23
00394 #define Exp_shift1 7
00395 #define Exp_msk1    0x80
00396 #define Exp_msk11   0x800000
00397 #define Exp_mask  0x7f80
00398 #define P 56
00399 #define Bias 129
00400 #define Exp_1  0x40800000
00401 #define Exp_11 0x4080
00402 #define Ebits 8
00403 #define Frac_mask  0x7fffff
00404 #define Frac_mask1 0xffff007f
00405 #define Ten_pmax 24
00406 #define Bletch 2
00407 #define Bndry_mask  0xffff007f
00408 #define Bndry_mask1 0xffff007f
00409 #define LSB 0x10000
00410 #define Sign_bit 0x8000
00411 #define Log2P 1
00412 #define Tiny0 0x80
00413 #define Tiny1 0
00414 #define Quick_max 15
00415 #define Int_max 15
00416 #endif /* IBM, VAX */
00417 #endif /* IEEE_Arith */
00418 
00419 #ifndef IEEE_Arith
00420 #define ROUND_BIASED
00421 #endif
00422 
00423 #ifdef RND_PRODQUOT
00424 #define rounded_product(a,b) a = rnd_prod(a, b)
00425 #define rounded_quotient(a,b) a = rnd_quot(a, b)
00426 extern double rnd_prod(double, double), rnd_quot(double, double);
00427 #else
00428 #define rounded_product(a,b) a *= b
00429 #define rounded_quotient(a,b) a /= b
00430 #endif
00431 
00432 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
00433 #define Big1 0xffffffff
00434 
00435 #ifndef Pack_32
00436 #define Pack_32
00437 #endif
00438 
00439 #define FFFFFFFF 0xffffffffUL
00440 
00441 #ifdef NO_LONG_LONG
00442 #undef ULLong
00443 #ifdef Just_16
00444 #undef Pack_32
00445 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
00446  * This makes some inner loops simpler and sometimes saves work
00447  * during multiplications, but it often seems to make things slightly
00448  * slower.  Hence the default is now to store 32 bits per Long.
00449  */
00450 #endif
00451 #else   /* long long available */
00452 #ifndef Llong
00453 #define Llong long long
00454 #endif
00455 #ifndef ULLong
00456 #define ULLong unsigned Llong
00457 #endif
00458 #endif /* NO_LONG_LONG */
00459 
00460 #ifndef MULTIPLE_THREADS
00461 #define ACQUIRE_DTOA_LOCK(n)    /*nothing*/
00462 #define FREE_DTOA_LOCK(n)   /*nothing*/
00463 #endif
00464 
00465 #define Kmax (sizeof(size_t) << 3)
00466 
00467  struct
00468 Bigint {
00469     struct Bigint *next;
00470     int k, maxwds, sign, wds;
00471     ULong x[1];
00472     };
00473 
00474  typedef struct Bigint Bigint;
00475 
00476  static Bigint *freelist[Kmax+1];
00477 
00478  static Bigint *
00479 Balloc
00480     (int k)
00481 {
00482     int x;
00483     Bigint *rv;
00484 #ifndef Omit_Private_Memory
00485     unsigned int len;
00486 #endif
00487 
00488     ACQUIRE_DTOA_LOCK(0);
00489     if ((rv = freelist[k])) {
00490         freelist[k] = rv->next;
00491         }
00492     else {
00493         x = 1 << k;
00494 #ifdef Omit_Private_Memory
00495         rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
00496 #else
00497         len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
00498             /sizeof(double);
00499         if (pmem_next - private_mem + len <= (unsigned)PRIVATE_mem) {
00500             rv = (Bigint*)pmem_next;
00501             pmem_next += len;
00502             }
00503         else
00504             rv = (Bigint*)MALLOC(len*sizeof(double));
00505 #endif
00506         rv->k = k;
00507         rv->maxwds = x;
00508         }
00509     FREE_DTOA_LOCK(0);
00510     rv->sign = rv->wds = 0;
00511     return rv;
00512     }
00513 
00514  static void
00515 Bfree
00516     (Bigint *v)
00517 {
00518     if (v) {
00519         ACQUIRE_DTOA_LOCK(0);
00520         v->next = freelist[v->k];
00521         freelist[v->k] = v;
00522         FREE_DTOA_LOCK(0);
00523         }
00524     }
00525 
00526 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
00527 y->wds*sizeof(Long) + 2*sizeof(int))
00528 
00529  static Bigint *
00530 multadd
00531     (Bigint *b, int m, int a)   /* multiply by m and add a */
00532 {
00533     int i, wds;
00534 #ifdef ULLong
00535     ULong *x;
00536     ULLong carry, y;
00537 #else
00538     ULong carry, *x, y;
00539 #ifdef Pack_32
00540     ULong xi, z;
00541 #endif
00542 #endif
00543     Bigint *b1;
00544 
00545     wds = b->wds;
00546     x = b->x;
00547     i = 0;
00548     carry = a;
00549     do {
00550 #ifdef ULLong
00551         y = *x * (ULLong)m + carry;
00552         carry = y >> 32;
00553         *x++ = (ULong)y & FFFFFFFF;
00554 #else
00555 #ifdef Pack_32
00556         xi = *x;
00557         y = (xi & 0xffff) * m + carry;
00558         z = (xi >> 16) * m + (y >> 16);
00559         carry = z >> 16;
00560         *x++ = (z << 16) + (y & 0xffff);
00561 #else
00562         y = *x * m + carry;
00563         carry = y >> 16;
00564         *x++ = y & 0xffff;
00565 #endif
00566 #endif
00567         }
00568         while(++i < wds);
00569     if (carry) {
00570         if (wds >= b->maxwds) {
00571             b1 = Balloc(b->k+1);
00572             Bcopy(b1, b);
00573             Bfree(b);
00574             b = b1;
00575             }
00576         b->x[wds++] = (ULong)carry;
00577         b->wds = wds;
00578         }
00579     return b;
00580     }
00581 
00582  static Bigint *
00583 s2b
00584     (CONST char *s, int nd0, int nd, ULong y9)
00585 {
00586     Bigint *b;
00587     int i, k;
00588     Long x, y;
00589 
00590     x = (nd + 8) / 9;
00591     for(k = 0, y = 1; x > y; y <<= 1, k++) ;
00592 #ifdef Pack_32
00593     b = Balloc(k);
00594     b->x[0] = y9;
00595     b->wds = 1;
00596 #else
00597     b = Balloc(k+1);
00598     b->x[0] = y9 & 0xffff;
00599     b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
00600 #endif
00601 
00602     i = 9;
00603     if (9 < nd0) {
00604         s += 9;
00605         do b = multadd(b, 10, *s++ - '0');
00606             while(++i < nd0);
00607         s++;
00608         }
00609     else
00610         s += 10;
00611     for(; i < nd; i++)
00612         b = multadd(b, 10, *s++ - '0');
00613     return b;
00614     }
00615 
00616  static int
00617 hi0bits
00618     (register ULong x)
00619 {
00620     register int k = 0;
00621 
00622     if (!(x & 0xffff0000)) {
00623         k = 16;
00624         x <<= 16;
00625         }
00626     if (!(x & 0xff000000)) {
00627         k += 8;
00628         x <<= 8;
00629         }
00630     if (!(x & 0xf0000000)) {
00631         k += 4;
00632         x <<= 4;
00633         }
00634     if (!(x & 0xc0000000)) {
00635         k += 2;
00636         x <<= 2;
00637         }
00638     if (!(x & 0x80000000)) {
00639         k++;
00640         if (!(x & 0x40000000))
00641             return 32;
00642         }
00643     return k;
00644     }
00645 
00646  static int
00647 lo0bits
00648     (ULong *y)
00649 {
00650     register int k;
00651     register ULong x = *y;
00652 
00653     if (x & 7) {
00654         if (x & 1)
00655             return 0;
00656         if (x & 2) {
00657             *y = x >> 1;
00658             return 1;
00659             }
00660         *y = x >> 2;
00661         return 2;
00662         }
00663     k = 0;
00664     if (!(x & 0xffff)) {
00665         k = 16;
00666         x >>= 16;
00667         }
00668     if (!(x & 0xff)) {
00669         k += 8;
00670         x >>= 8;
00671         }
00672     if (!(x & 0xf)) {
00673         k += 4;
00674         x >>= 4;
00675         }
00676     if (!(x & 0x3)) {
00677         k += 2;
00678         x >>= 2;
00679         }
00680     if (!(x & 1)) {
00681         k++;
00682         x >>= 1;
00683         if (!x & 1)
00684             return 32;
00685         }
00686     *y = x;
00687     return k;
00688     }
00689 
00690  static Bigint *
00691 i2b
00692     (int i)
00693 {
00694     Bigint *b;
00695 
00696     b = Balloc(1);
00697     b->x[0] = i;
00698     b->wds = 1;
00699     return b;
00700     }
00701 
00702  static Bigint *
00703 mult
00704     (Bigint *a, Bigint *b)
00705 {
00706     Bigint *c;
00707     int k, wa, wb, wc;
00708     ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00709     ULong y;
00710 #ifdef ULLong
00711     ULLong carry, z;
00712 #else
00713     ULong carry, z;
00714 #ifdef Pack_32
00715     ULong z2;
00716 #endif
00717 #endif
00718 
00719     if (a->wds < b->wds) {
00720         c = a;
00721         a = b;
00722         b = c;
00723         }
00724     k = a->k;
00725     wa = a->wds;
00726     wb = b->wds;
00727     wc = wa + wb;
00728     if (wc > a->maxwds)
00729         k++;
00730     c = Balloc(k);
00731     for(x = c->x, xa = x + wc; x < xa; x++)
00732         *x = 0;
00733     xa = a->x;
00734     xae = xa + wa;
00735     xb = b->x;
00736     xbe = xb + wb;
00737     xc0 = c->x;
00738 #ifdef ULLong
00739     for(; xb < xbe; xc0++) {
00740         if ((y = *xb++)) {
00741             x = xa;
00742             xc = xc0;
00743             carry = 0;
00744             do {
00745                 z = *x++ * (ULLong)y + *xc + carry;
00746                 carry = z >> 32;
00747                 *xc++ = (ULong)z & FFFFFFFF;
00748                 }
00749                 while(x < xae);
00750             *xc = (ULong)carry;
00751             }
00752         }
00753 #else
00754 #ifdef Pack_32
00755     for(; xb < xbe; xb++, xc0++) {
00756         if (y = *xb & 0xffff) {
00757             x = xa;
00758             xc = xc0;
00759             carry = 0;
00760             do {
00761                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
00762                 carry = z >> 16;
00763                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
00764                 carry = z2 >> 16;
00765                 Storeinc(xc, z2, z);
00766                 }
00767                 while(x < xae);
00768             *xc = carry;
00769             }
00770         if (y = *xb >> 16) {
00771             x = xa;
00772             xc = xc0;
00773             carry = 0;
00774             z2 = *xc;
00775             do {
00776                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
00777                 carry = z >> 16;
00778                 Storeinc(xc, z, z2);
00779                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
00780                 carry = z2 >> 16;
00781                 }
00782                 while(x < xae);
00783             *xc = z2;
00784             }
00785         }
00786 #else
00787     for(; xb < xbe; xc0++) {
00788         if (y = *xb++) {
00789             x = xa;
00790             xc = xc0;
00791             carry = 0;
00792             do {
00793                 z = *x++ * y + *xc + carry;
00794                 carry = z >> 16;
00795                 *xc++ = z & 0xffff;
00796                 }
00797                 while(x < xae);
00798             *xc = carry;
00799             }
00800         }
00801 #endif
00802 #endif
00803     for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
00804     c->wds = wc;
00805     return c;
00806     }
00807 
00808  static Bigint *p5s;
00809 
00810  static Bigint *
00811 pow5mult
00812     (Bigint *b, int k)
00813 {
00814     Bigint *b1, *p5, *p51;
00815     int i;
00816     static int p05[3] = { 5, 25, 125 };
00817 
00818     if ((i = k & 3))
00819         b = multadd(b, p05[i-1], 0);
00820 
00821     if (!(k >>= 2))
00822         return b;
00823     if (!(p5 = p5s)) {
00824         /* first time */
00825 #ifdef MULTIPLE_THREADS
00826         ACQUIRE_DTOA_LOCK(1);
00827         if (!(p5 = p5s)) {
00828             p5 = p5s = i2b(625);
00829             p5->next = 0;
00830             }
00831         FREE_DTOA_LOCK(1);
00832 #else
00833         p5 = p5s = i2b(625);
00834         p5->next = 0;
00835 #endif
00836         }
00837     for(;;) {
00838         if (k & 1) {
00839             b1 = mult(b, p5);
00840             Bfree(b);
00841             b = b1;
00842             }
00843         if (!(k >>= 1))
00844             break;
00845         if (!(p51 = p5->next)) {
00846 #ifdef MULTIPLE_THREADS
00847             ACQUIRE_DTOA_LOCK(1);
00848             if (!(p51 = p5->next)) {
00849                 p51 = p5->next = mult(p5,p5);
00850                 p51->next = 0;
00851                 }
00852             FREE_DTOA_LOCK(1);
00853 #else
00854             p51 = p5->next = mult(p5,p5);
00855             p51->next = 0;
00856 #endif
00857             }
00858         p5 = p51;
00859         }
00860     return b;
00861     }
00862 
00863  static Bigint *
00864 lshift
00865     (Bigint *b, int k)
00866 {
00867     int i, k1, n, n1;
00868     Bigint *b1;
00869     ULong *x, *x1, *xe, z;
00870 
00871 #ifdef Pack_32
00872     n = k >> 5;
00873 #else
00874     n = k >> 4;
00875 #endif
00876     k1 = b->k;
00877     n1 = n + b->wds + 1;
00878     for(i = b->maxwds; n1 > i; i <<= 1)
00879         k1++;
00880     b1 = Balloc(k1);
00881     x1 = b1->x;
00882     for(i = 0; i < n; i++)
00883         *x1++ = 0;
00884     x = b->x;
00885     xe = x + b->wds;
00886 #ifdef Pack_32
00887     if (k &= 0x1f) {
00888         k1 = 32 - k;
00889         z = 0;
00890         do {
00891             *x1++ = *x << k | z;
00892             z = *x++ >> k1;
00893             }
00894             while(x < xe);
00895         if ((*x1 = z))
00896             ++n1;
00897         }
00898 #else
00899     if (k &= 0xf) {
00900         k1 = 16 - k;
00901         z = 0;
00902         do {
00903             *x1++ = *x << k  & 0xffff | z;
00904             z = *x++ >> k1;
00905             }
00906             while(x < xe);
00907         if (*x1 = z)
00908             ++n1;
00909         }
00910 #endif
00911     else do
00912         *x1++ = *x++;
00913         while(x < xe);
00914     b1->wds = n1 - 1;
00915     Bfree(b);
00916     return b1;
00917     }
00918 
00919  static int
00920 cmp
00921     (Bigint *a, Bigint *b)
00922 {
00923     ULong *xa, *xa0, *xb, *xb0;
00924     int i, j;
00925 
00926     i = a->wds;
00927     j = b->wds;
00928 #ifdef DEBUG
00929     if (i > 1 && !a->x[i-1])
00930         Bug("cmp called with a->x[a->wds-1] == 0");
00931     if (j > 1 && !b->x[j-1])
00932         Bug("cmp called with b->x[b->wds-1] == 0");
00933 #endif
00934     if (i -= j)
00935         return i;
00936     xa0 = a->x;
00937     xa = xa0 + j;
00938     xb0 = b->x;
00939     xb = xb0 + j;
00940     for(;;) {
00941         if (*--xa != *--xb)
00942             return *xa < *xb ? -1 : 1;
00943         if (xa <= xa0)
00944             break;
00945         }
00946     return 0;
00947     }
00948 
00949  static Bigint *
00950 diff
00951     (Bigint *a, Bigint *b)
00952 {
00953     Bigint *c;
00954     int i, wa, wb;
00955     ULong *xa, *xae, *xb, *xbe, *xc;
00956 #ifdef ULLong
00957     ULLong borrow, y;
00958 #else
00959     ULong borrow, y;
00960 #ifdef Pack_32
00961     ULong z;
00962 #endif
00963 #endif
00964 
00965     i = cmp(a,b);
00966     if (!i) {
00967         c = Balloc(0);
00968         c->wds = 1;
00969         c->x[0] = 0;
00970         return c;
00971         }
00972     if (i < 0) {
00973         c = a;
00974         a = b;
00975         b = c;
00976         i = 1;
00977         }
00978     else
00979         i = 0;
00980     c = Balloc(a->k);
00981     c->sign = i;
00982     wa = a->wds;
00983     xa = a->x;
00984     xae = xa + wa;
00985     wb = b->wds;
00986     xb = b->x;
00987     xbe = xb + wb;
00988     xc = c->x;
00989     borrow = 0;
00990 #ifdef ULLong
00991     do {
00992         y = (ULLong)*xa++ - *xb++ - borrow;
00993         borrow = y >> 32 & (ULong)1;
00994         *xc++ = (ULong)y & FFFFFFFF;
00995         }
00996         while(xb < xbe);
00997     while(xa < xae) {
00998         y = *xa++ - borrow;
00999         borrow = y >> 32 & (ULong)1;
01000         *xc++ = (ULong)y & FFFFFFFF;
01001         }
01002 #else
01003 #ifdef Pack_32
01004     do {
01005         y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
01006         borrow = (y & 0x10000) >> 16;
01007         z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
01008         borrow = (z & 0x10000) >> 16;
01009         Storeinc(xc, z, y);
01010         }
01011         while(xb < xbe);
01012     while(xa < xae) {
01013         y = (*xa & 0xffff) - borrow;
01014         borrow = (y & 0x10000) >> 16;
01015         z = (*xa++ >> 16) - borrow;
01016         borrow = (z & 0x10000) >> 16;
01017         Storeinc(xc, z, y);
01018         }
01019 #else
01020     do {
01021         y = *xa++ - *xb++ - borrow;
01022         borrow = (y & 0x10000) >> 16;
01023         *xc++ = y & 0xffff;
01024         }
01025         while(xb < xbe);
01026     while(xa < xae) {
01027         y = *xa++ - borrow;
01028         borrow = (y & 0x10000) >> 16;
01029         *xc++ = y & 0xffff;
01030         }
01031 #endif
01032 #endif
01033     while(!*--xc)
01034         wa--;
01035     c->wds = wa;
01036     return c;
01037     }
01038 
01039  static double
01040 ulp
01041     (double dx)
01042 {
01043     register Long L;
01044     U x, a;
01045 
01046     dval(x) = dx;
01047     L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
01048 #ifndef Avoid_Underflow
01049 #ifndef Sudden_Underflow
01050     if (L > 0) {
01051 #endif
01052 #endif
01053 #ifdef IBM
01054         L |= Exp_msk1 >> 4;
01055 #endif
01056         word0(a) = L;
01057         word1(a) = 0;
01058 #ifndef Avoid_Underflow
01059 #ifndef Sudden_Underflow
01060         }
01061     else {
01062         L = -L >> Exp_shift;
01063         if (L < Exp_shift) {
01064             word0(a) = 0x80000 >> L;
01065             word1(a) = 0;
01066             }
01067         else {
01068             word0(a) = 0;
01069             L -= Exp_shift;
01070             word1(a) = L >= 31 ? 1 : 1 << 31 - L;
01071             }
01072         }
01073 #endif
01074 #endif
01075     return dval(a);
01076     }
01077 
01078  static double
01079 b2d
01080     (Bigint *a, int *e)
01081 {
01082     ULong *xa, *xa0, w, y, z;
01083     int k;
01084     U d;
01085 #ifdef VAX
01086     ULong d0, d1;
01087 #else
01088 #define d0 word0(d)
01089 #define d1 word1(d)
01090 #endif
01091 
01092     xa0 = a->x;
01093     xa = xa0 + a->wds;
01094     y = *--xa;
01095 #ifdef DEBUG
01096     if (!y) Bug("zero y in b2d");
01097 #endif
01098     k = hi0bits(y);
01099     *e = 32 - k;
01100 #ifdef Pack_32
01101     if (k < Ebits) {
01102         d0 = Exp_1 | y >> Ebits - k;
01103         w = xa > xa0 ? *--xa : 0;
01104         d1 = y << (32-Ebits) + k | w >> Ebits - k;
01105         goto ret_d;
01106         }
01107     z = xa > xa0 ? *--xa : 0;
01108     if (k -= Ebits) {
01109         d0 = Exp_1 | y << k | z >> 32 - k;
01110         y = xa > xa0 ? *--xa : 0;
01111         d1 = z << k | y >> 32 - k;
01112         }
01113     else {
01114         d0 = Exp_1 | y;
01115         d1 = z;
01116         }
01117 #else
01118     if (k < Ebits + 16) {
01119         z = xa > xa0 ? *--xa : 0;
01120         d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
01121         w = xa > xa0 ? *--xa : 0;
01122         y = xa > xa0 ? *--xa : 0;
01123         d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
01124         goto ret_d;
01125         }
01126     z = xa > xa0 ? *--xa : 0;
01127     w = xa > xa0 ? *--xa : 0;
01128     k -= Ebits + 16;
01129     d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
01130     y = xa > xa0 ? *--xa : 0;
01131     d1 = w << k + 16 | y << k;
01132 #endif
01133  ret_d:
01134 #ifdef VAX
01135     word0(d) = d0 >> 16 | d0 << 16;
01136     word1(d) = d1 >> 16 | d1 << 16;
01137 #else
01138 #undef d0
01139 #undef d1
01140 #endif
01141     return dval(d);
01142     }
01143 
01144  static Bigint *
01145 d2b
01146     (double dd, int *e, int *bits)
01147 {
01148     U d;
01149     Bigint *b;
01150     int de, k;
01151     ULong *x, y, z;
01152 #ifndef Sudden_Underflow
01153     int i;
01154 #endif
01155 #ifdef VAX
01156     ULong d0, d1;
01157 #endif
01158     dval(d) = dd;
01159 #ifdef VAX
01160     d0 = word0(d) >> 16 | word0(d) << 16;
01161     d1 = word1(d) >> 16 | word1(d) << 16;
01162 #else
01163 #define d0 word0(d)
01164 #define d1 word1(d)
01165 #endif
01166 
01167 #ifdef Pack_32
01168     b = Balloc(1);
01169 #else
01170     b = Balloc(2);
01171 #endif
01172     x = b->x;
01173 
01174     z = d0 & Frac_mask;
01175     d0 &= 0x7fffffff;   /* clear sign bit, which we ignore */
01176 #ifdef Sudden_Underflow
01177     de = (int)(d0 >> Exp_shift);
01178 #ifndef IBM
01179     z |= Exp_msk11;
01180 #endif
01181 #else
01182     if ((de = (int)(d0 >> Exp_shift)))
01183         z |= Exp_msk1;
01184 #endif
01185 #ifdef Pack_32
01186     if ((y = d1)) {
01187         if ((k = lo0bits(&y))) {
01188             x[0] = y | z << 32 - k;
01189             z >>= k;
01190             }
01191         else
01192             x[0] = y;
01193 #ifndef Sudden_Underflow
01194         i =
01195 #endif
01196             b->wds = (x[1] = z) ? 2 : 1;
01197         }
01198     else {
01199 #ifdef DEBUG
01200         if (!z)
01201             Bug("Zero passed to d2b");
01202 #endif
01203         k = lo0bits(&z);
01204         x[0] = z;
01205 #ifndef Sudden_Underflow
01206         i =
01207 #endif
01208             b->wds = 1;
01209         k += 32;
01210         }
01211 #else
01212     if (y = d1) {
01213         if (k = lo0bits(&y))
01214             if (k >= 16) {
01215                 x[0] = y | z << 32 - k & 0xffff;
01216                 x[1] = z >> k - 16 & 0xffff;
01217                 x[2] = z >> k;
01218                 i = 2;
01219                 }
01220             else {
01221                 x[0] = y & 0xffff;
01222                 x[1] = y >> 16 | z << 16 - k & 0xffff;
01223                 x[2] = z >> k & 0xffff;
01224                 x[3] = z >> k+16;
01225                 i = 3;
01226                 }
01227         else {
01228             x[0] = y & 0xffff;
01229             x[1] = y >> 16;
01230             x[2] = z & 0xffff;
01231             x[3] = z >> 16;
01232             i = 3;
01233             }
01234         }
01235     else {
01236 #ifdef DEBUG
01237         if (!z)
01238             Bug("Zero passed to d2b");
01239 #endif
01240         k = lo0bits(&z);
01241         if (k >= 16) {
01242             x[0] = z;
01243             i = 0;
01244             }
01245         else {
01246             x[0] = z & 0xffff;
01247             x[1] = z >> 16;
01248             i = 1;
01249             }
01250         k += 32;
01251         }
01252     while(!x[i])
01253         --i;
01254     b->wds = i + 1;
01255 #endif
01256 #ifndef Sudden_Underflow
01257     if (de) {
01258 #endif
01259 #ifdef IBM
01260         *e = (de - Bias - (P-1) << 2) + k;
01261         *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
01262 #else
01263         *e = de - Bias - (P-1) + k;
01264         *bits = P - k;
01265 #endif
01266 #ifndef Sudden_Underflow
01267         }
01268     else {
01269         *e = de - Bias - (P-1) + 1 + k;
01270 #ifdef Pack_32
01271         *bits = 32*i - hi0bits(x[i-1]);
01272 #else
01273         *bits = (i+2)*16 - hi0bits(x[i]);
01274 #endif
01275         }
01276 #endif
01277     return b;
01278     }
01279 #undef d0
01280 #undef d1
01281 
01282  static double
01283 ratio
01284     (Bigint *a, Bigint *b)
01285 {
01286     U da, db;
01287     int k, ka, kb;
01288 
01289     dval(da) = b2d(a, &ka);
01290     dval(db) = b2d(b, &kb);
01291 #ifdef Pack_32
01292     k = ka - kb + 32*(a->wds - b->wds);
01293 #else
01294     k = ka - kb + 16*(a->wds - b->wds);
01295 #endif
01296 #ifdef IBM
01297     if (k > 0) {
01298         word0(da) += (k >> 2)*Exp_msk1;
01299         if (k &= 3)
01300             dval(da) *= 1 << k;
01301         }
01302     else {
01303         k = -k;
01304         word0(db) += (k >> 2)*Exp_msk1;
01305         if (k &= 3)
01306             dval(db) *= 1 << k;
01307         }
01308 #else
01309     if (k > 0)
01310         word0(da) += k*Exp_msk1;
01311     else {
01312         k = -k;
01313         word0(db) += k*Exp_msk1;
01314         }
01315 #endif
01316     return dval(da) / dval(db);
01317     }
01318 
01319  static CONST double
01320 tens[] = {
01321         1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01322         1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01323         1e20, 1e21, 1e22
01324 #ifdef VAX
01325         , 1e23, 1e24
01326 #endif
01327         };
01328 
01329  static CONST double
01330 #ifdef IEEE_Arith
01331 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
01332 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
01333 #ifdef Avoid_Underflow
01334         9007199254740992.*9007199254740992.e-256
01335         /* = 2^106 * 1e-53 */
01336 #else
01337         1e-256
01338 #endif
01339         };
01340 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
01341 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
01342 #define Scale_Bit 0x10
01343 #define n_bigtens 5
01344 #else
01345 #ifdef IBM
01346 bigtens[] = { 1e16, 1e32, 1e64 };
01347 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
01348 #define n_bigtens 3
01349 #else
01350 bigtens[] = { 1e16, 1e32 };
01351 static CONST double tinytens[] = { 1e-16, 1e-32 };
01352 #define n_bigtens 2
01353 #endif
01354 #endif
01355 
01356 #ifndef IEEE_Arith
01357 #undef INFNAN_CHECK
01358 #endif
01359 
01360 #ifdef INFNAN_CHECK
01361 
01362 #ifndef NAN_WORD0
01363 #define NAN_WORD0 0x7ff80000
01364 #endif
01365 
01366 #ifndef NAN_WORD1
01367 #define NAN_WORD1 0
01368 #endif
01369 
01370  static int
01371 match
01372     (CONST char **sp, CONST char *t)
01373 {
01374     int c, d;
01375     CONST char *s = *sp;
01376 
01377     while((d = *t++)) {
01378         if ((c = *++s) >= 'A' && c <= 'Z')
01379             c += 'a' - 'A';
01380         if (c != d)
01381             return 0;
01382         }
01383     *sp = s + 1;
01384     return 1;
01385     }
01386 
01387 #ifndef No_Hex_NaN
01388  static void
01389 hexnan
01390     (U *rvp, CONST char **sp)
01391 {
01392     ULong c, x[2];
01393     CONST char *s;
01394     int havedig, udx0, xshift;
01395 
01396     x[0] = x[1] = 0;
01397     havedig = xshift = 0;
01398     udx0 = 1;
01399     s = *sp;
01400     while((c = *(CONST unsigned char*)++s)) {
01401         if (c >= '0' && c <= '9')
01402             c -= '0';
01403         else if (c >= 'a' && c <= 'f')
01404             c += 10 - 'a';
01405         else if (c >= 'A' && c <= 'F')
01406             c += 10 - 'A';
01407         else if (c <= ' ') {
01408             if (udx0 && havedig) {
01409                 udx0 = 0;
01410                 xshift = 1;
01411                 }
01412             continue;
01413             }
01414         else if (/*(*/ c == ')' && havedig) {
01415             *sp = s + 1;
01416             break;
01417             }
01418         else
01419             return; /* invalid form: don't change *sp */
01420         havedig = 1;
01421         if (xshift) {
01422             xshift = 0;
01423             x[0] = x[1];
01424             x[1] = 0;
01425             }
01426         if (udx0)
01427             x[0] = (x[0] << 4) | (x[1] >> 28);
01428         x[1] = (x[1] << 4) | c;
01429         }
01430     if ((x[0] &= 0xfffff) || x[1]) {
01431         word0(*rvp) = Exp_mask | x[0];
01432         word1(*rvp) = x[1];
01433         }
01434     }
01435 #endif /*No_Hex_NaN*/
01436 #endif /* INFNAN_CHECK */
01437 
01438  double
01439 strtod
01440     (CONST char *s00, char **se)
01441 {
01442 #ifdef Avoid_Underflow
01443     int scale;
01444 #endif
01445     int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01446          e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01447     CONST char *s, *s0, *s1;
01448     double aadj, aadj1, adj;
01449     U aadj2, rv, rv0;
01450     Long L;
01451     ULong y, z;
01452     Bigint *bb = NULL, *bb1 = NULL, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL;
01453 #ifdef SET_INEXACT
01454     int inexact, oldinexact;
01455 #endif
01456 #ifdef Honor_FLT_ROUNDS
01457     int rounding;
01458 #endif
01459 #ifdef USE_LOCALE
01460     CONST char *s2;
01461 #endif
01462 
01463     sign = nz0 = nz = 0;
01464     dval(rv) = 0.;
01465     for(s = s00;;s++) switch(*s) {
01466         case '-':
01467             sign = 1;
01468             /* no break */
01469         case '+':
01470             if (*++s)
01471                 goto break2;
01472             /* no break */
01473         case 0:
01474             goto ret0;
01475         case '\t':
01476         case '\n':
01477         case '\v':
01478         case '\f':
01479         case '\r':
01480         case ' ':
01481             continue;
01482         default:
01483             goto break2;
01484         }
01485  break2:
01486     if (*s == '0') {
01487         nz0 = 1;
01488         while(*++s == '0') ;
01489         if (!*s)
01490             goto ret;
01491         }
01492     s0 = s;
01493     y = z = 0;
01494     for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
01495         if (nd < 9)
01496             y = 10*y + c - '0';
01497         else if (nd < 16)
01498             z = 10*z + c - '0';
01499     nd0 = nd;
01500 #ifdef USE_LOCALE
01501     s1 = localeconv()->decimal_point;
01502     if (c == *s1) {
01503         c = '.';
01504         if (*++s1) {
01505             s2 = s;
01506             for(;;) {
01507                 if (*++s2 != *s1) {
01508                     c = 0;
01509                     break;
01510                     }
01511                 if (!*++s1) {
01512                     s = s2;
01513                     break;
01514                     }
01515                 }
01516             }
01517         }
01518 #endif
01519     if (c == '.') {
01520         c = *++s;
01521         if (!nd) {
01522             for(; c == '0'; c = *++s)
01523                 nz++;
01524             if (c > '0' && c <= '9') {
01525                 s0 = s;
01526                 nf += nz;
01527                 nz = 0;
01528                 goto have_dig;
01529                 }
01530             goto dig_done;
01531             }
01532         for(; c >= '0' && c <= '9'; c = *++s) {
01533  have_dig:
01534             nz++;
01535             if (c -= '0') {
01536                 nf += nz;
01537                 for(i = 1; i < nz; i++)
01538                     if (nd++ < 9)
01539                         y *= 10;
01540                     else if (nd <= DBL_DIG + 1)
01541                         z *= 10;
01542                 if (nd++ < 9)
01543                     y = 10*y + c;
01544                 else if (nd <= DBL_DIG + 1)
01545                     z = 10*z + c;
01546                 nz = 0;
01547                 }
01548             }
01549         }
01550  dig_done:
01551     e = 0;
01552     if (c == 'e' || c == 'E') {
01553         if (!nd && !nz && !nz0) {
01554             goto ret0;
01555             }
01556         s00 = s;
01557         esign = 0;
01558         switch(c = *++s) {
01559             case '-':
01560                 esign = 1;
01561             case '+':
01562                 c = *++s;
01563             }
01564         if (c >= '0' && c <= '9') {
01565             while(c == '0')
01566                 c = *++s;
01567             if (c > '0' && c <= '9') {
01568                 L = c - '0';
01569                 s1 = s;
01570                 while((c = *++s) >= '0' && c <= '9')
01571                     L = 10*L + c - '0';
01572                 if (s - s1 > 8 || L > 19999)
01573                     /* Avoid confusion from exponents
01574                      * so large that e might overflow.
01575                      */
01576                     e = 19999; /* safe for 16 bit ints */
01577                 else
01578                     e = (int)L;
01579                 if (esign)
01580                     e = -e;
01581                 }
01582             else
01583                 e = 0;
01584             }
01585         else
01586             s = s00;
01587         }
01588     if (!nd) {
01589         if (!nz && !nz0) {
01590 #ifdef INFNAN_CHECK
01591             /* Check for Nan and Infinity */
01592             switch(c) {
01593               case 'i':
01594               case 'I':
01595                 if (match(&s,"nf")) {
01596                     --s;
01597                     if (!match(&s,"inity"))
01598                         ++s;
01599                     word0(rv) = 0x7ff00000;
01600                     word1(rv) = 0;
01601                     goto ret;
01602                     }
01603                 break;
01604               case 'n':
01605               case 'N':
01606                 if (match(&s, "an")) {
01607                     word0(rv) = NAN_WORD0;
01608                     word1(rv) = NAN_WORD1;
01609 #ifndef No_Hex_NaN
01610                     if (*s == '(') /*)*/
01611                         hexnan(&rv, &s);
01612 #endif
01613                     goto ret;
01614                     }
01615               }
01616 #endif /* INFNAN_CHECK */
01617  ret0:
01618             s = s00;
01619             sign = 0;
01620             }
01621         goto ret;
01622         }
01623     e1 = e -= nf;
01624 
01625     /* Now we have nd0 digits, starting at s0, followed by a
01626      * decimal point, followed by nd-nd0 digits.  The number we're
01627      * after is the integer represented by those digits times
01628      * 10**e */
01629 
01630     if (!nd0)
01631         nd0 = nd;
01632     k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01633     dval(rv) = y;
01634     if (k > 9) {
01635 #ifdef SET_INEXACT
01636         if (k > DBL_DIG)
01637             oldinexact = get_inexact();
01638 #endif
01639         dval(rv) = tens[k - 9] * dval(rv) + z;
01640         }
01641     bd0 = 0;
01642     if (nd <= DBL_DIG
01643 #ifndef RND_PRODQUOT
01644 #ifndef Honor_FLT_ROUNDS
01645         && Flt_Rounds == 1
01646 #endif
01647 #endif
01648             ) {
01649         if (!e)
01650             goto ret;
01651         if (e > 0) {
01652             if (e <= Ten_pmax) {
01653 #ifdef VAX
01654                 goto vax_ovfl_check;
01655 #else
01656 #ifdef Honor_FLT_ROUNDS
01657                 /* round correctly FLT_ROUNDS = 2 or 3 */
01658                 if (sign) {
01659                     rv = -rv;
01660                     sign = 0;
01661                     }
01662 #endif
01663                 /* rv = */ rounded_product(dval(rv), tens[e]);
01664                 goto ret;
01665 #endif
01666                 }
01667             i = DBL_DIG - nd;
01668             if (e <= Ten_pmax + i) {
01669                 /* A fancier test would sometimes let us do
01670                  * this for larger i values.
01671                  */
01672 #ifdef Honor_FLT_ROUNDS
01673                 /* round correctly FLT_ROUNDS = 2 or 3 */
01674                 if (sign) {
01675                     rv = -rv;
01676                     sign = 0;
01677                     }
01678 #endif
01679                 e -= i;
01680                 dval(rv) *= tens[i];
01681 #ifdef VAX
01682                 /* VAX exponent range is so narrow we must
01683                  * worry about overflow here...
01684                  */
01685  vax_ovfl_check:
01686                 word0(rv) -= P*Exp_msk1;
01687                 /* rv = */ rounded_product(dval(rv), tens[e]);
01688                 if ((word0(rv) & Exp_mask)
01689                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
01690                     goto ovfl;
01691                 word0(rv) += P*Exp_msk1;
01692 #else
01693                 /* rv = */ rounded_product(dval(rv), tens[e]);
01694 #endif
01695                 goto ret;
01696                 }
01697             }
01698 #ifndef Inaccurate_Divide
01699         else if (e >= -Ten_pmax) {
01700 #ifdef Honor_FLT_ROUNDS
01701             /* round correctly FLT_ROUNDS = 2 or 3 */
01702             if (sign) {
01703                 rv = -rv;
01704                 sign = 0;
01705                 }
01706 #endif
01707             /* rv = */ rounded_quotient(dval(rv), tens[-e]);
01708             goto ret;
01709             }
01710 #endif
01711         }
01712     e1 += nd - k;
01713 
01714 #ifdef IEEE_Arith
01715 #ifdef SET_INEXACT
01716     inexact = 1;
01717     if (k <= DBL_DIG)
01718         oldinexact = get_inexact();
01719 #endif
01720 #ifdef Avoid_Underflow
01721     scale = 0;
01722 #endif
01723 #ifdef Honor_FLT_ROUNDS
01724     if ((rounding = Flt_Rounds) >= 2) {
01725         if (sign)
01726             rounding = rounding == 2 ? 0 : 2;
01727         else
01728             if (rounding != 2)
01729                 rounding = 0;
01730         }
01731 #endif
01732 #endif /*IEEE_Arith*/
01733 
01734     /* Get starting approximation = rv * 10**e1 */
01735 
01736     if (e1 > 0) {
01737         if ((i = e1 & 15))
01738             dval(rv) *= tens[i];
01739         if (e1 &= ~15) {
01740             if (e1 > DBL_MAX_10_EXP) {
01741  ovfl:
01742 #ifndef NO_ERRNO
01743                 errno = ERANGE;
01744 #endif
01745                 /* Can't trust HUGE_VAL */
01746 #ifdef IEEE_Arith
01747 #ifdef Honor_FLT_ROUNDS
01748                 switch(rounding) {
01749                   case 0: /* toward 0 */
01750                   case 3: /* toward -infinity */
01751                     word0(rv) = Big0;
01752                     word1(rv) = Big1;
01753                     break;
01754                   default:
01755                     word0(rv) = Exp_mask;
01756                     word1(rv) = 0;
01757                   }
01758 #else /*Honor_FLT_ROUNDS*/
01759                 word0(rv) = Exp_mask;
01760                 word1(rv) = 0;
01761 #endif /*Honor_FLT_ROUNDS*/
01762 #ifdef SET_INEXACT
01763                 /* set overflow bit */
01764                 dval(rv0) = 1e300;
01765                 dval(rv0) *= dval(rv0);
01766 #endif
01767 #else /*IEEE_Arith*/
01768                 word0(rv) = Big0;
01769                 word1(rv) = Big1;
01770 #endif /*IEEE_Arith*/
01771                 if (bd0)
01772                     goto retfree;
01773                 goto ret;
01774                 }
01775             e1 >>= 4;
01776             for(j = 0; e1 > 1; j++, e1 >>= 1)
01777                 if (e1 & 1)
01778                     dval(rv) *= bigtens[j];
01779         /* The last multiplication could overflow. */
01780             word0(rv) -= P*Exp_msk1;
01781             dval(rv) *= bigtens[j];
01782             if ((z = word0(rv) & Exp_mask)
01783              > Exp_msk1*(DBL_MAX_EXP+Bias-P))
01784                 goto ovfl;
01785             if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
01786                 /* set to largest number */
01787                 /* (Can't trust DBL_MAX) */
01788                 word0(rv) = Big0;
01789                 word1(rv) = Big1;
01790                 }
01791             else
01792                 word0(rv) += P*Exp_msk1;
01793             }
01794         }
01795     else if (e1 < 0) {
01796         e1 = -e1;
01797         if ((i = e1 & 15))
01798             dval(rv) /= tens[i];
01799         if (e1 >>= 4) {
01800             if (e1 >= 1 << n_bigtens)
01801                 goto undfl;
01802 #ifdef Avoid_Underflow
01803             if (e1 & Scale_Bit)
01804                 scale = 2*P;
01805             for(j = 0; e1 > 0; j++, e1 >>= 1)
01806                 if (e1 & 1)
01807                     dval(rv) *= tinytens[j];
01808             if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
01809                         >> Exp_shift)) > 0) {
01810                 /* scaled rv is denormal; zap j low bits */
01811                 if (j >= 32) {
01812                     word1(rv) = 0;
01813                     if (j >= 53)
01814                      word0(rv) = (P+2)*Exp_msk1;
01815                     else
01816                      word0(rv) &= 0xffffffff << j-32;
01817                     }
01818                 else
01819                     word1(rv) &= 0xffffffff << j;
01820                 }
01821 #else
01822             for(j = 0; e1 > 1; j++, e1 >>= 1)
01823                 if (e1 & 1)
01824                     dval(rv) *= tinytens[j];
01825             /* The last multiplication could underflow. */
01826             dval(rv0) = dval(rv);
01827             dval(rv) *= tinytens[j];
01828             if (!dval(rv)) {
01829                 dval(rv) = 2.*dval(rv0);
01830                 dval(rv) *= tinytens[j];
01831 #endif
01832                 if (!dval(rv)) {
01833  undfl:
01834                     dval(rv) = 0.;
01835 #ifndef NO_ERRNO
01836                     errno = ERANGE;
01837 #endif
01838                     if (bd0)
01839                         goto retfree;
01840                     goto ret;
01841                     }
01842 #ifndef Avoid_Underflow
01843                 word0(rv) = Tiny0;
01844                 word1(rv) = Tiny1;
01845                 /* The refinement below will clean
01846                  * this approximation up.
01847                  */
01848                 }
01849 #endif
01850             }
01851         }
01852 
01853     /* Now the hard part -- adjusting rv to the correct value.*/
01854 
01855     /* Put digits into bd: true value = bd * 10^e */
01856 
01857     bd0 = s2b(s0, nd0, nd, y);
01858 
01859     for(;;) {
01860         bd = Balloc(bd0->k);
01861         Bcopy(bd, bd0);
01862         bb = d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
01863         bs = i2b(1);
01864 
01865         if (e >= 0) {
01866             bb2 = bb5 = 0;
01867             bd2 = bd5 = e;
01868             }
01869         else {
01870             bb2 = bb5 = -e;
01871             bd2 = bd5 = 0;
01872             }
01873         if (bbe >= 0)
01874             bb2 += bbe;
01875         else
01876             bd2 -= bbe;
01877         bs2 = bb2;
01878 #ifdef Honor_FLT_ROUNDS
01879         if (rounding != 1)
01880             bs2++;
01881 #endif
01882 #ifdef Avoid_Underflow
01883         j = bbe - scale;
01884         i = j + bbbits - 1; /* logb(rv) */
01885         if (i < Emin)   /* denormal */
01886             j += P - Emin;
01887         else
01888             j = P + 1 - bbbits;
01889 #else /*Avoid_Underflow*/
01890 #ifdef Sudden_Underflow
01891 #ifdef IBM
01892         j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
01893 #else
01894         j = P + 1 - bbbits;
01895 #endif
01896 #else /*Sudden_Underflow*/
01897         j = bbe;
01898         i = j + bbbits - 1; /* logb(rv) */
01899         if (i < Emin)   /* denormal */
01900             j += P - Emin;
01901         else
01902             j = P + 1 - bbbits;
01903 #endif /*Sudden_Underflow*/
01904 #endif /*Avoid_Underflow*/
01905         bb2 += j;
01906         bd2 += j;
01907 #ifdef Avoid_Underflow
01908         bd2 += scale;
01909 #endif
01910         i = bb2 < bd2 ? bb2 : bd2;
01911         if (i > bs2)
01912             i = bs2;
01913         if (i > 0) {
01914             bb2 -= i;
01915             bd2 -= i;
01916             bs2 -= i;
01917             }
01918         if (bb5 > 0) {
01919             bs = pow5mult(bs, bb5);
01920             bb1 = mult(bs, bb);
01921             Bfree(bb);
01922             bb = bb1;
01923             }
01924         if (bb2 > 0)
01925             bb = lshift(bb, bb2);
01926         if (bd5 > 0)
01927             bd = pow5mult(bd, bd5);
01928         if (bd2 > 0)
01929             bd = lshift(bd, bd2);
01930         if (bs2 > 0)
01931             bs = lshift(bs, bs2);
01932         delta = diff(bb, bd);
01933         dsign = delta->sign;
01934         delta->sign = 0;
01935         i = cmp(delta, bs);
01936 #ifdef Honor_FLT_ROUNDS
01937         if (rounding != 1) {
01938             if (i < 0) {
01939                 /* Error is less than an ulp */
01940                 if (!delta->x[0] && delta->wds <= 1) {
01941                     /* exact */
01942 #ifdef SET_INEXACT
01943                     inexact = 0;
01944 #endif
01945                     break;
01946                     }
01947                 if (rounding) {
01948                     if (dsign) {
01949                         adj = 1.;
01950                         goto apply_adj;
01951                         }
01952                     }
01953                 else if (!dsign) {
01954                     adj = -1.;
01955                     if (!word1(rv)
01956                      && !(word0(rv) & Frac_mask)) {
01957                         y = word0(rv) & Exp_mask;
01958 #ifdef Avoid_Underflow
01959                         if (!scale || y > 2*P*Exp_msk1)
01960 #else
01961                         if (y)
01962 #endif
01963                           {
01964                           delta = lshift(delta,Log2P);
01965                           if (cmp(delta, bs) <= 0)
01966                             adj = -0.5;
01967                           }
01968                         }
01969  apply_adj:
01970 #ifdef Avoid_Underflow
01971                     if (scale && (y = word0(rv) & Exp_mask)
01972                         <= 2*P*Exp_msk1)
01973                       word0(adj) += (2*P+1)*Exp_msk1 - y;
01974 #else
01975 #ifdef Sudden_Underflow
01976                     if ((word0(rv) & Exp_mask) <=
01977                             P*Exp_msk1) {
01978                         word0(rv) += P*Exp_msk1;
01979                         dval(rv) += adj*ulp(dval(rv));
01980                         word0(rv) -= P*Exp_msk1;
01981                         }
01982                     else
01983 #endif /*Sudden_Underflow*/
01984 #endif /*Avoid_Underflow*/
01985                     dval(rv) += adj*ulp(dval(rv));
01986                     }
01987                 break;
01988                 }
01989             adj = ratio(delta, bs);
01990             if (adj < 1.)
01991                 adj = 1.;
01992             if (adj <= 0x7ffffffe) {
01993                 /* adj = rounding ? ceil(adj) : floor(adj); */
01994                 y = adj;
01995                 if (y != adj) {
01996                     if (!((rounding>>1) ^ dsign))
01997                         y++;
01998                     adj = y;
01999                     }
02000                 }
02001 #ifdef Avoid_Underflow
02002             if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02003                 word0(adj) += (2*P+1)*Exp_msk1 - y;
02004 #else
02005 #ifdef Sudden_Underflow
02006             if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02007                 word0(rv) += P*Exp_msk1;
02008                 adj *= ulp(dval(rv));
02009                 if (dsign)
02010                     dval(rv) += adj;
02011                 else
02012                     dval(rv) -= adj;
02013                 word0(rv) -= P*Exp_msk1;
02014                 goto cont;
02015                 }
02016 #endif /*Sudden_Underflow*/
02017 #endif /*Avoid_Underflow*/
02018             adj *= ulp(dval(rv));
02019             if (dsign)
02020                 dval(rv) += adj;
02021             else
02022                 dval(rv) -= adj;
02023             goto cont;
02024             }
02025 #endif /*Honor_FLT_ROUNDS*/
02026 
02027         if (i < 0) {
02028             /* Error is less than half an ulp -- check for
02029              * special case of mantissa a power of two.
02030              */
02031             if (dsign || word1(rv) || word0(rv) & Bndry_mask
02032 #ifdef IEEE_Arith
02033 #ifdef Avoid_Underflow
02034              || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02035 #else
02036              || (word0(rv) & Exp_mask) <= Exp_msk1
02037 #endif
02038 #endif
02039                 ) {
02040 #ifdef SET_INEXACT
02041                 if (!delta->x[0] && delta->wds <= 1)
02042                     inexact = 0;
02043 #endif
02044                 break;
02045                 }
02046             if (!delta->x[0] && delta->wds <= 1) {
02047                 /* exact result */
02048 #ifdef SET_INEXACT
02049                 inexact = 0;
02050 #endif
02051                 break;
02052                 }
02053             delta = lshift(delta,Log2P);
02054             if (cmp(delta, bs) > 0)
02055                 goto drop_down;
02056             break;
02057             }
02058         if (i == 0) {
02059             /* exactly half-way between */
02060             if (dsign) {
02061                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02062                  &&  word1(rv) == (
02063 #ifdef Avoid_Underflow
02064             (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02065         ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02066 #endif
02067                            0xffffffff)) {
02068                     /*boundary case -- increment exponent*/
02069                     word0(rv) = (word0(rv) & Exp_mask)
02070                         + Exp_msk1
02071 #ifdef IBM
02072                         | Exp_msk1 >> 4
02073 #endif
02074                         ;
02075                     word1(rv) = 0;
02076 #ifdef Avoid_Underflow
02077                     dsign = 0;
02078 #endif
02079                     break;
02080                     }
02081                 }
02082             else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02083  drop_down:
02084                 /* boundary case -- decrement exponent */
02085 #ifdef Sudden_Underflow /*{{*/
02086                 L = word0(rv) & Exp_mask;
02087 #ifdef IBM
02088                 if (L <  Exp_msk1)
02089 #else
02090 #ifdef Avoid_Underflow
02091                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02092 #else
02093                 if (L <= Exp_msk1)
02094 #endif /*Avoid_Underflow*/
02095 #endif /*IBM*/
02096                     goto undfl;
02097                 L -= Exp_msk1;
02098 #else /*Sudden_Underflow}{*/
02099 #ifdef Avoid_Underflow
02100                 if (scale) {
02101                     L = word0(rv) & Exp_mask;
02102                     if (L <= (2*P+1)*Exp_msk1) {
02103                         if (L > (P+2)*Exp_msk1)
02104                             /* round even ==> */
02105                             /* accept rv */
02106                             break;
02107                         /* rv = smallest denormal */
02108                         goto undfl;
02109                         }
02110                     }
02111 #endif /*Avoid_Underflow*/
02112                 L = (word0(rv) & Exp_mask) - Exp_msk1;
02113 #endif /*Sudden_Underflow}}*/
02114                 word0(rv) = L | Bndry_mask1;
02115                 word1(rv) = 0xffffffff;
02116 #ifdef IBM
02117                 goto cont;
02118 #else
02119                 break;
02120 #endif
02121                 }
02122 #ifndef ROUND_BIASED
02123             if (!(word1(rv) & LSB))
02124                 break;
02125 #endif
02126             if (dsign)
02127                 dval(rv) += ulp(dval(rv));
02128 #ifndef ROUND_BIASED
02129             else {
02130                 dval(rv) -= ulp(dval(rv));
02131 #ifndef Sudden_Underflow
02132                 if (!dval(rv))
02133                     goto undfl;
02134 #endif
02135                 }
02136 #ifdef Avoid_Underflow
02137             dsign = 1 - dsign;
02138 #endif
02139 #endif
02140             break;
02141             }
02142         if ((aadj = ratio(delta, bs)) <= 2.) {
02143             if (dsign)
02144                 aadj = aadj1 = 1.;
02145             else if (word1(rv) || word0(rv) & Bndry_mask) {
02146 #ifndef Sudden_Underflow
02147                 if (word1(rv) == Tiny1 && !word0(rv))
02148                     goto undfl;
02149 #endif
02150                 aadj = 1.;
02151                 aadj1 = -1.;
02152                 }
02153             else {
02154                 /* special case -- power of FLT_RADIX to be */
02155                 /* rounded down... */
02156 
02157                 if (aadj < 2./FLT_RADIX)
02158                     aadj = 1./FLT_RADIX;
02159                 else
02160                     aadj *= 0.5;
02161                 aadj1 = -aadj;
02162                 }
02163             }
02164         else {
02165             aadj *= 0.5;
02166             aadj1 = dsign ? aadj : -aadj;
02167 #ifdef Check_FLT_ROUNDS
02168             switch(Rounding) {
02169                 case 2: /* towards +infinity */
02170                     aadj1 -= 0.5;
02171                     break;
02172                 case 0: /* towards 0 */
02173                 case 3: /* towards -infinity */
02174                     aadj1 += 0.5;
02175                 }
02176 #else
02177             if (Flt_Rounds == 0)
02178                 aadj1 += 0.5;
02179 #endif /*Check_FLT_ROUNDS*/
02180             }
02181         y = word0(rv) & Exp_mask;
02182 
02183         /* Check for overflow */
02184 
02185         if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02186             dval(rv0) = dval(rv);
02187             word0(rv) -= P*Exp_msk1;
02188             adj = aadj1 * ulp(dval(rv));
02189             dval(rv) += adj;
02190             if ((word0(rv) & Exp_mask) >=
02191                     Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02192                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
02193                     goto ovfl;
02194                 word0(rv) = Big0;
02195                 word1(rv) = Big1;
02196                 goto cont;
02197                 }
02198             else
02199                 word0(rv) += P*Exp_msk1;
02200             }
02201         else {
02202 #ifdef Avoid_Underflow
02203             if (scale && y <= 2*P*Exp_msk1) {
02204                 if (aadj <= 0x7fffffff) {
02205                     if ((z = (ULong)aadj) <= 0)
02206                         z = 1;
02207                     aadj = z;
02208                     aadj1 = dsign ? aadj : -aadj;
02209                     }
02210                 dval(aadj2) = aadj1;
02211                 word0(aadj2) += (2*P+1)*Exp_msk1 - y;
02212                 aadj1 = dval(aadj2);
02213                 }
02214             adj = aadj1 * ulp(dval(rv));
02215             dval(rv) += adj;
02216 #else
02217 #ifdef Sudden_Underflow
02218             if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02219                 dval(rv0) = dval(rv);
02220                 word0(rv) += P*Exp_msk1;
02221                 adj = aadj1 * ulp(dval(rv));
02222                 dval(rv) += adj;
02223 #ifdef IBM
02224                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
02225 #else
02226                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02227 #endif
02228                     {
02229                     if (word0(rv0) == Tiny0
02230                      && word1(rv0) == Tiny1)
02231                         goto undfl;
02232                     word0(rv) = Tiny0;
02233                     word1(rv) = Tiny1;
02234                     goto cont;
02235                     }
02236                 else
02237                     word0(rv) -= P*Exp_msk1;
02238                 }
02239             else {
02240                 adj = aadj1 * ulp(dval(rv));
02241                 dval(rv) += adj;
02242                 }
02243 #else /*Sudden_Underflow*/
02244             /* Compute adj so that the IEEE rounding rules will
02245              * correctly round rv + adj in some half-way cases.
02246              * If rv * ulp(rv) is denormalized (i.e.,
02247              * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
02248              * trouble from bits lost to denormalization;
02249              * example: 1.2e-307 .
02250              */
02251             if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02252                 aadj1 = (double)(int)(aadj + 0.5);
02253                 if (!dsign)
02254                     aadj1 = -aadj1;
02255                 }
02256             adj = aadj1 * ulp(dval(rv));
02257             dval(rv) += adj;
02258 #endif /*Sudden_Underflow*/
02259 #endif /*Avoid_Underflow*/
02260             }
02261         z = word0(rv) & Exp_mask;
02262 #ifndef SET_INEXACT
02263 #ifdef Avoid_Underflow
02264         if (!scale)
02265 #endif
02266         if (y == z) {
02267             /* Can we stop now? */
02268             L = (Long)aadj;
02269             aadj -= L;
02270             /* The tolerances below are conservative. */
02271             if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02272                 if (aadj < .4999999 || aadj > .5000001)
02273                     break;
02274                 }
02275             else if (aadj < .4999999/FLT_RADIX)
02276                 break;
02277             }
02278 #endif
02279  cont:
02280         Bfree(bb);
02281         Bfree(bd);
02282         Bfree(bs);
02283         Bfree(delta);
02284         }
02285 #ifdef SET_INEXACT
02286     if (inexact) {
02287         if (!oldinexact) {
02288             word0(rv0) = Exp_1 + (70 << Exp_shift);
02289             word1(rv0) = 0;
02290             dval(rv0) += 1.;
02291             }
02292         }
02293     else if (!oldinexact)
02294         clear_inexact();
02295 #endif
02296 #ifdef Avoid_Underflow
02297     if (scale) {
02298         word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02299         word1(rv0) = 0;
02300         dval(rv) *= dval(rv0);
02301 #ifndef NO_ERRNO
02302         /* try to avoid the bug of testing an 8087 register value */
02303         if (word0(rv) == 0 && word1(rv) == 0)
02304             errno = ERANGE;
02305 #endif
02306         }
02307 #endif /* Avoid_Underflow */
02308 #ifdef SET_INEXACT
02309     if (inexact && !(word0(rv) & Exp_mask)) {
02310         /* set underflow bit */
02311         dval(rv0) = 1e-300;
02312         dval(rv0) *= dval(rv0);
02313         }
02314 #endif
02315  retfree:
02316     Bfree(bb);
02317     Bfree(bd);
02318     Bfree(bs);
02319     Bfree(bd0);
02320     Bfree(delta);
02321  ret:
02322     if (se)
02323         *se = (char *)s;
02324     return sign ? -dval(rv) : dval(rv);
02325     }
02326 
02327  static int
02328 quorem
02329     (Bigint *b, Bigint *S)
02330 {
02331     int n;
02332     ULong *bx, *bxe, q, *sx, *sxe;
02333 #ifdef ULLong
02334     ULLong borrow, carry, y, ys;
02335 #else
02336     ULong borrow, carry, y, ys;
02337 #ifdef Pack_32
02338     ULong si, z, zs;
02339 #endif
02340 #endif
02341 
02342     n = S->wds;
02343 #ifdef DEBUG
02344     /*debug*/ if (b->wds > n)
02345     /*debug*/   Bug("oversize b in quorem");
02346 #endif
02347     if (b->wds < n)
02348         return 0;
02349     sx = S->x;
02350     sxe = sx + --n;
02351     bx = b->x;
02352     bxe = bx + n;
02353     q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
02354 #ifdef DEBUG
02355     /*debug*/ if (q > 9)
02356     /*debug*/   Bug("oversized quotient in quorem");
02357 #endif
02358     if (q) {
02359         borrow = 0;
02360         carry = 0;
02361         do {
02362 #ifdef ULLong
02363             ys = *sx++ * (ULLong)q + carry;
02364             carry = ys >> 32;
02365             y = *bx - (ys & FFFFFFFF) - borrow;
02366             borrow = y >> 32 & (ULong)1;
02367             *bx++ = (ULong)y & FFFFFFFF;
02368 #else
02369 #ifdef Pack_32
02370             si = *sx++;
02371             ys = (si & 0xffff) * q + carry;
02372             zs = (si >> 16) * q + (ys >> 16);
02373             carry = zs >> 16;
02374             y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
02375             borrow = (y & 0x10000) >> 16;
02376             z = (*bx >> 16) - (zs & 0xffff) - borrow;
02377             borrow = (z & 0x10000) >> 16;
02378             Storeinc(bx, z, y);
02379 #else
02380             ys = *sx++ * q + carry;
02381             carry = ys >> 16;
02382             y = *bx - (ys & 0xffff) - borrow;
02383             borrow = (y & 0x10000) >> 16;
02384             *bx++ = y & 0xffff;
02385 #endif
02386 #endif
02387             }
02388             while(sx <= sxe);
02389         if (!*bxe) {
02390             bx = b->x;
02391             while(--bxe > bx && !*bxe)
02392                 --n;
02393             b->wds = n;
02394             }
02395         }
02396     if (cmp(b, S) >= 0) {
02397         q++;
02398         borrow = 0;
02399         carry = 0;
02400         bx = b->x;
02401         sx = S->x;
02402         do {
02403 #ifdef ULLong
02404             ys = *sx++ + carry;
02405             carry = ys >> 32;
02406             y = *bx - (ys & FFFFFFFF) - borrow;
02407             borrow = y >> 32 & (ULong)1;
02408             *bx++ = (ULong)y & FFFFFFFF;
02409 #else
02410 #ifdef Pack_32
02411             si = *sx++;
02412             ys = (si & 0xffff) + carry;
02413             zs = (si >> 16) + (ys >> 16);
02414             carry = zs >> 16;
02415             y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
02416             borrow = (y & 0x10000) >> 16;
02417             z = (*bx >> 16) - (zs & 0xffff) - borrow;
02418             borrow = (z & 0x10000) >> 16;
02419             Storeinc(bx, z, y);
02420 #else
02421             ys = *sx++ + carry;
02422             carry = ys >> 16;
02423             y = *bx - (ys & 0xffff) - borrow;
02424             borrow = (y & 0x10000) >> 16;
02425             *bx++ = y & 0xffff;
02426 #endif
02427 #endif
02428             }
02429             while(sx <= sxe);
02430         bx = b->x;
02431         bxe = bx + n;
02432         if (!*bxe) {
02433             while(--bxe > bx && !*bxe)
02434                 --n;
02435             b->wds = n;
02436             }
02437         }
02438     return q;
02439     }
02440 
02441 #ifndef MULTIPLE_THREADS
02442  static char *dtoa_result;
02443 #endif
02444 
02445  static char *
02446 rv_alloc(int i)
02447 {
02448     int j, k, *r;
02449 
02450     j = sizeof(ULong);
02451     for(k = 0;
02452         sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
02453         j <<= 1)
02454             k++;
02455     r = (int*)Balloc(k);
02456     *r = k;
02457     return
02458 #ifndef MULTIPLE_THREADS
02459     dtoa_result =
02460 #endif
02461         (char *)(r+1);
02462     }
02463 
02464  static char *
02465 nrv_alloc(CONST char *s, char **rve, int n)
02466 {
02467     char *rv, *t;
02468 
02469     t = rv = rv_alloc(n);
02470     while((*t = *s++)) t++;
02471     if (rve)
02472         *rve = t;
02473     return rv;
02474     }
02475 
02476 /* freedtoa(s) must be used to free values s returned by dtoa
02477  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
02478  * but for consistency with earlier versions of dtoa, it is optional
02479  * when MULTIPLE_THREADS is not defined.
02480  */
02481 
02482  void
02483 freedtoa(char *s)
02484 {
02485     Bigint *b = (Bigint *)((int *)s - 1);
02486     b->maxwds = 1 << (b->k = *(int*)b);
02487     Bfree(b);
02488 #ifndef MULTIPLE_THREADS
02489     if (s == dtoa_result)
02490         dtoa_result = 0;
02491 #endif
02492     }
02493 
02494 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
02495  *
02496  * Inspired by "How to Print Floating-Point Numbers Accurately" by
02497  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
02498  *
02499  * Modifications:
02500  *  1. Rather than iterating, we use a simple numeric overestimate
02501  *     to determine k = floor(log10(d)).  We scale relevant
02502  *     quantities using O(log2(k)) rather than O(k) multiplications.
02503  *  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
02504  *     try to generate digits strictly left to right.  Instead, we
02505  *     compute with fewer bits and propagate the carry if necessary
02506  *     when rounding the final digit up.  This is often faster.
02507  *  3. Under the assumption that input will be rounded nearest,
02508  *     mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
02509  *     That is, we allow equality in stopping tests when the
02510  *     round-nearest rule will give the same floating-point value
02511  *     as would satisfaction of the stopping test with strict
02512  *     inequality.
02513  *  4. We remove common factors of powers of 2 from relevant
02514  *     quantities.
02515  *  5. When converting floating-point integers less than 1e16,
02516  *     we use floating-point arithmetic rather than resorting
02517  *     to multiple-precision integers.
02518  *  6. When asked to produce fewer than 15 digits, we first try
02519  *     to get by with floating-point arithmetic; we resort to
02520  *     multiple-precision integer arithmetic only if we cannot
02521  *     guarantee that the floating-point calculation has given
02522  *     the correctly rounded result.  For k requested digits and
02523  *     "uniformly" distributed input, the probability is
02524  *     something like 10^(k-15) that we must resort to the Long
02525  *     calculation.
02526  */
02527 
02528  char *
02529 dtoa
02530     (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
02531 {
02532  /* Arguments ndigits, decpt, sign are similar to those
02533     of ecvt and fcvt; trailing zeros are suppressed from
02534     the returned string.  If not null, *rve is set to point
02535     to the end of the return value.  If d is +-Infinity or NaN,
02536     then *decpt is set to 9999.
02537 
02538     mode:
02539         0 ==> shortest string that yields d when read in
02540             and rounded to nearest.
02541         1 ==> like 0, but with Steele & White stopping rule;
02542             e.g. with IEEE P754 arithmetic , mode 0 gives
02543             1e23 whereas mode 1 gives 9.999999999999999e22.
02544         2 ==> max(1,ndigits) significant digits.  This gives a
02545             return value similar to that of ecvt, except
02546             that trailing zeros are suppressed.
02547         3 ==> through ndigits past the decimal point.  This
02548             gives a return value similar to that from fcvt,
02549             except that trailing zeros are suppressed, and
02550             ndigits can be negative.
02551         4,5 ==> similar to 2 and 3, respectively, but (in
02552             round-nearest mode) with the tests of mode 0 to
02553             possibly return a shorter string that rounds to d.
02554             With IEEE arithmetic and compilation with
02555             -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
02556             as modes 2 and 3 when FLT_ROUNDS != 1.
02557         6-9 ==> Debugging modes similar to mode - 4:  don't try
02558             fast floating-point estimate (if applicable).
02559 
02560         Values of mode other than 0-9 are treated as mode 0.
02561 
02562         Sufficient space is allocated to the return value
02563         to hold the suppressed trailing zeros.
02564     */
02565 
02566     int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0,
02567         j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
02568         spec_case, try_quick;
02569     Long L;
02570 #ifndef Sudden_Underflow
02571     int denorm;
02572     ULong x;
02573 #endif
02574     Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
02575     U d, d2, eps;
02576     double ds;
02577     char *s, *s0;
02578 #ifdef Honor_FLT_ROUNDS
02579     int rounding;
02580 #endif
02581 #ifdef SET_INEXACT
02582     int inexact, oldinexact;
02583 #endif
02584 
02585 #ifndef MULTIPLE_THREADS
02586     if (dtoa_result) {
02587         freedtoa(dtoa_result);
02588         dtoa_result = 0;
02589         }
02590 #endif
02591 
02592     dval(d) = dd;
02593     if (word0(d) & Sign_bit) {
02594         /* set sign for everything, including 0's and NaNs */
02595         *sign = 1;
02596         word0(d) &= ~Sign_bit;  /* clear sign bit */
02597         }
02598     else
02599         *sign = 0;
02600 
02601 #if defined(IEEE_Arith) + defined(VAX)
02602 #ifdef IEEE_Arith
02603     if ((word0(d) & Exp_mask) == Exp_mask)
02604 #else
02605     if (word0(d)  == 0x8000)
02606 #endif
02607         {
02608         /* Infinity or NaN */
02609         *decpt = 9999;
02610 #ifdef IEEE_Arith
02611         if (!word1(d) && !(word0(d) & 0xfffff))
02612             return nrv_alloc("Infinity", rve, 8);
02613 #endif
02614         return nrv_alloc("NaN", rve, 3);
02615         }
02616 #endif
02617 #ifdef IBM
02618     dval(d) += 0; /* normalize */
02619 #endif
02620     if (!dval(d)) {
02621         *decpt = 1;
02622         return nrv_alloc("0", rve, 1);
02623         }
02624 
02625 #ifdef SET_INEXACT
02626     try_quick = oldinexact = get_inexact();
02627     inexact = 1;
02628 #endif
02629 #ifdef Honor_FLT_ROUNDS
02630     if ((rounding = Flt_Rounds) >= 2) {
02631         if (*sign)
02632             rounding = rounding == 2 ? 0 : 2;
02633         else
02634             if (rounding != 2)
02635                 rounding = 0;
02636         }
02637 #endif
02638 
02639     b = d2b(dval(d), &be, &bbits);
02640 #ifdef Sudden_Underflow
02641     i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
02642 #else
02643     if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
02644 #endif
02645         dval(d2) = dval(d);
02646         word0(d2) &= Frac_mask1;
02647         word0(d2) |= Exp_11;
02648 #ifdef IBM
02649         if (j = 11 - hi0bits(word0(d2) & Frac_mask))
02650             dval(d2) /= 1 << j;
02651 #endif
02652 
02653         /* log(x)   ~=~ log(1.5) + (x-1.5)/1.5
02654          * log10(x)  =  log(x) / log(10)
02655          *      ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
02656          * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
02657          *
02658          * This suggests computing an approximation k to log10(d) by
02659          *
02660          * k = (i - Bias)*0.301029995663981
02661          *  + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
02662          *
02663          * We want k to be too large rather than too small.
02664          * The error in the first-order Taylor series approximation
02665          * is in our favor, so we just round up the constant enough
02666          * to compensate for any error in the multiplication of
02667          * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
02668          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
02669          * adding 1e-13 to the constant term more than suffices.
02670          * Hence we adjust the constant term to 0.1760912590558.
02671          * (We could get a more accurate k by invoking log10,
02672          *  but this is probably not worthwhile.)
02673          */
02674 
02675         i -= Bias;
02676 #ifdef IBM
02677         i <<= 2;
02678         i += j;
02679 #endif
02680 #ifndef Sudden_Underflow
02681         denorm = 0;
02682         }
02683     else {
02684         /* d is denormalized */
02685 
02686         i = bbits + be + (Bias + (P-1) - 1);
02687         x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
02688                 : word1(d) << 32 - i;
02689         dval(d2) = x;
02690         word0(d2) -= 31*Exp_msk1; /* adjust exponent */
02691         i -= (Bias + (P-1) - 1) + 1;
02692         denorm = 1;
02693         }
02694 #endif
02695     ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
02696     k = (int)ds;
02697     if (ds < 0. && ds != k)
02698         k--;    /* want k = floor(ds) */
02699     k_check = 1;
02700     if (k >= 0 && k <= Ten_pmax) {
02701         if (dval(d) < tens[k])
02702             k--;
02703         k_check = 0;
02704         }
02705     j = bbits - i - 1;
02706     if (j >= 0) {
02707         b2 = 0;
02708         s2 = j;
02709         }
02710     else {
02711         b2 = -j;
02712         s2 = 0;
02713         }
02714     if (k >= 0) {
02715         b5 = 0;
02716         s5 = k;
02717         s2 += k;
02718         }
02719     else {
02720         b2 -= k;
02721         b5 = -k;
02722         s5 = 0;
02723         }
02724     if (mode < 0 || mode > 9)
02725         mode = 0;
02726 
02727 #ifndef SET_INEXACT
02728 #ifdef Check_FLT_ROUNDS
02729     try_quick = Rounding == 1;
02730 #else
02731     try_quick = 1;
02732 #endif
02733 #endif /*SET_INEXACT*/
02734 
02735     if (mode > 5) {
02736         mode -= 4;
02737         try_quick = 0;
02738         }
02739     leftright = 1;
02740     switch(mode) {
02741         case 0:
02742         case 1:
02743             ilim = ilim1 = -1;
02744             i = 18;
02745             ndigits = 0;
02746             break;
02747         case 2:
02748             leftright = 0;
02749             /* no break */
02750         case 4:
02751             if (ndigits <= 0)
02752                 ndigits = 1;
02753             ilim = ilim1 = i = ndigits;
02754             break;
02755         case 3:
02756             leftright = 0;
02757             /* no break */
02758         case 5:
02759             i = ndigits + k + 1;
02760             ilim = i;
02761             ilim1 = i - 1;
02762             if (i <= 0)
02763                 i = 1;
02764         }
02765     s = s0 = rv_alloc(i);
02766 
02767 #ifdef Honor_FLT_ROUNDS
02768     if (mode > 1 && rounding != 1)
02769         leftright = 0;
02770 #endif
02771 
02772     if (ilim >= 0 && ilim <= Quick_max && try_quick) {
02773 
02774         /* Try to get by with floating-point arithmetic. */
02775 
02776         i = 0;
02777         dval(d2) = dval(d);
02778         k0 = k;
02779         ilim0 = ilim;
02780         ieps = 2; /* conservative */
02781         if (k > 0) {
02782             ds = tens[k&0xf];
02783             j = k >> 4;
02784             if (j & Bletch) {
02785                 /* prevent overflows */
02786                 j &= Bletch - 1;
02787                 dval(d) /= bigtens[n_bigtens-1];
02788                 ieps++;
02789                 }
02790             for(; j; j >>= 1, i++)
02791                 if (j & 1) {
02792                     ieps++;
02793                     ds *= bigtens[i];
02794                     }
02795             dval(d) /= ds;
02796             }
02797         else if ((j1 = -k)) {
02798             dval(d) *= tens[j1 & 0xf];
02799             for(j = j1 >> 4; j; j >>= 1, i++)
02800                 if (j & 1) {
02801                     ieps++;
02802                     dval(d) *= bigtens[i];
02803                     }
02804             }
02805         if (k_check && dval(d) < 1. && ilim > 0) {
02806             if (ilim1 <= 0)
02807                 goto fast_failed;
02808             ilim = ilim1;
02809             k--;
02810             dval(d) *= 10.;
02811             ieps++;
02812             }
02813         dval(eps) = ieps*dval(d) + 7.;
02814         word0(eps) -= (P-1)*Exp_msk1;
02815         if (ilim == 0) {
02816             S = mhi = 0;
02817             dval(d) -= 5.;
02818             if (dval(d) > dval(eps))
02819                 goto one_digit;
02820             if (dval(d) < -dval(eps))
02821                 goto no_digits;
02822             goto fast_failed;
02823             }
02824 #ifndef No_leftright
02825         if (leftright) {
02826             /* Use Steele & White method of only
02827              * generating digits needed.
02828              */
02829             dval(eps) = 0.5/tens[ilim-1] - dval(eps);
02830             for(i = 0;;) {
02831                 L = (long int)dval(d);
02832                 dval(d) -= L;
02833                 *s++ = '0' + (int)L;
02834                 if (dval(d) < dval(eps))
02835                     goto ret1;
02836                 if (1. - dval(d) < dval(eps))
02837                     goto bump_up;
02838                 if (++i >= ilim)
02839                     break;
02840                 dval(eps) *= 10.;
02841                 dval(d) *= 10.;
02842                 }
02843             }
02844         else {
02845 #endif
02846             /* Generate ilim digits, then fix them up. */
02847             dval(eps) *= tens[ilim-1];
02848             for(i = 1;; i++, dval(d) *= 10.) {
02849                 L = (Long)(dval(d));
02850                 if (!(dval(d) -= L))
02851                     ilim = i;
02852                 *s++ = '0' + (int)L;
02853                 if (i == ilim) {
02854                     if (dval(d) > 0.5 + dval(eps))
02855                         goto bump_up;
02856                     else if (dval(d) < 0.5 - dval(eps)) {
02857                         while(*--s == '0')
02858                             ;
02859                         s++;
02860                         goto ret1;
02861                         }
02862                     break;
02863                     }
02864                 }
02865 #ifndef No_leftright
02866             }
02867 #endif
02868  fast_failed:
02869         s = s0;
02870         dval(d) = dval(d2);
02871         k = k0;
02872         ilim = ilim0;
02873         }
02874 
02875     /* Do we have a "small" integer? */
02876 
02877     if (be >= 0 && k <= Int_max) {
02878         /* Yes. */
02879         ds = tens[k];
02880         if (ndigits < 0 && ilim <= 0) {
02881             S = mhi = 0;
02882             if (ilim < 0 || dval(d) <= 5*ds)
02883                 goto no_digits;
02884             goto one_digit;
02885             }
02886         for(i = 1;; i++, dval(d) *= 10.) {
02887             L = (Long)(dval(d) / ds);
02888             dval(d) -= L*ds;
02889 #ifdef Check_FLT_ROUNDS
02890             /* If FLT_ROUNDS == 2, L will usually be high by 1 */
02891             if (dval(d) < 0) {
02892                 L--;
02893                 dval(d) += ds;
02894                 }
02895 #endif
02896             *s++ = '0' + (int)L;
02897             if (!dval(d)) {
02898 #ifdef SET_INEXACT
02899                 inexact = 0;
02900 #endif
02901                 break;
02902                 }
02903             if (i == ilim) {
02904 #ifdef Honor_FLT_ROUNDS
02905                 if (mode > 1)
02906                 switch(rounding) {
02907                   case 0: goto ret1;
02908                   case 2: goto bump_up;
02909                   }
02910 #endif
02911                 dval(d) += dval(d);
02912                 if (dval(d) > ds || dval(d) == ds && L & 1) {
02913  bump_up:
02914                     while(*--s == '9')
02915                         if (s == s0) {
02916                             k++;
02917                             *s = '0';
02918                             break;
02919                             }
02920                     ++*s++;
02921                     }
02922                 break;
02923                 }
02924             }
02925         goto ret1;
02926         }
02927 
02928     m2 = b2;
02929     m5 = b5;
02930     mhi = mlo = 0;
02931     if (leftright) {
02932         i =
02933 #ifndef Sudden_Underflow
02934             denorm ? be + (Bias + (P-1) - 1 + 1) :
02935 #endif
02936 #ifdef IBM
02937             1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
02938 #else
02939             1 + P - bbits;
02940 #endif
02941         b2 += i;
02942         s2 += i;
02943         mhi = i2b(1);
02944         }
02945     if (m2 > 0 && s2 > 0) {
02946         i = m2 < s2 ? m2 : s2;
02947         b2 -= i;
02948         m2 -= i;
02949         s2 -= i;
02950         }
02951     if (b5 > 0) {
02952         if (leftright) {
02953             if (m5 > 0) {
02954                 mhi = pow5mult(mhi, m5);
02955                 b1 = mult(mhi, b);
02956                 Bfree(b);
02957                 b = b1;
02958                 }
02959             if ((j = b5 - m5))
02960                 b = pow5mult(b, j);
02961             }
02962         else
02963             b = pow5mult(b, b5);
02964         }
02965     S = i2b(1);
02966     if (s5 > 0)
02967         S = pow5mult(S, s5);
02968 
02969     /* Check for special case that d is a normalized power of 2. */
02970 
02971     spec_case = 0;
02972     if ((mode < 2 || leftright)
02973 #ifdef Honor_FLT_ROUNDS
02974             && rounding == 1
02975 #endif
02976                 ) {
02977         if (!word1(d) && !(word0(d) & Bndry_mask)
02978 #ifndef Sudden_Underflow
02979          && word0(d) & (Exp_mask & ~Exp_msk1)
02980 #endif
02981                 ) {
02982             /* The special case */
02983             b2 += Log2P;
02984             s2 += Log2P;
02985             spec_case = 1;
02986             }
02987         }
02988 
02989     /* Arrange for convenient computation of quotients:
02990      * shift left if necessary so divisor has 4 leading 0 bits.
02991      *
02992      * Perhaps we should just compute leading 28 bits of S once
02993      * and for all and pass them and a shift to quorem, so it
02994      * can do shifts and ors to compute the numerator for q.
02995      */
02996 #ifdef Pack_32
02997     if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
02998         i = 32 - i;
02999 #else
03000     if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
03001         i = 16 - i;
03002 #endif
03003     if (i > 4) {
03004         i -= 4;
03005         b2 += i;
03006         m2 += i;
03007         s2 += i;
03008         }
03009     else if (i < 4) {
03010         i += 28;
03011         b2 += i;
03012         m2 += i;
03013         s2 += i;
03014         }
03015     if (b2 > 0)
03016         b = lshift(b, b2);
03017     if (s2 > 0)
03018         S = lshift(S, s2);
03019     if (k_check) {
03020         if (cmp(b,S) < 0) {
03021             k--;
03022             b = multadd(b, 10, 0);  /* we botched the k estimate */
03023             if (leftright)
03024                 mhi = multadd(mhi, 10, 0);
03025             ilim = ilim1;
03026             }
03027         }
03028     if (ilim <= 0 && (mode == 3 || mode == 5)) {
03029         if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
03030             /* no digits, fcvt style */
03031  no_digits:
03032             k = -1 - ndigits;
03033             goto ret;
03034             }
03035  one_digit:
03036         *s++ = '1';
03037         k++;
03038         goto ret;
03039         }
03040     if (leftright) {
03041         if (m2 > 0)
03042             mhi = lshift(mhi, m2);
03043 
03044         /* Compute mlo -- check for special case
03045          * that d is a normalized power of 2.
03046          */
03047 
03048         mlo = mhi;
03049         if (spec_case) {
03050             mhi = Balloc(mhi->k);
03051             Bcopy(mhi, mlo);
03052             mhi = lshift(mhi, Log2P);
03053             }
03054 
03055         for(i = 1;;i++) {
03056             dig = quorem(b,S) + '0';
03057             /* Do we yet have the shortest decimal string
03058              * that will round to d?
03059              */
03060             j = cmp(b, mlo);
03061             delta = diff(S, mhi);
03062             j1 = delta->sign ? 1 : cmp(b, delta);
03063             Bfree(delta);
03064 #ifndef ROUND_BIASED
03065             if (j1 == 0 && mode != 1 && !(word1(d) & 1)
03066 #ifdef Honor_FLT_ROUNDS
03067                 && rounding >= 1
03068 #endif
03069                                    ) {
03070                 if (dig == '9')
03071                     goto round_9_up;
03072                 if (j > 0)
03073                     dig++;
03074 #ifdef SET_INEXACT
03075                 else if (!b->x[0] && b->wds <= 1)
03076                     inexact = 0;
03077 #endif
03078                 *s++ = dig;
03079                 goto ret;
03080                 }
03081 #endif
03082             if (j < 0 || j == 0 && mode != 1
03083 #ifndef ROUND_BIASED
03084                             && !(word1(d) & 1)
03085 #endif
03086                     ) {
03087                 if (!b->x[0] && b->wds <= 1) {
03088 #ifdef SET_INEXACT
03089                     inexact = 0;
03090 #endif
03091                     goto accept_dig;
03092                     }
03093 #ifdef Honor_FLT_ROUNDS
03094                 if (mode > 1)
03095                  switch(rounding) {
03096                   case 0: goto accept_dig;
03097                   case 2: goto keep_dig;
03098                   }
03099 #endif /*Honor_FLT_ROUNDS*/
03100                 if (j1 > 0) {
03101                     b = lshift(b, 1);
03102                     j1 = cmp(b, S);
03103                     if ((j1 > 0 || j1 == 0 && dig & 1)
03104                     && dig++ == '9')
03105                         goto round_9_up;
03106                     }
03107  accept_dig:
03108                 *s++ = dig;
03109                 goto ret;
03110                 }
03111             if (j1 > 0) {
03112 #ifdef Honor_FLT_ROUNDS
03113                 if (!rounding)
03114                     goto accept_dig;
03115 #endif
03116                 if (dig == '9') { /* possible if i == 1 */
03117  round_9_up:
03118                     *s++ = '9';
03119                     goto roundoff;
03120                     }
03121                 *s++ = dig + 1;
03122                 goto ret;
03123                 }
03124 #ifdef Honor_FLT_ROUNDS
03125  keep_dig:
03126 #endif
03127             *s++ = dig;
03128             if (i == ilim)
03129                 break;
03130             b = multadd(b, 10, 0);
03131             if (mlo == mhi)
03132                 mlo = mhi = multadd(mhi, 10, 0);
03133             else {
03134                 mlo = multadd(mlo, 10, 0);
03135                 mhi = multadd(mhi, 10, 0);
03136                 }
03137             }
03138         }
03139     else
03140         for(i = 1;; i++) {
03141             *s++ = dig = quorem(b,S) + '0';
03142             if (!b->x[0] && b->wds <= 1) {
03143 #ifdef SET_INEXACT
03144                 inexact = 0;
03145 #endif
03146                 goto ret;
03147                 }
03148             if (i >= ilim)
03149                 break;
03150             b = multadd(b, 10, 0);
03151             }
03152 
03153     /* Round off last digit */
03154 
03155 #ifdef Honor_FLT_ROUNDS
03156     switch(rounding) {
03157       case 0: goto trimzeros;
03158       case 2: goto roundoff;
03159       }
03160 #endif
03161     b = lshift(b, 1);
03162     j = cmp(b, S);
03163     if (j > 0 || j == 0 && dig & 1) {
03164  roundoff:
03165         while(*--s == '9')
03166             if (s == s0) {
03167                 k++;
03168                 *s++ = '1';
03169                 goto ret;
03170                 }
03171         ++*s++;
03172         }
03173     else {
03174 #ifdef Honor_FLT_ROUNDS
03175 trimzeros:
03176 #endif
03177         while(*--s == '0')
03178             ;
03179         s++;
03180         }
03181  ret:
03182     Bfree(S);
03183     if (mhi) {
03184         if (mlo && mlo != mhi)
03185             Bfree(mlo);
03186         Bfree(mhi);
03187         }
03188  ret1:
03189 #ifdef SET_INEXACT
03190     if (inexact) {
03191         if (!oldinexact) {
03192             word0(d) = Exp_1 + (70 << Exp_shift);
03193             word1(d) = 0;
03194             dval(d) += 1.;
03195             }
03196         }
03197     else if (!oldinexact)
03198         clear_inexact();
03199 #endif
03200     Bfree(b);
03201     *s = 0;
03202     *decpt = k + 1;
03203     if (rve)
03204         *rve = s;
03205     return s0;
03206     }
03207 #ifdef __cplusplus
03208 }
03209 #endif
KDE Home | KDE Accessibility Home | Description of Access Keys