diff options
Diffstat (limited to 'js/src/dtoa.c')
| -rw-r--r-- | js/src/dtoa.c | 1088 |
1 files changed, 1088 insertions, 0 deletions
diff --git a/js/src/dtoa.c b/js/src/dtoa.c new file mode 100644 index 0000000000..08dec3dd19 --- /dev/null +++ b/js/src/dtoa.c @@ -0,0 +1,1088 @@ +/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */ +/**************************************************************** + * + * The author of this software is David M. Gay. + * + * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. + * + * Permission to use, copy, modify, and distribute this software for any + * purpose without fee is hereby granted, provided that this entire notice + * is included in all copies of any software which is or includes a copy + * or modification of this software and in all copies of the supporting + * documentation for such software. + * + * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED + * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY + * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY + * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. + * + ***************************************************************/ + +/* Please send bug reports to David M. Gay (dmg at acm dot org, + * with " at " changed at "@" and " dot " changed to "."). */ + +/* On a machine with IEEE extended-precision registers, it is + * necessary to specify double-precision (53-bit) rounding precision + * before invoking strtod or dtoa. If the machine uses (the equivalent + * of) Intel 80x87 arithmetic, the call + * _control87(PC_53, MCW_PC); + * does this with many compilers. Whether this or another call is + * appropriate depends on the compiler; for this to work, it may be + * necessary to #include "float.h" or another system-dependent header + * file. + */ + +/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. + * + * This strtod returns a nearest machine number to the input decimal + * string (or sets errno to ERANGE). With IEEE arithmetic, ties are + * broken by the IEEE round-even rule. Otherwise ties are broken by + * biased rounding (add half and chop). + * + * Inspired loosely by William D. Clinger's paper "How to Read Floating + * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. + * + * Modifications: + * + * 1. We only require IEEE, IBM, or VAX double-precision + * arithmetic (not IEEE double-extended). + * 2. We get by with floating-point arithmetic in a case that + * Clinger missed -- when we're computing d * 10^n + * for a small integer d and the integer n is not too + * much larger than 22 (the maximum integer k for which + * we can represent 10^k exactly), we may be able to + * compute (d*10^k) * 10^(e-k) with just one roundoff. + * 3. Rather than a bit-at-a-time adjustment of the binary + * result in the hard case, we use floating-point + * arithmetic to determine the adjustment to within + * one bit; only in really hard cases do we need to + * compute a second residual. + * 4. Because of 3., we don't need a large table of powers of 10 + * for ten-to-e (just some small tables, e.g. of 10^k + * for 0 <= k <= 22). + */ + +/* + * #define IEEE_8087 for IEEE-arithmetic machines where the least + * significant byte has the lowest address. + * #define IEEE_MC68k for IEEE-arithmetic machines where the most + * significant byte has the lowest address. + * #define Long int on machines with 32-bit ints and 64-bit longs. + * #define IBM for IBM mainframe-style floating-point arithmetic. + * #define VAX for VAX-style floating-point arithmetic (D_floating). + * #define No_leftright to omit left-right logic in fast floating-point + * computation of dtoa. + * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 + * and strtod and dtoa should round accordingly. + * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 + * and Honor_FLT_ROUNDS is not #defined. + * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines + * that use extended-precision instructions to compute rounded + * products and quotients) with IBM. + * #define ROUND_BIASED for IEEE-format with biased rounding. + * #define Inaccurate_Divide for IEEE-format with correctly rounded + * products but inaccurate quotients, e.g., for Intel i860. + * #define NO_LONG_LONG on machines that do not have a "long long" + * integer type (of >= 64 bits). On such machines, you can + * #define Just_16 to store 16 bits per 32-bit Long when doing + * high-precision integer arithmetic. Whether this speeds things + * up or slows things down depends on the machine and the number + * being converted. If long long is available and the name is + * something other than "long long", #define Llong to be the name, + * and if "unsigned Llong" does not work as an unsigned version of + * Llong, #define #ULLong to be the corresponding unsigned type. + * #define KR_headers for old-style C function headers. + * #define Bad_float_h if your system lacks a float.h or if it does not + * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, + * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. + * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) + * if memory is available and otherwise does something you deem + * appropriate. If MALLOC is undefined, malloc will be invoked + * directly -- and assumed always to succeed. Similarly, if you + * want something other than the system's free() to be called to + * recycle memory acquired from MALLOC, #define FREE to be the + * name of the alternate routine. (Unless you #define + * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only + * called in pathological cases, e.g., in a dtoa call after a dtoa + * return in mode 3 with thousands of digits requested.) + * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making + * memory allocations from a private pool of memory when possible. + * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, + * unless #defined to be a different length. This default length + * suffices to get rid of MALLOC calls except for unusual cases, + * such as decimal-to-binary conversion of a very long string of + * digits. The longest string dtoa can return is about 751 bytes + * long. For conversions by strtod of strings of 800 digits and + * all dtoa conversions in single-threaded executions with 8-byte + * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte + * pointers, PRIVATE_MEM >= 7112 appears adequate. + * #define MULTIPLE_THREADS if the system offers preemptively scheduled + * multiple threads. In this case, you must provide (or suitably + * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed + * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed + * in pow5mult, ensures lazy evaluation of only one copy of high + * powers of 5; omitting this lock would introduce a small + * probability of wasting memory, but would otherwise be harmless.) + * You must also invoke freedtoa(s) to free the value s returned by + * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. + * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that + * avoids underflows on inputs whose result does not underflow. + * If you #define NO_IEEE_Scale on a machine that uses IEEE-format + * floating-point numbers and flushes underflows to zero rather + * than implementing gradual underflow, then you must also #define + * Sudden_Underflow. + * #define USE_LOCALE to use the current locale's decimal_point value. + * #define SET_INEXACT if IEEE arithmetic is being used and extra + * computation should be done to set the inexact flag when the + * result is inexact and avoid setting inexact when the result + * is exact. In this case, dtoa.c must be compiled in + * an environment, perhaps provided by #include "dtoa.c" in a + * suitable wrapper, that defines two functions, + * int get_inexact(void); + * void clear_inexact(void); + * such that get_inexact() returns a nonzero value if the + * inexact bit is already set, and clear_inexact() sets the + * inexact bit to 0. When SET_INEXACT is #defined, strtod + * also does extra computations to set the underflow and overflow + * flags when appropriate (i.e., when the result is tiny and + * inexact or when it is a numeric value rounded to +-infinity). + * #define NO_ERRNO if strtod should not assign errno = ERANGE when + * the result overflows to +-Infinity or underflows to 0. + * #define NO_GLOBAL_STATE to avoid defining any non-const global or + * static variables. Instead the necessary state is stored in an + * opaque struct, DtoaState, a pointer to which must be passed to + * every entry point. Two new functions are added to the API: + * DtoaState *newdtoa(void); + * void destroydtoa(DtoaState *); + */ + +#ifndef Long +#define Long long +#endif +#ifndef ULong +typedef unsigned Long ULong; +#endif + +#ifdef DEBUG +#include <stdio.h> +#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} +#endif + +#include <stdlib.h> +#include <string.h> + +#ifdef USE_LOCALE +#include <locale.h> +#endif + +#ifdef MALLOC +#ifdef KR_headers +extern char *MALLOC(); +#else +extern void *MALLOC(size_t); +#endif +#else +#define MALLOC malloc +#endif + +#ifndef FREE +#define FREE free +#endif + +#ifndef Omit_Private_Memory +#ifndef PRIVATE_MEM +#define PRIVATE_MEM 2304 +#endif +#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) +#endif + +#undef IEEE_Arith +#undef Avoid_Underflow +#ifdef IEEE_MC68k +#define IEEE_Arith +#endif +#ifdef IEEE_8087 +#define IEEE_Arith +#endif + +#include <errno.h> + +#ifdef Bad_float_h + +#ifdef IEEE_Arith +#define DBL_DIG 15 +#define DBL_MAX_10_EXP 308 +#define DBL_MAX_EXP 1024 +#define FLT_RADIX 2 +#endif /*IEEE_Arith*/ + +#ifdef IBM +#define DBL_DIG 16 +#define DBL_MAX_10_EXP 75 +#define DBL_MAX_EXP 63 +#define FLT_RADIX 16 +#define DBL_MAX 7.2370055773322621e+75 +#endif + +#ifdef VAX +#define DBL_DIG 16 +#define DBL_MAX_10_EXP 38 +#define DBL_MAX_EXP 127 +#define FLT_RADIX 2 +#define DBL_MAX 1.7014118346046923e+38 +#endif + +#ifndef LONG_MAX +#define LONG_MAX 2147483647 +#endif + +#else /* ifndef Bad_float_h */ +#include <float.h> +#endif /* Bad_float_h */ + +#ifndef __MATH_H__ +#include <math.h> +#endif + +#ifndef CONST +#ifdef KR_headers +#define CONST /* blank */ +#else +#define CONST const +#endif +#endif + +#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 +#error "Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined." +#endif + +typedef union { double d; ULong L[2]; } U; + +#define dval(x) ((x).d) +#ifdef IEEE_8087 +#define word0(x) ((x).L[1]) +#define word1(x) ((x).L[0]) +#else +#define word0(x) ((x).L[0]) +#define word1(x) ((x).L[1]) +#endif + +/* The following definition of Storeinc is appropriate for MIPS processors. + * An alternative that might be better on some machines is + * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) + */ +#if defined(IEEE_8087) + defined(VAX) +#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ +((unsigned short *)a)[0] = (unsigned short)c, a++) +#else +#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ +((unsigned short *)a)[1] = (unsigned short)c, a++) +#endif + +/* #define P DBL_MANT_DIG */ +/* Ten_pmax = floor(P*log(2)/log(5)) */ +/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ +/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ +/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ + +#ifdef IEEE_Arith +#define Exp_shift 20 +#define Exp_shift1 20 +#define Exp_msk1 0x100000 +#define Exp_msk11 0x100000 +#define Exp_mask 0x7ff00000 +#define P 53 +#define Bias 1023 +#define Emin (-1022) +#define Exp_1 0x3ff00000 +#define Exp_11 0x3ff00000 +#define Ebits 11 +#define Frac_mask 0xfffff +#define Frac_mask1 0xfffff +#define Ten_pmax 22 +#define Bletch 0x10 +#define Bndry_mask 0xfffff +#define Bndry_mask1 0xfffff +#define LSB 1 +#define Sign_bit 0x80000000 +#define Log2P 1 +#define Tiny0 0 +#define Tiny1 1 +#define Quick_max 14 +#define Int_max 14 +#ifndef NO_IEEE_Scale +#define Avoid_Underflow +#ifdef Flush_Denorm /* debugging option */ +#undef Sudden_Underflow +#endif +#endif + +#ifndef Flt_Rounds +#ifdef FLT_ROUNDS +#define Flt_Rounds FLT_ROUNDS +#else +#define Flt_Rounds 1 +#endif +#endif /*Flt_Rounds*/ + +#ifdef Honor_FLT_ROUNDS +#define Rounding rounding +#undef Check_FLT_ROUNDS +#define Check_FLT_ROUNDS +#else +#define Rounding Flt_Rounds +#endif + +#else /* ifndef IEEE_Arith */ +#undef Check_FLT_ROUNDS +#undef Honor_FLT_ROUNDS +#undef SET_INEXACT +#undef Sudden_Underflow +#define Sudden_Underflow +#ifdef IBM +#undef Flt_Rounds +#define Flt_Rounds 0 +#define Exp_shift 24 +#define Exp_shift1 24 +#define Exp_msk1 0x1000000 +#define Exp_msk11 0x1000000 +#define Exp_mask 0x7f000000 +#define P 14 +#define Bias 65 +#define Exp_1 0x41000000 +#define Exp_11 0x41000000 +#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ +#define Frac_mask 0xffffff +#define Frac_mask1 0xffffff +#define Bletch 4 +#define Ten_pmax 22 +#define Bndry_mask 0xefffff +#define Bndry_mask1 0xffffff +#define LSB 1 +#define Sign_bit 0x80000000 +#define Log2P 4 +#define Tiny0 0x100000 +#define Tiny1 0 +#define Quick_max 14 +#define Int_max 15 +#else /* VAX */ +#undef Flt_Rounds +#define Flt_Rounds 1 +#define Exp_shift 23 +#define Exp_shift1 7 +#define Exp_msk1 0x80 +#define Exp_msk11 0x800000 +#define Exp_mask 0x7f80 +#define P 56 +#define Bias 129 +#define Exp_1 0x40800000 +#define Exp_11 0x4080 +#define Ebits 8 +#define Frac_mask 0x7fffff +#define Frac_mask1 0xffff007f +#define Ten_pmax 24 +#define Bletch 2 +#define Bndry_mask 0xffff007f +#define Bndry_mask1 0xffff007f +#define LSB 0x10000 +#define Sign_bit 0x8000 +#define Log2P 1 +#define Tiny0 0x80 +#define Tiny1 0 +#define Quick_max 15 +#define Int_max 15 +#endif /* IBM, VAX */ +#endif /* IEEE_Arith */ + +#ifndef IEEE_Arith +#define ROUND_BIASED +#endif + +#ifdef RND_PRODQUOT +#define rounded_product(a,b) a = rnd_prod(a, b) +#define rounded_quotient(a,b) a = rnd_quot(a, b) +#ifdef KR_headers +extern double rnd_prod(), rnd_quot(); +#else +extern double rnd_prod(double, double), rnd_quot(double, double); +#endif +#else +#define rounded_product(a,b) a *= b +#define rounded_quotient(a,b) a /= b +#endif + +#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) +#define Big1 0xffffffff + +#ifndef Pack_32 +#define Pack_32 +#endif + +#ifdef KR_headers +#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff) +#else +#define FFFFFFFF 0xffffffffUL +#endif + +#ifdef NO_LONG_LONG +#undef ULLong +#ifdef Just_16 +#undef Pack_32 +/* When Pack_32 is not defined, we store 16 bits per 32-bit Long. + * This makes some inner loops simpler and sometimes saves work + * during multiplications, but it often seems to make things slightly + * slower. Hence the default is now to store 32 bits per Long. + */ +#endif +#else /* long long available */ +#ifndef Llong +#define Llong long long +#endif +#ifndef ULLong +#define ULLong unsigned Llong +#endif +#endif /* NO_LONG_LONG */ + +#ifndef MULTIPLE_THREADS +#define ACQUIRE_DTOA_LOCK(n) /*nothing*/ +#define FREE_DTOA_LOCK(n) /*nothing*/ +#endif + +#define Kmax 7 + + struct +Bigint { + struct Bigint *next; + int k, maxwds, sign, wds; + ULong x[1]; + }; + + typedef struct Bigint Bigint; + +#ifdef NO_GLOBAL_STATE +#ifdef MULTIPLE_THREADS +#error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS" +#endif + struct +DtoaState { +#define DECLARE_GLOBAL_STATE /* nothing */ +#else +#define DECLARE_GLOBAL_STATE static +#endif + + DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1]; + DECLARE_GLOBAL_STATE Bigint *p5s; +#ifndef Omit_Private_Memory + DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem]; + DECLARE_GLOBAL_STATE double *pmem_next +#ifndef NO_GLOBAL_STATE + = private_mem +#endif + ; +#endif +#ifdef NO_GLOBAL_STATE + }; + typedef struct DtoaState DtoaState; +#ifdef KR_headers +#define STATE_PARAM state, +#define STATE_PARAM_DECL DtoaState *state; +#else +#define STATE_PARAM DtoaState *state, +#endif +#define PASS_STATE state, +#define GET_STATE(field) (state->field) + + static DtoaState * +newdtoa(void) +{ + DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState)); + if (state) { + memset(state, 0, sizeof(DtoaState)); +#ifndef Omit_Private_Memory + state->pmem_next = state->private_mem; +#endif + } + return state; +} + + static void +destroydtoa +#ifdef KR_headers + (state) STATE_PARAM_DECL +#else + (DtoaState *state) +#endif +{ + int i; + Bigint *v, *next; + + for (i = 0; i <= Kmax; i++) { + for (v = GET_STATE(freelist)[i]; v; v = next) { + next = v->next; +#ifndef Omit_Private_Memory + if ((double*)v < GET_STATE(private_mem) || + (double*)v >= GET_STATE(private_mem) + PRIVATE_mem) +#endif + FREE((void*)v); + } + } +#ifdef Omit_Private_Memory + Bigint* p5 = GET_STATE(p5s); + while (p5) { + Bigint* tmp = p5; + p5 = p5->next; + FREE(tmp); + } +#endif + FREE((void *)state); +} + +#else +#define STATE_PARAM /* nothing */ +#define STATE_PARAM_DECL /* nothing */ +#define PASS_STATE /* nothing */ +#define GET_STATE(name) name +#endif + + static Bigint * +Balloc +#ifdef KR_headers + (STATE_PARAM k) STATE_PARAM_DECL int k; +#else + (STATE_PARAM int k) +#endif +{ + int x; + Bigint *rv; +#ifndef Omit_Private_Memory + size_t len; +#endif + + ACQUIRE_DTOA_LOCK(0); + /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ + /* but this case seems very unlikely. */ + if (k <= Kmax && (rv = GET_STATE(freelist)[k])) + GET_STATE(freelist)[k] = rv->next; + else { + x = 1 << k; +#ifdef Omit_Private_Memory + rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); +#else + len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) + /sizeof(double); + if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) { + rv = (Bigint*)GET_STATE(pmem_next); + GET_STATE(pmem_next) += len; + } + else + rv = (Bigint*)MALLOC(len*sizeof(double)); +#endif + rv->k = k; + rv->maxwds = x; + } + FREE_DTOA_LOCK(0); + rv->sign = rv->wds = 0; + return rv; + } + + static void +Bfree +#ifdef KR_headers + (STATE_PARAM v) STATE_PARAM_DECL Bigint *v; +#else + (STATE_PARAM Bigint *v) +#endif +{ + if (v) { + if (v->k > Kmax) + FREE((void*)v); + else { + ACQUIRE_DTOA_LOCK(0); + v->next = GET_STATE(freelist)[v->k]; + GET_STATE(freelist)[v->k] = v; + FREE_DTOA_LOCK(0); + } + } + } + +#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ +y->wds*sizeof(Long) + 2*sizeof(int)) + + static Bigint * +multadd +#ifdef KR_headers + (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a; +#else + (STATE_PARAM Bigint *b, int m, int a) /* multiply by m and add a */ +#endif +{ + int i, wds; +#ifdef ULLong + ULong *x; + ULLong carry, y; +#else + ULong carry, *x, y; +#ifdef Pack_32 + ULong xi, z; +#endif +#endif + Bigint *b1; + + wds = b->wds; + x = b->x; + i = 0; + carry = a; + do { +#ifdef ULLong + y = *x * (ULLong)m + carry; + carry = y >> 32; + *x++ = (ULong) y & FFFFFFFF; +#else +#ifdef Pack_32 + xi = *x; + y = (xi & 0xffff) * m + carry; + z = (xi >> 16) * m + (y >> 16); + carry = z >> 16; + *x++ = (z << 16) + (y & 0xffff); +#else + y = *x * m + carry; + carry = y >> 16; + *x++ = y & 0xffff; +#endif +#endif + } + while(++i < wds); + if (carry) { + if (wds >= b->maxwds) { + b1 = Balloc(PASS_STATE b->k+1); + Bcopy(b1, b); + Bfree(PASS_STATE b); + b = b1; + } + b->x[wds++] = (ULong) carry; + b->wds = wds; + } + return b; + } + + static int +hi0bits +#ifdef KR_headers + (x) ULong x; +#else + (ULong x) +#endif +{ + int k = 0; + + if (!(x & 0xffff0000)) { + k = 16; + x <<= 16; + } + if (!(x & 0xff000000)) { + k += 8; + x <<= 8; + } + if (!(x & 0xf0000000)) { + k += 4; + x <<= 4; + } + if (!(x & 0xc0000000)) { + k += 2; + x <<= 2; + } + if (!(x & 0x80000000)) { + k++; + if (!(x & 0x40000000)) + return 32; + } + return k; + } + + static int +lo0bits +#ifdef KR_headers + (y) ULong *y; +#else + (ULong *y) +#endif +{ + int k; + ULong x = *y; + + if (x & 7) { + if (x & 1) + return 0; + if (x & 2) { + *y = x >> 1; + return 1; + } + *y = x >> 2; + return 2; + } + k = 0; + if (!(x & 0xffff)) { + k = 16; + x >>= 16; + } + if (!(x & 0xff)) { + k += 8; + x >>= 8; + } + if (!(x & 0xf)) { + k += 4; + x >>= 4; + } + if (!(x & 0x3)) { + k += 2; + x >>= 2; + } + if (!(x & 1)) { + k++; + x >>= 1; + if (!x) + return 32; + } + *y = x; + return k; + } + + static Bigint * +i2b +#ifdef KR_headers + (STATE_PARAM i) STATE_PARAM_DECL int i; +#else + (STATE_PARAM int i) +#endif +{ + Bigint *b; + + b = Balloc(PASS_STATE 1); + b->x[0] = i; + b->wds = 1; + return b; + } + + static Bigint * +lshift +#ifdef KR_headers + (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k; +#else + (STATE_PARAM Bigint *b, int k) +#endif +{ + int i, k1, n, n1; + Bigint *b1; + ULong *x, *x1, *xe, z; + +#ifdef Pack_32 + n = k >> 5; +#else + n = k >> 4; +#endif + k1 = b->k; + n1 = n + b->wds + 1; + for(i = b->maxwds; n1 > i; i <<= 1) + k1++; + b1 = Balloc(PASS_STATE k1); + x1 = b1->x; + for(i = 0; i < n; i++) + *x1++ = 0; + x = b->x; + xe = x + b->wds; +#ifdef Pack_32 + if (k &= 0x1f) { + k1 = 32 - k; + z = 0; + do { + *x1++ = *x << k | z; + z = *x++ >> k1; + } + while(x < xe); + if ((*x1 = z)) + ++n1; + } +#else + if (k &= 0xf) { + k1 = 16 - k; + z = 0; + do { + *x1++ = *x << k & 0xffff | z; + z = *x++ >> k1; + } + while(x < xe); + if (*x1 = z) + ++n1; + } +#endif + else do + *x1++ = *x++; + while(x < xe); + b1->wds = n1 - 1; + Bfree(PASS_STATE b); + return b1; + } + + static int +cmp +#ifdef KR_headers + (a, b) Bigint *a, *b; +#else + (Bigint *a, Bigint *b) +#endif +{ + ULong *xa, *xa0, *xb, *xb0; + int i, j; + + i = a->wds; + j = b->wds; +#ifdef DEBUG + if (i > 1 && !a->x[i-1]) + Bug("cmp called with a->x[a->wds-1] == 0"); + if (j > 1 && !b->x[j-1]) + Bug("cmp called with b->x[b->wds-1] == 0"); +#endif + if (i -= j) + return i; + xa0 = a->x; + xa = xa0 + j; + xb0 = b->x; + xb = xb0 + j; + for(;;) { + if (*--xa != *--xb) + return *xa < *xb ? -1 : 1; + if (xa <= xa0) + break; + } + return 0; + } + + static Bigint * +diff +#ifdef KR_headers + (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b; +#else + (STATE_PARAM Bigint *a, Bigint *b) +#endif +{ + Bigint *c; + int i, wa, wb; + ULong *xa, *xae, *xb, *xbe, *xc; +#ifdef ULLong + ULLong borrow, y; +#else + ULong borrow, y; +#ifdef Pack_32 + ULong z; +#endif +#endif + + i = cmp(a,b); + if (!i) { + c = Balloc(PASS_STATE 0); + c->wds = 1; + c->x[0] = 0; + return c; + } + if (i < 0) { + c = a; + a = b; + b = c; + i = 1; + } + else + i = 0; + c = Balloc(PASS_STATE a->k); + c->sign = i; + wa = a->wds; + xa = a->x; + xae = xa + wa; + wb = b->wds; + xb = b->x; + xbe = xb + wb; + xc = c->x; + borrow = 0; +#ifdef ULLong + do { + y = (ULLong)*xa++ - *xb++ - borrow; + borrow = y >> 32 & (ULong)1; + *xc++ = (ULong) y & FFFFFFFF; + } + while(xb < xbe); + while(xa < xae) { + y = *xa++ - borrow; + borrow = y >> 32 & (ULong)1; + *xc++ = (ULong) y & FFFFFFFF; + } +#else +#ifdef Pack_32 + do { + y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; + borrow = (y & 0x10000) >> 16; + z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; + borrow = (z & 0x10000) >> 16; + Storeinc(xc, z, y); + } + while(xb < xbe); + while(xa < xae) { + y = (*xa & 0xffff) - borrow; + borrow = (y & 0x10000) >> 16; + z = (*xa++ >> 16) - borrow; + borrow = (z & 0x10000) >> 16; + Storeinc(xc, z, y); + } +#else + do { + y = *xa++ - *xb++ - borrow; + borrow = (y & 0x10000) >> 16; + *xc++ = y & 0xffff; + } + while(xb < xbe); + while(xa < xae) { + y = *xa++ - borrow; + borrow = (y & 0x10000) >> 16; + *xc++ = y & 0xffff; + } +#endif +#endif + while(!*--xc) + wa--; + c->wds = wa; + return c; + } + + static Bigint * +d2b +#ifdef KR_headers + (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits; +#else + (STATE_PARAM U d, int *e, int *bits) +#endif +{ + Bigint *b; + int de, k; + ULong *x, y, z; +#ifndef Sudden_Underflow + int i; +#endif +#ifdef VAX + ULong d0, d1; + d0 = word0(d) >> 16 | word0(d) << 16; + d1 = word1(d) >> 16 | word1(d) << 16; +#else +#define d0 word0(d) +#define d1 word1(d) +#endif + +#ifdef Pack_32 + b = Balloc(PASS_STATE 1); +#else + b = Balloc(PASS_STATE 2); +#endif + x = b->x; + + z = d0 & Frac_mask; + d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ +#ifdef Sudden_Underflow + de = (int)(d0 >> Exp_shift); +#ifndef IBM + z |= Exp_msk11; +#endif +#else + if ((de = (int)(d0 >> Exp_shift))) + z |= Exp_msk1; +#endif +#ifdef Pack_32 + if ((y = d1)) { + if ((k = lo0bits(&y))) { + x[0] = y | z << (32 - k); + z >>= k; + } + else + x[0] = y; +#ifndef Sudden_Underflow + i = +#endif + b->wds = (x[1] = z) ? 2 : 1; + } + else { + k = lo0bits(&z); + x[0] = z; +#ifndef Sudden_Underflow + i = +#endif + b->wds = 1; + k += 32; + } +#else + if (y = d1) { + if (k = lo0bits(&y)) + if (k >= 16) { + x[0] = y | z << 32 - k & 0xffff; + x[1] = z >> k - 16 & 0xffff; + x[2] = z >> k; + i = 2; + } + else { + x[0] = y & 0xffff; + x[1] = y >> 16 | z << 16 - k & 0xffff; + x[2] = z >> k & 0xffff; + x[3] = z >> k+16; + i = 3; + } + else { + x[0] = y & 0xffff; + x[1] = y >> 16; + x[2] = z & 0xffff; + x[3] = z >> 16; + i = 3; + } + } + else { +#ifdef DEBUG + if (!z) + Bug("Zero passed to d2b"); +#endif + k = lo0bits(&z); + if (k >= 16) { + x[0] = z; + i = 0; + } + else { + x[0] = z & 0xffff; + x[1] = z >> 16; + i = 1; + } + k += 32; + } + while(!x[i]) + --i; + b->wds = i + 1; +#endif +#ifndef Sudden_Underflow + if (de) { +#endif +#ifdef IBM + *e = (de - Bias - (P-1) << 2) + k; + *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); +#else + *e = de - Bias - (P-1) + k; + *bits = P - k; +#endif +#ifndef Sudden_Underflow + } + else { + *e = de - Bias - (P-1) + 1 + k; +#ifdef Pack_32 + *bits = 32*i - hi0bits(x[i-1]); +#else + *bits = (i+2)*16 - hi0bits(x[i]); +#endif + } +#endif + return b; + } +#undef d0 +#undef d1 + +#undef CONST |