js/src/dtoa.c

   1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
   2 /****************************************************************
   3  *
   4  * The author of this software is David M. Gay.
   5  *
   6  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
   7  *
   8  * Permission to use, copy, modify, and distribute this software for any
   9  * purpose without fee is hereby granted, provided that this entire notice
  10  * is included in all copies of any software which is or includes a copy
  11  * or modification of this software and in all copies of the supporting
  12  * documentation for such software.
  13  *
  14  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  15  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
  16  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  17  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  18  *
  19  ***************************************************************/
  20
  21 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  22  * with " at " changed at "@" and " dot " changed to ".").      */
  23
  24 /* On a machine with IEEE extended-precision registers, it is
  25  * necessary to specify double-precision (53-bit) rounding precision
  26  * before invoking strtod or dtoa.  If the machine uses (the equivalent
  27  * of) Intel 80x87 arithmetic, the call
  28  *      _control87(PC_53, MCW_PC);
  29  * does this with many compilers.  Whether this or another call is
  30  * appropriate depends on the compiler; for this to work, it may be
  31  * necessary to #include "float.h" or another system-dependent header
  32  * file.
  33  */
  34
  35 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  36  *
  37  * This strtod returns a nearest machine number to the input decimal
  38  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  39  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  40  * biased rounding (add half and chop).
  41  *
  42  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  43  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  44  *
  45  * Modifications:
  46  *
  47  *      1. We only require IEEE, IBM, or VAX double-precision
  48  *              arithmetic (not IEEE double-extended).
  49  *      2. We get by with floating-point arithmetic in a case that
  50  *              Clinger missed -- when we're computing d * 10^n
  51  *              for a small integer d and the integer n is not too
  52  *              much larger than 22 (the maximum integer k for which
  53  *              we can represent 10^k exactly), we may be able to
  54  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  55  *      3. Rather than a bit-at-a-time adjustment of the binary
  56  *              result in the hard case, we use floating-point
  57  *              arithmetic to determine the adjustment to within
  58  *              one bit; only in really hard cases do we need to
  59  *              compute a second residual.
  60  *      4. Because of 3., we don't need a large table of powers of 10
  61  *              for ten-to-e (just some small tables, e.g. of 10^k
  62  *              for 0 <= k <= 22).
  63  */
  64
  65 /*
  66  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  67  *      significant byte has the lowest address.
  68  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  69  *      significant byte has the lowest address.
  70  * #define Long int on machines with 32-bit ints and 64-bit longs.
  71  * #define IBM for IBM mainframe-style floating-point arithmetic.
  72  * #define VAX for VAX-style floating-point arithmetic (D_floating).
  73  * #define No_leftright to omit left-right logic in fast floating-point
  74  *      computation of dtoa.
  75  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  76  *      and strtod and dtoa should round accordingly.
  77  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  78  *      and Honor_FLT_ROUNDS is not #defined.
  79  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  80  *      that use extended-precision instructions to compute rounded
  81  *      products and quotients) with IBM.
  82  * #define ROUND_BIASED for IEEE-format with biased rounding.
  83  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  84  *      products but inaccurate quotients, e.g., for Intel i860.
  85  * #define NO_LONG_LONG on machines that do not have a "long long"
  86  *      integer type (of >= 64 bits).  On such machines, you can
  87  *      #define Just_16 to store 16 bits per 32-bit Long when doing
  88  *      high-precision integer arithmetic.  Whether this speeds things
  89  *      up or slows things down depends on the machine and the number
  90  *      being converted.  If long long is available and the name is
  91  *      something other than "long long", #define Llong to be the name,
  92  *      and if "unsigned Llong" does not work as an unsigned version of
  93  *      Llong, #define #ULLong to be the corresponding unsigned type.
  94  * #define KR_headers for old-style C function headers.
  95  * #define Bad_float_h if your system lacks a float.h or if it does not
  96  *      define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
  97  *      FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
  98  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
  99  *      if memory is available and otherwise does something you deem
 100  *      appropriate.  If MALLOC is undefined, malloc will be invoked
 101  *      directly -- and assumed always to succeed.  Similarly, if you
 102  *      want something other than the system's free() to be called to
 103  *      recycle memory acquired from MALLOC, #define FREE to be the
 104  *      name of the alternate routine.  (Unless you #define
 105  *      NO_GLOBAL_STATE and call destroydtoa, FREE or free is only
 106  *      called in pathological cases, e.g., in a dtoa call after a dtoa
 107  *      return in mode 3 with thousands of digits requested.)
 108  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
 109  *      memory allocations from a private pool of memory when possible.
 110  *      When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
 111  *      unless #defined to be a different length.  This default length
 112  *      suffices to get rid of MALLOC calls except for unusual cases,
 113  *      such as decimal-to-binary conversion of a very long string of
 114  *      digits.  The longest string dtoa can return is about 751 bytes
 115  *      long.  For conversions by strtod of strings of 800 digits and
 116  *      all dtoa conversions in single-threaded executions with 8-byte
 117  *      pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
 118  *      pointers, PRIVATE_MEM >= 7112 appears adequate.
 119  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
 120  *      multiple threads.  In this case, you must provide (or suitably
 121  *      #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
 122  *      by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
 123  *      in pow5mult, ensures lazy evaluation of only one copy of high
 124  *      powers of 5; omitting this lock would introduce a small
 125  *      probability of wasting memory, but would otherwise be harmless.)
 126  *      You must also invoke freedtoa(s) to free the value s returned by
 127  *      dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
 128  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
 129  *      avoids underflows on inputs whose result does not underflow.
 130  *      If you #define NO_IEEE_Scale on a machine that uses IEEE-format
 131  *      floating-point numbers and flushes underflows to zero rather
 132  *      than implementing gradual underflow, then you must also #define
 133  *      Sudden_Underflow.
 134  * #define USE_LOCALE to use the current locale's decimal_point value.
 135  * #define SET_INEXACT if IEEE arithmetic is being used and extra
 136  *      computation should be done to set the inexact flag when the
 137  *      result is inexact and avoid setting inexact when the result
 138  *      is exact.  In this case, dtoa.c must be compiled in
 139  *      an environment, perhaps provided by #include "dtoa.c" in a
 140  *      suitable wrapper, that defines two functions,
 141  *              int get_inexact(void);
 142  *              void clear_inexact(void);
 143  *      such that get_inexact() returns a nonzero value if the
 144  *      inexact bit is already set, and clear_inexact() sets the
 145  *      inexact bit to 0.  When SET_INEXACT is #defined, strtod
 146  *      also does extra computations to set the underflow and overflow
 147  *      flags when appropriate (i.e., when the result is tiny and
 148  *      inexact or when it is a numeric value rounded to +-infinity).
 149  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
 150  *      the result overflows to +-Infinity or underflows to 0.
 151  * #define NO_GLOBAL_STATE to avoid defining any non-const global or
 152  *      static variables. Instead the necessary state is stored in an
 153  *      opaque struct, DtoaState, a pointer to which must be passed to
 154  *      every entry point. Two new functions are added to the API:
 155  *              DtoaState *newdtoa(void);
 156  *              void destroydtoa(DtoaState *);
 157  */
 158
 159 #ifndef Long
 160 #define Long long
 161 #endif
 162 #ifndef ULong
 163 typedef unsigned Long ULong;
 164 #endif
 165
 166 #ifdef DEBUG
 167 #include <stdio.h>
 168 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
 169 #endif
 170
 171 #include <stdlib.h>
 172 #include <string.h>
 173
 174 #ifdef USE_LOCALE
 175 #include <locale.h>
 176 #endif
 177
 178 #ifdef MALLOC
 179 #ifdef KR_headers
 180 extern char *MALLOC();
 181 #else
 182 extern void *MALLOC(size_t);
 183 #endif
 184 #else
 185 #define MALLOC malloc
 186 #endif
 187
 188 #ifndef FREE
 189 #define FREE free
 190 #endif
 191
 192 #ifndef Omit_Private_Memory
 193 #ifndef PRIVATE_MEM
 194 #define PRIVATE_MEM 2304
 195 #endif
 196 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
 197 #endif
 198
 199 #undef IEEE_Arith
 200 #undef Avoid_Underflow
 201 #ifdef IEEE_MC68k
 202 #define IEEE_Arith
 203 #endif
 204 #ifdef IEEE_8087
 205 #define IEEE_Arith
 206 #endif
 207
 208 #include <errno.h>
 209
 210 #ifdef Bad_float_h
 211
 212 #ifdef IEEE_Arith
 213 #define DBL_DIG 15
 214 #define DBL_MAX_10_EXP 308
 215 #define DBL_MAX_EXP 1024
 216 #define FLT_RADIX 2
 217 #endif /*IEEE_Arith*/
 218
 219 #ifdef IBM
 220 #define DBL_DIG 16
 221 #define DBL_MAX_10_EXP 75
 222 #define DBL_MAX_EXP 63
 223 #define FLT_RADIX 16
 224 #define DBL_MAX 7.2370055773322621e+75
 225 #endif
 226
 227 #ifdef VAX
 228 #define DBL_DIG 16
 229 #define DBL_MAX_10_EXP 38
 230 #define DBL_MAX_EXP 127
 231 #define FLT_RADIX 2
 232 #define DBL_MAX 1.7014118346046923e+38
 233 #endif
 234
 235 #ifndef LONG_MAX
 236 #define LONG_MAX 2147483647
 237 #endif
 238
 239 #else /* ifndef Bad_float_h */
 240 #include <float.h>
 241 #endif /* Bad_float_h */
 242
 243 #ifndef __MATH_H__
 244 #include <math.h>
 245 #endif
 246
 247 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
 248 #error "Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined."
 249 #endif
 250
 251 typedef union { double d; ULong L[2]; } U;
 252
 253 #define dval(x) ((x).d)
 254 #ifdef IEEE_8087
 255 #define word0(x) ((x).L[1])
 256 #define word1(x) ((x).L[0])
 257 #else
 258 #define word0(x) ((x).L[0])
 259 #define word1(x) ((x).L[1])
 260 #endif
 261
 262 /* The following definition of Storeinc is appropriate for MIPS processors.
 263  * An alternative that might be better on some machines is
 264  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
 265  */
 266 #if defined(IEEE_8087) + defined(VAX)
 267 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
 268 ((unsigned short *)a)[0] = (unsigned short)c, a++)
 269 #else
 270 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
 271 ((unsigned short *)a)[1] = (unsigned short)c, a++)
 272 #endif
 273
 274 /* #define P DBL_MANT_DIG */
 275 /* Ten_pmax = floor(P*log(2)/log(5)) */
 276 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
 277 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
 278 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
 279
 280 #ifdef IEEE_Arith
 281 #define Exp_shift  20
 282 #define Exp_shift1 20
 283 #define Exp_msk1    0x100000
 284 #define Exp_msk11   0x100000
 285 #define Exp_mask  0x7ff00000
 286 #define P 53
 287 #define Bias 1023
 288 #define Emin (-1022)
 289 #define Exp_1  0x3ff00000
 290 #define Exp_11 0x3ff00000
 291 #define Ebits 11
 292 #define Frac_mask  0xfffff
 293 #define Frac_mask1 0xfffff
 294 #define Ten_pmax 22
 295 #define Bletch 0x10
 296 #define Bndry_mask  0xfffff
 297 #define Bndry_mask1 0xfffff
 298 #define LSB 1
 299 #define Sign_bit 0x80000000
 300 #define Log2P 1
 301 #define Tiny0 0
 302 #define Tiny1 1
 303 #define Quick_max 14
 304 #define Int_max 14
 305 #ifndef NO_IEEE_Scale
 306 #define Avoid_Underflow
 307 #ifdef Flush_Denorm     /* debugging option */
 308 #undef Sudden_Underflow
 309 #endif
 310 #endif
 311
 312 #ifndef Flt_Rounds
 313 #ifdef FLT_ROUNDS
 314 #define Flt_Rounds FLT_ROUNDS
 315 #else
 316 #define Flt_Rounds 1
 317 #endif
 318 #endif /*Flt_Rounds*/
 319
 320 #ifdef Honor_FLT_ROUNDS
 321 #define Rounding rounding
 322 #undef Check_FLT_ROUNDS
 323 #define Check_FLT_ROUNDS
 324 #else
 325 #define Rounding Flt_Rounds
 326 #endif
 327
 328 #else /* ifndef IEEE_Arith */
 329 #undef Check_FLT_ROUNDS
 330 #undef Honor_FLT_ROUNDS
 331 #undef SET_INEXACT
 332 #undef  Sudden_Underflow
 333 #define Sudden_Underflow
 334 #ifdef IBM
 335 #undef Flt_Rounds
 336 #define Flt_Rounds 0
 337 #define Exp_shift  24
 338 #define Exp_shift1 24
 339 #define Exp_msk1   0x1000000
 340 #define Exp_msk11  0x1000000
 341 #define Exp_mask  0x7f000000
 342 #define P 14
 343 #define Bias 65
 344 #define Exp_1  0x41000000
 345 #define Exp_11 0x41000000
 346 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
 347 #define Frac_mask  0xffffff
 348 #define Frac_mask1 0xffffff
 349 #define Bletch 4
 350 #define Ten_pmax 22
 351 #define Bndry_mask  0xefffff
 352 #define Bndry_mask1 0xffffff
 353 #define LSB 1
 354 #define Sign_bit 0x80000000
 355 #define Log2P 4
 356 #define Tiny0 0x100000
 357 #define Tiny1 0
 358 #define Quick_max 14
 359 #define Int_max 15
 360 #else /* VAX */
 361 #undef Flt_Rounds
 362 #define Flt_Rounds 1
 363 #define Exp_shift  23
 364 #define Exp_shift1 7
 365 #define Exp_msk1    0x80
 366 #define Exp_msk11   0x800000
 367 #define Exp_mask  0x7f80
 368 #define P 56
 369 #define Bias 129
 370 #define Exp_1  0x40800000
 371 #define Exp_11 0x4080
 372 #define Ebits 8
 373 #define Frac_mask  0x7fffff
 374 #define Frac_mask1 0xffff007f
 375 #define Ten_pmax 24
 376 #define Bletch 2
 377 #define Bndry_mask  0xffff007f
 378 #define Bndry_mask1 0xffff007f
 379 #define LSB 0x10000
 380 #define Sign_bit 0x8000
 381 #define Log2P 1
 382 #define Tiny0 0x80
 383 #define Tiny1 0
 384 #define Quick_max 15
 385 #define Int_max 15
 386 #endif /* IBM, VAX */
 387 #endif /* IEEE_Arith */
 388
 389 #ifndef IEEE_Arith
 390 #define ROUND_BIASED
 391 #endif
 392
 393 #ifdef RND_PRODQUOT
 394 #define rounded_product(a,b) a = rnd_prod(a, b)
 395 #define rounded_quotient(a,b) a = rnd_quot(a, b)
 396 #ifdef KR_headers
 397 extern double rnd_prod(), rnd_quot();
 398 #else
 399 extern double rnd_prod(double, double), rnd_quot(double, double);
 400 #endif
 401 #else
 402 #define rounded_product(a,b) a *= b
 403 #define rounded_quotient(a,b) a /= b
 404 #endif
 405
 406 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
 407 #define Big1 0xffffffff
 408
 409 #ifndef Pack_32
 410 #define Pack_32
 411 #endif
 412
 413 #ifdef KR_headers
 414 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
 415 #else
 416 #define FFFFFFFF 0xffffffffUL
 417 #endif
 418
 419 #ifdef NO_LONG_LONG
 420 #undef ULLong
 421 #ifdef Just_16
 422 #undef Pack_32
 423 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
 424  * This makes some inner loops simpler and sometimes saves work
 425  * during multiplications, but it often seems to make things slightly
 426  * slower.  Hence the default is now to store 32 bits per Long.
 427  */
 428 #endif
 429 #else   /* long long available */
 430 #ifndef Llong
 431 #define Llong long long
 432 #endif
 433 #ifndef ULLong
 434 #define ULLong unsigned Llong
 435 #endif
 436 #endif /* NO_LONG_LONG */
 437
 438 #ifndef MULTIPLE_THREADS
 439 #define ACQUIRE_DTOA_LOCK(n)    /*nothing*/
 440 #define FREE_DTOA_LOCK(n)       /*nothing*/
 441 #endif
 442
 443 #define Kmax 7
 444
 445  struct
 446 Bigint {
 447         struct Bigint *next;
 448         int k, maxwds, sign, wds;
 449         ULong x[1];
 450         };
 451
 452  typedef struct Bigint Bigint;
 453
 454 #ifdef NO_GLOBAL_STATE
 455 #ifdef MULTIPLE_THREADS
 456 #error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS"
 457 #endif
 458  struct
 459 DtoaState {
 460 #define DECLARE_GLOBAL_STATE  /* nothing */
 461 #else
 462 #define DECLARE_GLOBAL_STATE static
 463 #endif
 464
 465         DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1];
 466         DECLARE_GLOBAL_STATE Bigint *p5s;
 467 #ifndef Omit_Private_Memory
 468         DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem];
 469         DECLARE_GLOBAL_STATE double *pmem_next
 470 #ifndef NO_GLOBAL_STATE
 471                                                = private_mem
 472 #endif
 473                                                             ;
 474 #endif
 475 #ifdef NO_GLOBAL_STATE
 476         };
 477  typedef struct DtoaState DtoaState;
 478 #ifdef KR_headers
 479 #define STATE_PARAM state,
 480 #define STATE_PARAM_DECL DtoaState *state;
 481 #else
 482 #define STATE_PARAM DtoaState *state,
 483 #endif
 484 #define PASS_STATE state,
 485 #define GET_STATE(field) (state->field)
 486
 487  static DtoaState *
 488 newdtoa(void)
 489 {
 490         DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState));
 491         if (state) {
 492                 memset(state, 0, sizeof(DtoaState));
 493 #ifndef Omit_Private_Memory
 494                 state->pmem_next = state->private_mem;
 495 #endif
 496                 }
 497         return state;
 498 }
 499
 500  static void
 501 destroydtoa
 502 #ifdef KR_headers
 503         (state) STATE_PARAM_DECL
 504 #else
 505         (DtoaState *state)
 506 #endif
 507 {
 508         int i;
 509         Bigint *v, *next;
 510
 511         for (i = 0; i <= Kmax; i++) {
 512                 for (v = GET_STATE(freelist)[i]; v; v = next) {
 513                         next = v->next;
 514 #ifndef Omit_Private_Memory
 515                         if ((double*)v < GET_STATE(private_mem) ||
 516                             (double*)v >= GET_STATE(private_mem) + PRIVATE_mem)
 517 #endif
 518                                 FREE((void*)v);
 519                         }
 520                 }
 521 #ifdef Omit_Private_Memory
 522         Bigint* p5 = GET_STATE(p5s);
 523         while (p5) {
 524                 Bigint* tmp = p5;
 525                 p5 = p5->next;
 526                 FREE(tmp);
 527                 }
 528 #endif
 529         FREE((void *)state);
 530 }
 531
 532 #else
 533 #define STATE_PARAM      /* nothing */
 534 #define STATE_PARAM_DECL /* nothing */
 535 #define PASS_STATE       /* nothing */
 536 #define GET_STATE(name) name
 537 #endif
 538
 539  static Bigint *
 540 Balloc
 541 #ifdef KR_headers
 542         (STATE_PARAM k) STATE_PARAM_DECL int k;
 543 #else
 544         (STATE_PARAM int k)
 545 #endif
 546 {
 547         int x;
 548         Bigint *rv;
 549 #ifndef Omit_Private_Memory
 550         size_t len;
 551 #endif
 552
 553         ACQUIRE_DTOA_LOCK(0);
 554         /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
 555         /* but this case seems very unlikely. */
 556         if (k <= Kmax && (rv = GET_STATE(freelist)[k]))
 557                 GET_STATE(freelist)[k] = rv->next;
 558         else {
 559                 x = 1 << k;
 560 #ifdef Omit_Private_Memory
 561                 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
 562 #else
 563                 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
 564                         /sizeof(double);
 565                 if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) {
 566                         rv = (Bigint*)GET_STATE(pmem_next);
 567                         GET_STATE(pmem_next) += len;
 568                         }
 569                 else
 570                         rv = (Bigint*)MALLOC(len*sizeof(double));
 571 #endif
 572                 rv->k = k;
 573                 rv->maxwds = x;
 574                 }
 575         FREE_DTOA_LOCK(0);
 576         rv->sign = rv->wds = 0;
 577         return rv;
 578         }
 579
 580  static void
 581 Bfree
 582 #ifdef KR_headers
 583         (STATE_PARAM v) STATE_PARAM_DECL Bigint *v;
 584 #else
 585         (STATE_PARAM Bigint *v)
 586 #endif
 587 {
 588         if (v) {
 589                 if (v->k > Kmax)
 590                         FREE((void*)v);
 591                 else {
 592                         ACQUIRE_DTOA_LOCK(0);
 593                         v->next = GET_STATE(freelist)[v->k];
 594                         GET_STATE(freelist)[v->k] = v;
 595                         FREE_DTOA_LOCK(0);
 596                         }
 597                 }
 598         }
 599
 600 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
 601 y->wds*sizeof(Long) + 2*sizeof(int))
 602
 603  static Bigint *
 604 multadd
 605 #ifdef KR_headers
 606         (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a;
 607 #else
 608         (STATE_PARAM Bigint *b, int m, int a)   /* multiply by m and add a */
 609 #endif
 610 {
 611         int i, wds;
 612 #ifdef ULLong
 613         ULong *x;
 614         ULLong carry, y;
 615 #else
 616         ULong carry, *x, y;
 617 #ifdef Pack_32
 618         ULong xi, z;
 619 #endif
 620 #endif
 621         Bigint *b1;
 622
 623         wds = b->wds;
 624         x = b->x;
 625         i = 0;
 626         carry = a;
 627         do {
 628 #ifdef ULLong
 629                 y = *x * (ULLong)m + carry;
 630                 carry = y >> 32;
 631                 *x++ = (ULong) y & FFFFFFFF;
 632 #else
 633 #ifdef Pack_32
 634                 xi = *x;
 635                 y = (xi & 0xffff) * m + carry;
 636                 z = (xi >> 16) * m + (y >> 16);
 637                 carry = z >> 16;
 638                 *x++ = (z << 16) + (y & 0xffff);
 639 #else
 640                 y = *x * m + carry;
 641                 carry = y >> 16;
 642                 *x++ = y & 0xffff;
 643 #endif
 644 #endif
 645                 }
 646                 while(++i < wds);
 647         if (carry) {
 648                 if (wds >= b->maxwds) {
 649                         b1 = Balloc(PASS_STATE b->k+1);
 650                         Bcopy(b1, b);
 651                         Bfree(PASS_STATE b);
 652                         b = b1;
 653                         }
 654                 b->x[wds++] = (ULong) carry;
 655                 b->wds = wds;
 656                 }
 657         return b;
 658         }
 659
 660  static int
 661 hi0bits
 662 #ifdef KR_headers
 663         (x) ULong x;
 664 #else
 665         (ULong x)
 666 #endif
 667 {
 668         int k = 0;
 669
 670         if (!(x & 0xffff0000)) {
 671                 k = 16;
 672                 x <<= 16;
 673                 }
 674         if (!(x & 0xff000000)) {
 675                 k += 8;
 676                 x <<= 8;
 677                 }
 678         if (!(x & 0xf0000000)) {
 679                 k += 4;
 680                 x <<= 4;
 681                 }
 682         if (!(x & 0xc0000000)) {
 683                 k += 2;
 684                 x <<= 2;
 685                 }
 686         if (!(x & 0x80000000)) {
 687                 k++;
 688                 if (!(x & 0x40000000))
 689                         return 32;
 690                 }
 691         return k;
 692         }
 693
 694  static int
 695 lo0bits
 696 #ifdef KR_headers
 697         (y) ULong *y;
 698 #else
 699         (ULong *y)
 700 #endif
 701 {
 702         int k;
 703         ULong x = *y;
 704
 705         if (x & 7) {
 706                 if (x & 1)
 707                         return 0;
 708                 if (x & 2) {
 709                         *y = x >> 1;
 710                         return 1;
 711                         }
 712                 *y = x >> 2;
 713                 return 2;
 714                 }
 715         k = 0;
 716         if (!(x & 0xffff)) {
 717                 k = 16;
 718                 x >>= 16;
 719                 }
 720         if (!(x & 0xff)) {
 721                 k += 8;
 722                 x >>= 8;
 723                 }
 724         if (!(x & 0xf)) {
 725                 k += 4;
 726                 x >>= 4;
 727                 }
 728         if (!(x & 0x3)) {
 729                 k += 2;
 730                 x >>= 2;
 731                 }
 732         if (!(x & 1)) {
 733                 k++;
 734                 x >>= 1;
 735                 if (!x)
 736                         return 32;
 737                 }
 738         *y = x;
 739         return k;
 740         }
 741
 742  static Bigint *
 743 i2b
 744 #ifdef KR_headers
 745         (STATE_PARAM i) STATE_PARAM_DECL int i;
 746 #else
 747         (STATE_PARAM int i)
 748 #endif
 749 {
 750         Bigint *b;
 751
 752         b = Balloc(PASS_STATE 1);
 753         b->x[0] = i;
 754         b->wds = 1;
 755         return b;
 756         }
 757
 758  static Bigint *
 759 lshift
 760 #ifdef KR_headers
 761         (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
 762 #else
 763         (STATE_PARAM Bigint *b, int k)
 764 #endif
 765 {
 766         int i, k1, n, n1;
 767         Bigint *b1;
 768         ULong *x, *x1, *xe, z;
 769
 770 #ifdef Pack_32
 771         n = k >> 5;
 772 #else
 773         n = k >> 4;
 774 #endif
 775         k1 = b->k;
 776         n1 = n + b->wds + 1;
 777         for(i = b->maxwds; n1 > i; i <<= 1)
 778                 k1++;
 779         b1 = Balloc(PASS_STATE k1);
 780         x1 = b1->x;
 781         for(i = 0; i < n; i++)
 782                 *x1++ = 0;
 783         x = b->x;
 784         xe = x + b->wds;
 785 #ifdef Pack_32
 786         if (k &= 0x1f) {
 787                 k1 = 32 - k;
 788                 z = 0;
 789                 do {
 790                         *x1++ = *x << k | z;
 791                         z = *x++ >> k1;
 792                         }
 793                         while(x < xe);
 794                 if ((*x1 = z))
 795                         ++n1;
 796                 }
 797 #else
 798         if (k &= 0xf) {
 799                 k1 = 16 - k;
 800                 z = 0;
 801                 do {
 802                         *x1++ = *x << k  & 0xffff | z;
 803                         z = *x++ >> k1;
 804                         }
 805                         while(x < xe);
 806                 if (*x1 = z)
 807                         ++n1;
 808                 }
 809 #endif
 810         else do
 811                 *x1++ = *x++;
 812                 while(x < xe);
 813         b1->wds = n1 - 1;
 814         Bfree(PASS_STATE b);
 815         return b1;
 816         }
 817
 818  static int
 819 cmp
 820 #ifdef KR_headers
 821         (a, b) Bigint *a, *b;
 822 #else
 823         (Bigint *a, Bigint *b)
 824 #endif
 825 {
 826         ULong *xa, *xa0, *xb, *xb0;
 827         int i, j;
 828
 829         i = a->wds;
 830         j = b->wds;
 831 #ifdef DEBUG
 832         if (i > 1 && !a->x[i-1])
 833                 Bug("cmp called with a->x[a->wds-1] == 0");
 834         if (j > 1 && !b->x[j-1])
 835                 Bug("cmp called with b->x[b->wds-1] == 0");
 836 #endif
 837         if (i -= j)
 838                 return i;
 839         xa0 = a->x;
 840         xa = xa0 + j;
 841         xb0 = b->x;
 842         xb = xb0 + j;
 843         for(;;) {
 844                 if (*--xa != *--xb)
 845                         return *xa < *xb ? -1 : 1;
 846                 if (xa <= xa0)
 847                         break;
 848                 }
 849         return 0;
 850         }
 851
 852  static Bigint *
 853 diff
 854 #ifdef KR_headers
 855         (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
 856 #else
 857         (STATE_PARAM Bigint *a, Bigint *b)
 858 #endif
 859 {
 860         Bigint *c;
 861         int i, wa, wb;
 862         ULong *xa, *xae, *xb, *xbe, *xc;
 863 #ifdef ULLong
 864         ULLong borrow, y;
 865 #else
 866         ULong borrow, y;
 867 #ifdef Pack_32
 868         ULong z;
 869 #endif
 870 #endif
 871
 872         i = cmp(a,b);
 873         if (!i) {
 874                 c = Balloc(PASS_STATE 0);
 875                 c->wds = 1;
 876                 c->x[0] = 0;
 877                 return c;
 878                 }
 879         if (i < 0) {
 880                 c = a;
 881                 a = b;
 882                 b = c;
 883                 i = 1;
 884                 }
 885         else
 886                 i = 0;
 887         c = Balloc(PASS_STATE a->k);
 888         c->sign = i;
 889         wa = a->wds;
 890         xa = a->x;
 891         xae = xa + wa;
 892         wb = b->wds;
 893         xb = b->x;
 894         xbe = xb + wb;
 895         xc = c->x;
 896         borrow = 0;
 897 #ifdef ULLong
 898         do {
 899                 y = (ULLong)*xa++ - *xb++ - borrow;
 900                 borrow = y >> 32 & (ULong)1;
 901                 *xc++ = (ULong) y & FFFFFFFF;
 902                 }
 903                 while(xb < xbe);
 904         while(xa < xae) {
 905                 y = *xa++ - borrow;
 906                 borrow = y >> 32 & (ULong)1;
 907                 *xc++ = (ULong) y & FFFFFFFF;
 908                 }
 909 #else
 910 #ifdef Pack_32
 911         do {
 912                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
 913                 borrow = (y & 0x10000) >> 16;
 914                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
 915                 borrow = (z & 0x10000) >> 16;
 916                 Storeinc(xc, z, y);
 917                 }
 918                 while(xb < xbe);
 919         while(xa < xae) {
 920                 y = (*xa & 0xffff) - borrow;
 921                 borrow = (y & 0x10000) >> 16;
 922                 z = (*xa++ >> 16) - borrow;
 923                 borrow = (z & 0x10000) >> 16;
 924                 Storeinc(xc, z, y);
 925                 }
 926 #else
 927         do {
 928                 y = *xa++ - *xb++ - borrow;
 929                 borrow = (y & 0x10000) >> 16;
 930                 *xc++ = y & 0xffff;
 931                 }
 932                 while(xb < xbe);
 933         while(xa < xae) {
 934                 y = *xa++ - borrow;
 935                 borrow = (y & 0x10000) >> 16;
 936                 *xc++ = y & 0xffff;
 937                 }
 938 #endif
 939 #endif
 940         while(!*--xc)
 941                 wa--;
 942         c->wds = wa;
 943         return c;
 944         }
 945
 946  static Bigint *
 947 d2b
 948 #ifdef KR_headers
 949         (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits;
 950 #else
 951         (STATE_PARAM U d, int *e, int *bits)
 952 #endif
 953 {
 954         Bigint *b;
 955         int de, k;
 956         ULong *x, y, z;
 957 #ifndef Sudden_Underflow
 958         int i;
 959 #endif
 960 #ifdef VAX
 961         ULong d0, d1;
 962         d0 = word0(d) >> 16 | word0(d) << 16;
 963         d1 = word1(d) >> 16 | word1(d) << 16;
 964 #else
 965 #define d0 word0(d)
 966 #define d1 word1(d)
 967 #endif
 968
 969 #ifdef Pack_32
 970         b = Balloc(PASS_STATE 1);
 971 #else
 972         b = Balloc(PASS_STATE 2);
 973 #endif
 974         x = b->x;
 975
 976         z = d0 & Frac_mask;
 977         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
 978 #ifdef Sudden_Underflow
 979         de = (int)(d0 >> Exp_shift);
 980 #ifndef IBM
 981         z |= Exp_msk11;
 982 #endif
 983 #else
 984         if ((de = (int)(d0 >> Exp_shift)))
 985                 z |= Exp_msk1;
 986 #endif
 987 #ifdef Pack_32
 988         if ((y = d1)) {
 989                 if ((k = lo0bits(&y))) {
 990                         x[0] = y | z << (32 - k);
 991                         z >>= k;
 992                         }
 993                 else
 994                         x[0] = y;
 995 #ifndef Sudden_Underflow
 996                 i =
 997 #endif
 998                     b->wds = (x[1] = z) ? 2 : 1;
 999                 }
1000         else {
1001                 k = lo0bits(&z);
1002                 x[0] = z;
1003 #ifndef Sudden_Underflow
1004                 i =
1005 #endif
1006                     b->wds = 1;
1007                 k += 32;
1008                 }
1009 #else
1010         if (y = d1) {
1011                 if (k = lo0bits(&y))
1012                         if (k >= 16) {
1013                                 x[0] = y | z << 32 - k & 0xffff;
1014                                 x[1] = z >> k - 16 & 0xffff;
1015                                 x[2] = z >> k;
1016                                 i = 2;
1017                                 }
1018                         else {
1019                                 x[0] = y & 0xffff;
1020                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
1021                                 x[2] = z >> k & 0xffff;
1022                                 x[3] = z >> k+16;
1023                                 i = 3;
1024                                 }
1025                 else {
1026                         x[0] = y & 0xffff;
1027                         x[1] = y >> 16;
1028                         x[2] = z & 0xffff;
1029                         x[3] = z >> 16;
1030                         i = 3;
1031                         }
1032                 }
1033         else {
1034 #ifdef DEBUG
1035                 if (!z)
1036                         Bug("Zero passed to d2b");
1037 #endif
1038                 k = lo0bits(&z);
1039                 if (k >= 16) {
1040                         x[0] = z;
1041                         i = 0;
1042                         }
1043                 else {
1044                         x[0] = z & 0xffff;
1045                         x[1] = z >> 16;
1046                         i = 1;
1047                         }
1048                 k += 32;
1049                 }
1050         while(!x[i])
1051                 --i;
1052         b->wds = i + 1;
1053 #endif
1054 #ifndef Sudden_Underflow
1055         if (de) {
1056 #endif
1057 #ifdef IBM
1058                 *e = (de - Bias - (P-1) << 2) + k;
1059                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1060 #else
1061                 *e = de - Bias - (P-1) + k;
1062                 *bits = P - k;
1063 #endif
1064 #ifndef Sudden_Underflow
1065                 }
1066         else {
1067                 *e = de - Bias - (P-1) + 1 + k;
1068 #ifdef Pack_32
1069                 *bits = 32*i - hi0bits(x[i-1]);
1070 #else
1071                 *bits = (i+2)*16 - hi0bits(x[i]);
1072 #endif
1073                 }
1074 #endif
1075         return b;
1076         }
1077 #undef d0
1078 #undef d1