• Skip to content
  • Skip to link menu
Trinity API Reference
  • Trinity API Reference
  • kjs
 

kjs

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

kjs

Skip menu "kjs"
  • Main Page
  • Class Hierarchy
  • Alphabetical List
  • Class List
  • File List
  • Class Members
  • Related Pages

kjs

Skip menu "kjs"
  • arts
  • dcop
  • dnssd
  • interfaces
  •   kspeech
  •     interface
  •     library
  •   tdetexteditor
  • kate
  • kded
  • kdoctools
  • kimgio
  • kjs
  • libtdemid
  • libtdescreensaver
  • tdeabc
  • tdecmshell
  • tdecore
  • tdefx
  • tdehtml
  • tdeinit
  • tdeio
  •   bookmarks
  •   httpfilter
  •   kpasswdserver
  •   kssl
  •   tdefile
  •   tdeio
  •   tdeioexec
  • tdeioslave
  •   http
  • tdemdi
  •   tdemdi
  • tdenewstuff
  • tdeparts
  • tdeprint
  • tderandr
  • tderesources
  • tdespell2
  • tdesu
  • tdeui
  • tdeunittest
  • tdeutils
  • tdewallet
Generated for kjs by doxygen 1.6.3
This website is maintained by Timothy Pearson.