Quelle  dtoa.c   Sprache: C

/* -*- 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