mono/utils/strtod.c

   1 /****************************************************************
   2  *
   3  * The author of this software is David M. Gay.
   4  *
   5  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
   6  *
   7  * Permission to use, copy, modify, and distribute this software for any
   8  * purpose without fee is hereby granted, provided that this entire notice
   9  * is included in all copies of any software which is or includes a copy
  10  * or modification of this software and in all copies of the supporting
  11  * documentation for such software.
  12  *
  13  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
  15  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17  *
  18  ***************************************************************/
  19 #include "strtod.h"
  20 #include <glib.h>
  21 #define freedtoa __freedtoa
  22 #define dtoa __dtoa
  23
  24 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  25  * with " at " changed at "@" and " dot " changed to ".").      */
  26
  27 /* On a machine with IEEE extended-precision registers, it is
  28  * necessary to specify double-precision (53-bit) rounding precision
  29  * before invoking strtod or dtoa.  If the machine uses (the equivalent
  30  * of) Intel 80x87 arithmetic, the call
  31  *      _control87(PC_53, MCW_PC);
  32  * does this with many compilers.  Whether this or another call is
  33  * appropriate depends on the compiler; for this to work, it may be
  34  * necessary to #include "float.h" or another system-dependent header
  35  * file.
  36  */
  37
  38 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  39  *
  40  * This strtod returns a nearest machine number to the input decimal
  41  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  42  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  43  * biased rounding (add half and chop).
  44  *
  45  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  46  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  47  *
  48  * Modifications:
  49  *
  50  *      1. We only require IEEE, IBM, or VAX double-precision
  51  *              arithmetic (not IEEE double-extended).
  52  *      2. We get by with floating-point arithmetic in a case that
  53  *              Clinger missed -- when we're computing d * 10^n
  54  *              for a small integer d and the integer n is not too
  55  *              much larger than 22 (the maximum integer k for which
  56  *              we can represent 10^k exactly), we may be able to
  57  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  58  *      3. Rather than a bit-at-a-time adjustment of the binary
  59  *              result in the hard case, we use floating-point
  60  *              arithmetic to determine the adjustment to within
  61  *              one bit; only in really hard cases do we need to
  62  *              compute a second residual.
  63  *      4. Because of 3., we don't need a large table of powers of 10
  64  *              for ten-to-e (just some small tables, e.g. of 10^k
  65  *              for 0 <= k <= 22).
  66  */
  67
  68 /*
  69  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  70  *      significant byte has the lowest address.
  71  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  72  *      significant byte has the lowest address.
  73  * #define Long int on machines with 32-bit ints and 64-bit longs.
  74  * #define IBM for IBM mainframe-style floating-point arithmetic.
  75  * #define VAX for VAX-style floating-point arithmetic (D_floating).
  76  * #define No_leftright to omit left-right logic in fast floating-point
  77  *      computation of dtoa.
  78  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  79  *      and strtod and dtoa should round accordingly.
  80  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  81  *      and Honor_FLT_ROUNDS is not #defined.
  82  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  83  *      that use extended-precision instructions to compute rounded
  84  *      products and quotients) with IBM.
  85  * #define ROUND_BIASED for IEEE-format with biased rounding.
  86  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  87  *      products but inaccurate quotients, e.g., for Intel i860.
  88  * #define NO_LONG_LONG on machines that do not have a "long long"
  89  *      integer type (of >= 64 bits).  On such machines, you can
  90  *      #define Just_16 to store 16 bits per 32-bit Long when doing
  91  *      high-precision integer arithmetic.  Whether this speeds things
  92  *      up or slows things down depends on the machine and the number
  93  *      being converted.  If long long is available and the name is
  94  *      something other than "long long", #define Llong to be the name,
  95  *      and if "unsigned Llong" does not work as an unsigned version of
  96  *      Llong, #define #ULLong to be the corresponding unsigned type.
  97  * #define KR_headers for old-style C function headers.
  98  * #define Bad_float_h if your system lacks a float.h or if it does not
  99  *      define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
 100  *      FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
 101  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
 102  *      if memory is available and otherwise does something you deem
 103  *      appropriate.  If MALLOC is undefined, malloc will be invoked
 104  *      directly -- and assumed always to succeed.
 105  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
 106  *      memory allocations from a private pool of memory when possible.
 107  *      When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
 108  *      unless #defined to be a different length.  This default length
 109  *      suffices to get rid of MALLOC calls except for unusual cases,
 110  *      such as decimal-to-binary conversion of a very long string of
 111  *      digits.  The longest string dtoa can return is about 751 bytes
 112  *      long.  For conversions by strtod of strings of 800 digits and
 113  *      all dtoa conversions in single-threaded executions with 8-byte
 114  *      pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
 115  *      pointers, PRIVATE_MEM >= 7112 appears adequate.
 116  * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
 117  *      Infinity and NaN (case insensitively).  On some systems (e.g.,
 118  *      some HP systems), it may be necessary to #define NAN_WORD0
 119  *      appropriately -- to the most significant word of a quiet NaN.
 120  *      (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
 121  *      When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
 122  *      strtod also accepts (case insensitively) strings of the form
 123  *      NaN(x), where x is a string of hexadecimal digits and spaces;
 124  *      if there is only one string of hexadecimal digits, it is taken
 125  *      for the 52 fraction bits of the resulting NaN; if there are two
 126  *      or more strings of hex digits, the first is for the high 20 bits,
 127  *      the second and subsequent for the low 32 bits, with intervening
 128  *      white space ignored; but if this results in none of the 52
 129  *      fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
 130  *      and NAN_WORD1 are used instead.
 131  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
 132  *      multiple threads.  In this case, you must provide (or suitably
 133  *      #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
 134  *      by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
 135  *      in pow5mult, ensures lazy evaluation of only one copy of high
 136  *      powers of 5; omitting this lock would introduce a small
 137  *      probability of wasting memory, but would otherwise be harmless.)
 138  *      You must also invoke freedtoa(s) to free the value s returned by
 139  *      dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
 140  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
 141  *      avoids underflows on inputs whose result does not underflow.
 142  *      If you #define NO_IEEE_Scale on a machine that uses IEEE-format
 143  *      floating-point numbers and flushes underflows to zero rather
 144  *      than implementing gradual underflow, then you must also #define
 145  *      Sudden_Underflow.
 146  * #define YES_ALIAS to permit aliasing certain double values with
 147  *      arrays of ULongs.  This leads to slightly better code with
 148  *      some compilers and was always used prior to 19990916, but it
 149  *      is not strictly legal and can cause trouble with aggressively
 150  *      optimizing compilers (e.g., gcc 2.95.1 under -O2).
 151  * #define USE_LOCALE to use the current locale's decimal_point value.
 152  * #define SET_INEXACT if IEEE arithmetic is being used and extra
 153  *      computation should be done to set the inexact flag when the
 154  *      result is inexact and avoid setting inexact when the result
 155  *      is exact.  In this case, dtoa.c must be compiled in
 156  *      an environment, perhaps provided by #include "dtoa.c" in a
 157  *      suitable wrapper, that defines two functions,
 158  *              int get_inexact(void);
 159  *              void clear_inexact(void);
 160  *      such that get_inexact() returns a nonzero value if the
 161  *      inexact bit is already set, and clear_inexact() sets the
 162  *      inexact bit to 0.  When SET_INEXACT is #defined, strtod
 163  *      also does extra computations to set the underflow and overflow
 164  *      flags when appropriate (i.e., when the result is tiny and
 165  *      inexact or when it is a numeric value rounded to +-infinity).
 166  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
 167  *      the result overflows to +-Infinity or underflows to 0.
 168  */
 169 #if defined(i386) || defined(mips) && defined(MIPSEL) || defined (__arm__)
 170
 171 #   define IEEE_8087
 172
 173 #elif defined(__x86_64__) || defined(__alpha__)
 174
 175 #   define IEEE_8087
 176
 177 #elif defined(__ia64)
 178
 179 #   ifdef __hpux
 180 #       define IEEE_MC68k
 181 #   else
 182 #       define IEEE_8087
 183 #   endif
 184
 185 #elif defined(__hppa)
 186
 187 #   define IEEE_MC68k
 188
 189 #else
 190 #define IEEE_MC68k
 191 #endif
 192
 193 #define Long gint32
 194 #define ULong guint32
 195
 196 #ifdef DEBUG
 197 #include "stdio.h"
 198 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
 199 #endif
 200
 201 #include "stdlib.h"
 202 #include "string.h"
 203
 204 #undef USE_LOCALE
 205 #ifdef USE_LOCALE
 206 #include "locale.h"
 207 #endif
 208
 209 #ifdef MALLOC
 210 #ifdef KR_headers
 211 extern char *MALLOC();
 212 #else
 213 extern void *MALLOC(size_t);
 214 #endif
 215 #else
 216 #define MALLOC malloc
 217 #endif
 218
 219 #ifndef Omit_Private_Memory
 220 #ifndef PRIVATE_MEM
 221 #define PRIVATE_MEM 2304
 222 #endif
 223 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
 224 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
 225 #endif
 226
 227 #undef IEEE_Arith
 228 #undef Avoid_Underflow
 229 #ifdef IEEE_MC68k
 230 #define IEEE_Arith
 231 #endif
 232 #ifdef IEEE_8087
 233 #define IEEE_Arith
 234 #endif
 235
 236 #include "errno.h"
 237
 238 #ifdef Bad_float_h
 239
 240 #ifdef IEEE_Arith
 241 #define DBL_DIG 15
 242 #define DBL_MAX_10_EXP 308
 243 #define DBL_MAX_EXP 1024
 244 #define FLT_RADIX 2
 245 #endif /*IEEE_Arith*/
 246
 247 #ifdef IBM
 248 #define DBL_DIG 16
 249 #define DBL_MAX_10_EXP 75
 250 #define DBL_MAX_EXP 63
 251 #define FLT_RADIX 16
 252 #define DBL_MAX 7.2370055773322621e+75
 253 #endif
 254
 255 #ifdef VAX
 256 #define DBL_DIG 16
 257 #define DBL_MAX_10_EXP 38
 258 #define DBL_MAX_EXP 127
 259 #define FLT_RADIX 2
 260 #define DBL_MAX 1.7014118346046923e+38
 261 #endif
 262
 263 #ifndef LONG_MAX
 264 #define LONG_MAX 2147483647
 265 #endif
 266
 267 #else /* ifndef Bad_float_h */
 268 #include "float.h"
 269 #endif /* Bad_float_h */
 270
 271 #ifndef __MATH_H__
 272 #include "math.h"
 273 #endif
 274
 275 #ifdef __cplusplus
 276 extern "C" {
 277 #endif
 278
 279 #ifndef CONST
 280 #ifdef KR_headers
 281 #define CONST /* blank */
 282 #else
 283 #define CONST const
 284 #endif
 285 #endif
 286
 287 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
 288 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
 289 #endif
 290
 291 typedef union { double d; ULong L[2]; } U;
 292
 293 #ifdef YES_ALIAS
 294 #define dval(x) x
 295 #ifdef IEEE_8087
 296 #define word0(x) ((ULong *)&x)[1]
 297 #define word1(x) ((ULong *)&x)[0]
 298 #else
 299 #define word0(x) ((ULong *)&x)[0]
 300 #define word1(x) ((ULong *)&x)[1]
 301 #endif
 302 #else
 303 #ifdef IEEE_8087
 304 #define word0(x) ((U*)&x)->L[1]
 305 #define word1(x) ((U*)&x)->L[0]
 306 #else
 307 #define word0(x) ((U*)&x)->L[0]
 308 #define word1(x) ((U*)&x)->L[1]
 309 #endif
 310 #define dval(x) ((U*)&x)->d
 311 #endif
 312
 313 /* The following definition of Storeinc is appropriate for MIPS processors.
 314  * An alternative that might be better on some machines is
 315  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
 316  */
 317 #if defined(IEEE_8087) + defined(VAX)
 318 #define Storeinc(a,b,c) do { (((unsigned short *)a)[1] = (unsigned short)b, \
 319 ((unsigned short *)a)[0] = (unsigned short)c, a++) } while (0)
 320 #else
 321 #define Storeinc(a,b,c) do { (((unsigned short *)a)[0] = (unsigned short)b, \
 322 ((unsigned short *)a)[1] = (unsigned short)c, a++) } while (0)
 323 #endif
 324
 325 /* #define P DBL_MANT_DIG */
 326 /* Ten_pmax = floor(P*log(2)/log(5)) */
 327 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
 328 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
 329 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
 330
 331 #ifdef IEEE_Arith
 332 #define Exp_shift  20
 333 #define Exp_shift1 20
 334 #define Exp_msk1    0x100000
 335 #define Exp_msk11   0x100000
 336 #define Exp_mask  0x7ff00000
 337 #define P 53
 338 #define Bias 1023
 339 #define Emin (-1022)
 340 #define Exp_1  0x3ff00000
 341 #define Exp_11 0x3ff00000
 342 #define Ebits 11
 343 #define Frac_mask  0xfffff
 344 #define Frac_mask1 0xfffff
 345 #define Ten_pmax 22
 346 #define Bletch 0x10
 347 #define Bndry_mask  0xfffff
 348 #define Bndry_mask1 0xfffff
 349 #define LSB 1
 350 #define Sign_bit 0x80000000
 351 #define Log2P 1
 352 #define Tiny0 0
 353 #define Tiny1 1
 354 #define Quick_max 14
 355 #define Int_max 14
 356 #ifndef NO_IEEE_Scale
 357 #define Avoid_Underflow
 358 #ifdef Flush_Denorm     /* debugging option */
 359 #undef Sudden_Underflow
 360 #endif
 361 #endif
 362
 363 #ifndef Flt_Rounds
 364 #ifdef FLT_ROUNDS
 365 #define Flt_Rounds FLT_ROUNDS
 366 #else
 367 #define Flt_Rounds 1
 368 #endif
 369 #endif /*Flt_Rounds*/
 370
 371 #ifdef Honor_FLT_ROUNDS
 372 #define Rounding rounding
 373 #undef Check_FLT_ROUNDS
 374 #define Check_FLT_ROUNDS
 375 #else
 376 #define Rounding Flt_Rounds
 377 #endif
 378
 379 #else /* ifndef IEEE_Arith */
 380 #undef Check_FLT_ROUNDS
 381 #undef Honor_FLT_ROUNDS
 382 #undef SET_INEXACT
 383 #undef  Sudden_Underflow
 384 #define Sudden_Underflow
 385 #ifdef IBM
 386 #undef Flt_Rounds
 387 #define Flt_Rounds 0
 388 #define Exp_shift  24
 389 #define Exp_shift1 24
 390 #define Exp_msk1   0x1000000
 391 #define Exp_msk11  0x1000000
 392 #define Exp_mask  0x7f000000
 393 #define P 14
 394 #define Bias 65
 395 #define Exp_1  0x41000000
 396 #define Exp_11 0x41000000
 397 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
 398 #define Frac_mask  0xffffff
 399 #define Frac_mask1 0xffffff
 400 #define Bletch 4
 401 #define Ten_pmax 22
 402 #define Bndry_mask  0xefffff
 403 #define Bndry_mask1 0xffffff
 404 #define LSB 1
 405 #define Sign_bit 0x80000000
 406 #define Log2P 4
 407 #define Tiny0 0x100000
 408 #define Tiny1 0
 409 #define Quick_max 14
 410 #define Int_max 15
 411 #else /* VAX */
 412 #undef Flt_Rounds
 413 #define Flt_Rounds 1
 414 #define Exp_shift  23
 415 #define Exp_shift1 7
 416 #define Exp_msk1    0x80
 417 #define Exp_msk11   0x800000
 418 #define Exp_mask  0x7f80
 419 #define P 56
 420 #define Bias 129
 421 #define Exp_1  0x40800000
 422 #define Exp_11 0x4080
 423 #define Ebits 8
 424 #define Frac_mask  0x7fffff
 425 #define Frac_mask1 0xffff007f
 426 #define Ten_pmax 24
 427 #define Bletch 2
 428 #define Bndry_mask  0xffff007f
 429 #define Bndry_mask1 0xffff007f
 430 #define LSB 0x10000
 431 #define Sign_bit 0x8000
 432 #define Log2P 1
 433 #define Tiny0 0x80
 434 #define Tiny1 0
 435 #define Quick_max 15
 436 #define Int_max 15
 437 #endif /* IBM, VAX */
 438 #endif /* IEEE_Arith */
 439
 440 #ifndef IEEE_Arith
 441 #define ROUND_BIASED
 442 #endif
 443
 444 #ifdef RND_PRODQUOT
 445 #define rounded_product(a,b) a = rnd_prod(a, b)
 446 #define rounded_quotient(a,b) a = rnd_quot(a, b)
 447 #ifdef KR_headers
 448 extern double rnd_prod(), rnd_quot();
 449 #else
 450 extern double rnd_prod(double, double), rnd_quot(double, double);
 451 #endif
 452 #else
 453 #define rounded_product(a,b) a *= b
 454 #define rounded_quotient(a,b) a /= b
 455 #endif
 456
 457 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
 458 #define Big1 0xffffffff
 459
 460 #ifndef Pack_32
 461 #define Pack_32
 462 #endif
 463
 464 #ifdef KR_headers
 465 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
 466 #else
 467 #define FFFFFFFF 0xffffffffUL
 468 #endif
 469
 470 #ifdef NO_LONG_LONG
 471 #undef ULLong
 472 #ifdef Just_16
 473 #undef Pack_32
 474 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
 475  * This makes some inner loops simpler and sometimes saves work
 476  * during multiplications, but it often seems to make things slightly
 477  * slower.  Hence the default is now to store 32 bits per Long.
 478  */
 479 #endif
 480 #else   /* long long available */
 481 #ifndef Llong
 482 #define Llong long long
 483 #endif
 484 #ifndef ULLong
 485 #define ULLong unsigned Llong
 486 #endif
 487 #endif /* NO_LONG_LONG */
 488
 489 #ifndef MULTIPLE_THREADS
 490 #define ACQUIRE_DTOA_LOCK(n)    /*nothing*/
 491 #define FREE_DTOA_LOCK(n)       /*nothing*/
 492 #endif
 493
 494 #define Kmax 15
 495
 496 #ifdef __cplusplus
 497 extern "C" double strtod(const char *s00, char **se);
 498 extern "C" char *dtoa(double d, int mode, int ndigits,
 499                         int *decpt, int *sign, char **rve);
 500 #endif
 501
 502  struct
 503 Bigint {
 504         struct Bigint *next;
 505         int k, maxwds, sign, wds;
 506         ULong x[1];
 507         };
 508
 509  typedef struct Bigint Bigint;
 510
 511  static Bigint *freelist[Kmax+1];
 512
 513  static Bigint *
 514 Balloc
 515 #ifdef KR_headers
 516         (k) int k;
 517 #else
 518         (int k)
 519 #endif
 520 {
 521         int x;
 522         Bigint *rv;
 523 #ifndef Omit_Private_Memory
 524         unsigned int len;
 525 #endif
 526
 527         ACQUIRE_DTOA_LOCK(0);
 528         if ((rv = freelist[k])) {
 529                 freelist[k] = rv->next;
 530                 }
 531         else {
 532                 x = 1 << k;
 533 #ifdef Omit_Private_Memory
 534                 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
 535 #else
 536                 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
 537                         /sizeof(double);
 538                 if (pmem_next - private_mem + len <= PRIVATE_mem) {
 539                         rv = (Bigint*)pmem_next;
 540                         pmem_next += len;
 541                         }
 542                 else
 543                         rv = (Bigint*)MALLOC(len*sizeof(double));
 544 #endif
 545                 rv->k = k;
 546                 rv->maxwds = x;
 547                 }
 548         FREE_DTOA_LOCK(0);
 549         rv->sign = rv->wds = 0;
 550         return rv;
 551         }
 552
 553  static void
 554 Bfree
 555 #ifdef KR_headers
 556         (v) Bigint *v;
 557 #else
 558         (Bigint *v)
 559 #endif
 560 {
 561         if (v) {
 562                 ACQUIRE_DTOA_LOCK(0);
 563                 v->next = freelist[v->k];
 564                 freelist[v->k] = v;
 565                 FREE_DTOA_LOCK(0);
 566                 }
 567         }
 568
 569 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
 570 y->wds*sizeof(Long) + 2*sizeof(int))
 571
 572  static Bigint *
 573 multadd
 574 #ifdef KR_headers
 575         (b, m, a) Bigint *b; int m, a;
 576 #else
 577         (Bigint *b, int m, int a)       /* multiply by m and add a */
 578 #endif
 579 {
 580         int i, wds;
 581 #ifdef ULLong
 582         ULong *x;
 583         ULLong carry, y;
 584 #else
 585         ULong carry, *x, y;
 586 #ifdef Pack_32
 587         ULong xi, z;
 588 #endif
 589 #endif
 590         Bigint *b1;
 591
 592         wds = b->wds;
 593         x = b->x;
 594         i = 0;
 595         carry = a;
 596         do {
 597 #ifdef ULLong
 598                 y = *x * (ULLong)m + carry;
 599                 carry = y >> 32;
 600                 *x++ = y & FFFFFFFF;
 601 #else
 602 #ifdef Pack_32
 603                 xi = *x;
 604                 y = (xi & 0xffff) * m + carry;
 605                 z = (xi >> 16) * m + (y >> 16);
 606                 carry = z >> 16;
 607                 *x++ = (z << 16) + (y & 0xffff);
 608 #else
 609                 y = *x * m + carry;
 610                 carry = y >> 16;
 611                 *x++ = y & 0xffff;
 612 #endif
 613 #endif
 614                 }
 615                 while(++i < wds);
 616         if (carry) {
 617                 if (wds >= b->maxwds) {
 618                         b1 = Balloc(b->k+1);
 619                         Bcopy(b1, b);
 620                         Bfree(b);
 621                         b = b1;
 622                         }
 623                 b->x[wds++] = carry;
 624                 b->wds = wds;
 625                 }
 626         return b;
 627         }
 628
 629  static Bigint *
 630 s2b
 631 #ifdef KR_headers
 632         (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
 633 #else
 634         (CONST char *s, int nd0, int nd, ULong y9)
 635 #endif
 636 {
 637         Bigint *b;
 638         int i, k;
 639         Long x, y;
 640
 641         x = (nd + 8) / 9;
 642         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
 643 #ifdef Pack_32
 644         b = Balloc(k);
 645         b->x[0] = y9;
 646         b->wds = 1;
 647 #else
 648         b = Balloc(k+1);
 649         b->x[0] = y9 & 0xffff;
 650         b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
 651 #endif
 652
 653         i = 9;
 654         if (9 < nd0) {
 655                 s += 9;
 656                 do b = multadd(b, 10, *s++ - '0');
 657                         while(++i < nd0);
 658                 s++;
 659                 }
 660         else
 661                 s += 10;
 662         for(; i < nd; i++)
 663                 b = multadd(b, 10, *s++ - '0');
 664         return b;
 665         }
 666
 667  static int
 668 hi0bits
 669 #ifdef KR_headers
 670         (x) register ULong x;
 671 #else
 672         (register ULong x)
 673 #endif
 674 {
 675         register int k = 0;
 676
 677         if (!(x & 0xffff0000)) {
 678                 k = 16;
 679                 x <<= 16;
 680                 }
 681         if (!(x & 0xff000000)) {
 682                 k += 8;
 683                 x <<= 8;
 684                 }
 685         if (!(x & 0xf0000000)) {
 686                 k += 4;
 687                 x <<= 4;
 688                 }
 689         if (!(x & 0xc0000000)) {
 690                 k += 2;
 691                 x <<= 2;
 692                 }
 693         if (!(x & 0x80000000)) {
 694                 k++;
 695                 if (!(x & 0x40000000))
 696                         return 32;
 697                 }
 698         return k;
 699         }
 700
 701  static int
 702 lo0bits
 703 #ifdef KR_headers
 704         (y) ULong *y;
 705 #else
 706         (ULong *y)
 707 #endif
 708 {
 709         register int k;
 710         register ULong x = *y;
 711
 712         if (x & 7) {
 713                 if (x & 1)
 714                         return 0;
 715                 if (x & 2) {
 716                         *y = x >> 1;
 717                         return 1;
 718                         }
 719                 *y = x >> 2;
 720                 return 2;
 721                 }
 722         k = 0;
 723         if (!(x & 0xffff)) {
 724                 k = 16;
 725                 x >>= 16;
 726                 }
 727         if (!(x & 0xff)) {
 728                 k += 8;
 729                 x >>= 8;
 730                 }
 731         if (!(x & 0xf)) {
 732                 k += 4;
 733                 x >>= 4;
 734                 }
 735         if (!(x & 0x3)) {
 736                 k += 2;
 737                 x >>= 2;
 738                 }
 739         if (!(x & 1)) {
 740                 k++;
 741                 x >>= 1;
 742                 if (!x)
 743                         return 32;
 744                 }
 745         *y = x;
 746         return k;
 747         }
 748
 749  static Bigint *
 750 i2b
 751 #ifdef KR_headers
 752         (i) int i;
 753 #else
 754         (int i)
 755 #endif
 756 {
 757         Bigint *b;
 758
 759         b = Balloc(1);
 760         b->x[0] = i;
 761         b->wds = 1;
 762         return b;
 763         }
 764
 765  static Bigint *
 766 mult
 767 #ifdef KR_headers
 768         (a, b) Bigint *a, *b;
 769 #else
 770         (Bigint *a, Bigint *b)
 771 #endif
 772 {
 773         Bigint *c;
 774         int k, wa, wb, wc;
 775         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 776         ULong y;
 777 #ifdef ULLong
 778         ULLong carry, z;
 779 #else
 780         ULong carry, z;
 781 #ifdef Pack_32
 782         ULong z2;
 783 #endif
 784 #endif
 785
 786         if (a->wds < b->wds) {
 787                 c = a;
 788                 a = b;
 789                 b = c;
 790                 }
 791         k = a->k;
 792         wa = a->wds;
 793         wb = b->wds;
 794         wc = wa + wb;
 795         if (wc > a->maxwds)
 796                 k++;
 797         c = Balloc(k);
 798         for(x = c->x, xa = x + wc; x < xa; x++)
 799                 *x = 0;
 800         xa = a->x;
 801         xae = xa + wa;
 802         xb = b->x;
 803         xbe = xb + wb;
 804         xc0 = c->x;
 805 #ifdef ULLong
 806         for(; xb < xbe; xc0++) {
 807                 if ((y = *xb++)) {
 808                         x = xa;
 809                         xc = xc0;
 810                         carry = 0;
 811                         do {
 812                                 z = *x++ * (ULLong)y + *xc + carry;
 813                                 carry = z >> 32;
 814                                 *xc++ = z & FFFFFFFF;
 815                                 }
 816                                 while(x < xae);
 817                         *xc = carry;
 818                         }
 819                 }
 820 #else
 821 #ifdef Pack_32
 822         for(; xb < xbe; xb++, xc0++) {
 823                 if (y = *xb & 0xffff) {
 824                         x = xa;
 825                         xc = xc0;
 826                         carry = 0;
 827                         do {
 828                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 829                                 carry = z >> 16;
 830                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 831                                 carry = z2 >> 16;
 832                                 Storeinc(xc, z2, z);
 833                                 }
 834                                 while(x < xae);
 835                         *xc = carry;
 836                         }
 837                 if (y = *xb >> 16) {
 838                         x = xa;
 839                         xc = xc0;
 840                         carry = 0;
 841                         z2 = *xc;
 842                         do {
 843                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 844                                 carry = z >> 16;
 845                                 Storeinc(xc, z, z2);
 846                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 847                                 carry = z2 >> 16;
 848                                 }
 849                                 while(x < xae);
 850                         *xc = z2;
 851                         }
 852                 }
 853 #else
 854         for(; xb < xbe; xc0++) {
 855                 if (y = *xb++) {
 856                         x = xa;
 857                         xc = xc0;
 858                         carry = 0;
 859                         do {
 860                                 z = *x++ * y + *xc + carry;
 861                                 carry = z >> 16;
 862                                 *xc++ = z & 0xffff;
 863                                 }
 864                                 while(x < xae);
 865                         *xc = carry;
 866                         }
 867                 }
 868 #endif
 869 #endif
 870         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
 871         c->wds = wc;
 872         return c;
 873         }
 874
 875  static Bigint *p5s;
 876
 877  static Bigint *
 878 pow5mult
 879 #ifdef KR_headers
 880         (b, k) Bigint *b; int k;
 881 #else
 882         (Bigint *b, int k)
 883 #endif
 884 {
 885         Bigint *b1, *p5, *p51;
 886         int i;
 887         static int p05[3] = { 5, 25, 125 };
 888
 889         if ((i = k & 3))
 890                 b = multadd(b, p05[i-1], 0);
 891
 892         if (!(k >>= 2))
 893                 return b;
 894         if (!(p5 = p5s)) {
 895                 /* first time */
 896 #ifdef MULTIPLE_THREADS
 897                 ACQUIRE_DTOA_LOCK(1);
 898                 if (!(p5 = p5s)) {
 899                         p5 = p5s = i2b(625);
 900                         p5->next = 0;
 901                         }
 902                 FREE_DTOA_LOCK(1);
 903 #else
 904                 p5 = p5s = i2b(625);
 905                 p5->next = 0;
 906 #endif
 907                 }
 908         for(;;) {
 909                 if (k & 1) {
 910                         b1 = mult(b, p5);
 911                         Bfree(b);
 912                         b = b1;
 913                         }
 914                 if (!(k >>= 1))
 915                         break;
 916                 if (!(p51 = p5->next)) {
 917 #ifdef MULTIPLE_THREADS
 918                         ACQUIRE_DTOA_LOCK(1);
 919                         if (!(p51 = p5->next)) {
 920                                 p51 = p5->next = mult(p5,p5);
 921                                 p51->next = 0;
 922                                 }
 923                         FREE_DTOA_LOCK(1);
 924 #else
 925                         p51 = p5->next = mult(p5,p5);
 926                         p51->next = 0;
 927 #endif
 928                         }
 929                 p5 = p51;
 930                 }
 931         return b;
 932         }
 933
 934  static Bigint *
 935 lshift
 936 #ifdef KR_headers
 937         (b, k) Bigint *b; int k;
 938 #else
 939         (Bigint *b, int k)
 940 #endif
 941 {
 942         int i, k1, n, n1;
 943         Bigint *b1;
 944         ULong *x, *x1, *xe, z;
 945
 946 #ifdef Pack_32
 947         n = k >> 5;
 948 #else
 949         n = k >> 4;
 950 #endif
 951         k1 = b->k;
 952         n1 = n + b->wds + 1;
 953         for(i = b->maxwds; n1 > i; i <<= 1)
 954                 k1++;
 955         b1 = Balloc(k1);
 956         x1 = b1->x;
 957         for(i = 0; i < n; i++)
 958                 *x1++ = 0;
 959         x = b->x;
 960         xe = x + b->wds;
 961 #ifdef Pack_32
 962         if (k &= 0x1f) {
 963                 k1 = 32 - k;
 964                 z = 0;
 965                 do {
 966                         *x1++ = *x << k | z;
 967                         z = *x++ >> k1;
 968                         }
 969                         while(x < xe);
 970                 if ((*x1 = z))
 971                         ++n1;
 972                 }
 973 #else
 974         if (k &= 0xf) {
 975                 k1 = 16 - k;
 976                 z = 0;
 977                 do {
 978                         *x1++ = *x << k  & 0xffff | z;
 979                         z = *x++ >> k1;
 980                         }
 981                         while(x < xe);
 982                 if (*x1 = z)
 983                         ++n1;
 984                 }
 985 #endif
 986         else do
 987                 *x1++ = *x++;
 988                 while(x < xe);
 989         b1->wds = n1 - 1;
 990         Bfree(b);
 991         return b1;
 992         }
 993
 994  static int
 995 cmp
 996 #ifdef KR_headers
 997         (a, b) Bigint *a, *b;
 998 #else
 999         (Bigint *a, Bigint *b)
1000 #endif
1001 {
1002         ULong *xa, *xa0, *xb, *xb0;
1003         int i, j;
1004
1005         i = a->wds;
1006         j = b->wds;
1007 #ifdef DEBUG
1008         if (i > 1 && !a->x[i-1])
1009                 Bug("cmp called with a->x[a->wds-1] == 0");
1010         if (j > 1 && !b->x[j-1])
1011                 Bug("cmp called with b->x[b->wds-1] == 0");
1012 #endif
1013         if (i -= j)
1014                 return i;
1015         xa0 = a->x;
1016         xa = xa0 + j;
1017         xb0 = b->x;
1018         xb = xb0 + j;
1019         for(;;) {
1020                 if (*--xa != *--xb)
1021                         return *xa < *xb ? -1 : 1;
1022                 if (xa <= xa0)
1023                         break;
1024                 }
1025         return 0;
1026         }
1027
1028  static Bigint *
1029 diff
1030 #ifdef KR_headers
1031         (a, b) Bigint *a, *b;
1032 #else
1033         (Bigint *a, Bigint *b)
1034 #endif
1035 {
1036         Bigint *c;
1037         int i, wa, wb;
1038         ULong *xa, *xae, *xb, *xbe, *xc;
1039 #ifdef ULLong
1040         ULLong borrow, y;
1041 #else
1042         ULong borrow, y;
1043 #ifdef Pack_32
1044         ULong z;
1045 #endif
1046 #endif
1047
1048         i = cmp(a,b);
1049         if (!i) {
1050                 c = Balloc(0);
1051                 c->wds = 1;
1052                 c->x[0] = 0;
1053                 return c;
1054                 }
1055         if (i < 0) {
1056                 c = a;
1057                 a = b;
1058                 b = c;
1059                 i = 1;
1060                 }
1061         else
1062                 i = 0;
1063         c = Balloc(a->k);
1064         c->sign = i;
1065         wa = a->wds;
1066         xa = a->x;
1067         xae = xa + wa;
1068         wb = b->wds;
1069         xb = b->x;
1070         xbe = xb + wb;
1071         xc = c->x;
1072         borrow = 0;
1073 #ifdef ULLong
1074         do {
1075                 y = (ULLong)*xa++ - *xb++ - borrow;
1076                 borrow = y >> 32 & (ULong)1;
1077                 *xc++ = y & FFFFFFFF;
1078                 }
1079                 while(xb < xbe);
1080         while(xa < xae) {
1081                 y = *xa++ - borrow;
1082                 borrow = y >> 32 & (ULong)1;
1083                 *xc++ = y & FFFFFFFF;
1084                 }
1085 #else
1086 #ifdef Pack_32
1087         do {
1088                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1089                 borrow = (y & 0x10000) >> 16;
1090                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1091                 borrow = (z & 0x10000) >> 16;
1092                 Storeinc(xc, z, y);
1093                 }
1094                 while(xb < xbe);
1095         while(xa < xae) {
1096                 y = (*xa & 0xffff) - borrow;
1097                 borrow = (y & 0x10000) >> 16;
1098                 z = (*xa++ >> 16) - borrow;
1099                 borrow = (z & 0x10000) >> 16;
1100                 Storeinc(xc, z, y);
1101                 }
1102 #else
1103         do {
1104                 y = *xa++ - *xb++ - borrow;
1105                 borrow = (y & 0x10000) >> 16;
1106                 *xc++ = y & 0xffff;
1107                 }
1108                 while(xb < xbe);
1109         while(xa < xae) {
1110                 y = *xa++ - borrow;
1111                 borrow = (y & 0x10000) >> 16;
1112                 *xc++ = y & 0xffff;
1113                 }
1114 #endif
1115 #endif
1116         while(!*--xc)
1117                 wa--;
1118         c->wds = wa;
1119         return c;
1120         }
1121
1122  static double
1123 ulp
1124 #ifdef KR_headers
1125         (x) double x;
1126 #else
1127         (double x)
1128 #endif
1129 {
1130         register Long L;
1131         double a;
1132
1133         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1134 #ifndef Avoid_Underflow
1135 #ifndef Sudden_Underflow
1136         if (L > 0) {
1137 #endif
1138 #endif
1139 #ifdef IBM
1140                 L |= Exp_msk1 >> 4;
1141 #endif
1142                 word0(a) = L;
1143                 word1(a) = 0;
1144 #ifndef Avoid_Underflow
1145 #ifndef Sudden_Underflow
1146                 }
1147         else {
1148                 L = -L >> Exp_shift;
1149                 if (L < Exp_shift) {
1150                         word0(a) = 0x80000 >> L;
1151                         word1(a) = 0;
1152                         }
1153                 else {
1154                         word0(a) = 0;
1155                         L -= Exp_shift;
1156                         word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1157                         }
1158                 }
1159 #endif
1160 #endif
1161         return dval(a);
1162         }
1163
1164  static double
1165 b2d
1166 #ifdef KR_headers
1167         (a, e) Bigint *a; int *e;
1168 #else
1169         (Bigint *a, int *e)
1170 #endif
1171 {
1172         ULong *xa, *xa0, w, y, z;
1173         int k;
1174         double d;
1175 #ifdef VAX
1176         ULong d0, d1;
1177 #else
1178 #define d0 word0(d)
1179 #define d1 word1(d)
1180 #endif
1181
1182         xa0 = a->x;
1183         xa = xa0 + a->wds;
1184         y = *--xa;
1185 #ifdef DEBUG
1186         if (!y) Bug("zero y in b2d");
1187 #endif
1188         k = hi0bits(y);
1189         *e = 32 - k;
1190 #ifdef Pack_32
1191         if (k < Ebits) {
1192                 d0 = Exp_1 | y >> Ebits - k;
1193                 w = xa > xa0 ? *--xa : 0;
1194                 d1 = y << (32-Ebits) + k | w >> Ebits - k;
1195                 goto ret_d;
1196                 }
1197         z = xa > xa0 ? *--xa : 0;
1198         if (k -= Ebits) {
1199                 d0 = Exp_1 | y << k | z >> 32 - k;
1200                 y = xa > xa0 ? *--xa : 0;
1201                 d1 = z << k | y >> 32 - k;
1202                 }
1203         else {
1204                 d0 = Exp_1 | y;
1205                 d1 = z;
1206                 }
1207 #else
1208         if (k < Ebits + 16) {
1209                 z = xa > xa0 ? *--xa : 0;
1210                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1211                 w = xa > xa0 ? *--xa : 0;
1212                 y = xa > xa0 ? *--xa : 0;
1213                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1214                 goto ret_d;
1215                 }
1216         z = xa > xa0 ? *--xa : 0;
1217         w = xa > xa0 ? *--xa : 0;
1218         k -= Ebits + 16;
1219         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1220         y = xa > xa0 ? *--xa : 0;
1221         d1 = w << k + 16 | y << k;
1222 #endif
1223  ret_d:
1224 #ifdef VAX
1225         word0(d) = d0 >> 16 | d0 << 16;
1226         word1(d) = d1 >> 16 | d1 << 16;
1227 #else
1228 #undef d0
1229 #undef d1
1230 #endif
1231         return dval(d);
1232         }
1233
1234  static Bigint *
1235 d2b
1236 #ifdef KR_headers
1237         (d, e, bits) double d; int *e, *bits;
1238 #else
1239         (double d, int *e, int *bits)
1240 #endif
1241 {
1242         Bigint *b;
1243         int de, k;
1244         ULong *x, y, z;
1245 #ifndef Sudden_Underflow
1246         int i;
1247 #endif
1248 #ifdef VAX
1249         ULong d0, d1;
1250         d0 = word0(d) >> 16 | word0(d) << 16;
1251         d1 = word1(d) >> 16 | word1(d) << 16;
1252 #else
1253 #define d0 word0(d)
1254 #define d1 word1(d)
1255 #endif
1256
1257 #ifdef Pack_32
1258         b = Balloc(1);
1259 #else
1260         b = Balloc(2);
1261 #endif
1262         x = b->x;
1263
1264         z = d0 & Frac_mask;
1265         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
1266 #ifdef Sudden_Underflow
1267         de = (int)(d0 >> Exp_shift);
1268 #ifndef IBM
1269         z |= Exp_msk11;
1270 #endif
1271 #else
1272         if ((de = (int)(d0 >> Exp_shift)))
1273                 z |= Exp_msk1;
1274 #endif
1275 #ifdef Pack_32
1276         if ((y = d1)) {
1277                 if ((k = lo0bits(&y))) {
1278                         x[0] = y | z << 32 - k;
1279                         z >>= k;
1280                         }
1281                 else
1282                         x[0] = y;
1283 #ifndef Sudden_Underflow
1284                 i =
1285 #endif
1286                     b->wds = (x[1] = z) ? 2 : 1;
1287                 }
1288         else {
1289 #ifdef DEBUG
1290                 if (!z)
1291                         Bug("Zero passed to d2b");
1292 #endif
1293                 k = lo0bits(&z);
1294                 x[0] = z;
1295 #ifndef Sudden_Underflow
1296                 i =
1297 #endif
1298                     b->wds = 1;
1299                 k += 32;
1300                 }
1301 #else
1302         if (y = d1) {
1303                 if (k = lo0bits(&y))
1304                         if (k >= 16) {
1305                                 x[0] = y | z << 32 - k & 0xffff;
1306                                 x[1] = z >> k - 16 & 0xffff;
1307                                 x[2] = z >> k;
1308                                 i = 2;
1309                                 }
1310                         else {
1311                                 x[0] = y & 0xffff;
1312                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
1313                                 x[2] = z >> k & 0xffff;
1314                                 x[3] = z >> k+16;
1315                                 i = 3;
1316                                 }
1317                 else {
1318                         x[0] = y & 0xffff;
1319                         x[1] = y >> 16;
1320                         x[2] = z & 0xffff;
1321                         x[3] = z >> 16;
1322                         i = 3;
1323                         }
1324                 }
1325         else {
1326 #ifdef DEBUG
1327                 if (!z)
1328                         Bug("Zero passed to d2b");
1329 #endif
1330                 k = lo0bits(&z);
1331                 if (k >= 16) {
1332                         x[0] = z;
1333                         i = 0;
1334                         }
1335                 else {
1336                         x[0] = z & 0xffff;
1337                         x[1] = z >> 16;
1338                         i = 1;
1339                         }
1340                 k += 32;
1341                 }
1342         while(!x[i])
1343                 --i;
1344         b->wds = i + 1;
1345 #endif
1346 #ifndef Sudden_Underflow
1347         if (de) {
1348 #endif
1349 #ifdef IBM
1350                 *e = (de - Bias - (P-1) << 2) + k;
1351                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1352 #else
1353                 *e = de - Bias - (P-1) + k;
1354                 *bits = P - k;
1355 #endif
1356 #ifndef Sudden_Underflow
1357                 }
1358         else {
1359                 *e = de - Bias - (P-1) + 1 + k;
1360 #ifdef Pack_32
1361                 *bits = 32*i - hi0bits(x[i-1]);
1362 #else
1363                 *bits = (i+2)*16 - hi0bits(x[i]);
1364 #endif
1365                 }
1366 #endif
1367         return b;
1368         }
1369 #undef d0
1370 #undef d1
1371
1372  static double
1373 ratio
1374 #ifdef KR_headers
1375         (a, b) Bigint *a, *b;
1376 #else
1377         (Bigint *a, Bigint *b)
1378 #endif
1379 {
1380         double da, db;
1381         int k, ka, kb;
1382
1383         dval(da) = b2d(a, &ka);
1384         dval(db) = b2d(b, &kb);
1385 #ifdef Pack_32
1386         k = ka - kb + 32*(a->wds - b->wds);
1387 #else
1388         k = ka - kb + 16*(a->wds - b->wds);
1389 #endif
1390 #ifdef IBM
1391         if (k > 0) {
1392                 word0(da) += (k >> 2)*Exp_msk1;
1393                 if (k &= 3)
1394                         dval(da) *= 1 << k;
1395                 }
1396         else {
1397                 k = -k;
1398                 word0(db) += (k >> 2)*Exp_msk1;
1399                 if (k &= 3)
1400                         dval(db) *= 1 << k;
1401                 }
1402 #else
1403         if (k > 0)
1404                 word0(da) += k*Exp_msk1;
1405         else {
1406                 k = -k;
1407                 word0(db) += k*Exp_msk1;
1408                 }
1409 #endif
1410         return dval(da) / dval(db);
1411         }
1412
1413  static CONST double
1414 tens[] = {
1415                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1416                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1417                 1e20, 1e21, 1e22
1418 #ifdef VAX
1419                 , 1e23, 1e24
1420 #endif
1421                 };
1422
1423  static CONST double
1424 #ifdef IEEE_Arith
1425 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1426 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1427 #ifdef Avoid_Underflow
1428                 9007199254740992.*9007199254740992.e-256
1429                 /* = 2^106 * 1e-53 */
1430 #else
1431                 1e-256
1432 #endif
1433                 };
1434 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1435 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
1436 #define Scale_Bit 0x10
1437 #define n_bigtens 5
1438 #else
1439 #ifdef IBM
1440 bigtens[] = { 1e16, 1e32, 1e64 };
1441 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1442 #define n_bigtens 3
1443 #else
1444 bigtens[] = { 1e16, 1e32 };
1445 static CONST double tinytens[] = { 1e-16, 1e-32 };
1446 #define n_bigtens 2
1447 #endif
1448 #endif
1449
1450 #ifndef IEEE_Arith
1451 #undef INFNAN_CHECK
1452 #endif
1453
1454 #ifdef INFNAN_CHECK
1455
1456 #ifndef NAN_WORD0
1457 #define NAN_WORD0 0x7ff80000
1458 #endif
1459
1460 #ifndef NAN_WORD1
1461 #define NAN_WORD1 0
1462 #endif
1463
1464  static int
1465 match
1466 #ifdef KR_headers
1467         (sp, t) char **sp, *t;
1468 #else
1469         (CONST char **sp, char *t)
1470 #endif
1471 {
1472         int c, d;
1473         CONST char *s = *sp;
1474
1475         while(d = *t++) {
1476                 if ((c = *++s) >= 'A' && c <= 'Z')
1477                         c += 'a' - 'A';
1478                 if (c != d)
1479                         return 0;
1480                 }
1481         *sp = s + 1;
1482         return 1;
1483         }
1484
1485 #ifndef No_Hex_NaN
1486  static void
1487 hexnan
1488 #ifdef KR_headers
1489         (rvp, sp) double *rvp; CONST char **sp;
1490 #else
1491         (double *rvp, CONST char **sp)
1492 #endif
1493 {
1494         ULong c, x[2];
1495         CONST char *s;
1496         int havedig, udx0, xshift;
1497
1498         x[0] = x[1] = 0;
1499         havedig = xshift = 0;
1500         udx0 = 1;
1501         s = *sp;
1502         while(c = *(CONST unsigned char*)++s) {
1503                 if (c >= '0' && c <= '9')
1504                         c -= '0';
1505                 else if (c >= 'a' && c <= 'f')
1506                         c += 10 - 'a';
1507                 else if (c >= 'A' && c <= 'F')
1508                         c += 10 - 'A';
1509                 else if (c <= ' ') {
1510                         if (udx0 && havedig) {
1511                                 udx0 = 0;
1512                                 xshift = 1;
1513                                 }
1514                         continue;
1515                         }
1516                 else if (/*(*/ c == ')' && havedig) {
1517                         *sp = s + 1;
1518                         break;
1519                         }
1520                 else
1521                         return; /* invalid form: don't change *sp */
1522                 havedig = 1;
1523                 if (xshift) {
1524                         xshift = 0;
1525                         x[0] = x[1];
1526                         x[1] = 0;
1527                         }
1528                 if (udx0)
1529                         x[0] = (x[0] << 4) | (x[1] >> 28);
1530                 x[1] = (x[1] << 4) | c;
1531                 }
1532         if ((x[0] &= 0xfffff) || x[1]) {
1533                 word0(*rvp) = Exp_mask | x[0];
1534                 word1(*rvp) = x[1];
1535                 }
1536         }
1537 #endif /*No_Hex_NaN*/
1538 #endif /* INFNAN_CHECK */
1539
1540  double
1541 bsd_strtod
1542 #ifdef KR_headers
1543         (s00, se) CONST char *s00; char **se;
1544 #else
1545         (CONST char *s00, char **se)
1546 #endif
1547 {
1548 #ifdef Avoid_Underflow
1549         int scale;
1550 #endif
1551         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1552                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1553         CONST char *s, *s0, *s1;
1554         double aadj, aadj1, adj, rv, rv0;
1555         Long L;
1556         ULong y, z;
1557         Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1558 #ifdef SET_INEXACT
1559         int inexact, oldinexact;
1560 #endif
1561 #ifdef Honor_FLT_ROUNDS
1562         int rounding;
1563 #endif
1564 #ifdef USE_LOCALE
1565         CONST char *s2;
1566 #endif
1567
1568         sign = nz0 = nz = 0;
1569         dval(rv) = 0.;
1570         for(s = s00;;s++) switch(*s) {
1571                 case '-':
1572                         sign = 1;
1573                         /* no break */
1574                 case '+':
1575                         if (*++s)
1576                                 goto break2;
1577                         /* no break */
1578                 case 0:
1579                         goto ret0;
1580                 case '\t':
1581                 case '\n':
1582                 case '\v':
1583                 case '\f':
1584                 case '\r':
1585                 case ' ':
1586                         continue;
1587                 default:
1588                         goto break2;
1589                 }
1590  break2:
1591         if (*s == '0') {
1592                 nz0 = 1;
1593                 while(*++s == '0') ;
1594                 if (!*s)
1595                         goto ret;
1596                 }
1597         s0 = s;
1598         y = z = 0;
1599         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1600                 if (nd < 9)
1601                         y = 10*y + c - '0';
1602                 else if (nd < 16)
1603                         z = 10*z + c - '0';
1604         nd0 = nd;
1605 #ifdef USE_LOCALE
1606         s1 = localeconv()->decimal_point;
1607         if (c == *s1) {
1608                 c = '.';
1609                 if (*++s1) {
1610                         s2 = s;
1611                         for(;;) {
1612                                 if (*++s2 != *s1) {
1613                                         c = 0;
1614                                         break;
1615                                         }
1616                                 if (!*++s1) {
1617                                         s = s2;
1618                                         break;
1619                                         }
1620                                 }
1621                         }
1622                 }
1623 #endif
1624         if (c == '.') {
1625                 c = *++s;
1626                 if (!nd) {
1627                         for(; c == '0'; c = *++s)
1628                                 nz++;
1629                         if (c > '0' && c <= '9') {
1630                                 s0 = s;
1631                                 nf += nz;
1632                                 nz = 0;
1633                                 goto have_dig;
1634                                 }
1635                         goto dig_done;
1636                         }
1637                 for(; c >= '0' && c <= '9'; c = *++s) {
1638  have_dig:
1639                         nz++;
1640                         if (c -= '0') {
1641                                 nf += nz;
1642                                 for(i = 1; i < nz; i++)
1643                                         if (nd++ < 9)
1644                                                 y *= 10;
1645                                         else if (nd <= DBL_DIG + 1)
1646                                                 z *= 10;
1647                                 if (nd++ < 9)
1648                                         y = 10*y + c;
1649                                 else if (nd <= DBL_DIG + 1)
1650                                         z = 10*z + c;
1651                                 nz = 0;
1652                                 }
1653                         }
1654                 }
1655  dig_done:
1656         e = 0;
1657         if (c == 'e' || c == 'E') {
1658                 if (!nd && !nz && !nz0) {
1659                         goto ret0;
1660                         }
1661                 s00 = s;
1662                 esign = 0;
1663                 switch(c = *++s) {
1664                         case '-':
1665                                 esign = 1;
1666                         case '+':
1667                                 c = *++s;
1668                         }
1669                 if (c >= '0' && c <= '9') {
1670                         while(c == '0')
1671                                 c = *++s;
1672                         if (c > '0' && c <= '9') {
1673                                 L = c - '0';
1674                                 s1 = s;
1675                                 while((c = *++s) >= '0' && c <= '9')
1676                                         L = 10*L + c - '0';
1677                                 if (s - s1 > 8 || L > 19999)
1678                                         /* Avoid confusion from exponents
1679                                          * so large that e might overflow.
1680                                          */
1681                                         e = 19999; /* safe for 16 bit ints */
1682                                 else
1683                                         e = (int)L;
1684                                 if (esign)
1685                                         e = -e;
1686                                 }
1687                         else
1688                                 e = 0;
1689                         }
1690                 else
1691                         s = s00;
1692                 }
1693         if (!nd) {
1694                 if (!nz && !nz0) {
1695 #ifdef INFNAN_CHECK
1696                         /* Check for Nan and Infinity */
1697                         switch(c) {
1698                           case 'i':
1699                           case 'I':
1700                                 if (match(&s,"nf")) {
1701                                         --s;
1702                                         if (!match(&s,"inity"))
1703                                                 ++s;
1704                                         word0(rv) = 0x7ff00000;
1705                                         word1(rv) = 0;
1706                                         goto ret;
1707                                         }
1708                                 break;
1709                           case 'n':
1710                           case 'N':
1711                                 if (match(&s, "an")) {
1712                                         word0(rv) = NAN_WORD0;
1713                                         word1(rv) = NAN_WORD1;
1714 #ifndef No_Hex_NaN
1715                                         if (*s == '(') /*)*/
1716                                                 hexnan(&rv, &s);
1717 #endif
1718                                         goto ret;
1719                                         }
1720                           }
1721 #endif /* INFNAN_CHECK */
1722  ret0:
1723                         s = s00;
1724                         sign = 0;
1725                         }
1726                 goto ret;
1727                 }
1728         e1 = e -= nf;
1729
1730         /* Now we have nd0 digits, starting at s0, followed by a
1731          * decimal point, followed by nd-nd0 digits.  The number we're
1732          * after is the integer represented by those digits times
1733          * 10**e */
1734
1735         if (!nd0)
1736                 nd0 = nd;
1737         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1738         dval(rv) = y;
1739         if (k > 9) {
1740 #ifdef SET_INEXACT
1741                 if (k > DBL_DIG)
1742                         oldinexact = get_inexact();
1743 #endif
1744                 dval(rv) = tens[k - 9] * dval(rv) + z;
1745                 }
1746         bd0 = 0;
1747         if (nd <= DBL_DIG
1748 #ifndef RND_PRODQUOT
1749 #ifndef Honor_FLT_ROUNDS
1750                 && Flt_Rounds == 1
1751 #endif
1752 #endif
1753                         ) {
1754                 if (!e)
1755                         goto ret;
1756                 if (e > 0) {
1757                         if (e <= Ten_pmax) {
1758 #ifdef VAX
1759                                 goto vax_ovfl_check;
1760 #else
1761 #ifdef Honor_FLT_ROUNDS
1762                                 /* round correctly FLT_ROUNDS = 2 or 3 */
1763                                 if (sign) {
1764                                         rv = -rv;
1765                                         sign = 0;
1766                                         }
1767 #endif
1768                                 /* rv = */ rounded_product(dval(rv), tens[e]);
1769                                 goto ret;
1770 #endif
1771                                 }
1772                         i = DBL_DIG - nd;
1773                         if (e <= Ten_pmax + i) {
1774                                 /* A fancier test would sometimes let us do
1775                                  * this for larger i values.
1776                                  */
1777 #ifdef Honor_FLT_ROUNDS
1778                                 /* round correctly FLT_ROUNDS = 2 or 3 */
1779                                 if (sign) {
1780                                         rv = -rv;
1781                                         sign = 0;
1782                                         }
1783 #endif
1784                                 e -= i;
1785                                 dval(rv) *= tens[i];
1786 #ifdef VAX
1787                                 /* VAX exponent range is so narrow we must
1788                                  * worry about overflow here...
1789                                  */
1790  vax_ovfl_check:
1791                                 word0(rv) -= P*Exp_msk1;
1792                                 /* rv = */ rounded_product(dval(rv), tens[e]);
1793                                 if ((word0(rv) & Exp_mask)
1794                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1795                                         goto ovfl;
1796                                 word0(rv) += P*Exp_msk1;
1797 #else
1798                                 /* rv = */ rounded_product(dval(rv), tens[e]);
1799 #endif
1800                                 goto ret;
1801                                 }
1802                         }
1803 #ifndef Inaccurate_Divide
1804                 else if (e >= -Ten_pmax) {
1805 #ifdef Honor_FLT_ROUNDS
1806                         /* round correctly FLT_ROUNDS = 2 or 3 */
1807                         if (sign) {
1808                                 rv = -rv;
1809                                 sign = 0;
1810                                 }
1811 #endif
1812                         /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1813                         goto ret;
1814                         }
1815 #endif
1816                 }
1817         e1 += nd - k;
1818
1819 #ifdef IEEE_Arith
1820 #ifdef SET_INEXACT
1821         inexact = 1;
1822         if (k <= DBL_DIG)
1823                 oldinexact = get_inexact();
1824 #endif
1825 #ifdef Avoid_Underflow
1826         scale = 0;
1827 #endif
1828 #ifdef Honor_FLT_ROUNDS
1829         if ((rounding = Flt_Rounds) >= 2) {
1830                 if (sign)
1831                         rounding = rounding == 2 ? 0 : 2;
1832                 else
1833                         if (rounding != 2)
1834                                 rounding = 0;
1835                 }
1836 #endif
1837 #endif /*IEEE_Arith*/
1838
1839         /* Get starting approximation = rv * 10**e1 */
1840
1841         if (e1 > 0) {
1842                 if ((i = e1 & 15))
1843                         dval(rv) *= tens[i];
1844                 if (e1 &= ~15) {
1845                         if (e1 > DBL_MAX_10_EXP) {
1846  ovfl:
1847 #ifndef NO_ERRNO
1848                                 errno = ERANGE;
1849 #endif
1850                                 /* Can't trust HUGE_VAL */
1851 #ifdef IEEE_Arith
1852 #ifdef Honor_FLT_ROUNDS
1853                                 switch(rounding) {
1854                                   case 0: /* toward 0 */
1855                                   case 3: /* toward -infinity */
1856                                         word0(rv) = Big0;
1857                                         word1(rv) = Big1;
1858                                         break;
1859                                   default:
1860                                         word0(rv) = Exp_mask;
1861                                         word1(rv) = 0;
1862                                   }
1863 #else /*Honor_FLT_ROUNDS*/
1864                                 word0(rv) = Exp_mask;
1865                                 word1(rv) = 0;
1866 #endif /*Honor_FLT_ROUNDS*/
1867 #ifdef SET_INEXACT
1868                                 /* set overflow bit */
1869                                 dval(rv0) = 1e300;
1870                                 dval(rv0) *= dval(rv0);
1871 #endif
1872 #else /*IEEE_Arith*/
1873                                 word0(rv) = Big0;
1874                                 word1(rv) = Big1;
1875 #endif /*IEEE_Arith*/
1876                                 if (bd0)
1877                                         goto retfree;
1878                                 goto ret;
1879                                 }
1880                         e1 >>= 4;
1881                         for(j = 0; e1 > 1; j++, e1 >>= 1)
1882                                 if (e1 & 1)
1883                                         dval(rv) *= bigtens[j];
1884                 /* The last multiplication could overflow. */
1885                         word0(rv) -= P*Exp_msk1;
1886                         dval(rv) *= bigtens[j];
1887                         if ((z = word0(rv) & Exp_mask)
1888                          > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1889                                 goto ovfl;
1890                         if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1891                                 /* set to largest number */
1892                                 /* (Can't trust DBL_MAX) */
1893                                 word0(rv) = Big0;
1894                                 word1(rv) = Big1;
1895                                 }
1896                         else
1897                                 word0(rv) += P*Exp_msk1;
1898                         }
1899                 }
1900         else if (e1 < 0) {
1901                 e1 = -e1;
1902                 if ((i = e1 & 15))
1903                         dval(rv) /= tens[i];
1904                 if (e1 >>= 4) {
1905                         if (e1 >= 1 << n_bigtens)
1906                                 goto undfl;
1907 #ifdef Avoid_Underflow
1908                         if (e1 & Scale_Bit)
1909                                 scale = 2*P;
1910                         for(j = 0; e1 > 0; j++, e1 >>= 1)
1911                                 if (e1 & 1)
1912                                         dval(rv) *= tinytens[j];
1913                         if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1914                                                 >> Exp_shift)) > 0) {
1915                                 /* scaled rv is denormal; zap j low bits */
1916                                 if (j >= 32) {
1917                                         word1(rv) = 0;
1918                                         if (j >= 53)
1919                                          word0(rv) = (P+2)*Exp_msk1;
1920                                         else
1921                                          word0(rv) &= 0xffffffff << j-32;
1922                                         }
1923                                 else
1924                                         word1(rv) &= 0xffffffff << j;
1925                                 }
1926 #else
1927                         for(j = 0; e1 > 1; j++, e1 >>= 1)
1928                                 if (e1 & 1)
1929                                         dval(rv) *= tinytens[j];
1930                         /* The last multiplication could underflow. */
1931                         dval(rv0) = dval(rv);
1932                         dval(rv) *= tinytens[j];
1933                         if (!dval(rv)) {
1934                                 dval(rv) = 2.*dval(rv0);
1935                                 dval(rv) *= tinytens[j];
1936 #endif
1937                                 if (!dval(rv)) {
1938  undfl:
1939                                         dval(rv) = 0.;
1940 #ifndef NO_ERRNO
1941                                         errno = ERANGE;
1942 #endif
1943                                         if (bd0)
1944                                                 goto retfree;
1945                                         goto ret;
1946                                         }
1947 #ifndef Avoid_Underflow
1948                                 word0(rv) = Tiny0;
1949                                 word1(rv) = Tiny1;
1950                                 /* The refinement below will clean
1951                                  * this approximation up.
1952                                  */
1953                                 }
1954 #endif
1955                         }
1956                 }
1957
1958         /* Now the hard part -- adjusting rv to the correct value.*/
1959
1960         /* Put digits into bd: true value = bd * 10^e */
1961
1962         bd0 = s2b(s0, nd0, nd, y);
1963
1964         for(;;) {
1965                 bd = Balloc(bd0->k);
1966                 Bcopy(bd, bd0);
1967                 bb = d2b(dval(rv), &bbe, &bbbits);      /* rv = bb * 2^bbe */
1968                 bs = i2b(1);
1969
1970                 if (e >= 0) {
1971                         bb2 = bb5 = 0;
1972                         bd2 = bd5 = e;
1973                         }
1974                 else {
1975                         bb2 = bb5 = -e;
1976                         bd2 = bd5 = 0;
1977                         }
1978                 if (bbe >= 0)
1979                         bb2 += bbe;
1980                 else
1981                         bd2 -= bbe;
1982                 bs2 = bb2;
1983 #ifdef Honor_FLT_ROUNDS
1984                 if (rounding != 1)
1985                         bs2++;
1986 #endif
1987 #ifdef Avoid_Underflow
1988                 j = bbe - scale;
1989                 i = j + bbbits - 1;     /* logb(rv) */
1990                 if (i < Emin)   /* denormal */
1991                         j += P - Emin;
1992                 else
1993                         j = P + 1 - bbbits;
1994 #else /*Avoid_Underflow*/
1995 #ifdef Sudden_Underflow
1996 #ifdef IBM
1997                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1998 #else
1999                 j = P + 1 - bbbits;
2000 #endif
2001 #else /*Sudden_Underflow*/
2002                 j = bbe;
2003                 i = j + bbbits - 1;     /* logb(rv) */
2004                 if (i < Emin)   /* denormal */
2005                         j += P - Emin;
2006                 else
2007                         j = P + 1 - bbbits;
2008 #endif /*Sudden_Underflow*/
2009 #endif /*Avoid_Underflow*/
2010                 bb2 += j;
2011                 bd2 += j;
2012 #ifdef Avoid_Underflow
2013                 bd2 += scale;
2014 #endif
2015                 i = bb2 < bd2 ? bb2 : bd2;
2016                 if (i > bs2)
2017                         i = bs2;
2018                 if (i > 0) {
2019                         bb2 -= i;
2020                         bd2 -= i;
2021                         bs2 -= i;
2022                         }
2023                 if (bb5 > 0) {
2024                         bs = pow5mult(bs, bb5);
2025                         bb1 = mult(bs, bb);
2026                         Bfree(bb);
2027                         bb = bb1;
2028                         }
2029                 if (bb2 > 0)
2030                         bb = lshift(bb, bb2);
2031                 if (bd5 > 0)
2032                         bd = pow5mult(bd, bd5);
2033                 if (bd2 > 0)
2034                         bd = lshift(bd, bd2);
2035                 if (bs2 > 0)
2036                         bs = lshift(bs, bs2);
2037                 delta = diff(bb, bd);
2038                 dsign = delta->sign;
2039                 delta->sign = 0;
2040                 i = cmp(delta, bs);
2041 #ifdef Honor_FLT_ROUNDS
2042                 if (rounding != 1) {
2043                         if (i < 0) {
2044                                 /* Error is less than an ulp */
2045                                 if (!delta->x[0] && delta->wds <= 1) {
2046                                         /* exact */
2047 #ifdef SET_INEXACT
2048                                         inexact = 0;
2049 #endif
2050                                         break;
2051                                         }
2052                                 if (rounding) {
2053                                         if (dsign) {
2054                                                 adj = 1.;
2055                                                 goto apply_adj;
2056                                                 }
2057                                         }
2058                                 else if (!dsign) {
2059                                         adj = -1.;
2060                                         if (!word1(rv)
2061                                          && !(word0(rv) & Frac_mask)) {
2062                                                 y = word0(rv) & Exp_mask;
2063 #ifdef Avoid_Underflow
2064                                                 if (!scale || y > 2*P*Exp_msk1)
2065 #else
2066                                                 if (y)
2067 #endif
2068                                                   {
2069                                                   delta = lshift(delta,Log2P);
2070                                                   if (cmp(delta, bs) <= 0)
2071                                                         adj = -0.5;
2072                                                   }
2073                                                 }
2074  apply_adj:
2075 #ifdef Avoid_Underflow
2076                                         if (scale && (y = word0(rv) & Exp_mask)
2077                                                 <= 2*P*Exp_msk1)
2078                                           word0(adj) += (2*P+1)*Exp_msk1 - y;
2079 #else
2080 #ifdef Sudden_Underflow
2081                                         if ((word0(rv) & Exp_mask) <=
2082                                                         P*Exp_msk1) {
2083                                                 word0(rv) += P*Exp_msk1;
2084                                                 dval(rv) += adj*ulp(dval(rv));
2085                                                 word0(rv) -= P*Exp_msk1;
2086                                                 }
2087                                         else
2088 #endif /*Sudden_Underflow*/
2089 #endif /*Avoid_Underflow*/
2090                                         dval(rv) += adj*ulp(dval(rv));
2091                                         }
2092                                 break;
2093                                 }
2094                         adj = ratio(delta, bs);
2095                         if (adj < 1.)
2096                                 adj = 1.;
2097                         if (adj <= 0x7ffffffe) {
2098                                 /* adj = rounding ? ceil(adj) : floor(adj); */
2099                                 y = adj;
2100                                 if (y != adj) {
2101                                         if (!((rounding>>1) ^ dsign))
2102                                                 y++;
2103                                         adj = y;
2104                                         }
2105                                 }
2106 #ifdef Avoid_Underflow
2107                         if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2108                                 word0(adj) += (2*P+1)*Exp_msk1 - y;
2109 #else
2110 #ifdef Sudden_Underflow
2111                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2112                                 word0(rv) += P*Exp_msk1;
2113                                 adj *= ulp(dval(rv));
2114                                 if (dsign)
2115                                         dval(rv) += adj;
2116                                 else
2117                                         dval(rv) -= adj;
2118                                 word0(rv) -= P*Exp_msk1;
2119                                 goto cont;
2120                                 }
2121 #endif /*Sudden_Underflow*/
2122 #endif /*Avoid_Underflow*/
2123                         adj *= ulp(dval(rv));
2124                         if (dsign)
2125                                 dval(rv) += adj;
2126                         else
2127                                 dval(rv) -= adj;
2128                         goto cont;
2129                         }
2130 #endif /*Honor_FLT_ROUNDS*/
2131
2132                 if (i < 0) {
2133                         /* Error is less than half an ulp -- check for
2134                          * special case of mantissa a power of two.
2135                          */
2136                         if (dsign || word1(rv) || word0(rv) & Bndry_mask
2137 #ifdef IEEE_Arith
2138 #ifdef Avoid_Underflow
2139                          || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2140 #else
2141                          || (word0(rv) & Exp_mask) <= Exp_msk1
2142 #endif
2143 #endif
2144                                 ) {
2145 #ifdef SET_INEXACT
2146                                 if (!delta->x[0] && delta->wds <= 1)
2147                                         inexact = 0;
2148 #endif
2149                                 break;
2150                                 }
2151                         if (!delta->x[0] && delta->wds <= 1) {
2152                                 /* exact result */
2153 #ifdef SET_INEXACT
2154                                 inexact = 0;
2155 #endif
2156                                 break;
2157                                 }
2158                         delta = lshift(delta,Log2P);
2159                         if (cmp(delta, bs) > 0)
2160                                 goto drop_down;
2161                         break;
2162                         }
2163                 if (i == 0) {
2164                         /* exactly half-way between */
2165                         if (dsign) {
2166                                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2167                                  &&  word1(rv) == (
2168 #ifdef Avoid_Underflow
2169                         (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2170                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2171 #endif
2172                                                    0xffffffff)) {
2173                                         /*boundary case -- increment exponent*/
2174                                         word0(rv) = (word0(rv) & Exp_mask)
2175                                                 + Exp_msk1
2176 #ifdef IBM
2177                                                 | Exp_msk1 >> 4
2178 #endif
2179                                                 ;
2180                                         word1(rv) = 0;
2181 #ifdef Avoid_Underflow
2182                                         dsign = 0;
2183 #endif
2184                                         break;
2185                                         }
2186                                 }
2187                         else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2188  drop_down:
2189                                 /* boundary case -- decrement exponent */
2190 #ifdef Sudden_Underflow /*{{*/
2191                                 L = word0(rv) & Exp_mask;
2192 #ifdef IBM
2193                                 if (L <  Exp_msk1)
2194 #else
2195 #ifdef Avoid_Underflow
2196                                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2197 #else
2198                                 if (L <= Exp_msk1)
2199 #endif /*Avoid_Underflow*/
2200 #endif /*IBM*/
2201                                         goto undfl;
2202                                 L -= Exp_msk1;
2203 #else /*Sudden_Underflow}{*/
2204 #ifdef Avoid_Underflow
2205                                 if (scale) {
2206                                         L = word0(rv) & Exp_mask;
2207                                         if (L <= (2*P+1)*Exp_msk1) {
2208                                                 if (L > (P+2)*Exp_msk1)
2209                                                         /* round even ==> */
2210                                                         /* accept rv */
2211                                                         break;
2212                                                 /* rv = smallest denormal */
2213                                                 goto undfl;
2214                                                 }
2215                                         }
2216 #endif /*Avoid_Underflow*/
2217                                 L = (word0(rv) & Exp_mask) - Exp_msk1;
2218 #endif /*Sudden_Underflow}}*/
2219                                 word0(rv) = L | Bndry_mask1;
2220                                 word1(rv) = 0xffffffff;
2221 #ifdef IBM
2222                                 goto cont;
2223 #else
2224                                 break;
2225 #endif
2226                                 }
2227 #ifndef ROUND_BIASED
2228                         if (!(word1(rv) & LSB))
2229                                 break;
2230 #endif
2231                         if (dsign)
2232                                 dval(rv) += ulp(dval(rv));
2233 #ifndef ROUND_BIASED
2234                         else {
2235                                 dval(rv) -= ulp(dval(rv));
2236 #ifndef Sudden_Underflow
2237                                 if (!dval(rv))
2238                                         goto undfl;
2239 #endif
2240                                 }
2241 #ifdef Avoid_Underflow
2242                         dsign = 1 - dsign;
2243 #endif
2244 #endif
2245                         break;
2246                         }
2247                 if ((aadj = ratio(delta, bs)) <= 2.) {
2248                         if (dsign)
2249                                 aadj = aadj1 = 1.;
2250                         else if (word1(rv) || word0(rv) & Bndry_mask) {
2251 #ifndef Sudden_Underflow
2252                                 if (word1(rv) == Tiny1 && !word0(rv))
2253                                         goto undfl;
2254 #endif
2255                                 aadj = 1.;
2256                                 aadj1 = -1.;
2257                                 }
2258                         else {
2259                                 /* special case -- power of FLT_RADIX to be */
2260                                 /* rounded down... */
2261
2262                                 if (aadj < 2./FLT_RADIX)
2263                                         aadj = 1./FLT_RADIX;
2264                                 else
2265                                         aadj *= 0.5;
2266                                 aadj1 = -aadj;
2267                                 }
2268                         }
2269                 else {
2270                         aadj *= 0.5;
2271                         aadj1 = dsign ? aadj : -aadj;
2272 #ifdef Check_FLT_ROUNDS
2273                         switch(Rounding) {
2274                                 case 2: /* towards +infinity */
2275                                         aadj1 -= 0.5;
2276                                         break;
2277                                 case 0: /* towards 0 */
2278                                 case 3: /* towards -infinity */
2279                                         aadj1 += 0.5;
2280                                 }
2281 #else
2282                         if (Flt_Rounds == 0)
2283                                 aadj1 += 0.5;
2284 #endif /*Check_FLT_ROUNDS*/
2285                         }
2286                 y = word0(rv) & Exp_mask;
2287
2288                 /* Check for overflow */
2289
2290                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2291                         dval(rv0) = dval(rv);
2292                         word0(rv) -= P*Exp_msk1;
2293                         adj = aadj1 * ulp(dval(rv));
2294                         dval(rv) += adj;
2295                         if ((word0(rv) & Exp_mask) >=
2296                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2297                                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2298                                         goto ovfl;
2299                                 word0(rv) = Big0;
2300                                 word1(rv) = Big1;
2301                                 goto cont;
2302                                 }
2303                         else
2304                                 word0(rv) += P*Exp_msk1;
2305                         }
2306                 else {
2307 #ifdef Avoid_Underflow
2308                         if (scale && y <= 2*P*Exp_msk1) {
2309                                 if (aadj <= 0x7fffffff) {
2310                                         if ((z = aadj) <= 0)
2311                                                 z = 1;
2312                                         aadj = z;
2313                                         aadj1 = dsign ? aadj : -aadj;
2314                                         }
2315                                 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2316                                 }
2317                         adj = aadj1 * ulp(dval(rv));
2318                         dval(rv) += adj;
2319 #else
2320 #ifdef Sudden_Underflow
2321                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2322                                 dval(rv0) = dval(rv);
2323                                 word0(rv) += P*Exp_msk1;
2324                                 adj = aadj1 * ulp(dval(rv));
2325                                 dval(rv) += adj;
2326 #ifdef IBM
2327                                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
2328 #else
2329                                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2330 #endif
2331                                         {
2332                                         if (word0(rv0) == Tiny0
2333                                          && word1(rv0) == Tiny1)
2334                                                 goto undfl;
2335                                         word0(rv) = Tiny0;
2336                                         word1(rv) = Tiny1;
2337                                         goto cont;
2338                                         }
2339                                 else
2340                                         word0(rv) -= P*Exp_msk1;
2341                                 }
2342                         else {
2343                                 adj = aadj1 * ulp(dval(rv));
2344                                 dval(rv) += adj;
2345                                 }
2346 #else /*Sudden_Underflow*/
2347                         /* Compute adj so that the IEEE rounding rules will
2348                          * correctly round rv + adj in some half-way cases.
2349                          * If rv * ulp(rv) is denormalized (i.e.,
2350                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2351                          * trouble from bits lost to denormalization;
2352                          * example: 1.2e-307 .
2353                          */
2354                         if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2355                                 aadj1 = (double)(int)(aadj + 0.5);
2356                                 if (!dsign)
2357                                         aadj1 = -aadj1;
2358                                 }
2359                         adj = aadj1 * ulp(dval(rv));
2360                         dval(rv) += adj;
2361 #endif /*Sudden_Underflow*/
2362 #endif /*Avoid_Underflow*/
2363                         }
2364                 z = word0(rv) & Exp_mask;
2365 #ifndef SET_INEXACT
2366 #ifdef Avoid_Underflow
2367                 if (!scale)
2368 #endif
2369                 if (y == z) {
2370                         /* Can we stop now? */
2371                         L = (Long)aadj;
2372                         aadj -= L;
2373                         /* The tolerances below are conservative. */
2374                         if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2375                                 if (aadj < .4999999 || aadj > .5000001)
2376                                         break;
2377                                 }
2378                         else if (aadj < .4999999/FLT_RADIX)
2379                                 break;
2380                         }
2381 #endif
2382  cont:
2383                 Bfree(bb);
2384                 Bfree(bd);
2385                 Bfree(bs);
2386                 Bfree(delta);
2387                 }
2388 #ifdef SET_INEXACT
2389         if (inexact) {
2390                 if (!oldinexact) {
2391                         word0(rv0) = Exp_1 + (70 << Exp_shift);
2392                         word1(rv0) = 0;
2393                         dval(rv0) += 1.;
2394                         }
2395                 }
2396         else if (!oldinexact)
2397                 clear_inexact();
2398 #endif
2399 #ifdef Avoid_Underflow
2400         if (scale) {
2401                 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2402                 word1(rv0) = 0;
2403                 dval(rv) *= dval(rv0);
2404 #ifndef NO_ERRNO
2405                 /* try to avoid the bug of testing an 8087 register value */
2406                 if (word0(rv) == 0 && word1(rv) == 0)
2407                         errno = ERANGE;
2408 #endif
2409                 }
2410 #endif /* Avoid_Underflow */
2411 #ifdef SET_INEXACT
2412         if (inexact && !(word0(rv) & Exp_mask)) {
2413                 /* set underflow bit */
2414                 dval(rv0) = 1e-300;
2415                 dval(rv0) *= dval(rv0);
2416                 }
2417 #endif
2418  retfree:
2419         Bfree(bb);
2420         Bfree(bd);
2421         Bfree(bs);
2422         Bfree(bd0);
2423         Bfree(delta);
2424  ret:
2425         if (se)
2426                 *se = (char *)s;
2427         return sign ? -dval(rv) : dval(rv);
2428         }
2429
2430  static int
2431 quorem
2432 #ifdef KR_headers
2433         (b, S) Bigint *b, *S;
2434 #else
2435         (Bigint *b, Bigint *S)
2436 #endif
2437 {
2438         int n;
2439         ULong *bx, *bxe, q, *sx, *sxe;
2440 #ifdef ULLong
2441         ULLong borrow, carry, y, ys;
2442 #else
2443         ULong borrow, carry, y, ys;
2444 #ifdef Pack_32
2445         ULong si, z, zs;
2446 #endif
2447 #endif
2448
2449         n = S->wds;
2450 #ifdef DEBUG
2451         /*debug*/ if (b->wds > n)
2452         /*debug*/       Bug("oversize b in quorem");
2453 #endif
2454         if (b->wds < n)
2455                 return 0;
2456         sx = S->x;
2457         sxe = sx + --n;
2458         bx = b->x;
2459         bxe = bx + n;
2460         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
2461 #ifdef DEBUG
2462         /*debug*/ if (q > 9)
2463         /*debug*/       Bug("oversized quotient in quorem");
2464 #endif
2465         if (q) {
2466                 borrow = 0;
2467                 carry = 0;
2468                 do {
2469 #ifdef ULLong
2470                         ys = *sx++ * (ULLong)q + carry;
2471                         carry = ys >> 32;
2472                         y = *bx - (ys & FFFFFFFF) - borrow;
2473                         borrow = y >> 32 & (ULong)1;
2474                         *bx++ = y & FFFFFFFF;
2475 #else
2476 #ifdef Pack_32
2477                         si = *sx++;
2478                         ys = (si & 0xffff) * q + carry;
2479                         zs = (si >> 16) * q + (ys >> 16);
2480                         carry = zs >> 16;
2481                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2482                         borrow = (y & 0x10000) >> 16;
2483                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
2484                         borrow = (z & 0x10000) >> 16;
2485                         Storeinc(bx, z, y);
2486 #else
2487                         ys = *sx++ * q + carry;
2488                         carry = ys >> 16;
2489                         y = *bx - (ys & 0xffff) - borrow;
2490                         borrow = (y & 0x10000) >> 16;
2491                         *bx++ = y & 0xffff;
2492 #endif
2493 #endif
2494                         }
2495                         while(sx <= sxe);
2496                 if (!*bxe) {
2497                         bx = b->x;
2498                         while(--bxe > bx && !*bxe)
2499                                 --n;
2500                         b->wds = n;
2501                         }
2502                 }
2503         if (cmp(b, S) >= 0) {
2504                 q++;
2505                 borrow = 0;
2506                 carry = 0;
2507                 bx = b->x;
2508                 sx = S->x;
2509                 do {
2510 #ifdef ULLong
2511                         ys = *sx++ + carry;
2512                         carry = ys >> 32;
2513                         y = *bx - (ys & FFFFFFFF) - borrow;
2514                         borrow = y >> 32 & (ULong)1;
2515                         *bx++ = y & FFFFFFFF;
2516 #else
2517 #ifdef Pack_32
2518                         si = *sx++;
2519                         ys = (si & 0xffff) + carry;
2520                         zs = (si >> 16) + (ys >> 16);
2521                         carry = zs >> 16;
2522                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2523                         borrow = (y & 0x10000) >> 16;
2524                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
2525                         borrow = (z & 0x10000) >> 16;
2526                         Storeinc(bx, z, y);
2527 #else
2528                         ys = *sx++ + carry;
2529                         carry = ys >> 16;
2530                         y = *bx - (ys & 0xffff) - borrow;
2531                         borrow = (y & 0x10000) >> 16;
2532                         *bx++ = y & 0xffff;
2533 #endif
2534 #endif
2535                         }
2536                         while(sx <= sxe);
2537                 bx = b->x;
2538                 bxe = bx + n;
2539                 if (!*bxe) {
2540                         while(--bxe > bx && !*bxe)
2541                                 --n;
2542                         b->wds = n;
2543                         }
2544                 }
2545         return q;
2546         }
2547
2548 #ifndef MULTIPLE_THREADS
2549  static char *dtoa_result;
2550 #endif
2551
2552  static char *
2553 #ifdef KR_headers
2554 rv_alloc(i) int i;
2555 #else
2556 rv_alloc(int i)
2557 #endif
2558 {
2559         int j, k, *r;
2560
2561         j = sizeof(ULong);
2562         for(k = 0;
2563                 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
2564                 j <<= 1)
2565                         k++;
2566         r = (int*)Balloc(k);
2567         *r = k;
2568         return
2569 #ifndef MULTIPLE_THREADS
2570         dtoa_result =
2571 #endif
2572                 (char *)(r+1);
2573         }
2574
2575  static char *
2576 #ifdef KR_headers
2577 nrv_alloc(s, rve, n) char *s, **rve; int n;
2578 #else
2579 nrv_alloc(char *s, char **rve, int n)
2580 #endif
2581 {
2582         char *rv, *t;
2583
2584         t = rv = rv_alloc(n);
2585         while((*t = *s++)) t++;
2586         if (rve)
2587                 *rve = t;
2588         return rv;
2589         }
2590
2591 /* freedtoa(s) must be used to free values s returned by dtoa
2592  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
2593  * but for consistency with earlier versions of dtoa, it is optional
2594  * when MULTIPLE_THREADS is not defined.
2595  */
2596
2597 void freedtoa (char *s);
2598
2599  void
2600 #ifdef KR_headers
2601 freedtoa(s) char *s;
2602 #else
2603 freedtoa(char *s)
2604 #endif
2605 {
2606         Bigint *b = (Bigint *)((int *)s - 1);
2607         b->maxwds = 1 << (b->k = *(int*)b);
2608         Bfree(b);
2609 #ifndef MULTIPLE_THREADS
2610         if (s == dtoa_result)
2611                 dtoa_result = 0;
2612 #endif
2613         }
2614
2615 #if 0
2616 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2617  *
2618  * Inspired by "How to Print Floating-Point Numbers Accurately" by
2619  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2620  *
2621  * Modifications:
2622  *      1. Rather than iterating, we use a simple numeric overestimate
2623  *         to determine k = floor(log10(d)).  We scale relevant
2624  *         quantities using O(log2(k)) rather than O(k) multiplications.
2625  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2626  *         try to generate digits strictly left to right.  Instead, we
2627  *         compute with fewer bits and propagate the carry if necessary
2628  *         when rounding the final digit up.  This is often faster.
2629  *      3. Under the assumption that input will be rounded nearest,
2630  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2631  *         That is, we allow equality in stopping tests when the
2632  *         round-nearest rule will give the same floating-point value
2633  *         as would satisfaction of the stopping test with strict
2634  *         inequality.
2635  *      4. We remove common factors of powers of 2 from relevant
2636  *         quantities.
2637  *      5. When converting floating-point integers less than 1e16,
2638  *         we use floating-point arithmetic rather than resorting
2639  *         to multiple-precision integers.
2640  *      6. When asked to produce fewer than 15 digits, we first try
2641  *         to get by with floating-point arithmetic; we resort to
2642  *         multiple-precision integer arithmetic only if we cannot
2643  *         guarantee that the floating-point calculation has given
2644  *         the correctly rounded result.  For k requested digits and
2645  *         "uniformly" distributed input, the probability is
2646  *         something like 10^(k-15) that we must resort to the Long
2647  *         calculation.
2648  */
2649
2650  char *
2651 dtoa
2652 #ifdef KR_headers
2653         (d, mode, ndigits, decpt, sign, rve)
2654         double d; int mode, ndigits, *decpt, *sign; char **rve;
2655 #else
2656         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2657 #endif
2658 {
2659  /*     Arguments ndigits, decpt, sign are similar to those
2660         of ecvt and fcvt; trailing zeros are suppressed from
2661         the returned string.  If not null, *rve is set to point
2662         to the end of the return value.  If d is +-Infinity or NaN,
2663         then *decpt is set to 9999.
2664
2665         mode:
2666                 0 ==> shortest string that yields d when read in
2667                         and rounded to nearest.
2668                 1 ==> like 0, but with Steele & White stopping rule;
2669                         e.g. with IEEE P754 arithmetic , mode 0 gives
2670                         1e23 whereas mode 1 gives 9.999999999999999e22.
2671                 2 ==> max(1,ndigits) significant digits.  This gives a
2672                         return value similar to that of ecvt, except
2673                         that trailing zeros are suppressed.
2674                 3 ==> through ndigits past the decimal point.  This
2675                         gives a return value similar to that from fcvt,
2676                         except that trailing zeros are suppressed, and
2677                         ndigits can be negative.
2678                 4,5 ==> similar to 2 and 3, respectively, but (in
2679                         round-nearest mode) with the tests of mode 0 to
2680                         possibly return a shorter string that rounds to d.
2681                         With IEEE arithmetic and compilation with
2682                         -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2683                         as modes 2 and 3 when FLT_ROUNDS != 1.
2684                 6-9 ==> Debugging modes similar to mode - 4:  don't try
2685                         fast floating-point estimate (if applicable).
2686
2687                 Values of mode other than 0-9 are treated as mode 0.
2688
2689                 Sufficient space is allocated to the return value
2690                 to hold the suppressed trailing zeros.
2691         */
2692
2693         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2694                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2695                 spec_case, try_quick;
2696         Long L;
2697 #ifndef Sudden_Underflow
2698         int denorm;
2699         ULong x;
2700 #endif
2701         Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2702         double d2, ds, eps;
2703         char *s, *s0;
2704 #ifdef Honor_FLT_ROUNDS
2705         int rounding;
2706 #endif
2707 #ifdef SET_INEXACT
2708         int inexact, oldinexact;
2709 #endif
2710
2711 #ifndef MULTIPLE_THREADS
2712         if (dtoa_result) {
2713                 freedtoa(dtoa_result);
2714                 dtoa_result = 0;
2715                 }
2716 #endif
2717
2718         if (word0(d) & Sign_bit) {
2719                 /* set sign for everything, including 0's and NaNs */
2720                 *sign = 1;
2721                 word0(d) &= ~Sign_bit;  /* clear sign bit */
2722                 }
2723         else
2724                 *sign = 0;
2725
2726 #if defined(IEEE_Arith) + defined(VAX)
2727 #ifdef IEEE_Arith
2728         if ((word0(d) & Exp_mask) == Exp_mask)
2729 #else
2730         if (word0(d)  == 0x8000)
2731 #endif
2732                 {
2733                 /* Infinity or NaN */
2734                 *decpt = 9999;
2735 #ifdef IEEE_Arith
2736                 if (!word1(d) && !(word0(d) & 0xfffff))
2737                         return nrv_alloc("Infinity", rve, 8);
2738 #endif
2739                 return nrv_alloc("NaN", rve, 3);
2740                 }
2741 #endif
2742 #ifdef IBM
2743         dval(d) += 0; /* normalize */
2744 #endif
2745         if (!dval(d)) {
2746                 *decpt = 1;
2747                 return nrv_alloc("0", rve, 1);
2748                 }
2749
2750 #ifdef SET_INEXACT
2751         try_quick = oldinexact = get_inexact();
2752         inexact = 1;
2753 #endif
2754 #ifdef Honor_FLT_ROUNDS
2755         if ((rounding = Flt_Rounds) >= 2) {
2756                 if (*sign)
2757                         rounding = rounding == 2 ? 0 : 2;
2758                 else
2759                         if (rounding != 2)
2760                                 rounding = 0;
2761                 }
2762 #endif
2763
2764         b = d2b(dval(d), &be, &bbits);
2765 #ifdef Sudden_Underflow
2766         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2767 #else
2768         if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
2769 #endif
2770                 dval(d2) = dval(d);
2771                 word0(d2) &= Frac_mask1;
2772                 word0(d2) |= Exp_11;
2773 #ifdef IBM
2774                 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2775                         dval(d2) /= 1 << j;
2776 #endif
2777
2778                 /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
2779                  * log10(x)      =  log(x) / log(10)
2780                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2781                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2782                  *
2783                  * This suggests computing an approximation k to log10(d) by
2784                  *
2785                  * k = (i - Bias)*0.301029995663981
2786                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2787                  *
2788                  * We want k to be too large rather than too small.
2789                  * The error in the first-order Taylor series approximation
2790                  * is in our favor, so we just round up the constant enough
2791                  * to compensate for any error in the multiplication of
2792                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2793                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2794                  * adding 1e-13 to the constant term more than suffices.
2795                  * Hence we adjust the constant term to 0.1760912590558.
2796                  * (We could get a more accurate k by invoking log10,
2797                  *  but this is probably not worthwhile.)
2798                  */
2799
2800                 i -= Bias;
2801 #ifdef IBM
2802                 i <<= 2;
2803                 i += j;
2804 #endif
2805 #ifndef Sudden_Underflow
2806                 denorm = 0;
2807                 }
2808         else {
2809                 /* d is denormalized */
2810
2811                 i = bbits + be + (Bias + (P-1) - 1);
2812                 x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
2813                             : word1(d) << 32 - i;
2814                 dval(d2) = x;
2815                 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2816                 i -= (Bias + (P-1) - 1) + 1;
2817                 denorm = 1;
2818                 }
2819 #endif
2820         ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2821         k = (int)ds;
2822         if (ds < 0. && ds != k)
2823                 k--;    /* want k = floor(ds) */
2824         k_check = 1;
2825         if (k >= 0 && k <= Ten_pmax) {
2826                 if (dval(d) < tens[k])
2827                         k--;
2828                 k_check = 0;
2829                 }
2830         j = bbits - i - 1;
2831         if (j >= 0) {
2832                 b2 = 0;
2833                 s2 = j;
2834                 }
2835         else {
2836                 b2 = -j;
2837                 s2 = 0;
2838                 }
2839         if (k >= 0) {
2840                 b5 = 0;
2841                 s5 = k;
2842                 s2 += k;
2843                 }
2844         else {
2845                 b2 -= k;
2846                 b5 = -k;
2847                 s5 = 0;
2848                 }
2849         if (mode < 0 || mode > 9)
2850                 mode = 0;
2851
2852 #ifndef SET_INEXACT
2853 #ifdef Check_FLT_ROUNDS
2854         try_quick = Rounding == 1;
2855 #else
2856         try_quick = 1;
2857 #endif
2858 #endif /*SET_INEXACT*/
2859
2860         if (mode > 5) {
2861                 mode -= 4;
2862                 try_quick = 0;
2863                 }
2864         leftright = 1;
2865         switch(mode) {
2866                 case 0:
2867                 case 1:
2868                         ilim = ilim1 = -1;
2869                         i = 18;
2870                         ndigits = 0;
2871                         break;
2872                 case 2:
2873                         leftright = 0;
2874                         /* no break */
2875                 case 4:
2876                         if (ndigits <= 0)
2877                                 ndigits = 1;
2878                         ilim = ilim1 = i = ndigits;
2879                         break;
2880                 case 3:
2881                         leftright = 0;
2882                         /* no break */
2883                 case 5:
2884                         i = ndigits + k + 1;
2885                         ilim = i;
2886                         ilim1 = i - 1;
2887                         if (i <= 0)
2888                                 i = 1;
2889                 }
2890         s = s0 = rv_alloc(i);
2891
2892 #ifdef Honor_FLT_ROUNDS
2893         if (mode > 1 && rounding != 1)
2894                 leftright = 0;
2895 #endif
2896
2897         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2898
2899                 /* Try to get by with floating-point arithmetic. */
2900
2901                 i = 0;
2902                 dval(d2) = dval(d);
2903                 k0 = k;
2904                 ilim0 = ilim;
2905                 ieps = 2; /* conservative */
2906                 if (k > 0) {
2907                         ds = tens[k&0xf];
2908                         j = k >> 4;
2909                         if (j & Bletch) {
2910                                 /* prevent overflows */
2911                                 j &= Bletch - 1;
2912                                 dval(d) /= bigtens[n_bigtens-1];
2913                                 ieps++;
2914                                 }
2915                         for(; j; j >>= 1, i++)
2916                                 if (j & 1) {
2917                                         ieps++;
2918                                         ds *= bigtens[i];
2919                                         }
2920                         dval(d) /= ds;
2921                         }
2922                 else if (j1 = -k) {
2923                         dval(d) *= tens[j1 & 0xf];
2924                         for(j = j1 >> 4; j; j >>= 1, i++)
2925                                 if (j & 1) {
2926                                         ieps++;
2927                                         dval(d) *= bigtens[i];
2928                                         }
2929                         }
2930                 if (k_check && dval(d) < 1. && ilim > 0) {
2931                         if (ilim1 <= 0)
2932                                 goto fast_failed;
2933                         ilim = ilim1;
2934                         k--;
2935                         dval(d) *= 10.;
2936                         ieps++;
2937                         }
2938                 dval(eps) = ieps*dval(d) + 7.;
2939                 word0(eps) -= (P-1)*Exp_msk1;
2940                 if (ilim == 0) {
2941                         S = mhi = 0;
2942                         dval(d) -= 5.;
2943                         if (dval(d) > dval(eps))
2944                                 goto one_digit;
2945                         if (dval(d) < -dval(eps))
2946                                 goto no_digits;
2947                         goto fast_failed;
2948                         }
2949 #ifndef No_leftright
2950                 if (leftright) {
2951                         /* Use Steele & White method of only
2952                          * generating digits needed.
2953                          */
2954                         dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2955                         for(i = 0;;) {
2956                                 L = dval(d);
2957                                 dval(d) -= L;
2958                                 *s++ = '0' + (int)L;
2959                                 if (dval(d) < dval(eps))
2960                                         goto ret1;
2961                                 if (1. - dval(d) < dval(eps))
2962                                         goto bump_up;
2963                                 if (++i >= ilim)
2964                                         break;
2965                                 dval(eps) *= 10.;
2966                                 dval(d) *= 10.;
2967                                 }
2968                         }
2969                 else {
2970 #endif
2971                         /* Generate ilim digits, then fix them up. */
2972                         dval(eps) *= tens[ilim-1];
2973                         for(i = 1;; i++, dval(d) *= 10.) {
2974                                 L = (Long)(dval(d));
2975                                 if (!(dval(d) -= L))
2976                                         ilim = i;
2977                                 *s++ = '0' + (int)L;
2978                                 if (i == ilim) {
2979                                         if (dval(d) > 0.5 + dval(eps))
2980                                                 goto bump_up;
2981                                         else if (dval(d) < 0.5 - dval(eps)) {
2982                                                 while(*--s == '0');
2983                                                 s++;
2984                                                 goto ret1;
2985                                                 }
2986                                         break;
2987                                         }
2988                                 }
2989 #ifndef No_leftright
2990                         }
2991 #endif
2992  fast_failed:
2993                 s = s0;
2994                 dval(d) = dval(d2);
2995                 k = k0;
2996                 ilim = ilim0;
2997                 }
2998
2999         /* Do we have a "small" integer? */
3000
3001         if (be >= 0 && k <= Int_max) {
3002                 /* Yes. */
3003                 ds = tens[k];
3004                 if (ndigits < 0 && ilim <= 0) {
3005                         S = mhi = 0;
3006                         if (ilim < 0 || dval(d) <= 5*ds)
3007                                 goto no_digits;
3008                         goto one_digit;
3009                         }
3010                 for(i = 1;; i++, dval(d) *= 10.) {
3011                         L = (Long)(dval(d) / ds);
3012                         dval(d) -= L*ds;
3013 #ifdef Check_FLT_ROUNDS
3014                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3015                         if (dval(d) < 0) {
3016                                 L--;
3017                                 dval(d) += ds;
3018                                 }
3019 #endif
3020                         *s++ = '0' + (int)L;
3021                         if (!dval(d)) {
3022 #ifdef SET_INEXACT
3023                                 inexact = 0;
3024 #endif
3025                                 break;
3026                                 }
3027                         if (i == ilim) {
3028 #ifdef Honor_FLT_ROUNDS
3029                                 if (mode > 1)
3030                                 switch(rounding) {
3031                                   case 0: goto ret1;
3032                                   case 2: goto bump_up;
3033                                   }
3034 #endif
3035                                 dval(d) += dval(d);
3036                                 if (dval(d) > ds || dval(d) == ds && L & 1) {
3037  bump_up:
3038                                         while(*--s == '9')
3039                                                 if (s == s0) {
3040                                                         k++;
3041                                                         *s = '0';
3042                                                         break;
3043                                                         }
3044                                         ++*s++;
3045                                         }
3046                                 break;
3047                                 }
3048                         }
3049                 goto ret1;
3050                 }
3051
3052         m2 = b2;
3053         m5 = b5;
3054         mhi = mlo = 0;
3055         if (leftright) {
3056                 i =
3057 #ifndef Sudden_Underflow
3058                         denorm ? be + (Bias + (P-1) - 1 + 1) :
3059 #endif
3060 #ifdef IBM
3061                         1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3062 #else
3063                         1 + P - bbits;
3064 #endif
3065                 b2 += i;
3066                 s2 += i;
3067                 mhi = i2b(1);
3068                 }
3069         if (m2 > 0 && s2 > 0) {
3070                 i = m2 < s2 ? m2 : s2;
3071                 b2 -= i;
3072                 m2 -= i;
3073                 s2 -= i;
3074                 }
3075         if (b5 > 0) {
3076                 if (leftright) {
3077                         if (m5 > 0) {
3078                                 mhi = pow5mult(mhi, m5);
3079                                 b1 = mult(mhi, b);
3080                                 Bfree(b);
3081                                 b = b1;
3082                                 }
3083                         if (j = b5 - m5)
3084                                 b = pow5mult(b, j);
3085                         }
3086                 else
3087                         b = pow5mult(b, b5);
3088                 }
3089         S = i2b(1);
3090         if (s5 > 0)
3091                 S = pow5mult(S, s5);
3092
3093         /* Check for special case that d is a normalized power of 2. */
3094
3095         spec_case = 0;
3096         if ((mode < 2 || leftright)
3097 #ifdef Honor_FLT_ROUNDS
3098                         && rounding == 1
3099 #endif
3100                                 ) {
3101                 if (!word1(d) && !(word0(d) & Bndry_mask)
3102 #ifndef Sudden_Underflow
3103                  && word0(d) & (Exp_mask & ~Exp_msk1)
3104 #endif
3105                                 ) {
3106                         /* The special case */
3107                         b2 += Log2P;
3108                         s2 += Log2P;
3109                         spec_case = 1;
3110                         }
3111                 }
3112
3113         /* Arrange for convenient computation of quotients:
3114          * shift left if necessary so divisor has 4 leading 0 bits.
3115          *
3116          * Perhaps we should just compute leading 28 bits of S once
3117          * and for all and pass them and a shift to quorem, so it
3118          * can do shifts and ors to compute the numerator for q.
3119          */
3120 #ifdef Pack_32
3121         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
3122                 i = 32 - i;
3123 #else
3124         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
3125                 i = 16 - i;
3126 #endif
3127         if (i > 4) {
3128                 i -= 4;
3129                 b2 += i;
3130                 m2 += i;
3131                 s2 += i;
3132                 }
3133         else if (i < 4) {
3134                 i += 28;
3135                 b2 += i;
3136                 m2 += i;
3137                 s2 += i;
3138                 }
3139         if (b2 > 0)
3140                 b = lshift(b, b2);
3141         if (s2 > 0)
3142                 S = lshift(S, s2);
3143         if (k_check) {
3144                 if (cmp(b,S) < 0) {
3145                         k--;
3146                         b = multadd(b, 10, 0);  /* we botched the k estimate */
3147                         if (leftright)
3148                                 mhi = multadd(mhi, 10, 0);
3149                         ilim = ilim1;
3150                         }
3151                 }
3152         if (ilim <= 0 && (mode == 3 || mode == 5)) {
3153                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3154                         /* no digits, fcvt style */
3155  no_digits:
3156                         k = -1 - ndigits;
3157                         goto ret;
3158                         }
3159  one_digit:
3160                 *s++ = '1';
3161                 k++;
3162                 goto ret;
3163                 }
3164         if (leftright) {
3165                 if (m2 > 0)
3166                         mhi = lshift(mhi, m2);
3167
3168                 /* Compute mlo -- check for special case
3169                  * that d is a normalized power of 2.
3170                  */
3171
3172                 mlo = mhi;
3173                 if (spec_case) {
3174                         mhi = Balloc(mhi->k);
3175                         Bcopy(mhi, mlo);
3176                         mhi = lshift(mhi, Log2P);
3177                         }
3178
3179                 for(i = 1;;i++) {
3180                         dig = quorem(b,S) + '0';
3181                         /* Do we yet have the shortest decimal string
3182                          * that will round to d?
3183                          */
3184                         j = cmp(b, mlo);
3185                         delta = diff(S, mhi);
3186                         j1 = delta->sign ? 1 : cmp(b, delta);
3187                         Bfree(delta);
3188 #ifndef ROUND_BIASED
3189                         if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3190 #ifdef Honor_FLT_ROUNDS
3191                                 && rounding >= 1
3192 #endif
3193                                                                    ) {
3194                                 if (dig == '9')
3195                                         goto round_9_up;
3196                                 if (j > 0)
3197                                         dig++;
3198 #ifdef SET_INEXACT
3199                                 else if (!b->x[0] && b->wds <= 1)
3200                                         inexact = 0;
3201 #endif
3202                                 *s++ = dig;
3203                                 goto ret;
3204                                 }
3205 #endif
3206                         if (j < 0 || j == 0 && mode != 1
3207 #ifndef ROUND_BIASED
3208                                                         && !(word1(d) & 1)
3209 #endif
3210                                         ) {
3211                                 if (!b->x[0] && b->wds <= 1) {
3212 #ifdef SET_INEXACT
3213                                         inexact = 0;
3214 #endif
3215                                         goto accept_dig;
3216                                         }
3217 #ifdef Honor_FLT_ROUNDS
3218                                 if (mode > 1)
3219                                  switch(rounding) {
3220                                   case 0: goto accept_dig;
3221                                   case 2: goto keep_dig;
3222                                   }
3223 #endif /*Honor_FLT_ROUNDS*/
3224                                 if (j1 > 0) {
3225                                         b = lshift(b, 1);
3226                                         j1 = cmp(b, S);
3227                                         if ((j1 > 0 || j1 == 0 && dig & 1)
3228                                         && dig++ == '9')
3229                                                 goto round_9_up;
3230                                         }
3231  accept_dig:
3232                                 *s++ = dig;
3233                                 goto ret;
3234                                 }
3235                         if (j1 > 0) {
3236 #ifdef Honor_FLT_ROUNDS
3237                                 if (!rounding)
3238                                         goto accept_dig;
3239 #endif
3240                                 if (dig == '9') { /* possible if i == 1 */
3241  round_9_up:
3242                                         *s++ = '9';
3243                                         goto roundoff;
3244                                         }
3245                                 *s++ = dig + 1;
3246                                 goto ret;
3247                                 }
3248 #ifdef Honor_FLT_ROUNDS
3249  keep_dig:
3250 #endif
3251                         *s++ = dig;
3252                         if (i == ilim)
3253                                 break;
3254                         b = multadd(b, 10, 0);
3255                         if (mlo == mhi)
3256                                 mlo = mhi = multadd(mhi, 10, 0);
3257                         else {
3258                                 mlo = multadd(mlo, 10, 0);
3259                                 mhi = multadd(mhi, 10, 0);
3260                                 }
3261                         }
3262                 }
3263         else
3264                 for(i = 1;; i++) {
3265                         *s++ = dig = quorem(b,S) + '0';
3266                         if (!b->x[0] && b->wds <= 1) {
3267 #ifdef SET_INEXACT
3268                                 inexact = 0;
3269 #endif
3270                                 goto ret;
3271                                 }
3272                         if (i >= ilim)
3273                                 break;
3274                         b = multadd(b, 10, 0);
3275                         }
3276
3277         /* Round off last digit */
3278
3279 #ifdef Honor_FLT_ROUNDS
3280         switch(rounding) {
3281           case 0: goto trimzeros;
3282           case 2: goto roundoff;
3283           }
3284 #endif
3285         b = lshift(b, 1);
3286         j = cmp(b, S);
3287         if (j > 0 || j == 0 && dig & 1) {
3288  roundoff:
3289                 while(*--s == '9')
3290                         if (s == s0) {
3291                                 k++;
3292                                 *s++ = '1';
3293                                 goto ret;
3294                                 }
3295                 ++*s++;
3296                 }
3297         else {
3298  trimzeros:
3299                 while(*--s == '0');
3300                 s++;
3301                 }
3302  ret:
3303         Bfree(S);
3304         if (mhi) {
3305                 if (mlo && mlo != mhi)
3306                         Bfree(mlo);
3307                 Bfree(mhi);
3308                 }
3309  ret1:
3310 #ifdef SET_INEXACT
3311         if (inexact) {
3312                 if (!oldinexact) {
3313                         word0(d) = Exp_1 + (70 << Exp_shift);
3314                         word1(d) = 0;
3315                         dval(d) += 1.;
3316                         }
3317                 }
3318         else if (!oldinexact)
3319                 clear_inexact();
3320 #endif
3321         Bfree(b);
3322         *s = 0;
3323         *decpt = k + 1;
3324         if (rve)
3325                 *rve = s;
3326         return s0;
3327         }
3328 #endif
3329
3330 #ifdef __cplusplus
3331 }
3332 #endif