Source code
Revision control
Copy as Markdown
Other Tools
/* -*- 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
#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