/* -*-C-*-

$Id: bignum.c,v 9.49 2000/12/05 21:23:43 cph Exp $

Copyright (c) 1989-2000 Massachusetts Institute of Technology

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/* Implementation of Bignums (unlimited precision integers) */

#ifdef MIT_SCHEME
#include "scheme.h"
#undef CHAR_BIT /* redefined in "limits.h" */
#else
#include "bignum.h"
#endif

#include "bignmint.h"
#include "limits.h"
#include <math.h>

#ifndef MIT_SCHEME

static bignum_type
DEFUN (bignum_malloc, (length), bignum_length_type length)
{
  extern char * malloc ();
  char * result = (malloc ((length + 1) * (sizeof (bignum_digit_type))));
  BIGNUM_ASSERT (result != ((char *) 0));
  return ((bignum_type) result);
}

static bignum_type
DEFUN (bignum_realloc, (bignum, length),
       bignum_type bignum AND bignum_length_type length)
{
  extern char * realloc ();
  char * result =
    (realloc (((char *) bignum),
	      ((length + 1) * (sizeof (bignum_digit_type)))));
  BIGNUM_ASSERT (result != ((char *) 0));
  return ((bignum_type) result);
}

#endif /* not MIT_SCHEME */

/* Forward references */
static int EXFUN (bignum_equal_p_unsigned,
		  (bignum_type, bignum_type));
static enum bignum_comparison EXFUN (bignum_compare_unsigned,
				     (bignum_type, bignum_type));
static bignum_type EXFUN (bignum_add_unsigned,
			  (bignum_type, bignum_type, int));
static bignum_type EXFUN (bignum_subtract_unsigned,
			  (bignum_type, bignum_type));
static bignum_type EXFUN (bignum_multiply_unsigned,
			  (bignum_type, bignum_type, int));
static bignum_type EXFUN (bignum_multiply_unsigned_small_factor,
			  (bignum_type, bignum_digit_type, int));
static void EXFUN (bignum_destructive_scale_up,
		   (bignum_type, bignum_digit_type));
static void EXFUN (bignum_destructive_add,
		   (bignum_type, bignum_digit_type));
static void EXFUN (bignum_divide_unsigned_large_denominator,
		   (bignum_type, bignum_type, bignum_type *, bignum_type *,
		    int, int));
static void EXFUN (bignum_destructive_normalization,
		   (bignum_type, bignum_type, int));
static void EXFUN (bignum_destructive_unnormalization,
		   (bignum_type, int));
static void EXFUN (bignum_divide_unsigned_normalized,
		   (bignum_type, bignum_type, bignum_type));
static bignum_digit_type EXFUN (bignum_divide_subtract,
				(bignum_digit_type *, bignum_digit_type *,
				 bignum_digit_type, bignum_digit_type *));
static void EXFUN (bignum_divide_unsigned_medium_denominator,
		   (bignum_type, bignum_digit_type, bignum_type *,
		    bignum_type *, int, int));
static bignum_digit_type EXFUN (bignum_digit_divide,
				(bignum_digit_type, bignum_digit_type,
				 bignum_digit_type, bignum_digit_type *));
static bignum_digit_type EXFUN (bignum_digit_divide_subtract,
				(bignum_digit_type, bignum_digit_type,
				 bignum_digit_type, bignum_digit_type *));
static void EXFUN (bignum_divide_unsigned_small_denominator,
		   (bignum_type, bignum_digit_type, bignum_type *,
		    bignum_type *, int, int));
static bignum_digit_type EXFUN (bignum_destructive_scale_down,
				(bignum_type, bignum_digit_type));
static bignum_type EXFUN (bignum_remainder_unsigned_small_denominator,
			  (bignum_type, bignum_digit_type, int));
static bignum_type EXFUN (bignum_digit_to_bignum,
			  (bignum_digit_type, int));
static bignum_type EXFUN (bignum_allocate,
			  (bignum_length_type, int));
static bignum_type EXFUN (bignum_allocate_zeroed,
			  (bignum_length_type, int));
static bignum_type EXFUN (bignum_shorten_length,
			  (bignum_type, bignum_length_type));
static bignum_type EXFUN (bignum_trim,
			  (bignum_type));
static bignum_type EXFUN (bignum_copy,
			  (bignum_type));
static bignum_type EXFUN (bignum_new_sign,
			  (bignum_type, int));
static bignum_type EXFUN (bignum_maybe_new_sign,
			  (bignum_type, int));
static void EXFUN (bignum_destructive_copy,
		   (bignum_type, bignum_type));

#define ULONG_LENGTH_IN_BITS(digit, len)		\
do {							\
  fast unsigned long w = digit;				\
  len = 0;						\
  while (w > 0xff) { len += 8; w >>= 8; }		\
  while (w > 0)    { len += 1; w >>= 1; }		\
} while (0)

/* Exports */

bignum_type
DEFUN_VOID (bignum_make_zero)
{
  fast bignum_type result = (BIGNUM_ALLOCATE (0));
  BIGNUM_SET_HEADER (result, 0, 0);
  return (result);
}

bignum_type
DEFUN (bignum_make_one, (negative_p), int negative_p)
{
  fast bignum_type result = (BIGNUM_ALLOCATE (1));
  BIGNUM_SET_HEADER (result, 1, negative_p);
  (BIGNUM_REF (result, 0)) = 1;
  return (result);
}

int
DEFUN (bignum_equal_p, (x, y),
       fast bignum_type x AND fast bignum_type y)
{
  return
    ((BIGNUM_ZERO_P (x))
     ? (BIGNUM_ZERO_P (y))
     : ((! (BIGNUM_ZERO_P (y)))
	&& ((BIGNUM_NEGATIVE_P (x))
	    ? (BIGNUM_NEGATIVE_P (y))
	    : (! (BIGNUM_NEGATIVE_P (y))))
	&& (bignum_equal_p_unsigned (x, y))));
}

enum bignum_comparison
DEFUN (bignum_test, (bignum), fast bignum_type bignum)
{
  return
    ((BIGNUM_ZERO_P (bignum))
     ? bignum_comparison_equal
     : (BIGNUM_NEGATIVE_P (bignum))
     ? bignum_comparison_less
     : bignum_comparison_greater);
}

enum bignum_comparison
DEFUN (bignum_compare, (x, y),
       fast bignum_type x AND fast bignum_type y)
{
  return
    ((BIGNUM_ZERO_P (x))
     ? ((BIGNUM_ZERO_P (y))
	? bignum_comparison_equal
	: (BIGNUM_NEGATIVE_P (y))
	? bignum_comparison_greater
	: bignum_comparison_less)
     : (BIGNUM_ZERO_P (y))
     ? ((BIGNUM_NEGATIVE_P (x))
	? bignum_comparison_less
	: bignum_comparison_greater)
     : (BIGNUM_NEGATIVE_P (x))
     ? ((BIGNUM_NEGATIVE_P (y))
	? (bignum_compare_unsigned (y, x))
	: (bignum_comparison_less))
     : ((BIGNUM_NEGATIVE_P (y))
	? (bignum_comparison_greater)
	: (bignum_compare_unsigned (x, y))));
}

bignum_type
DEFUN (bignum_add, (x, y),
       fast bignum_type x AND fast bignum_type y)
{
  return
    ((BIGNUM_ZERO_P (x))
     ? (BIGNUM_MAYBE_COPY (y))
     : (BIGNUM_ZERO_P (y))
     ? (BIGNUM_MAYBE_COPY (x))
     : ((BIGNUM_NEGATIVE_P (x))
	? ((BIGNUM_NEGATIVE_P (y))
	   ? (bignum_add_unsigned (x, y, 1))
	   : (bignum_subtract_unsigned (y, x)))
	: ((BIGNUM_NEGATIVE_P (y))
	   ? (bignum_subtract_unsigned (x, y))
	   : (bignum_add_unsigned (x, y, 0)))));
}

bignum_type
DEFUN (bignum_subtract, (x, y),
       fast bignum_type x AND fast bignum_type y)
{
  return
    ((BIGNUM_ZERO_P (x))
     ? ((BIGNUM_ZERO_P (y))
	? (BIGNUM_MAYBE_COPY (y))
	: (bignum_new_sign (y, (! (BIGNUM_NEGATIVE_P (y))))))
     : ((BIGNUM_ZERO_P (y))
	? (BIGNUM_MAYBE_COPY (x))
	: ((BIGNUM_NEGATIVE_P (x))
	   ? ((BIGNUM_NEGATIVE_P (y))
	      ? (bignum_subtract_unsigned (y, x))
	      : (bignum_add_unsigned (x, y, 1)))
	   : ((BIGNUM_NEGATIVE_P (y))
	      ? (bignum_add_unsigned (x, y, 0))
	      : (bignum_subtract_unsigned (x, y))))));
}

bignum_type
DEFUN (bignum_negate, (x), fast bignum_type x)
{
  return
    ((BIGNUM_ZERO_P (x))
     ? (BIGNUM_MAYBE_COPY (x))
     : (bignum_new_sign (x, (! (BIGNUM_NEGATIVE_P (x))))));
}

bignum_type
DEFUN (bignum_multiply, (x, y),
       fast bignum_type x AND fast bignum_type y)
{
  fast bignum_length_type x_length = (BIGNUM_LENGTH (x));
  fast bignum_length_type y_length = (BIGNUM_LENGTH (y));
  fast int negative_p =
    ((BIGNUM_NEGATIVE_P (x))
     ? (! (BIGNUM_NEGATIVE_P (y)))
     : (BIGNUM_NEGATIVE_P (y)));
  if (BIGNUM_ZERO_P (x))
    return (BIGNUM_MAYBE_COPY (x));
  if (BIGNUM_ZERO_P (y))
    return (BIGNUM_MAYBE_COPY (y));
  if (x_length == 1)
    {
      bignum_digit_type digit = (BIGNUM_REF (x, 0));
      if (digit == 1)
	return (bignum_maybe_new_sign (y, negative_p));
      if (digit < BIGNUM_RADIX_ROOT)
	return (bignum_multiply_unsigned_small_factor (y, digit, negative_p));
    }
  if (y_length == 1)
    {
      bignum_digit_type digit = (BIGNUM_REF (y, 0));
      if (digit == 1)
	return (bignum_maybe_new_sign (x, negative_p));
      if (digit < BIGNUM_RADIX_ROOT)
	return (bignum_multiply_unsigned_small_factor (x, digit, negative_p));
    }
  return (bignum_multiply_unsigned (x, y, negative_p));
}

int
DEFUN (bignum_divide, (numerator, denominator, quotient, remainder),
       bignum_type numerator AND bignum_type denominator
       AND bignum_type * quotient AND bignum_type * remainder)
{
  if (BIGNUM_ZERO_P (denominator))
    return (1);
  if (BIGNUM_ZERO_P (numerator))
    {
      (*quotient) = (BIGNUM_MAYBE_COPY (numerator));
      (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
    }
  else
    {
      int r_negative_p = (BIGNUM_NEGATIVE_P (numerator));
      int q_negative_p =
	((BIGNUM_NEGATIVE_P (denominator)) ? (! r_negative_p) : r_negative_p);
      switch (bignum_compare_unsigned (numerator, denominator))
	{
	case bignum_comparison_equal:
	  {
	    (*quotient) = (BIGNUM_ONE (q_negative_p));
	    (*remainder) = (BIGNUM_ZERO ());
	    break;
	  }
	case bignum_comparison_less:
	  {
	    (*quotient) = (BIGNUM_ZERO ());
	    (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
	    break;
	  }
	case bignum_comparison_greater:
	  {
	    if ((BIGNUM_LENGTH (denominator)) == 1)
	      {
		bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
		if (digit == 1)
		  {
		    (*quotient) =
		      (bignum_maybe_new_sign (numerator, q_negative_p));
		    (*remainder) = (BIGNUM_ZERO ());
		    break;
		  }
		else if (digit < BIGNUM_RADIX_ROOT)
		  {
		    bignum_divide_unsigned_small_denominator
		      (numerator, digit,
		       quotient, remainder,
		       q_negative_p, r_negative_p);
		    break;
		  }
		else
		  {
		    bignum_divide_unsigned_medium_denominator
		      (numerator, digit,
		       quotient, remainder,
		       q_negative_p, r_negative_p);
		    break;
		  }
	      }
	    bignum_divide_unsigned_large_denominator
	      (numerator, denominator,
	       quotient, remainder,
	       q_negative_p, r_negative_p);
	    break;
	  }
	}
    }
  return (0);
}

bignum_type
DEFUN (bignum_quotient, (numerator, denominator),
       bignum_type numerator AND bignum_type denominator)
{
  if (BIGNUM_ZERO_P (denominator))
    return (BIGNUM_OUT_OF_BAND);
  if (BIGNUM_ZERO_P (numerator))
    return (BIGNUM_MAYBE_COPY (numerator));
  {
    int q_negative_p =
      ((BIGNUM_NEGATIVE_P (denominator))
       ? (! (BIGNUM_NEGATIVE_P (numerator)))
       : (BIGNUM_NEGATIVE_P (numerator)));
    switch (bignum_compare_unsigned (numerator, denominator))
      {
      case bignum_comparison_equal:
	return (BIGNUM_ONE (q_negative_p));
      case bignum_comparison_less:
	return (BIGNUM_ZERO ());
      case bignum_comparison_greater:
	{
	  bignum_type quotient;
	  if ((BIGNUM_LENGTH (denominator)) == 1)
	    {
	      bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
	      if (digit == 1)
		return (bignum_maybe_new_sign (numerator, q_negative_p));
	      if (digit < BIGNUM_RADIX_ROOT)
		bignum_divide_unsigned_small_denominator
		  (numerator, digit,
		   (&quotient), ((bignum_type *) 0),
		   q_negative_p, 0);
	      else
		bignum_divide_unsigned_medium_denominator
		  (numerator, digit,
		   (&quotient), ((bignum_type *) 0),
		   q_negative_p, 0);
	    }
	  else
	    bignum_divide_unsigned_large_denominator
	      (numerator, denominator,
	       (&quotient), ((bignum_type *) 0),
	       q_negative_p, 0);
	  return (quotient);
	}
      default:
	/*NOTREACHED*/
	return (0);
      }
  }
}

bignum_type
DEFUN (bignum_remainder, (numerator, denominator),
       bignum_type numerator AND bignum_type denominator)
{
  if (BIGNUM_ZERO_P (denominator))
    return (BIGNUM_OUT_OF_BAND);
  if (BIGNUM_ZERO_P (numerator))
    return (BIGNUM_MAYBE_COPY (numerator));
  switch (bignum_compare_unsigned (numerator, denominator))
    {
    case bignum_comparison_equal:
      return (BIGNUM_ZERO ());
    case bignum_comparison_less:
      return (BIGNUM_MAYBE_COPY (numerator));
    case bignum_comparison_greater:
      {
	bignum_type remainder;
	if ((BIGNUM_LENGTH (denominator)) == 1)
	  {
	    bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
	    if ((digit & (digit-1)) == 0) /* i.e. digit = 2^k, including 1 */
	      {
		bignum_digit_type unsigned_remainder;
		if (BIGNUM_LENGTH (numerator) == 0)
		  return (BIGNUM_ZERO ());
		unsigned_remainder = (digit-1) & (BIGNUM_REF (numerator, 0));
		if (unsigned_remainder == 0)
		  return (BIGNUM_ZERO ());
		return (bignum_digit_to_bignum
			(unsigned_remainder, BIGNUM_NEGATIVE_P (numerator)));
	      }
	    if (digit < BIGNUM_RADIX_ROOT)
	      return
		(bignum_remainder_unsigned_small_denominator
		 (numerator, digit, (BIGNUM_NEGATIVE_P (numerator))));
	    bignum_divide_unsigned_medium_denominator
	      (numerator, digit,
	       ((bignum_type *) 0), (&remainder),
	       0, (BIGNUM_NEGATIVE_P (numerator)));
	  }
	else
	  bignum_divide_unsigned_large_denominator
	    (numerator, denominator,
	     ((bignum_type *) 0), (&remainder),
	     0, (BIGNUM_NEGATIVE_P (numerator)));
	return (remainder);
      }
    default:
      /*NOTREACHED*/
      return (0);
    }
}

/* These procedures depend on the non-portable type `unsigned long'.
   If your compiler doesn't support this type, either define the
   switch `BIGNUM_NO_ULONG' to disable them (in "bignum.h"), or write
   new versions that don't use this type. */

#ifndef BIGNUM_NO_ULONG

bignum_type
DEFUN (long_to_bignum, (n), long n)
{
  int negative_p;
  bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG];
  fast bignum_digit_type * end_digits = result_digits;
  /* Special cases win when these small constants are cached. */
  if (n == 0) return (BIGNUM_ZERO ());
  if (n == 1) return (BIGNUM_ONE (0));
  if (n == -1) return (BIGNUM_ONE (1));
  {
    fast unsigned long accumulator = ((negative_p = (n < 0)) ? (-n) : n);
    do
      {
	(*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);
	accumulator >>= BIGNUM_DIGIT_LENGTH;
      }
    while (accumulator != 0);
  }
  {
    bignum_type result =
      (bignum_allocate ((end_digits - result_digits), negative_p));
    fast bignum_digit_type * scan_digits = result_digits;
    fast bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
    while (scan_digits < end_digits)
      (*scan_result++) = (*scan_digits++);
    return (result);
  }
}

long
DEFUN (bignum_to_long, (bignum), bignum_type bignum)
{
  if (BIGNUM_ZERO_P (bignum))
    return (0);
  {
    fast unsigned long accumulator = 0;
    fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
    fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
    while (start < scan)
      accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));
    return
      ((BIGNUM_NEGATIVE_P (bignum))
       ? (- ((long) accumulator))
       : ((long) accumulator));
  }
}

bignum_type
DEFUN (ulong_to_bignum, (n), unsigned long n)
{
  bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG];
  fast bignum_digit_type * end_digits = result_digits;
  /* Special cases win when these small constants are cached. */
  if (n == 0) return (BIGNUM_ZERO ());
  if (n == 1) return (BIGNUM_ONE (0));
  {
    fast unsigned long accumulator = n;
    do
      {
	(*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);
	accumulator >>= BIGNUM_DIGIT_LENGTH;
      }
    while (accumulator != 0);
  }
  {
    bignum_type result =
      (bignum_allocate ((end_digits - result_digits), 0));
    fast bignum_digit_type * scan_digits = result_digits;
    fast bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
    while (scan_digits < end_digits)
      (*scan_result++) = (*scan_digits++);
    return (result);
  }
}

unsigned long
DEFUN (bignum_to_ulong, (bignum), bignum_type bignum)
{
  if (BIGNUM_ZERO_P (bignum))
    return (0);
  {
    fast unsigned long accumulator = 0;
    fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
    fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
    while (start < scan)
      accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));
    return (accumulator);
  }
}

#endif /* not BIGNUM_NO_ULONG */

#define DTB_WRITE_DIGIT(factor)						\
{									\
  significand *= (factor);						\
  digit = ((bignum_digit_type) significand);				\
  (*--scan) = digit;							\
  significand -= ((double) digit);					\
}

bignum_type
DEFUN (double_to_bignum, (x), double x)
{
  extern double frexp ();
  int exponent;
  fast double significand = (frexp (x, (&exponent)));
  if (exponent <= 0) return (BIGNUM_ZERO ());
  if (exponent == 1) return (BIGNUM_ONE (x < 0));
  if (significand < 0) significand = (-significand);
  {
    bignum_length_type length = (BIGNUM_BITS_TO_DIGITS (exponent));
    bignum_type result = (bignum_allocate (length, (x < 0)));
    bignum_digit_type * start = (BIGNUM_START_PTR (result));
    fast bignum_digit_type * scan = (start + length);
    fast bignum_digit_type digit;
    int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
    if (odd_bits > 0)
      DTB_WRITE_DIGIT (1L << odd_bits);
    while (start < scan)
      {
	if (significand == 0)
	  {
	    while (start < scan)
	      (*--scan) = 0;
	    break;
	  }
	DTB_WRITE_DIGIT (BIGNUM_RADIX);
      }
    return (result);
  }
}

#undef DTB_WRITE_DIGIT

/*  This version sometimes gets the answer wrong due to rounding errors.
    It would be slightly better if the digits were accumulated lsb to msb.
 */
/*
double
DEFUN (bignum_to_double, (bignum), bignum_type bignum)
{
  if (BIGNUM_ZERO_P (bignum))
    return (0);
  {
    fast double accumulator = 0;
    fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
    fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
    while (start < scan)
      accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
    return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator);
  }
}
*/

double
DEFUN (bignum_to_double, (bignum), bignum_type bignum)
{
  if (BIGNUM_ZERO_P (bignum))
    return (0.0);

  {
    bignum_length_type length = (BIGNUM_LENGTH (bignum));
    bignum_length_type index = length - 1;
    bignum_length_type scale_words = length - 1;
    bignum_digit_type msd = (BIGNUM_REF (bignum, (index)));
#if (FLT_RADIX == 2)
    int bits_to_get = DBL_MANT_DIG; /* includes implicit 1 */
#else
#include "error: must have FLT_RADIX==2"
#endif
    double value = 0;
    bignum_digit_type mask = 0;
    bignum_digit_type guard_bit_mask = BIGNUM_RADIX>>1;
    bignum_digit_type rounding_correction = 0;
    int current_digit_bit_count = 0;

    ULONG_LENGTH_IN_BITS (msd, current_digit_bit_count);
    mask = (1 << (current_digit_bit_count)) - 1;

    while (1) {
      if (current_digit_bit_count > bits_to_get) {
	guard_bit_mask = (1 << (current_digit_bit_count - bits_to_get - 1));
	mask &= ~((guard_bit_mask << 1) - 1);
	current_digit_bit_count = bits_to_get;
	bits_to_get = 0;
      } else {
	bits_to_get -= current_digit_bit_count;
      }

      value = (value * BIGNUM_RADIX) + ((BIGNUM_REF (bignum, index)) & mask);

      if (bits_to_get == 0) {
	scale_words = index;
	if (current_digit_bit_count == BIGNUM_DIGIT_LENGTH) {
	  if (index == 0) /* there is no guard bit */
	    goto finished;
	  guard_bit_mask = (1 << (BIGNUM_DIGIT_LENGTH - 1));
	  rounding_correction = 1;
	  index -= 1;
	} else {
	  rounding_correction = (guard_bit_mask << 1);
	}
	break;
      }
      if (index == 0)  /* fewer than DBL_MANT_DIG bits */
	goto finished;

      index -= 1;
      current_digit_bit_count = BIGNUM_DIGIT_LENGTH;
      mask = BIGNUM_DIGIT_MASK;
    }

    /* round-to-even depending on lsb, guard and following bits: lgfffff */

    if ((BIGNUM_REF(bignum,index) & guard_bit_mask) == 0) /* case x0xxxx */
      goto round_down;

    if ((BIGNUM_REF(bignum,index) & (guard_bit_mask-1)) != 0) /* case x1xx1x */
      goto round_up;

    /* cases 110000 and 1101xx: test "odd?", i.e. round-to-even rounds up */
    if ((guard_bit_mask << 1) == BIGNUM_RADIX) {
      if (((BIGNUM_REF (bignum, index+1)) & 1) != 0)  /* "odd?" */
	goto round_up;
    } else {
      if (((BIGNUM_REF (bignum, index)) & (guard_bit_mask << 1)) != 0)
	goto round_up;
    }

    if (index==0)   /* case 010000, no more words of following bits */
      goto finished;

    { /* distinguish between cases 0100...00 and 0100..1xx, multiple words */
      bignum_length_type index2 = index - 1;
      while (index2 >= 0) {
	if ((BIGNUM_REF (bignum, index2)) != 0)
	  goto round_up;
	index2--;
      }
      goto round_down;
    }

  round_up:
    value += rounding_correction;
  round_down:
    /* note, ldexp `sticks' at the maximal non-infinity value, which
       is a reasonable compromise for numbers with DBL_MAX_EXP bits
       that round up */
    if (scale_words > 0)
      value = ldexp (value, scale_words * BIGNUM_DIGIT_LENGTH);

  finished:
    return ((BIGNUM_NEGATIVE_P (bignum)) ? (-value) : value);
  }
}

/*
;; Code to test bignum_to_double

(declare (usual-integrations))

(define integer->flonum (make-primitive-procedure 'integer->flonum 2))

(define (check n)
  (let ((f1 (integer->flonum n #b10))
	(f2 (exact->inexact n)))
    (if (and f1 (= f1 f2))
	(number->string f1 2)
	(begin
	  (pp n)
	  (pp (number->string f1 2))
	  (pp (number->string f2 2))))))

(define (test)
  (define n 0)
  (do ((i 0 (+ i 1)))			; guard bit zone patterns
      ((= i 256))
    (do ((j 0 (+ j 1)))			; general word alignment
	((= j 30))
      (do ((e 0 (+ e 1)))		; random `insignificant' bit
	  ((= e 2))
	(do ((up 0 (+ up 1)))		; test potential for carry
	    ((= up 2))
	  (do ((end-pad 0 (+ end-pad 100)))
	      ((> end-pad 100))
	    (set! n (+ n 1)) (if (= (remainder n 1000) 0) (pp `(test ,n)))
	    (let ((s (string-append "1"
				   (make-string 48 (vector-ref '#(#\0 #\1) up))
				   (string-pad-left (number->string i 2) 8 #\0)
				    (make-string (* j 23) #\0) ; gcd(23,30)=1
				    (number->string e)
				    (make-string end-pad #\0))))
	      (check (string->number s 2)))))))))
*/

int
DEFUN (bignum_fits_in_word_p, (bignum, word_length, twos_complement_p),
       bignum_type bignum AND long word_length AND int twos_complement_p)
{
  unsigned int n_bits = (twos_complement_p ? (word_length - 1) : word_length);
  BIGNUM_ASSERT (n_bits > 0);
  {
    fast bignum_length_type length = (BIGNUM_LENGTH (bignum));
    fast unsigned int max_digits = (BIGNUM_BITS_TO_DIGITS (n_bits));
    if (((unsigned int) length) < max_digits)
      return (1);
    if (((unsigned int) length) > max_digits)
      return (0);
    {
      bignum_digit_type msd = (BIGNUM_REF (bignum, (length - 1)));
      bignum_digit_type max
	= (1L << (n_bits - ((length - 1) * BIGNUM_DIGIT_LENGTH)));
      return
	(((msd < max)
	  || (twos_complement_p
	      && (msd == max)
	      && (BIGNUM_NEGATIVE_P (bignum)))));
    }
  }
}

bignum_type
DEFUN (bignum_length_in_bits, (bignum), bignum_type bignum)
{
  if (BIGNUM_ZERO_P (bignum))
    return (BIGNUM_ZERO ());
  {
    bignum_length_type index = ((BIGNUM_LENGTH (bignum)) - 1);
    fast bignum_digit_type digit = (BIGNUM_REF (bignum, index));
    fast bignum_digit_type delta = 0;
    fast bignum_type result = (bignum_allocate (2, 0));
    (BIGNUM_REF (result, 0)) = index;
    (BIGNUM_REF (result, 1)) = 0;
    bignum_destructive_scale_up (result, BIGNUM_DIGIT_LENGTH);
    ULONG_LENGTH_IN_BITS (digit, delta);
    bignum_destructive_add (result, ((bignum_digit_type) delta));
    return (bignum_trim (result));
  }
}

bignum_type
DEFUN_VOID (bignum_length_upper_limit)
{
  fast bignum_type result = (bignum_allocate (2, 0));
  (BIGNUM_REF (result, 0)) = 0;
  (BIGNUM_REF (result, 1)) = BIGNUM_DIGIT_LENGTH;
  return (result);
}

bignum_type
DEFUN (bignum_shift_left, (n, m), bignum_type n AND unsigned long m)
{
  unsigned long ln = (BIGNUM_LENGTH (n));
  unsigned long delta = 0;
  if (m == 0)
    return (n);

  ULONG_LENGTH_IN_BITS ((BIGNUM_REF (n, (ln - 1))), delta);

  {
    unsigned long zeroes = (m / BIGNUM_DIGIT_LENGTH);
    unsigned long shift = (m % BIGNUM_DIGIT_LENGTH);
    unsigned long ln2
      = (((ln - 1) + ((delta + m) / BIGNUM_DIGIT_LENGTH))
	 + (((delta + m) % BIGNUM_DIGIT_LENGTH) != 0));
    bignum_type result = (bignum_allocate (ln2, (BIGNUM_NEGATIVE_P (n))));
    bignum_digit_type * scan_n = (BIGNUM_START_PTR (n));
    bignum_digit_type * end_n = (scan_n + ln);
    bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
    while ((zeroes--) > 0)
      (*scan_result++) = 0;
    if (shift == 0)
      while (scan_n < end_n)
	(*scan_result++) = (*scan_n++);
    else
      {
	unsigned long temp = 0;
	while (scan_n < end_n)
	  {
	    bignum_digit_type digit = (*scan_n++);
	    (*scan_result++) = (((digit << shift) & BIGNUM_DIGIT_MASK) | temp);
	    temp = (digit >> (BIGNUM_DIGIT_LENGTH - shift));
	  }
	if (temp != 0)
	  (*scan_result) = temp;
      }
    return (result);
  }
}

bignum_type
DEFUN (unsigned_long_to_shifted_bignum, (n, m, sign),
       unsigned long n AND
       unsigned long m AND
       int sign)
{
  unsigned long delta = 0;
  if (n == 0)
    return (BIGNUM_ZERO ());

  ULONG_LENGTH_IN_BITS (n, delta);

  {
    unsigned long zeroes = (m / BIGNUM_DIGIT_LENGTH);
    unsigned long shift = (m % BIGNUM_DIGIT_LENGTH);
    unsigned long ln
      = (((delta + m) / BIGNUM_DIGIT_LENGTH)
	 + (((delta + m) % BIGNUM_DIGIT_LENGTH) != 0));
    bignum_type result = (bignum_allocate (ln, sign));
    bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
    while ((zeroes--) > 0)
      (*scan_result++) = 0;
    (*scan_result++) = ((n << shift) & BIGNUM_DIGIT_MASK);
    n >>= (BIGNUM_DIGIT_LENGTH - shift);
    while (n > 0)
      {
	(*scan_result++) = (n & BIGNUM_DIGIT_MASK);
	n >>= BIGNUM_DIGIT_LENGTH;
      }
    return (result);
  }
}

bignum_type
DEFUN (digit_stream_to_bignum,
       (n_digits, producer, context, radix, negative_p),
       fast unsigned int n_digits
       AND unsigned int EXFUN ((*producer), (bignum_procedure_context))
       AND bignum_procedure_context context
       AND fast unsigned int radix
       AND int negative_p)
{
  BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  if (n_digits == 0)
    return (BIGNUM_ZERO ());
  if (n_digits == 1)
    {
      long digit = ((long) ((*producer) (context)));
      return (long_to_bignum (negative_p ? (- digit) : digit));
    }
  {
    bignum_length_type length;
    {
      fast unsigned int log_radix = 0;
      ULONG_LENGTH_IN_BITS (radix, log_radix);
      /* This length will be at least as large as needed. */
      length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix));
    }
    {
      fast bignum_type result = (bignum_allocate_zeroed (length, negative_p));
      while ((n_digits--) > 0)
	{
	  bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
	  bignum_destructive_add
	    (result, ((bignum_digit_type) ((*producer) (context))));
	}
      return (bignum_trim (result));
    }
  }
}

void
DEFUN (bignum_to_digit_stream, (bignum, radix, consumer, context),
       bignum_type bignum
       AND unsigned int radix
       AND void EXFUN ((*consumer), (bignum_procedure_context, long))
       AND bignum_procedure_context context)
{
  BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  if (! (BIGNUM_ZERO_P (bignum)))
    {
      fast bignum_type working_copy = (bignum_copy (bignum));
      fast bignum_digit_type * start = (BIGNUM_START_PTR (working_copy));
      fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (working_copy)));
      while (start < scan)
	{
	  if ((scan[-1]) == 0)
	    scan -= 1;
	  else
	    (*consumer)
	      (context, (bignum_destructive_scale_down (working_copy, radix)));
	}
      BIGNUM_DEALLOCATE (working_copy);
    }
  return;
}

long
DEFUN_VOID (bignum_max_digit_stream_radix)
{
  return (BIGNUM_RADIX_ROOT);
}

/* Comparisons */

static int
DEFUN (bignum_equal_p_unsigned, (x, y),
       bignum_type x AND bignum_type y)
{
  bignum_length_type length = (BIGNUM_LENGTH (x));
  if (length != ((bignum_length_type) (BIGNUM_LENGTH (y))))
    return (0);
  else
    {
      fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
      fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
      fast bignum_digit_type * end_x = (scan_x + length);
      while (scan_x < end_x)
	if ((*scan_x++) != (*scan_y++))
	  return (0);
      return (1);
    }
}

static enum bignum_comparison
DEFUN (bignum_compare_unsigned, (x, y),
       bignum_type x AND bignum_type y)
{
  bignum_length_type x_length = (BIGNUM_LENGTH (x));
  bignum_length_type y_length = (BIGNUM_LENGTH (y));
  if (x_length < y_length)
    return (bignum_comparison_less);
  if (x_length > y_length)
    return (bignum_comparison_greater);
  {
    fast bignum_digit_type * start_x = (BIGNUM_START_PTR (x));
    fast bignum_digit_type * scan_x = (start_x + x_length);
    fast bignum_digit_type * scan_y = ((BIGNUM_START_PTR (y)) + y_length);
    while (start_x < scan_x)
      {
	fast bignum_digit_type digit_x = (*--scan_x);
	fast bignum_digit_type digit_y = (*--scan_y);
	if (digit_x < digit_y)
	  return (bignum_comparison_less);
	if (digit_x > digit_y)
	  return (bignum_comparison_greater);
      }
  }
  return (bignum_comparison_equal);
}

/* Addition */

static bignum_type
DEFUN (bignum_add_unsigned, (x, y, negative_p),
       bignum_type x AND bignum_type y AND int negative_p)
{
  if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
    {
      bignum_type z = x;
      x = y;
      y = z;
    }
  {
    bignum_length_type x_length = (BIGNUM_LENGTH (x));
    bignum_type r = (bignum_allocate ((x_length + 1), negative_p));
    fast bignum_digit_type sum;
    fast bignum_digit_type carry = 0;
    fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
    fast bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
    {
      fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
      fast bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
      while (scan_y < end_y)
	{
	  sum = ((*scan_x++) + (*scan_y++) + carry);
	  if (sum < BIGNUM_RADIX)
	    {
	      (*scan_r++) = sum;
	      carry = 0;
	    }
	  else
	    {
	      (*scan_r++) = (sum - BIGNUM_RADIX);
	      carry = 1;
	    }
	}
    }
    {
      fast bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
      if (carry != 0)
	while (scan_x < end_x)
	  {
	    sum = ((*scan_x++) + 1);
	    if (sum < BIGNUM_RADIX)
	      {
		(*scan_r++) = sum;
		carry = 0;
		break;
	      }
	    else
	      (*scan_r++) = (sum - BIGNUM_RADIX);
	  }
      while (scan_x < end_x)
	(*scan_r++) = (*scan_x++);
    }
    if (carry != 0)
      {
	(*scan_r) = 1;
	return (r);
      }
    return (bignum_shorten_length (r, x_length));
  }
}

/* Subtraction */

static bignum_type
DEFUN (bignum_subtract_unsigned, (x, y),
       bignum_type x AND bignum_type y)
{
  int negative_p = 0;
  switch (bignum_compare_unsigned (x, y))
    {
    case bignum_comparison_equal:
      return (BIGNUM_ZERO ());
    case bignum_comparison_less:
      {
	bignum_type z = x;
	x = y;
	y = z;
      }
      negative_p = 1;
      break;
    case bignum_comparison_greater:
      negative_p = 0;
      break;
    }
  {
    bignum_length_type x_length = (BIGNUM_LENGTH (x));
    bignum_type r = (bignum_allocate (x_length, negative_p));
    fast bignum_digit_type difference;
    fast bignum_digit_type borrow = 0;
    fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
    fast bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
    {
      fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
      fast bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
      while (scan_y < end_y)
	{
	  difference = (((*scan_x++) - (*scan_y++)) - borrow);
	  if (difference < 0)
	    {
	      (*scan_r++) = (difference + BIGNUM_RADIX);
	      borrow = 1;
	    }
	  else
	    {
	      (*scan_r++) = difference;
	      borrow = 0;
	    }
	}
    }
    {
      fast bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
      if (borrow != 0)
	while (scan_x < end_x)
	  {
	    difference = ((*scan_x++) - borrow);
	    if (difference < 0)
	      (*scan_r++) = (difference + BIGNUM_RADIX);
	    else
	      {
		(*scan_r++) = difference;
		borrow = 0;
		break;
	      }
	  }
      BIGNUM_ASSERT (borrow == 0);
      while (scan_x < end_x)
	(*scan_r++) = (*scan_x++);
    }
    return (bignum_trim (r));
  }
}

/* Multiplication
   Maximum value for product_low or product_high:
	((R * R) + (R * (R - 2)) + (R - 1))
   Maximum value for carry: ((R * (R - 1)) + (R - 1))
	where R == BIGNUM_RADIX_ROOT */

static bignum_type
DEFUN (bignum_multiply_unsigned, (x, y, negative_p),
       bignum_type x AND bignum_type y AND int negative_p)
{
  if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
    {
      bignum_type z = x;
      x = y;
      y = z;
    }
  {
    fast bignum_digit_type carry;
    fast bignum_digit_type y_digit_low;
    fast bignum_digit_type y_digit_high;
    fast bignum_digit_type x_digit_low;
    fast bignum_digit_type x_digit_high;
    bignum_digit_type product_low;
    fast bignum_digit_type * scan_r;
    fast bignum_digit_type * scan_y;
    bignum_length_type x_length = (BIGNUM_LENGTH (x));
    bignum_length_type y_length = (BIGNUM_LENGTH (y));
    bignum_type r =
      (bignum_allocate_zeroed ((x_length + y_length), negative_p));
    bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
    bignum_digit_type * end_x = (scan_x + x_length);
    bignum_digit_type * start_y = (BIGNUM_START_PTR (y));
    bignum_digit_type * end_y = (start_y + y_length);
    bignum_digit_type * start_r = (BIGNUM_START_PTR (r));
#define x_digit x_digit_high
#define y_digit y_digit_high
#define product_high carry
    while (scan_x < end_x)
      {
	x_digit = (*scan_x++);
	x_digit_low = (HD_LOW (x_digit));
	x_digit_high = (HD_HIGH (x_digit));
	carry = 0;
	scan_y = start_y;
	scan_r = (start_r++);
	while (scan_y < end_y)
	  {
	    y_digit = (*scan_y++);
	    y_digit_low = (HD_LOW (y_digit));
	    y_digit_high = (HD_HIGH (y_digit));
	    product_low =
	      ((*scan_r) +
	       (x_digit_low * y_digit_low) +
	       (HD_LOW (carry)));
	    product_high =
	      ((x_digit_high * y_digit_low) +
	       (x_digit_low * y_digit_high) +
	       (HD_HIGH (product_low)) +
	       (HD_HIGH (carry)));
	    (*scan_r++) =
	      (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
	    carry =
	      ((x_digit_high * y_digit_high) +
	       (HD_HIGH (product_high)));
	  }
	(*scan_r) += carry;
      }
    return (bignum_trim (r));
#undef x_digit
#undef y_digit
#undef product_high
  }
}

static bignum_type
DEFUN (bignum_multiply_unsigned_small_factor, (x, y, negative_p),
       bignum_type x AND bignum_digit_type y AND int negative_p)
{
  bignum_length_type length_x = (BIGNUM_LENGTH (x));
  bignum_type p = (bignum_allocate ((length_x + 1), negative_p));
  bignum_destructive_copy (x, p);
  (BIGNUM_REF (p, length_x)) = 0;
  bignum_destructive_scale_up (p, y);
  return (bignum_trim (p));
}

static void
DEFUN (bignum_destructive_scale_up, (bignum, factor),
       bignum_type bignum AND bignum_digit_type factor)
{
  fast bignum_digit_type carry = 0;
  fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  fast bignum_digit_type two_digits;
  fast bignum_digit_type product_low;
#define product_high carry
  bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
  BIGNUM_ASSERT ((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
  while (scan < end)
    {
      two_digits = (*scan);
      product_low = ((factor * (HD_LOW (two_digits))) + (HD_LOW (carry)));
      product_high =
	((factor * (HD_HIGH (two_digits))) +
	 (HD_HIGH (product_low)) +
	 (HD_HIGH (carry)));
      (*scan++) = (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
      carry = (HD_HIGH (product_high));
    }
  /* A carry here would be an overflow, i.e. it would not fit.
     Hopefully the callers allocate enough space that this will
     never happen.
   */
  BIGNUM_ASSERT (carry == 0);
  return;
#undef product_high
}

static void
DEFUN (bignum_destructive_add, (bignum, n),
       bignum_type bignum AND bignum_digit_type n)
{
  fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  fast bignum_digit_type digit;
  digit = ((*scan) + n);
  if (digit < BIGNUM_RADIX)
    {
      (*scan) = digit;
      return;
    }
  (*scan++) = (digit - BIGNUM_RADIX);
  while (1)
    {
      digit = ((*scan) + 1);
      if (digit < BIGNUM_RADIX)
	{
	  (*scan) = digit;
	  return;
	}
      (*scan++) = (digit - BIGNUM_RADIX);
    }
}

/* Division */

/* For help understanding this algorithm, see:
   Knuth, Donald E., "The Art of Computer Programming",
   volume 2, "Seminumerical Algorithms"
   section 4.3.1, "Multiple-Precision Arithmetic". */

static void
DEFUN (bignum_divide_unsigned_large_denominator, (numerator, denominator,
						  quotient, remainder,
						  q_negative_p, r_negative_p),
       bignum_type numerator
       AND bignum_type denominator
       AND bignum_type * quotient
       AND bignum_type * remainder
       AND int q_negative_p
       AND int r_negative_p)
{
  bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1);
  bignum_length_type length_d = (BIGNUM_LENGTH (denominator));
  bignum_type q =
    ((quotient != ((bignum_type *) 0))
     ? (bignum_allocate ((length_n - length_d), q_negative_p))
     : BIGNUM_OUT_OF_BAND);
  bignum_type u = (bignum_allocate (length_n, r_negative_p));
  int shift = 0;
  BIGNUM_ASSERT (length_d > 1);
  {
    fast bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1)));
    while (v1 < (BIGNUM_RADIX / 2))
      {
	v1 <<= 1;
	shift += 1;
      }
  }
  if (shift == 0)
    {
      bignum_destructive_copy (numerator, u);
      (BIGNUM_REF (u, (length_n - 1))) = 0;
      bignum_divide_unsigned_normalized (u, denominator, q);
    }
  else
    {
      bignum_type v = (bignum_allocate (length_d, 0));
      bignum_destructive_normalization (numerator, u, shift);
      bignum_destructive_normalization (denominator, v, shift);
      bignum_divide_unsigned_normalized (u, v, q);
      BIGNUM_DEALLOCATE (v);
      if (remainder != ((bignum_type *) 0))
	bignum_destructive_unnormalization (u, shift);
    }
  if (quotient != ((bignum_type *) 0))
    (*quotient) = (bignum_trim (q));
  if (remainder != ((bignum_type *) 0))
    (*remainder) = (bignum_trim (u));
  else
    BIGNUM_DEALLOCATE (u);
  return;
}

static void
DEFUN (bignum_divide_unsigned_normalized, (u, v, q),
       bignum_type u AND bignum_type v AND bignum_type q)
{
  bignum_length_type u_length = (BIGNUM_LENGTH (u));
  bignum_length_type v_length = (BIGNUM_LENGTH (v));
  bignum_digit_type * u_start = (BIGNUM_START_PTR (u));
  bignum_digit_type * u_scan = (u_start + u_length);
  bignum_digit_type * u_scan_limit = (u_start + v_length);
  bignum_digit_type * u_scan_start = (u_scan - v_length);
  bignum_digit_type * v_start = (BIGNUM_START_PTR (v));
  bignum_digit_type * v_end = (v_start + v_length);
  bignum_digit_type * q_scan = 0;
  bignum_digit_type v1 = (v_end[-1]);
  bignum_digit_type v2 = (v_end[-2]);
  fast bignum_digit_type ph;	/* high half of double-digit product */
  fast bignum_digit_type pl;	/* low half of double-digit product */
  fast bignum_digit_type guess;
  fast bignum_digit_type gh;	/* high half-digit of guess */
  fast bignum_digit_type ch;	/* high half of double-digit comparand */
  fast bignum_digit_type v2l = (HD_LOW (v2));
  fast bignum_digit_type v2h = (HD_HIGH (v2));
  fast bignum_digit_type cl;	/* low half of double-digit comparand */
#define gl ph			/* low half-digit of guess */
#define uj pl
#define qj ph
  bignum_digit_type gm;		/* memory loc for reference parameter */
  if (q != BIGNUM_OUT_OF_BAND)
    q_scan = ((BIGNUM_START_PTR (q)) + (BIGNUM_LENGTH (q)));
  while (u_scan_limit < u_scan)
    {
      uj = (*--u_scan);
      if (uj != v1)
	{
	  /* comparand =
	     (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
	     guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
	  cl = (u_scan[-2]);
	  ch = (bignum_digit_divide (uj, (u_scan[-1]), v1, (&gm)));
	  guess = gm;
	}
      else
	{
	  cl = (u_scan[-2]);
	  ch = ((u_scan[-1]) + v1);
	  guess = (BIGNUM_RADIX - 1);
	}
      while (1)
	{
	  /* product = (guess * v2); */
	  gl = (HD_LOW (guess));
	  gh = (HD_HIGH (guess));
	  pl = (v2l * gl);
	  ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH (pl)));
	  pl = (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl))));
	  ph = ((v2h * gh) + (HD_HIGH (ph)));
	  /* if (comparand >= product) */
	  if ((ch > ph) || ((ch == ph) && (cl >= pl)))
	    break;
	  guess -= 1;
	  /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
	  ch += v1;
	  /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
	  if (ch >= BIGNUM_RADIX)
	    break;
	}
      qj = (bignum_divide_subtract (v_start, v_end, guess, (--u_scan_start)));
      if (q != BIGNUM_OUT_OF_BAND)
	(*--q_scan) = qj;
    }
  return;
#undef gl
#undef uj
#undef qj
}

static bignum_digit_type
DEFUN (bignum_divide_subtract, (v_start, v_end, guess, u_start),
       bignum_digit_type * v_start
       AND bignum_digit_type * v_end
       AND bignum_digit_type guess
       AND bignum_digit_type * u_start)
{
  bignum_digit_type * v_scan = v_start;
  bignum_digit_type * u_scan = u_start;
  fast bignum_digit_type carry = 0;
  if (guess == 0) return (0);
  {
    bignum_digit_type gl = (HD_LOW (guess));
    bignum_digit_type gh = (HD_HIGH (guess));
    fast bignum_digit_type v;
    fast bignum_digit_type pl;
    fast bignum_digit_type vl;
#define vh v
#define ph carry
#define diff pl
    while (v_scan < v_end)
      {
	v = (*v_scan++);
	vl = (HD_LOW (v));
	vh = (HD_HIGH (v));
	pl = ((vl * gl) + (HD_LOW (carry)));
	ph = ((vl * gh) + (vh * gl) + (HD_HIGH (pl)) + (HD_HIGH (carry)));
	diff = ((*u_scan) - (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl)))));
	if (diff < 0)
	  {
	    (*u_scan++) = (diff + BIGNUM_RADIX);
	    carry = ((vh * gh) + (HD_HIGH (ph)) + 1);
	  }
	else
	  {
	    (*u_scan++) = diff;
	    carry = ((vh * gh) + (HD_HIGH (ph)));
	  }
      }
    if (carry == 0)
      return (guess);
    diff = ((*u_scan) - carry);
    if (diff < 0)
      (*u_scan) = (diff + BIGNUM_RADIX);
    else
      {
	(*u_scan) = diff;
	return (guess);
      }
#undef vh
#undef ph
#undef diff
  }
  /* Subtraction generated carry, implying guess is one too large.
     Add v back in to bring it back down. */
  v_scan = v_start;
  u_scan = u_start;
  carry = 0;
  while (v_scan < v_end)
    {
      bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
      if (sum < BIGNUM_RADIX)
	{
	  (*u_scan++) = sum;
	  carry = 0;
	}
      else
	{
	  (*u_scan++) = (sum - BIGNUM_RADIX);
	  carry = 1;
	}
    }
  if (carry == 1)
    {
      bignum_digit_type sum = ((*u_scan) + carry);
      (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
    }
  return (guess - 1);
}

static void
DEFUN (bignum_divide_unsigned_medium_denominator, (numerator, denominator,
						   quotient, remainder,
						   q_negative_p, r_negative_p),
       bignum_type numerator
       AND bignum_digit_type denominator
       AND bignum_type * quotient
       AND bignum_type * remainder
       AND int q_negative_p
       AND int r_negative_p)
{
  bignum_length_type length_n = (BIGNUM_LENGTH (numerator));
  bignum_length_type length_q;
  bignum_type q;
  int shift = 0;
  /* Because `bignum_digit_divide' requires a normalized denominator. */
  while (denominator < (BIGNUM_RADIX / 2))
    {
      denominator <<= 1;
      shift += 1;
    }
  if (shift == 0)
    {
      length_q = length_n;
      q = (bignum_allocate (length_q, q_negative_p));
      bignum_destructive_copy (numerator, q);
    }
  else
    {
      length_q = (length_n + 1);
      q = (bignum_allocate (length_q, q_negative_p));
      bignum_destructive_normalization (numerator, q, shift);
    }
  {
    fast bignum_digit_type r = 0;
    fast bignum_digit_type * start = (BIGNUM_START_PTR (q));
    fast bignum_digit_type * scan = (start + length_q);
    bignum_digit_type qj;
    if (quotient != ((bignum_type *) 0))
      {
	while (start < scan)
	  {
	    r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
	    (*scan) = qj;
	  }
	(*quotient) = (bignum_trim (q));
      }
    else
      {
	while (start < scan)
	  r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
	BIGNUM_DEALLOCATE (q);
      }
    if (remainder != ((bignum_type *) 0))
      {
	if (shift != 0)
	  r >>= shift;
	(*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
      }
  }
  return;
}

static void
DEFUN (bignum_destructive_normalization, (source, target, shift_left),
       bignum_type source AND bignum_type target AND int shift_left)
{
  fast bignum_digit_type digit;
  fast bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
  fast bignum_digit_type carry = 0;
  fast bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
  bignum_digit_type * end_source = (scan_source + (BIGNUM_LENGTH (source)));
  bignum_digit_type * end_target = (scan_target + (BIGNUM_LENGTH (target)));
  int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
  bignum_digit_type mask = ((1L << shift_right) - 1);
  while (scan_source < end_source)
    {
      digit = (*scan_source++);
      (*scan_target++) = (((digit & mask) << shift_left) | carry);
      carry = (digit >> shift_right);
    }
  if (scan_target < end_target)
    (*scan_target) = carry;
  else
    BIGNUM_ASSERT (carry == 0);
  return;
}

static void
DEFUN (bignum_destructive_unnormalization, (bignum, shift_right),
       bignum_type bignum AND int shift_right)
{
  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  fast bignum_digit_type digit;
  fast bignum_digit_type carry = 0;
  int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
  bignum_digit_type mask = ((1L << shift_right) - 1);
  while (start < scan)
    {
      digit = (*--scan);
      (*scan) = ((digit >> shift_right) | carry);
      carry = ((digit & mask) << shift_left);
    }
  BIGNUM_ASSERT (carry == 0);
  return;
}

/* This is a reduced version of the division algorithm, applied to the
   case of dividing two bignum digits by one bignum digit.  It is
   assumed that the numerator and denominator are normalized. */

#define BDD_STEP(qn, j)							\
{									\
  uj = (u[j]);								\
  if (uj != v1)								\
    {									\
      uj_uj1 = (HD_CONS (uj, (u[j + 1])));				\
      guess = (uj_uj1 / v1);						\
      comparand = (HD_CONS ((uj_uj1 % v1), (u[j + 2])));		\
    }									\
  else									\
    {									\
      guess = (BIGNUM_RADIX_ROOT - 1);					\
      comparand = (HD_CONS (((u[j + 1]) + v1), (u[j + 2])));		\
    }									\
  while ((guess * v2) > comparand)					\
    {									\
      guess -= 1;							\
      comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);			\
      if (comparand >= BIGNUM_RADIX)					\
	break;								\
    }									\
  qn = (bignum_digit_divide_subtract (v1, v2, guess, (&u[j])));		\
}

static bignum_digit_type
DEFUN (bignum_digit_divide, (uh, ul, v, q),
       bignum_digit_type uh AND bignum_digit_type ul
       AND bignum_digit_type v AND bignum_digit_type * q) /* return value */
{
  fast bignum_digit_type guess;
  fast bignum_digit_type comparand;
  fast bignum_digit_type v1 = (HD_HIGH (v));
  fast bignum_digit_type v2 = (HD_LOW (v));
  fast bignum_digit_type uj;
  fast bignum_digit_type uj_uj1;
  bignum_digit_type q1;
  bignum_digit_type q2;
  bignum_digit_type u [4];
  if (uh == 0)
    {
      if (ul < v)
	{
	  (*q) = 0;
	  return (ul);
	}
      else if (ul == v)
	{
	  (*q) = 1;
	  return (0);
	}
    }
  (u[0]) = (HD_HIGH (uh));
  (u[1]) = (HD_LOW (uh));
  (u[2]) = (HD_HIGH (ul));
  (u[3]) = (HD_LOW (ul));
  v1 = (HD_HIGH (v));
  v2 = (HD_LOW (v));
  BDD_STEP (q1, 0);
  BDD_STEP (q2, 1);
  (*q) = (HD_CONS (q1, q2));
  return (HD_CONS ((u[2]), (u[3])));
}

#undef BDD_STEP

#define BDDS_MULSUB(vn, un, carry_in)					\
{									\
  product = ((vn * guess) + carry_in);					\
  diff = (un - (HD_LOW (product)));					\
  if (diff < 0)								\
    {									\
      un = (diff + BIGNUM_RADIX_ROOT);					\
      carry = ((HD_HIGH (product)) + 1);				\
    }									\
  else									\
    {									\
      un = diff;							\
      carry = (HD_HIGH (product));					\
    }									\
}

#define BDDS_ADD(vn, un, carry_in)					\
{									\
  sum = (vn + un + carry_in);						\
  if (sum < BIGNUM_RADIX_ROOT)						\
    {									\
      un = sum;								\
      carry = 0;							\
    }									\
  else									\
    {									\
      un = (sum - BIGNUM_RADIX_ROOT);					\
      carry = 1;							\
    }									\
}

static bignum_digit_type
DEFUN (bignum_digit_divide_subtract, (v1, v2, guess, u),
       bignum_digit_type v1 AND bignum_digit_type v2
       AND bignum_digit_type guess AND bignum_digit_type * u)
{
  {
    fast bignum_digit_type product;
    fast bignum_digit_type diff;
    fast bignum_digit_type carry;
    BDDS_MULSUB (v2, (u[2]), 0);
    BDDS_MULSUB (v1, (u[1]), carry);
    if (carry == 0)
      return (guess);
    diff = ((u[0]) - carry);
    if (diff < 0)
      (u[0]) = (diff + BIGNUM_RADIX);
    else
      {
	(u[0]) = diff;
	return (guess);
      }
  }
  {
    fast bignum_digit_type sum;
    fast bignum_digit_type carry;
    BDDS_ADD(v2, (u[2]), 0);
    BDDS_ADD(v1, (u[1]), carry);
    if (carry == 1)
      (u[0]) += 1;
  }
  return (guess - 1);
}

#undef BDDS_MULSUB
#undef BDDS_ADD

static void
DEFUN (bignum_divide_unsigned_small_denominator, (numerator, denominator,
						  quotient, remainder,
						  q_negative_p, r_negative_p),
       bignum_type numerator
       AND bignum_digit_type denominator
       AND bignum_type * quotient
       AND bignum_type * remainder
       AND int q_negative_p
       AND int r_negative_p)
{
  bignum_type q = (bignum_new_sign (numerator, q_negative_p));
  bignum_digit_type r = (bignum_destructive_scale_down (q, denominator));
  (*quotient) = (bignum_trim (q));
  if (remainder != ((bignum_type *) 0))
    (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
  return;
}

/* Given (denominator > 1), it is fairly easy to show that
   (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
   that all digits are < BIGNUM_RADIX. */

static bignum_digit_type
DEFUN (bignum_destructive_scale_down, (bignum, denominator),
       bignum_type bignum AND fast bignum_digit_type denominator)
{
  fast bignum_digit_type numerator;
  fast bignum_digit_type remainder = 0;
  fast bignum_digit_type two_digits;
#define quotient_high remainder
  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  BIGNUM_ASSERT ((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
  while (start < scan)
    {
      two_digits = (*--scan);
      numerator = (HD_CONS (remainder, (HD_HIGH (two_digits))));
      quotient_high = (numerator / denominator);
      numerator = (HD_CONS ((numerator % denominator), (HD_LOW (two_digits))));
      (*scan) = (HD_CONS (quotient_high, (numerator / denominator)));
      remainder = (numerator % denominator);
    }
  return (remainder);
#undef quotient_high
}

static bignum_type
DEFUN (bignum_remainder_unsigned_small_denominator, (n, d, negative_p),
       bignum_type n AND bignum_digit_type d AND int negative_p)
{
  fast bignum_digit_type two_digits;
  bignum_digit_type * start = (BIGNUM_START_PTR (n));
  fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (n)));
  fast bignum_digit_type r = 0;
  BIGNUM_ASSERT ((d > 1) && (d < BIGNUM_RADIX_ROOT));
  while (start < scan)
    {
      two_digits = (*--scan);
      r =
	((HD_CONS (((HD_CONS (r, (HD_HIGH (two_digits)))) % d),
		   (HD_LOW (two_digits))))
	 % d);
    }
  return (bignum_digit_to_bignum (r, negative_p));
}

static bignum_type
DEFUN (bignum_digit_to_bignum, (digit, negative_p),
       fast bignum_digit_type digit AND int negative_p)
{
  if (digit == 0)
    return (BIGNUM_ZERO ());
  else
    {
      fast bignum_type result = (bignum_allocate (1, negative_p));
      (BIGNUM_REF (result, 0)) = digit;
      return (result);
    }
}

/* Allocation */

static bignum_type
DEFUN (bignum_allocate, (length, negative_p),
       fast bignum_length_type length AND int negative_p)
{
  BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
  {
    fast bignum_type result = (BIGNUM_ALLOCATE (length));
    BIGNUM_SET_HEADER (result, length, negative_p);
    return (result);
  }
}

static bignum_type
DEFUN (bignum_allocate_zeroed, (length, negative_p),
       fast bignum_length_type length AND int negative_p)
{
  BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
  {
    fast bignum_type result = (BIGNUM_ALLOCATE (length));
    fast bignum_digit_type * scan = (BIGNUM_START_PTR (result));
    fast bignum_digit_type * end = (scan + length);
    BIGNUM_SET_HEADER (result, length, negative_p);
    while (scan < end)
      (*scan++) = 0;
    return (result);
  }
}

static bignum_type
DEFUN (bignum_shorten_length, (bignum, length),
       fast bignum_type bignum AND fast bignum_length_type length)
{
  fast bignum_length_type current_length = (BIGNUM_LENGTH (bignum));
  BIGNUM_ASSERT ((length >= 0) || (length <= current_length));
  if (length < current_length)
    {
      BIGNUM_SET_HEADER
	(bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
      BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
    }
  return (bignum);
}

static bignum_type
DEFUN (bignum_trim, (bignum), bignum_type bignum)
{
  fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  fast bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
  fast bignum_digit_type * scan = end;
  while ((start <= scan) && ((*--scan) == 0))
    ;
  scan += 1;
  if (scan < end)
    {
      fast bignum_length_type length = (scan - start);
      BIGNUM_SET_HEADER
	(bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
      BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
    }
  return (bignum);
}

/* Copying */

static bignum_type
DEFUN (bignum_copy, (source), fast bignum_type source)
{
  fast bignum_type target =
    (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source))));
  bignum_destructive_copy (source, target);
  return (target);
}

static bignum_type
DEFUN (bignum_new_sign, (bignum, negative_p),
       fast bignum_type bignum AND int negative_p)
{
  fast bignum_type result =
    (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
  bignum_destructive_copy (bignum, result);
  return (result);
}

static bignum_type
DEFUN (bignum_maybe_new_sign, (bignum, negative_p),
       fast bignum_type bignum AND int negative_p)
{
#ifndef BIGNUM_FORCE_NEW_RESULTS
  if ((BIGNUM_NEGATIVE_P (bignum)) ? negative_p : (! negative_p))
    return (bignum);
  else
#endif /* not BIGNUM_FORCE_NEW_RESULTS */
    {
      fast bignum_type result =
	(bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
      bignum_destructive_copy (bignum, result);
      return (result);
    }
}

static void
DEFUN (bignum_destructive_copy, (source, target),
       bignum_type source AND bignum_type target)
{
  fast bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
  fast bignum_digit_type * end_source =
    (scan_source + (BIGNUM_LENGTH (source)));
  fast bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
  while (scan_source < end_source)
    (*scan_target++) = (*scan_source++);
  return;
}