factor/vm/bignum.cpp

/*
   Copyright (C) 1989-94 Massachusetts Institute of Technology
   Portions copyright (C) 2004-2008 Slava Pestov

   This material was developed by the Scheme project at the Massachusetts
   Institute of Technology, Department of Electrical Engineering and
   Computer Science.  Permission to copy and modify this software, to
   redistribute either the original software or a modified version, and
   to use this software for any purpose is granted, subject to the
   following restrictions and understandings.

   1. Any copy made of this software must include this copyright notice
   in full.

   2. Users of this software agree to make their best efforts (a) to
   return to the MIT Scheme project any improvements or extensions that
   they make, so that these may be included in future releases; and (b)
   to inform MIT of noteworthy uses of this software.

   3. All materials developed as a consequence of the use of this
   software shall duly acknowledge such use, in accordance with the usual
   standards of acknowledging credit in academic research.

   4. MIT has made no warrantee or representation that the operation of
   this software will be error-free, and MIT is under no obligation to
   provide any services, by way of maintenance, update, or otherwise.

   5. In conjunction with products arising from the use of this material,
   there shall be no use of the name of the Massachusetts Institute of
   Technology nor of any adaptation thereof in any advertising,
   promotional, or sales literature without prior written consent from
   MIT in each case. */

/* Changes for Scheme 48:
 *  - Converted to ANSI.
 *  - Added bitwise operations.
 *  - Added s48 to the beginning of all externally visible names.
 *  - Cached the bignum representations of -1, 0, and 1.
 */

/* Changes for Factor:
 *  - Adapt bignumint.h for Factor memory manager
 *  - Add more bignum <-> C type conversions
 *  - Remove unused functions
 *  - Add local variable GC root recording
 *  - Remove s48 prefix from function names
 *  - Various fixes for Win64
 *  - Port to C++
 *  - Added bignum_gcd implementation
 */

#include "master.hpp"

namespace factor {

/* Exports */

int factor_vm::bignum_equal_p(bignum* x, bignum* 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 factor_vm::bignum_compare(bignum* x, bignum* 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))));
}

/* Allocates memory */
bignum* factor_vm::bignum_add(bignum* x, bignum* y) {
  return (
      (BIGNUM_ZERO_P(x)) ? (y) : (BIGNUM_ZERO_P(y))
          ? (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)))));
}

/* Allocates memory */
bignum* factor_vm::bignum_subtract(bignum* x, bignum* y) {
  return ((BIGNUM_ZERO_P(x))
              ? ((BIGNUM_ZERO_P(y)) ? (y) : (bignum_new_sign(
                                                y, (!(BIGNUM_NEGATIVE_P(y))))))
              : ((BIGNUM_ZERO_P(y))
                     ? (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))))));
}

#ifdef _WIN64
bignum *factor_vm::bignum_square(bignum* x_)
{
    return bignum_multiply(x_, x_);
}
#else
/* Allocates memory */
bignum *factor_vm::bignum_square(bignum* x_)
{
    data_root<bignum> x(x_, this);

    bignum_length_type length = (BIGNUM_LENGTH (x));
    bignum * z = (allot_bignum_zeroed ((length + length), 0));

    bignum_digit_type * scan_z = BIGNUM_START_PTR (z);
    bignum_digit_type * scan_x = BIGNUM_START_PTR (x);
    bignum_digit_type * end_x = scan_x + length;

    for (int i = 0; i < length; ++i) {
        bignum_twodigit_type carry;
        bignum_twodigit_type f = BIGNUM_REF (x, i);
        bignum_digit_type *pz = scan_z + (i << 1);
        bignum_digit_type *px = scan_x + i + 1;

        carry = *pz + f * f;
        *pz++ = carry & BIGNUM_DIGIT_MASK;
        carry >>= BIGNUM_DIGIT_LENGTH;
        BIGNUM_ASSERT (carry <= BIGNUM_DIGIT_MASK);

        f <<= 1;
        while (px < end_x)
        {
            carry += *pz + *px++ * f;
            *pz++ = carry & BIGNUM_DIGIT_MASK;
            carry >>= BIGNUM_DIGIT_LENGTH;
            BIGNUM_ASSERT (carry <= (BIGNUM_DIGIT_MASK << 1));
        }
        if (carry) {
            carry += *pz;
            *pz++ = carry & BIGNUM_DIGIT_MASK;
            carry >>= BIGNUM_DIGIT_LENGTH;
        }
        if (carry)
            *pz += carry & BIGNUM_DIGIT_MASK;
        BIGNUM_ASSERT ((carry >> BIGNUM_DIGIT_LENGTH) == 0);
    }
    return (bignum_trim (z));
}
#endif

/* Allocates memory */
bignum* factor_vm::bignum_multiply(bignum* x, bignum* y) {

#ifndef _WIN64
  if (x == y) {
    return bignum_square(x);
  }
#endif

  bignum_length_type x_length = (BIGNUM_LENGTH(x));
  bignum_length_type y_length = (BIGNUM_LENGTH(y));
  int negative_p = ((BIGNUM_NEGATIVE_P(x)) ? (!(BIGNUM_NEGATIVE_P(y)))
                                           : (BIGNUM_NEGATIVE_P(y)));
  if (BIGNUM_ZERO_P(x))
    return (x);
  if (BIGNUM_ZERO_P(y))
    return (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));
}

/* Allocates memory */
void factor_vm::bignum_divide(bignum* numerator, bignum* denominator,
                              bignum** quotient, bignum** remainder) {
  if (BIGNUM_ZERO_P(denominator)) {
    divide_by_zero_error();
    return;
  }
  if (BIGNUM_ZERO_P(numerator)) {
    (*quotient) = numerator;
    (*remainder) = 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) = 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;
      }
    }
  }
}

/* Allocates memory */
bignum* factor_vm::bignum_quotient(bignum* numerator, bignum* denominator) {
  if (BIGNUM_ZERO_P(denominator)) {
    divide_by_zero_error();
    return (BIGNUM_OUT_OF_BAND);
  }
  if (BIGNUM_ZERO_P(numerator))
    return 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:
      default: /* to appease gcc -Wall */
               {
        bignum* 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**)0), q_negative_p, 0);
          else
            bignum_divide_unsigned_medium_denominator(
                numerator, digit, (&quotient), ((bignum**)0), q_negative_p, 0);
        } else
          bignum_divide_unsigned_large_denominator(
              numerator, denominator, (&quotient), ((bignum**)0), q_negative_p,
              0);
        return (quotient);
      }
    }
  }
}

/* Allocates memory */
bignum* factor_vm::bignum_remainder(bignum* numerator, bignum* denominator) {
  if (BIGNUM_ZERO_P(denominator)) {
    divide_by_zero_error();
    return (BIGNUM_OUT_OF_BAND);
  }
  if (BIGNUM_ZERO_P(numerator))
    return numerator;
  switch (bignum_compare_unsigned(numerator, denominator)) {
    case bignum_comparison_equal:
      return (BIGNUM_ZERO());
    case bignum_comparison_less:
      return numerator;
    case bignum_comparison_greater:
    default: /* to appease gcc -Wall */
             {
      bignum* remainder;
      if ((BIGNUM_LENGTH(denominator)) == 1) {
        bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
        if (digit == 1)
          return (BIGNUM_ZERO());
        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**)0), (&remainder), 0,
            (BIGNUM_NEGATIVE_P(numerator)));
      } else
        bignum_divide_unsigned_large_denominator(
            numerator, denominator, ((bignum**)0), (&remainder), 0,
            (BIGNUM_NEGATIVE_P(numerator)));
      return (remainder);
    }
  }
}

/* cell_to_bignum, fixnum_to_bignum, long_long_to_bignum, ulong_long_to_bignum
 */
/* Allocates memory */
#define FOO_TO_BIGNUM(name, type, stype, utype)                       \
  bignum* factor_vm::name##_to_bignum(type n) {                       \
    int negative_p;                                                   \
    /* Special cases win when these small constants are cached. */    \
    if (n == 0)                                                       \
      return (BIGNUM_ZERO());                                         \
    if (n == 1)                                                       \
      return (BIGNUM_ONE(0));                                         \
    if (n < (type) 0 && n == (type) - 1)                              \
      return (BIGNUM_ONE(1));                                         \
    {                                                                 \
      utype accumulator =                                             \
          ((negative_p = (n < (type) 0)) ? ((type)(-(stype) n)) : n); \
      if (accumulator < BIGNUM_RADIX)                                 \
      {                                                               \
        bignum* result = allot_bignum(1, negative_p);                 \
        bignum_digit_type* scan = (BIGNUM_START_PTR(result));         \
        *scan = (accumulator & BIGNUM_DIGIT_MASK);                    \
        return (result);                                              \
      } else {                                                        \
        bignum_digit_type result_digits[BIGNUM_DIGITS_FOR(type)];     \
        bignum_digit_type* end_digits = result_digits;                \
        do {                                                          \
          (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);        \
          accumulator >>= BIGNUM_DIGIT_LENGTH;                        \
        } while (accumulator != 0);                                   \
        bignum* result =                                              \
           (allot_bignum((end_digits - result_digits), negative_p));  \
        bignum_digit_type* scan_digits = result_digits;               \
        bignum_digit_type* scan_result = (BIGNUM_START_PTR(result));  \
        while (scan_digits < end_digits)                              \
          (*scan_result++) = (*scan_digits++);                        \
        return (result);                                              \
      }                                                               \
    }                                                                 \
  }

FOO_TO_BIGNUM(cell, cell, fixnum, cell)
FOO_TO_BIGNUM(fixnum, fixnum, fixnum, cell)
FOO_TO_BIGNUM(long_long, int64_t, int64_t, uint64_t)
FOO_TO_BIGNUM(ulong_long, uint64_t, int64_t, uint64_t)

/* cannot allocate memory */
/* bignum_to_cell, fixnum_to_cell, long_long_to_cell, ulong_long_to_cell */
#define BIGNUM_TO_FOO(name, type, stype, utype)                            \
  type factor_vm::bignum_to_##name(bignum* bn) {                           \
    if (BIGNUM_ZERO_P(bn))                                                 \
      return (0);                                                          \
    {                                                                      \
      utype accumulator = 0;                                               \
      bignum_digit_type* start = (BIGNUM_START_PTR(bn));                   \
      bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));             \
      while (start < scan)                                                 \
        accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));  \
      return ((BIGNUM_NEGATIVE_P(bn)) ? ((type)(-(stype) accumulator))     \
                                      : accumulator);                      \
    }                                                                      \
  }

BIGNUM_TO_FOO(cell, cell, fixnum, cell)
BIGNUM_TO_FOO(fixnum, fixnum, fixnum, cell)
BIGNUM_TO_FOO(long_long, int64_t, int64_t, uint64_t)
BIGNUM_TO_FOO(ulong_long, uint64_t, int64_t, uint64_t)

/* cannot allocate memory */
fixnum factor_vm::bignum_to_fixnum_strict(bignum* bn) {
  fixnum len = BIGNUM_LENGTH(bn);
  bignum_digit_type *digits = BIGNUM_START_PTR(bn);
  if ((len == 1 && digits[0] > fixnum_max) || (len > 1)) {
    general_error(ERROR_OUT_OF_FIXNUM_RANGE, tag<bignum>(bn), false_object);
  }
  fixnum fix = bignum_to_fixnum(bn);
  FACTOR_ASSERT(fix <= fixnum_max && fix >= fixnum_min);
  return fix;
}

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

#define inf std::numeric_limits<double>::infinity()

/* Allocates memory */
bignum* factor_vm::double_to_bignum(double x) {
  if (x == inf || x == -inf || x != x)
    return (BIGNUM_ZERO());
  int exponent;
  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* result = (allot_bignum(length, (x < 0)));
    bignum_digit_type* start = (BIGNUM_START_PTR(result));
    bignum_digit_type* scan = (start + length);
    bignum_digit_type digit;
    int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
    if (odd_bits > 0)
      DTB_WRITE_DIGIT((fixnum)1 << 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

/* Comparisons */

int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) {
  bignum_length_type length = (BIGNUM_LENGTH(x));
  if (length != (BIGNUM_LENGTH(y)))
    return (0);
  else {
    bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
    bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
    bignum_digit_type* end_x = (scan_x + length);
    while (scan_x < end_x)
      if ((*scan_x++) != (*scan_y++))
        return (0);
    return (1);
  }
}

enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x,
                                                          bignum* 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);
  {
    bignum_digit_type* start_x = (BIGNUM_START_PTR(x));
    bignum_digit_type* scan_x = (start_x + x_length);
    bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length);
    while (start_x < scan_x) {
      bignum_digit_type digit_x = (*--scan_x);
      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 */

/* Allocates memory */
bignum* factor_vm::bignum_add_unsigned(bignum* x_, bignum* y_, int negative_p) {

  data_root<bignum> x(x_, this);
  data_root<bignum> y(y_, this);
  if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
    swap(x, y);
  }
  {
    bignum_length_type x_length = (BIGNUM_LENGTH(x));

    bignum* r = (allot_bignum((x_length + 1), negative_p));

    bignum_digit_type sum;
    bignum_digit_type carry = 0;
    bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
    bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
    {
      bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
      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;
        }
      }
    }
    {
      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 */

/* Allocates memory */
bignum* factor_vm::bignum_subtract_unsigned(bignum* x_, bignum* y_) {

  data_root<bignum> x(x_, this);
  data_root<bignum> y(y_, this);

  int negative_p = 0;
  switch (bignum_compare_unsigned(x.untagged(), y.untagged())) {
    case bignum_comparison_equal:
      return (BIGNUM_ZERO());
    case bignum_comparison_less:
      swap(x, y);
      negative_p = 1;
      break;
    case bignum_comparison_greater:
      negative_p = 0;
      break;
  }
  {
    bignum_length_type x_length = (BIGNUM_LENGTH(x));

    bignum* r = (allot_bignum(x_length, negative_p));

    bignum_digit_type difference;
    bignum_digit_type borrow = 0;
    bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
    bignum_digit_type* scan_r = BIGNUM_START_PTR(r);
    {
      bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
      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;
        }
      }
    }
    {
      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 */

/* Allocates memory */
bignum* factor_vm::bignum_multiply_unsigned(bignum* x_, bignum* y_,
                                            int negative_p) {

  data_root<bignum> x(x_, this);
  data_root<bignum> y(y_, this);

  if (BIGNUM_LENGTH(y) > BIGNUM_LENGTH(x)) {
    swap(x, y);
  }
  {
    bignum_digit_type carry;
    bignum_digit_type y_digit_low;
    bignum_digit_type y_digit_high;
    bignum_digit_type x_digit_low;
    bignum_digit_type x_digit_high;
    bignum_digit_type product_low;
    bignum_digit_type* scan_r;
    bignum_digit_type* scan_y;
    bignum_length_type x_length = BIGNUM_LENGTH(x);
    bignum_length_type y_length = BIGNUM_LENGTH(y);

    bignum* r = (allot_bignum_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
  }
}

/* Allocates memory */
bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x_,
                                                         bignum_digit_type y,
                                                         int negative_p) {
  data_root<bignum> x(x_, this);

  bignum_length_type length_x = (BIGNUM_LENGTH(x));

  bignum* p = (allot_bignum((length_x + 1), negative_p));

  bignum_destructive_copy(x.untagged(), p);
  (BIGNUM_REF(p, length_x)) = 0;
  bignum_destructive_scale_up(p, y);
  return (bignum_trim(p));
}

void factor_vm::bignum_destructive_add(bignum* bn, bignum_digit_type n) {
  bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
  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);
  }
}

void factor_vm::bignum_destructive_scale_up(bignum* bn,
                                            bignum_digit_type factor) {
  bignum_digit_type carry = 0;
  bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
  bignum_digit_type two_digits;
  bignum_digit_type product_low;
#define product_high carry
  bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bn)));
  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
}

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

/* Allocates memory */
void factor_vm::bignum_divide_unsigned_large_denominator(
    bignum* numerator_, bignum* denominator_,
    bignum** quotient, bignum** remainder,
    int q_negative_p, int r_negative_p) {

  data_root<bignum> numerator(numerator_, this);
  data_root<bignum> denominator(denominator_, this);

  bignum_length_type length_n = BIGNUM_LENGTH(numerator) + 1;
  bignum_length_type length_d = BIGNUM_LENGTH(denominator);

  data_root<bignum> u(allot_bignum(length_n, r_negative_p), this);

  int shift = 0;
  BIGNUM_ASSERT(length_d > 1);
  {
    bignum_digit_type v1 = BIGNUM_REF(denominator.untagged(), length_d - 1);
    while (v1 < (BIGNUM_RADIX / 2)) {
      v1 <<= 1;
      shift += 1;
    }
  }

  if (quotient != NULL) {
    bignum *q_ = allot_bignum(length_n - length_d, q_negative_p);
    data_root<bignum> q(q_, this);

    if (shift == 0) {
      bignum_destructive_copy(numerator.untagged(), u.untagged());
      (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0;
      bignum_divide_unsigned_normalized(u.untagged(),
                                        denominator.untagged(),
                                        q.untagged());
    } else {
      bignum* v = allot_bignum(length_d, 0);
      bignum_destructive_normalization(numerator.untagged(),
                                       u.untagged(),
                                       shift);
      bignum_destructive_normalization(denominator.untagged(), v, shift);
      bignum_divide_unsigned_normalized(u.untagged(), v, q.untagged());
      if (remainder != NULL)
        bignum_destructive_unnormalization(u.untagged(), shift);
    }

    q = bignum_trim(q.untagged());
    *quotient = q.untagged();
  } else {

    if (shift == 0) {
        bignum_destructive_copy(numerator.untagged(), u.untagged());
        (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0;
        bignum_divide_unsigned_normalized(u.untagged(),
                                          denominator.untagged(),
                                          NULL);
      } else {
        bignum* v = allot_bignum(length_d, 0);
        bignum_destructive_normalization(numerator.untagged(),
                                         u.untagged(),
                                         shift);
        bignum_destructive_normalization(denominator.untagged(),
                                         v,
                                         shift);
        bignum_divide_unsigned_normalized(u.untagged(), v, NULL);
        if (remainder != NULL)
          bignum_destructive_unnormalization(u.untagged(), shift);
      }
  }

  u = bignum_trim(u.untagged());
  if (remainder != NULL)
    *remainder = u.untagged();
}

void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v,
                                                  bignum* 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 = NULL;
  bignum_digit_type v1 = (v_end[-1]);
  bignum_digit_type v2 = (v_end[-2]);
  bignum_digit_type ph; /* high half of double-digit product */
  bignum_digit_type pl; /* low half of double-digit product */
  bignum_digit_type guess;
  bignum_digit_type gh; /* high half-digit of guess */
  bignum_digit_type ch; /* high half of double-digit comparand */
  bignum_digit_type v2l = (HD_LOW(v2));
  bignum_digit_type v2h = (HD_HIGH(v2));
  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
}

bignum_digit_type factor_vm::bignum_divide_subtract(
    bignum_digit_type* v_start, bignum_digit_type* v_end,
    bignum_digit_type guess, bignum_digit_type* u_start) {
  bignum_digit_type* v_scan = v_start;
  bignum_digit_type* u_scan = u_start;
  bignum_digit_type carry = 0;
  if (guess == 0)
    return (0);
  {
    bignum_digit_type gl = (HD_LOW(guess));
    bignum_digit_type gh = (HD_HIGH(guess));
    bignum_digit_type v;
    bignum_digit_type pl;
    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);
}

/* Allocates memory */
void factor_vm::bignum_divide_unsigned_medium_denominator(
    bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
    bignum** remainder, int q_negative_p, int r_negative_p) {

  data_root<bignum> numerator(numerator_, this);

  bignum_length_type length_n = (BIGNUM_LENGTH(numerator));

  int shift = 0;
  /* Because `bignum_digit_divide' requires a normalized denominator. */
  while (denominator < (BIGNUM_RADIX / 2)) {
    denominator <<= 1;
    shift += 1;
  }

  bignum_length_type length_q = (shift == 0) ? length_n : length_n + 1;
  data_root<bignum> q(allot_bignum(length_q, q_negative_p), this);
  if (shift == 0) {
    bignum_destructive_copy(numerator.untagged(), q.untagged());
  } else {
    bignum_destructive_normalization(numerator.untagged(), q.untagged(), shift);
  }
  {
    bignum_digit_type r = 0;
    bignum_digit_type* start = (BIGNUM_START_PTR(q));
    bignum_digit_type* scan = (start + length_q);
    bignum_digit_type qj;

    while (start < scan) {
      r = (bignum_digit_divide(r, (*--scan), denominator, (&qj)));
      (*scan) = qj;
    }

    q = bignum_trim(q.untagged());

    if (remainder != ((bignum**)0)) {
      if (shift != 0)
        r >>= shift;

      (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
    }

    if (quotient != ((bignum**)0))
      (*quotient) = q.untagged();
  }
  return;
}

void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target,
                                                 int shift_left) {
  bignum_digit_type digit;
  bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
  bignum_digit_type carry = 0;
  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 = (((cell)1 << 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;
}

void factor_vm::bignum_destructive_unnormalization(bignum* bn,
                                                   int shift_right) {
  bignum_digit_type* start = (BIGNUM_START_PTR(bn));
  bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
  bignum_digit_type digit;
  bignum_digit_type carry = 0;
  int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
  bignum_digit_type mask = (((fixnum)1 << 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, 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]))); \
  }

bignum_digit_type factor_vm::bignum_digit_divide(
    bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v,
    bignum_digit_type* q) /* return value */
    {
  bignum_digit_type guess;
  bignum_digit_type comparand;
  bignum_digit_type v1 = (HD_HIGH(v));
  bignum_digit_type v2 = (HD_LOW(v));
  bignum_digit_type uj;
  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;                      \
    }                                 \
  }

bignum_digit_type factor_vm::bignum_digit_divide_subtract(
    bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess,
    bignum_digit_type* u) {
  {
    bignum_digit_type product;
    bignum_digit_type diff;
    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);
    }
  }
  {
    bignum_digit_type sum;
    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

/* Allocates memory */
void factor_vm::bignum_divide_unsigned_small_denominator(
    bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
    bignum** remainder, int q_negative_p, int r_negative_p) {
  data_root<bignum> numerator(numerator_, this);

  bignum* q_ = bignum_new_sign(numerator.untagged(), q_negative_p);
  data_root<bignum> q(q_, this);

  bignum_digit_type r = bignum_destructive_scale_down(q.untagged(), denominator);

  q = bignum_trim(q.untagged());

  if (remainder != ((bignum**)0))
    (*remainder) = bignum_digit_to_bignum(r, r_negative_p);

  (*quotient) = q.untagged();

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

bignum_digit_type factor_vm::bignum_destructive_scale_down(
    bignum* bn, bignum_digit_type denominator) {
  bignum_digit_type numerator;
  bignum_digit_type remainder = 0;
  bignum_digit_type two_digits;
#define quotient_high remainder
  bignum_digit_type* start = (BIGNUM_START_PTR(bn));
  bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
  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
}

/* Allocates memory */
bignum* factor_vm::bignum_remainder_unsigned_small_denominator(
    bignum* n, bignum_digit_type d, int negative_p) {
  bignum_digit_type two_digits;
  bignum_digit_type* start = (BIGNUM_START_PTR(n));
  bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n)));
  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));
}

/* Allocates memory */
bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit,
                                          int negative_p) {
  if (digit == 0)
    return (BIGNUM_ZERO());
  else {
    bignum* result = (allot_bignum(1, negative_p));
    (BIGNUM_REF(result, 0)) = digit;
    return (result);
  }
}

/* Allocates memory */
bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) {
  BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX));
  bignum* result = allot_uninitialized_array<bignum>(length + 1);
  BIGNUM_SET_NEGATIVE_P(result, negative_p);
  return (result);
}

/* Allocates memory */
bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length,
                                       int negative_p) {
  bignum* result = allot_bignum(length, negative_p);
  bignum_digit_type* scan = (BIGNUM_START_PTR(result));
  bignum_digit_type* end = (scan + length);
  while (scan < end)
    (*scan++) = 0;
  return (result);
}

/* Allocates memory */
bignum* factor_vm::bignum_shorten_length(bignum* bn,
                                         bignum_length_type length) {
  bignum_length_type current_length = (BIGNUM_LENGTH(bn));
  BIGNUM_ASSERT((length >= 0) || (length <= current_length));
  if (length < current_length) {
    bn = reallot_array(bn, length + 1);
    BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
  }
  return (bn);
}

/* Allocates memory */
bignum* factor_vm::bignum_trim(bignum* bn) {
  bignum_digit_type* start = (BIGNUM_START_PTR(bn));
  bignum_digit_type* end = (start + (BIGNUM_LENGTH(bn)));
  bignum_digit_type* scan = end;
  while ((start <= scan) && ((*--scan) == 0))
    ;
  scan += 1;
  if (scan < end) {
    bignum_length_type length = (scan - start);
    bn = reallot_array(bn, length + 1);
    BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
  }
  return (bn);
}

/* Copying */

/* Allocates memory */
bignum* factor_vm::bignum_new_sign(bignum* x_, int negative_p) {
  data_root<bignum> x(x_, this);
  bignum* result = allot_bignum(BIGNUM_LENGTH(x), negative_p);
  bignum_destructive_copy(x.untagged(), result);
  return result;
}

/* Allocates memory */
bignum* factor_vm::bignum_maybe_new_sign(bignum* x_, int negative_p) {
  if ((BIGNUM_NEGATIVE_P(x_)) ? negative_p : (!negative_p))
    return x_;
  else {
    return bignum_new_sign(x_, negative_p);
  }
}

void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) {
  bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
  bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
  bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
  while (scan_source < end_source)
    (*scan_target++) = (*scan_source++);
  return;
}

/*
 * Added bitwise operations (and oddp).
 */

/* Allocates memory */
bignum* factor_vm::bignum_bitwise_not(bignum* x_) {

  int carry = 1;
  bignum_length_type size = BIGNUM_LENGTH(x_);
  int is_negative = BIGNUM_NEGATIVE_P(x_);
  data_root<bignum> x(x_, this);
  data_root<bignum> y(allot_bignum(size, is_negative ? 0 : 1), this);

  bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
  bignum_digit_type* end_x = scan_x + size;
  bignum_digit_type* scan_y = BIGNUM_START_PTR(y);

  if (is_negative) {
    while (scan_x < end_x) {
      if (*scan_x == 0) {
        *scan_y++ = BIGNUM_RADIX - 1;
        scan_x++;
      } else {
        *scan_y++ = *scan_x++ - 1;
        carry = 0;
        break;
      }
    }
  } else {
    while (scan_x < end_x) {
      if (*scan_x == (BIGNUM_RADIX - 1)) {
        *scan_y++ = 0;
        scan_x++;
      } else {
        *scan_y++ = *scan_x++ + 1;
        carry = 0;
        break;
      }
    }
  }

  while (scan_x < end_x) {
    *scan_y++ = *scan_x++;
  }

  if (carry) {
    bignum* ret = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y));
    bignum_destructive_copy(y.untagged(), ret);
    bignum_digit_type* ret_start = BIGNUM_START_PTR(ret);
    *(ret_start + size) = 1;
    return ret;
  } else {
    return bignum_trim(y.untagged());
  }
}

/* Allocates memory */
bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) {
  if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
    return bignum_bitwise_not(
        bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
  else
    return bignum_magnitude_ash(arg1, n);
}

#define AND_OP 0
#define IOR_OP 1
#define XOR_OP 2

/* Allocates memory */
bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) {
  return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
              ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
              : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
              : (BIGNUM_NEGATIVE_P(arg2))
              ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
              : bignum_pospos_bitwise_op(AND_OP, arg1, arg2));
}

/* Allocates memory */
bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) {
  return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
              ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
              : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
              : (BIGNUM_NEGATIVE_P(arg2))
              ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
              : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2));
}

/* Allocates memory */
bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) {
  return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
              ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
              : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
              : (BIGNUM_NEGATIVE_P(arg2))
              ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
              : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2));
}

/* Allocates memory */
/* ash for the magnitude */
/* assume arg1 is a big number, n is a long */
bignum* factor_vm::bignum_magnitude_ash(bignum* arg1_, fixnum n) {

  data_root<bignum> arg1(arg1_, this);

  bignum* result = NULL;
  bignum_digit_type* scan1;
  bignum_digit_type* scanr;
  bignum_digit_type* end;

  fixnum digit_offset, bit_offset;

  if (BIGNUM_ZERO_P(arg1))
    return arg1.untagged();

  if (n > 0) {
    digit_offset = n / BIGNUM_DIGIT_LENGTH;
    bit_offset = n % BIGNUM_DIGIT_LENGTH;

    result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1,
                                 BIGNUM_NEGATIVE_P(arg1));

    scanr = BIGNUM_START_PTR(result) + digit_offset;
    scan1 = BIGNUM_START_PTR(arg1);
    end = scan1 + BIGNUM_LENGTH(arg1);

    while (scan1 < end) {
      *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
      *scanr = *scanr & BIGNUM_DIGIT_MASK;
      scanr++;
      *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
      *scanr = *scanr & BIGNUM_DIGIT_MASK;
    }
  } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type)
                              BIGNUM_DIGIT_LENGTH))) {
    result = BIGNUM_ZERO();
  } else if (n < 0) {
    digit_offset = -n / BIGNUM_DIGIT_LENGTH;
    bit_offset = -n % BIGNUM_DIGIT_LENGTH;

    result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset,
                                 BIGNUM_NEGATIVE_P(arg1));

    scanr = BIGNUM_START_PTR(result);
    scan1 = BIGNUM_START_PTR(arg1) + digit_offset;
    end = scanr + BIGNUM_LENGTH(result) - 1;

    while (scanr < end) {
      *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
      *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) &
               BIGNUM_DIGIT_MASK;
      scanr++;
    }
    *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
  } else if (n == 0) {
    result = arg1.untagged();
  }

  return bignum_trim(result);
}

/* Allocates memory */
bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1_,
                                            bignum* arg2_) {
  data_root<bignum> arg1(arg1_, this);
  data_root<bignum> arg2(arg2_, this);

  bignum_length_type max_length;

  bignum_digit_type* scan1, *end1, digit1;
  bignum_digit_type* scan2, *end2, digit2;
  bignum_digit_type* scanr, *endr;

  max_length =
      (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1)
                                                  : BIGNUM_LENGTH(arg2);

  bignum* result = allot_bignum(max_length, 0);

  scanr = BIGNUM_START_PTR(result);
  scan1 = BIGNUM_START_PTR(arg1);
  scan2 = BIGNUM_START_PTR(arg2);
  endr = scanr + max_length;
  end1 = scan1 + BIGNUM_LENGTH(arg1);
  end2 = scan2 + BIGNUM_LENGTH(arg2);

  while (scanr < endr) {
    digit1 = (scan1 < end1) ? *scan1++ : 0;
    digit2 = (scan2 < end2) ? *scan2++ : 0;
    *scanr++ =
        (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
                                                          : digit1 ^ digit2;
  }
  return bignum_trim(result);
}

/* Allocates memory */
bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1_,
                                            bignum* arg2_) {
  data_root<bignum> arg1(arg1_, this);
  data_root<bignum> arg2(arg2_, this);

  bignum_length_type max_length;

  bignum_digit_type* scan1, *end1, digit1;
  bignum_digit_type* scan2, *end2, digit2, carry2;
  bignum_digit_type* scanr, *endr;

  char neg_p = op == IOR_OP || op == XOR_OP;

  max_length =
      (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1)
                                                      : BIGNUM_LENGTH(arg2) + 1;

  bignum* result = allot_bignum(max_length, neg_p);

  scanr = BIGNUM_START_PTR(result);
  scan1 = BIGNUM_START_PTR(arg1);
  scan2 = BIGNUM_START_PTR(arg2);
  endr = scanr + max_length;
  end1 = scan1 + BIGNUM_LENGTH(arg1);
  end2 = scan2 + BIGNUM_LENGTH(arg2);

  carry2 = 1;

  while (scanr < endr) {
    digit1 = (scan1 < end1) ? *scan1++ : 0;
    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;

    if (digit2 < BIGNUM_RADIX)
      carry2 = 0;
    else {
      digit2 = (digit2 - BIGNUM_RADIX);
      carry2 = 1;
    }

    *scanr++ =
        (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
                                                          : digit1 ^ digit2;
  }

  if (neg_p)
    bignum_negate_magnitude(result);

  return bignum_trim(result);
}

/* Allocates memory */
bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1_,
                                            bignum* arg2_) {
  data_root<bignum> arg1(arg1_, this);
  data_root<bignum> arg2(arg2_, this);

  bignum_length_type max_length;

  bignum_digit_type* scan1, *end1, digit1, carry1;
  bignum_digit_type* scan2, *end2, digit2, carry2;
  bignum_digit_type* scanr, *endr;

  char neg_p = op == AND_OP || op == IOR_OP;

  max_length =
      (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1
                                                  : BIGNUM_LENGTH(arg2) + 1;

  bignum* result = allot_bignum(max_length, neg_p);

  scanr = BIGNUM_START_PTR(result);
  scan1 = BIGNUM_START_PTR(arg1);
  scan2 = BIGNUM_START_PTR(arg2);
  endr = scanr + max_length;
  end1 = scan1 + BIGNUM_LENGTH(arg1);
  end2 = scan2 + BIGNUM_LENGTH(arg2);

  carry1 = 1;
  carry2 = 1;

  while (scanr < endr) {
    digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
    digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;

    if (digit1 < BIGNUM_RADIX)
      carry1 = 0;
    else {
      digit1 = (digit1 - BIGNUM_RADIX);
      carry1 = 1;
    }

    if (digit2 < BIGNUM_RADIX)
      carry2 = 0;
    else {
      digit2 = (digit2 - BIGNUM_RADIX);
      carry2 = 1;
    }

    *scanr++ =
        (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
                                                          : digit1 ^ digit2;
  }

  if (neg_p)
    bignum_negate_magnitude(result);

  return bignum_trim(result);
}

void factor_vm::bignum_negate_magnitude(bignum* arg) {
  bignum_digit_type* scan;
  bignum_digit_type* end;
  bignum_digit_type digit;
  bignum_digit_type carry;

  scan = BIGNUM_START_PTR(arg);
  end = scan + BIGNUM_LENGTH(arg);

  carry = 1;

  while (scan < end) {
    digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry;

    if (digit < BIGNUM_RADIX)
      carry = 0;
    else {
      digit = (digit - BIGNUM_RADIX);
      carry = 1;
    }

    *scan++ = digit;
  }
}

/* Allocates memory */
bignum* factor_vm::bignum_integer_length(bignum* x_) {
  data_root<bignum> x(x_, this);
  bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1);
  bignum_digit_type digit = (BIGNUM_REF(x, index));
  bignum_digit_type carry = 0;
  bignum* result;

  while (digit > 1) {
    carry += 1;
    digit >>= 1;
  }

  if (index < BIGNUM_RADIX_ROOT) {
     result = allot_bignum(1, 0);
     (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry;
  } else {
     result = allot_bignum(2, 0);
     (BIGNUM_REF(result, 0)) = index;
     (BIGNUM_REF(result, 1)) = 0;
     bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH);
     bignum_destructive_add(result, carry);
  }
  return (bignum_trim(result));
}

/* Allocates memory */
int factor_vm::bignum_logbitp(int shift, bignum* arg) {
  return ((BIGNUM_NEGATIVE_P(arg))
              ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg))
              : bignum_unsigned_logbitp(shift, arg));
}

int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bn) {
  bignum_length_type len = (BIGNUM_LENGTH(bn));
  int index = shift / BIGNUM_DIGIT_LENGTH;
  if (index >= len)
    return 0;
  bignum_digit_type digit = (BIGNUM_REF(bn, index));
  int p = shift % BIGNUM_DIGIT_LENGTH;
  bignum_digit_type mask = ((fixnum)1) << p;
  return (digit & mask) ? 1 : 0;
}

#ifdef _WIN64
/* Allocates memory. */
bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {

  data_root<bignum> a(a_, this);
  data_root<bignum> b(b_, this);

  /* Copies of a and b with that are both positive. */
  data_root<bignum> ac(bignum_maybe_new_sign(a.untagged(), 0), this);
  data_root<bignum> bc(bignum_maybe_new_sign(b.untagged(), 0), this);

  if (bignum_compare(ac.untagged(), bc.untagged()) == bignum_comparison_less) {
    swap(ac, bc);
  }

  while (BIGNUM_LENGTH(bc) != 0) {
    data_root<bignum> d(bignum_remainder(ac.untagged(), bc.untagged()), this);
    if (d.untagged() == BIGNUM_OUT_OF_BAND) {
      return d.untagged();
    }
    ac = bc;
    bc = d;
  }
  return ac.untagged();
}
#else
/* Allocates memory */
bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
  data_root<bignum> a(a_, this);
  data_root<bignum> b(b_, this);
  bignum_twodigit_type x, y, q, s, t, A, B, C, D;
  int nbits, k;
  bignum_length_type size_a, size_b, size_c;
  bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d;
  bignum_digit_type* a_end, *b_end, *c_end;

  /* clone the bignums so we can modify them in-place */
  size_a = BIGNUM_LENGTH(a);
  data_root<bignum> c(allot_bignum(size_a, 0), this);
  // c = allot_bignum(size_a, 0);
  scan_a = BIGNUM_START_PTR(a);
  a_end = scan_a + size_a;
  scan_c = BIGNUM_START_PTR(c);
  while (scan_a < a_end)
    (*scan_c++) = (*scan_a++);
  a = c;
  size_b = BIGNUM_LENGTH(b);
  data_root<bignum> d(allot_bignum(size_b, 0), this);
  scan_b = BIGNUM_START_PTR(b);
  b_end = scan_b + size_b;
  scan_d = BIGNUM_START_PTR(d);
  while (scan_b < b_end)
    (*scan_d++) = (*scan_b++);
  b = d;

  /* Initial reduction: make sure that 0 <= b <= a. */
  if (bignum_compare(a.untagged(), b.untagged()) == bignum_comparison_less) {
    swap(a, b);
    std::swap(size_a, size_b);
  }

  while (size_a > 1) {
    nbits = log2(BIGNUM_REF(a, size_a - 1));
    x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) |
         (BIGNUM_REF(a, size_a - 2) >> nbits));
    y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) |
         (size_b >= size_a
              ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)
              : 0));

    /* inner loop of Lehmer's algorithm; */
    A = 1;
    B = 0;
    C = 0;
    D = 1;
    for (k = 0;; k++) {
      if (y - C == 0)
        break;

      q = (x + (A - 1)) / (y - C);

      s = B + (q * D);
      t = x - (q * y);

      if (s > t)
        break;

      x = y;
      y = t;

      t = A + (q * C);

      A = D;
      B = C;
      C = s;
      D = t;
    }

    if (k == 0) {
      /* no progress; do a Euclidean step */
      if (size_b == 0) {
        return bignum_trim(a.untagged());
      }
      data_root<bignum> e(bignum_trim(a.untagged()), this);
      data_root<bignum> f(bignum_trim(b.untagged()), this);
      c = bignum_remainder(e.untagged(), f.untagged());
      if (c.untagged() == BIGNUM_OUT_OF_BAND) {
        return c.untagged();
      }

      // copy 'b' to 'a'
      scan_a = BIGNUM_START_PTR(a);
      scan_b = BIGNUM_START_PTR(b);
      a_end = scan_a + size_a;
      b_end = scan_b + size_b;
      while (scan_b < b_end)
        *(scan_a++) = *(scan_b++);
      while (scan_a < a_end)
        *(scan_a++) = 0;
      size_a = size_b;

      // copy 'c' to 'b'
      scan_b = BIGNUM_START_PTR(b);
      scan_c = BIGNUM_START_PTR(c);
      size_c = BIGNUM_LENGTH(c);
      c_end = scan_c + size_c;
      while (scan_c < c_end)
        *(scan_b++) = *(scan_c++);
      while (scan_b < b_end)
        *(scan_b++) = 0;
      size_b = size_c;

      continue;
    }

    /*
      a, b = A*b - B*a, D*a - C*b if k is odd
      a, b = A*a - B*b, D*b - C*a if k is even
    */
    scan_a = BIGNUM_START_PTR(a);
    scan_b = BIGNUM_START_PTR(b);
    scan_c = scan_a;
    scan_d = scan_b;
    a_end = scan_a + size_a;
    b_end = scan_b + size_b;
    s = 0;
    t = 0;
    if (k & 1) {
      while (scan_b < b_end) {
        s += (A * *scan_b) - (B * *scan_a);
        t += (D * *scan_a++) - (C * *scan_b++);
        *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
        *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
        s >>= BIGNUM_DIGIT_LENGTH;
        t >>= BIGNUM_DIGIT_LENGTH;
      }
      while (scan_a < a_end) {
        s -= (B * *scan_a);
        t += (D * *scan_a++);
        *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
        //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
        s >>= BIGNUM_DIGIT_LENGTH;
        t >>= BIGNUM_DIGIT_LENGTH;
      }
    } else {
      while (scan_b < b_end) {
        s += (A * *scan_a) - (B * *scan_b);
        t += (D * *scan_b++) - (C * *scan_a++);
        *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
        *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
        s >>= BIGNUM_DIGIT_LENGTH;
        t >>= BIGNUM_DIGIT_LENGTH;
      }
      while (scan_a < a_end) {
        s += (A * *scan_a);
        t -= (C * *scan_a++);
        *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
        //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
        s >>= BIGNUM_DIGIT_LENGTH;
        t >>= BIGNUM_DIGIT_LENGTH;
      }
    }
    BIGNUM_ASSERT(s == 0);
    BIGNUM_ASSERT(t == 0);

    // update size_a and size_b to remove any zeros at end
    while (size_a > 0 && *(--scan_a) == 0)
      size_a--;
    while (size_b > 0 && *(--scan_b) == 0)
      size_b--;

    BIGNUM_ASSERT(size_a >= size_b);
  }

  /* a fits into a fixnum, so b must too */
  fixnum xx = bignum_to_fixnum(a.untagged());
  fixnum yy = bignum_to_fixnum(b.untagged());
  fixnum tt;

  /* usual Euclidean algorithm for longs */
  while (yy != 0) {
    tt = yy;
    yy = xx % yy;
    xx = tt;
  }

  return fixnum_to_bignum(xx);
}
#endif

}