source: project/release/4/numbers/trunk/numbers-c.c @ 31291

/* numbers-c.c
 *
 * Copyright 2010-2014 The CHICKEN Team
 *
 * This contains a barely recognizable version of c/bignum.c from Scheme48 1.8:
 * Copyright (c) 1993-2008 Richard Kelsey and Jonathan Rees
 * Copyright 1986,1987,1988,1989,1990,1991 Massachusetts Institute of Technology
 * Copyright 1992,1993,1994,2004 Massachusetts Institute of Technology
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 *
 *    1. Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 * 
 *    2. Redistributions in binary form must reproduce the above
 *       copyright notice, this list of conditions and the following
 *       disclaimer in the documentation and/or other materials provided
 *       with the distribution.
 * 
 *    3. The name of the author may not be used to endorse or promote
 *       products derived from this software without specific prior
 *       written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
 * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include <assert.h>
#include <errno.h>

static void *tags;

#include "numbers-c.h"

#define fix_to_flo(p, n, f)       C_flonum(p, C_unfix(f))
#define big_of(v)                 ((bignum_type)C_data_pointer(C_block_item(v, 1)))

static C_word
init_tags(___scheme_value tagvec)
{
  tags = CHICKEN_new_gc_root();
  CHICKEN_gc_root_set(tags, tagvec);
  return C_SCHEME_UNDEFINED;
}

/*
 * This is an odd one out. It doesn't accept a continuation.
 * I've put in the (more or less) verbatim code for s48_bignum_to_double.
 */
static C_word
big_to_flo(C_word **p, C_word n, C_word value)
{
  bignum_type bignum = big_of(value);
  double accumulator = 0;
  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  while (start < scan)
    accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
  return C_flonum(p, (BIGNUM_NEGATIVE_P(bignum) ? -accumulator : accumulator));
}


static C_word
check_number(C_word x)
{
  C_word tagvec;

  if((x & C_FIXNUM_BIT) != 0) return C_fix(FIX);
  else if(C_immediatep(x)) return C_fix(NONE);
  else if(C_header_bits(x) == C_FLONUM_TYPE) return C_fix(FLO);
  else if(C_header_bits(x) == C_STRUCTURE_TYPE) {
    tagvec = CHICKEN_gc_root_ref(tags);
    
    if (C_block_item(x, 0) == C_block_item(tagvec, BIG_TAG))
      return C_fix(BIG);
    else if (C_block_item(x, 0) == C_block_item(tagvec, COMP_TAG))
      return C_fix(COMP);
    else if (C_block_item(x, 0) == C_block_item(tagvec, RAT_TAG))
      return C_fix(RAT);
    else
      return C_fix(NONE);
  } else
      return C_fix(NONE);
}

static void
C_wrap_bignum(C_word c, C_word closure, C_word bigvec)
{
  C_word ab[3], *a = ab;
  bignum_type src = (bignum_type)C_block_item(C_block_item(closure, 2), 0);
  bignum_digit_type *scan = BIGNUM_TO_POINTER(src);
  bignum_digit_type *end = BIGNUM_START_PTR(src) + BIGNUM_LENGTH(src);
  bignum_digit_type *tgt = C_data_pointer(bigvec);
  C_word tagvec = CHICKEN_gc_root_ref(tags);
  C_word result = C_structure(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
  while (scan < end)
    (*tgt++) = (*scan++);
  BIGNUM_DEALLOCATE(src);
  C_kontinue(C_block_item(closure, 1), result);
}

static void
C_return_bignum(C_word k, bignum_type b)
{
  C_word result;
  
  assert(b != BIGNUM_OUT_OF_BAND);

  switch(BIGNUM_LENGTH(b)) {
  case 0:
    BIGNUM_DEALLOCATE(b);
    C_kontinue(k, C_fix(0));
  /* This code "knows" that bignums have 2 "reserved" bits, like fixnums */
  case 1:
    result = C_fix(BIGNUM_NEGATIVE_P(b) ? -BIGNUM_REF(b, 0) : BIGNUM_REF(b, 0));
    BIGNUM_DEALLOCATE(b);
    C_kontinue(k, result);
  case 2:
    /* Edge case: most negative fixnum */
    if (BIGNUM_NEGATIVE_P(b) && BIGNUM_REF(b, 0) == 0 && BIGNUM_REF(b, 1) == 1) {
      BIGNUM_DEALLOCATE(b);
      C_kontinue(k, C_fix(C_MOST_NEGATIVE_FIXNUM));
    }
    /* FALLTHROUGH */
  default:
    {
      C_word pab[2], *pa = pab, kab[4], *ka = kab, k2;
      /* Make result a wrapped pointer because C_closure wants scheme objects */
      result = C_mpointer(&pa, b);
      k2 = C_closure(&ka, 3, (C_word)C_wrap_bignum, k, result);
      /* Here we assume bignum digits are C words.. */
      C_allocate_vector(6, (C_word)NULL, k2,
                        C_fix(sizeof(C_word) * (BIGNUM_LENGTH(b) + 1)),
                        /* Byte vec, no initialization, align at 8 bytes */
                        C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
    }
  }
}

static void
C_bignum_wrapped_return_bigobj(C_word c, C_word closure, C_word wrapped_big)
{
  C_word obj = C_block_item(closure, 2);
  C_values(4, C_SCHEME_UNDEFINED, C_block_item(closure, 1), wrapped_big, obj);
}

static void
C_return_big_fix(C_word k, bignum_type big, C_word fix)
{
  C_word kab[4], *ka = kab, k2;
  k2 = C_closure(&ka, 3, (C_word)C_bignum_wrapped_return_bigobj, k, fix);
  C_return_bignum(k2, big);
}

static void
C_b2_wrapped(C_word c, C_word closure, C_word wrapped_b2)
{
  C_word k = C_block_item(closure, 1); /* Original closure */
  C_word kab[4], *ka = kab, k2; /* Next continuation */
  k2 = C_closure(&ka, 3, (C_word)C_bignum_wrapped_return_bigobj, k, wrapped_b2);
  C_return_bignum(k2, (bignum_type)C_block_item(C_block_item(closure, 2), 0));
}

static void
C_return_2_bignums(C_word k, bignum_type b1, bignum_type b2)
{
  C_word bab[2], *ba = bab, kab[4], *ka = kab, k2, b1_ptr;
  /* Make b1 a wrapped pointer because C_closure wants scheme objects */
  b1_ptr = C_mpointer(&ba, b1);
  /* Wrap b2 first, then b1. Return them to k in the original (b1,b2) order */
  k2 = C_closure(&ka, 3, (C_word)C_b2_wrapped, k, b1_ptr);
  C_return_bignum(k2, b2);
}


static bignum_type
bignum_allocate(bignum_length_type length, int negative_p)
{
  bignum_type result;
  bignum_digit_type *digits;

  digits = (bignum_digit_type *)C_malloc(sizeof(bignum_digit_type)*(length + 1));
  
  if(digits == NULL) {
    fprintf(stderr, "out of memory - can not allocate number");
    exit(EXIT_FAILURE);
  }

  result = (bignum_type)digits;
  BIGNUM_SET_HEADER(result, length, negative_p);
  return result;
}

static bignum_type
bignum_allocate_zeroed(bignum_length_type length, int negative_p)
{
  assert((length >= 0) || (length < BIGNUM_RADIX));
  {
    bignum_type result;
    bignum_digit_type *digits;
    bignum_digit_type *scan;
    bignum_digit_type *end;
    digits=(bignum_digit_type *)C_malloc(sizeof(bignum_digit_type)*(length+1));
  
    if(digits == NULL) {
      fprintf(stderr, "out of memory - can not allocate number");
      exit(EXIT_FAILURE);
    }

    result = (bignum_type)digits;
    BIGNUM_SET_HEADER(result, length, negative_p);
    scan = BIGNUM_START_PTR(result);
    end = scan + length;
    while (scan < end)
      (*scan++) = 0;
    return (result);
  }
}

static bignum_type
bignum_allocate_from_fixnum(C_word fix)
{
  bignum_type ret;

  if (fix == C_fix(0)) {
    ret = bignum_allocate(0, 0);
  } else if (fix == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
    ret = bignum_allocate(2, 1);
    BIGNUM_REF(ret, 0) = 0;
    BIGNUM_REF(ret, 1) = 1;
  } else {
    bignum_digit_type digit = C_unfix(fix);
    ret = bignum_allocate(1, digit < 0);
    BIGNUM_REF(ret, 0) = ((digit < 0) ? -digit : digit);
  }
  return ret;
}

static bignum_type
bignum_digit_to_bignum(bignum_digit_type digit, int neg_p)
{
  bignum_type ret;
  if (digit == 0)
    return bignum_allocate(0, 0);

  ret = bignum_allocate(1, neg_p);
  BIGNUM_REF(ret, 0) = digit;
  return ret;
}

static bignum_type
shorten_bignum(bignum_type big, bignum_length_type newlength)
{
  bignum_digit_type *digits, *newdigits;
  digits = BIGNUM_TO_POINTER(big);
  newdigits = (bignum_digit_type *)C_realloc(digits, sizeof(bignum_digit_type)*(newlength + 1));
  if (newdigits == NULL) {
    fprintf(stderr, "out of memory - can not reallocate number");
    exit(EXIT_FAILURE);
  }
  return (bignum_type)newdigits;
}

static bignum_type
bignum_trim(bignum_type bignum)
{
  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
  bignum_digit_type * scan = end;
  while ((start <= scan) && ((*--scan) == 0))
    ;
  scan += 1;
  if (scan < end)
    {
      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
bignum_copy(bignum_type source)
{
  bignum_type target =
    (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source))));
  bignum_destructive_copy (source, target);
  return (target);
}

static void
bignum_destructive_copy(bignum_type source, bignum_type 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;
}

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

#define fix_randomize  C_randomize
#define fix_random C_random_fixnum

/*
 * This random number generator is very simple. Probably too simple...
 */
static C_word
big_randomize(C_word bignum)
{
  bignum_type big;
  C_word seed = 0;
  bignum_digit_type * scan, * end;

  big = big_of(bignum);
  scan = (BIGNUM_START_PTR (big));
  end = scan + BIGNUM_LENGTH(big);
  /* What a cheap way to initialize the random generator. I feel dirty! */
  while (scan < end)
    seed ^= *scan++;

  srand(C_unfix(seed));
  return C_SCHEME_UNDEFINED;
}

static void
big_random(C_word c, C_word self, C_word k, C_word max)
{
  bignum_type bigmax;
  bignum_type result;
  bignum_length_type randlen, max_len, max_bits;
  bignum_digit_type max_top_digit, d;
  bignum_digit_type * scan, * end;
  
  bigmax = big_of(max);

  max_len = BIGNUM_LENGTH(bigmax);
  max_top_digit = d = BIGNUM_REF(bigmax, max_len - 1);
  
  max_bits = (max_len - 1) * BIGNUM_DIGIT_LENGTH;
  while(d) {
    max_bits++;
    d >>= 1;
  }
  /* Subtract/add one because we don't want zero to be over-represented */
  randlen = (((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1));
  randlen = BIGNUM_BITS_TO_DIGITS(randlen + 1);
  
  result = bignum_allocate(randlen, BIGNUM_NEGATIVE_P(bigmax));
  scan = BIGNUM_START_PTR(result);
  end = scan + randlen - 1; /* Go to just before the end */
  while(scan < end)
    *scan++ = (((double)rand())/(RAND_MAX + 1.0) * (double)BIGNUM_RADIX);
  /*
   * Last word is special when length is max_len: It must be less than
   * max's most significant digit, instead of BIGNUM_RADIX.
   */
  if (max_len == randlen)
    *scan = (((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit);
  else
    *scan = (((double)rand())/(RAND_MAX + 1.0) * (double)BIGNUM_RADIX);

  result = bignum_trim(result);  
  C_return_bignum(k, result);
}

static void
fix_plus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx, result;
  bigx = bignum_allocate_from_fixnum(x);
  result = bignum_add(bigx, big_of(y));
  BIGNUM_DEALLOCATE(bigx);
  C_return_bignum(k, result);
}


static void
fix_minus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type tmp, result;
  tmp = bignum_allocate_from_fixnum(x);

  result = bignum_subtract(tmp, big_of(y));
  BIGNUM_DEALLOCATE(tmp);
  C_return_bignum(k, result);
}


static void
big_minus_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type tmp, result;
  tmp = bignum_allocate_from_fixnum(y);

  result = bignum_subtract(big_of(x), tmp);
  BIGNUM_DEALLOCATE(tmp);
  C_return_bignum(k, result);
}

static void
big_divrem_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type numerator = big_of(x);
  bignum_type denominator = big_of(y);
  int q_neg_p = ((BIGNUM_NEGATIVE_P (denominator))
                 ? (! (BIGNUM_NEGATIVE_P (numerator)))
                 : (BIGNUM_NEGATIVE_P (numerator)));

  switch (bignum_compare_unsigned (numerator, denominator))
    {
    case bignum_comparison_equal:
      C_values(4, C_SCHEME_UNDEFINED, k,
               q_neg_p ? C_fix(-1) : C_fix(1), C_fix(0));
    case bignum_comparison_less:
      C_values(4, C_SCHEME_UNDEFINED, k, C_fix(0), x);
    case bignum_comparison_greater:
    default:                                  /* to appease gcc -Wall */
      {
        bignum_type quotient, remainder;
        int r_neg_p = BIGNUM_NEGATIVE_P (numerator) ? 1 : 0;
        
        assert(BIGNUM_LENGTH(denominator) > 1);
        bignum_divide_unsigned_large_denominator
          (numerator, denominator, (&quotient), (&remainder), q_neg_p, r_neg_p);
        C_return_2_bignums(k, quotient, remainder);
      }
    }
}

static void
big_divrem_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx = big_of(x);
  y = C_unfix(y);

  if (y == 1)
    C_values(4, C_SCHEME_UNDEFINED, k, x, C_fix(0));
  else if (y == -1)
    C_return_big_fix(k, bignum_new_sign(bigx, !(BIGNUM_NEGATIVE_P(bigx))), C_fix(0));

  /* Too bad, we really need to do some work... */
  {
    int q_neg_p = (y < 0) ? !(BIGNUM_NEGATIVE_P(bigx)) : BIGNUM_NEGATIVE_P(bigx);
    int r_neg_p = BIGNUM_NEGATIVE_P(bigx);
    bignum_digit_type abs_y = (y < 0) ? -y : y;
    bignum_type quotient, remainder;
    
    if (y == C_MOST_NEGATIVE_FIXNUM) {
      if (!BIGNUM_NEGATIVE_P(bigx) && BIGNUM_LENGTH(bigx) == 1
          && BIGNUM_REF(bigx, 1) == 1 && BIGNUM_REF(bigx, 0) == 0) {
        /*
         * Very very special case:
         * quotient(MOST_NEGATIVE_FIXNUM, -(MOST_NEGATIVE_FIXNUM)) => -1
         */
        C_values(4, C_SCHEME_UNDEFINED, k, C_fix(-1), C_fix(0));
      } else {
        /* This is the only case we need to go allocate a bignum for */
        bignum_type bigy =
          bignum_allocate_from_fixnum(C_fix(C_MOST_NEGATIVE_FIXNUM));

        bignum_divide_unsigned_large_denominator
          (bigx, bigy, (&quotient), (&remainder), q_neg_p, r_neg_p);
        BIGNUM_DEALLOCATE(bigy);
        C_return_2_bignums(k, quotient, remainder);
      }
    } else if (abs_y < BIGNUM_RADIX_ROOT) {
      bignum_divide_unsigned_small_denominator
        (bigx, abs_y, (&quotient), (&remainder), q_neg_p, r_neg_p);
      C_return_2_bignums(k, quotient, remainder);
    } else {
      bignum_divide_unsigned_medium_denominator
        (bigx, abs_y, (&quotient), (&remainder), q_neg_p, r_neg_p);
      C_return_2_bignums(k, quotient, remainder);
    }
  }
}

/*
 * Big_gcd is a huge function and it sucks that it needs to be in C.
 *
 * Why does it need to be in C?  Because if you have a very big bignum
 * that doesn't cleanly divide another big bignum, you end up calling
 * the remainder procedure a lot in Scheme.  This produces tons of
 * intermediate bignums, which means a lot of copies into GC'able memory
 * need to be made (and the GC will be triggered more often).
 * That's a major slowdown.  Doing the loop in C means the intermediate
 * results can be cleaned up right away each loop step, and returning
 * just one result to Scheme.
 * Once (if?) we find a way to avoid copying bignums (instead allocating
 * directly in GCable memory) this function can be cut out and replaced by
 * a recursive call to gcd-0 in Scheme.
 */
static void
big_gcd(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx = big_of(x), bigy = big_of(y), bigr;

  assert(BIGNUM_LENGTH(bigx) > 1);
  assert(BIGNUM_LENGTH(bigy) > 1);
  
  switch(bignum_compare_unsigned (bigx, bigy)) {
  case bignum_comparison_equal:
    C_kontinue(k, x);
  case bignum_comparison_less:
    /* Swap, since remainder of bigx, bigy would be bigx, causing an extra loop */
    bigr = bigy;
    bigy = bigx;
    bigx = bigr;

    /* Make x and y match for the special case where gcd(x, y) = y */
    {
      C_word tmp = y;
      y = x;
      x = tmp;
    }
    /* FALLTHROUGH */
  default: /* Continue below */
    break;
  }
  
  /*
   * Be careful! Don't deallocate live objects. We could start with a copy
   * or compare pointers with big_of(x) or y every time but that seems wasteful.
   */
  bignum_divide_unsigned_large_denominator
   (bigx, bigy, ((bignum_type *) 0), (&bigr), 0, 0);
  bigx = bigy;
  bigy = bigr;
  /* Original bigx is forgotten now */
  assert(bigx != big_of(x));
  assert(bigy != big_of(x));
  /* Only bigx points to y */
  assert(bigy != big_of(y));
  assert(bigx == big_of(y));

  switch (BIGNUM_LENGTH(bigy)) {
  case 0:  /* bigy = 0 */
    /* remainder(x, y) = 0 => y  |  x < y */
    BIGNUM_DEALLOCATE(bigy); /* Got allocated in previous step */
    C_kontinue(k, y);
  case 1:
    if (BIGNUM_REF(bigy, 0) == 1) { /* y = 1? Then 1 is the result */
      BIGNUM_DEALLOCATE(bigy); /* Got allocated in previous step */
      C_kontinue(k, C_fix(1));
    } else if (BIGNUM_REF(bigy, 0) < BIGNUM_RADIX_ROOT)
      bigr = bignum_remainder_unsigned_small_denominator
              (bigx, BIGNUM_REF(bigy, 0), 0);
    else
      bignum_divide_unsigned_medium_denominator
       (bigx, BIGNUM_REF(bigy, 0), (bignum_type *)0, (&bigr), 0, 0);
    break;
  default:
    bignum_divide_unsigned_large_denominator
     (bigx, bigy, ((bignum_type *) 0), (&bigr), 0, 0);
  }
  /* Swap, but don't deallocate x since it holds the original value of y */
  bigx = bigy;
  bigy = bigr;
  
  /* Original bigy is forgotten now, we can safely always deallocate bigx */

  /* Assume that bignums coming from outside are never length 1 */
  assert(bigx != big_of(y));
  
  while(BIGNUM_LENGTH(bigy) > 1) {
    bignum_divide_unsigned_large_denominator
     (bigx, bigy, ((bignum_type *) 0), (&bigr), 0, 0);
    BIGNUM_DEALLOCATE(bigx);
    bigx = bigy;
    bigy = bigr;
  }

  /* Finish up with a faster loop until y = 0 (ie, length(bigy) = 0) */
  while (BIGNUM_LENGTH(bigy) == 1) {
    if (BIGNUM_REF(bigy, 0) == 1) {
      BIGNUM_DEALLOCATE(bigx);
      BIGNUM_DEALLOCATE(bigy);
      C_kontinue(k, C_fix(1));
      break;
    } else if (BIGNUM_REF(bigy, 0) < BIGNUM_RADIX_ROOT) {
      bigr = bignum_remainder_unsigned_small_denominator
              (bigx, BIGNUM_REF(bigy, 0), 0);
    } else {
      bignum_divide_unsigned_medium_denominator
       (bigx, BIGNUM_REF(bigy, 0), (bignum_type *)0, (&bigr), 0, 0);
    }
    BIGNUM_DEALLOCATE(bigx);
    bigx = bigy;
    bigy = bigr;
  }
  BIGNUM_DEALLOCATE(bigy);
  C_return_bignum(k, bigx);
}

static C_word
C_a_i_flonum_gcd(C_word **p, C_word n, C_word x, C_word y)
{
   double xub, yub, r;
   
   xub = C_flonum_magnitude(x);
   yub = C_flonum_magnitude(y);

   if (xub < 0.0) xub = -xub;
   if (yub < 0.0) yub = -yub;
   
   while(yub != 0.0) {
     r = fmod(xub, yub);
     xub = yub;
     yub = r;
   }
   C_return(C_flonum(p, xub));
}


static void
big_plus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type result;
  result = bignum_add(big_of(x), big_of(y));
  C_return_bignum(k, result);
}

/*
 * This now makes the assumption it is never passed a bignum of LENGTH 0.
 * This should always be valid in Chicken.
 */
static bignum_type
bignum_add(bignum_type x, bignum_type y)
{
  return
    (((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)))));
}

static bignum_type
bignum_add_unsigned(bignum_type x, bignum_type y, 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));
    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;
    } else { /* r is one word too big (to hold a possible carry), readjust */
      BIGNUM_SET_HEADER(r, x_length, negative_p);
      BIGNUM_REDUCE_LENGTH(r, r, x_length);
    }
    return (r);
  }
}

static void
big_minus_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type result;
  result = bignum_subtract(big_of(x), big_of(y));
  C_return_bignum(k, result);
}

/*
 * This now makes the assumption it is never passed a bignum of LENGTH 0.
 * This should always be valid in Chicken.
 */
static bignum_type
bignum_subtract(bignum_type x, bignum_type y)
{
  return
    (((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)))));
}

static bignum_type
bignum_subtract_unsigned(bignum_type x, bignum_type y)
{
  int negative_p;
  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));
    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;
              }
          }
      assert(borrow == 0);
      while (scan_x < end_x)
        (*scan_r++) = (*scan_x++);
    }
    return (bignum_trim (r));
  }
}


static void
big_times_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx = big_of(x), bigy = big_of(y);
  int neg_p = ((BIGNUM_NEGATIVE_P (bigx))
               ? (! (BIGNUM_NEGATIVE_P (bigy)))
               : (BIGNUM_NEGATIVE_P (bigy)));
  /* If length 1 or 0, it should be a fixnum */
  assert(BIGNUM_LENGTH(bigx) > 1);
  assert(BIGNUM_LENGTH(bigy) > 1);
  C_return_bignum(k, bignum_multiply_unsigned(bigx, bigy, neg_p));
}

static void
fix_times_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx, bigy, result;
  bignum_digit_type absx, absy;
  C_word neg_p;

  absx = C_unfix(x);
  absx = absx < 0 ? -absx : absx;
  absy = C_unfix(y);
  absy = absy < 0 ? -absy : absy;
  neg_p = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT));

  if (absx < BIGNUM_RADIX_ROOT) {
     if (absx == 0 || absx == 1 || absy < BIGNUM_RADIX_ROOT) {
       C_kontinue(k, C_fix(neg_p ? -(absx * absy) : (absx * absy)));
     } else {
       bigy = bignum_allocate_from_fixnum(y);
       result = bignum_multiply_unsigned_small_factor(bigy, absx, neg_p ? 1 : 0);
       BIGNUM_DEALLOCATE(bigy);
       C_return_bignum(k, result);
     }
  } else if (absy < BIGNUM_RADIX_ROOT) {
     if (absy == 0 || absy == 1 /*|| absx < BIGNUM_RADIX_ROOT */) {
       C_kontinue(k, C_fix(neg_p ? -(absx * absy) : (absx * absy)));
     } else {
       bigx = bignum_allocate_from_fixnum(x);
       result = bignum_multiply_unsigned_small_factor(bigx, absy, neg_p ? 1 : 0);
       BIGNUM_DEALLOCATE(bigx);
       C_return_bignum(k, result);
     }
  } else {
    bigx = bignum_allocate_from_fixnum(x);
    bigy = bignum_allocate_from_fixnum(y);
    result = bignum_multiply_unsigned (bigx, bigy, neg_p ? 1 : 0);
    BIGNUM_DEALLOCATE(bigy);
    BIGNUM_DEALLOCATE(bigx);
    C_return_bignum(k, result);
  }
}


static void
fix_times_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx, result, bigy;
  bignum_digit_type absx;
  C_word neg_p;
  
  bigy = big_of(y);
  
  if (x == C_fix(0))
    C_kontinue(k, C_fix(0));
  else if (x == C_fix(1))
    C_kontinue(k, y);
  else if (x == C_fix(-1))
    C_return_bignum(k, bignum_new_sign(bigy, !(BIGNUM_NEGATIVE_P(bigy))));

  absx = C_unfix(x);
  absx = absx < 0 ? -absx : absx;  
  neg_p = ((x & C_INT_SIGN_BIT)
           ? !(BIGNUM_NEGATIVE_P(bigy)) : (BIGNUM_NEGATIVE_P(bigy)));
  
  if (absx < BIGNUM_RADIX_ROOT) {
     C_return_bignum(k, bignum_multiply_unsigned_small_factor(bigy, absx, neg_p));
  } else {
    bigx = bignum_allocate_from_fixnum(x);
    result = bignum_multiply_unsigned (bigx, bigy, neg_p);
    BIGNUM_DEALLOCATE(bigx);
    C_return_bignum(k, result);
  }
}

/* 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
bignum_multiply_unsigned(bignum_type x, bignum_type y, int negative_p)
{
  if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
    {
      bignum_type z = x;
      x = y;
      y = z;
    }
  {
    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_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
bignum_multiply_unsigned_small_factor(bignum_type x, bignum_digit_type y,
                                      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
bignum_destructive_scale_up(bignum_type bignum, bignum_digit_type factor)
{
  bignum_digit_type carry = 0;
  bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  bignum_digit_type two_digits;
  bignum_digit_type product_low;
#define product_high carry
  bignum_digit_type * end = (scan + (BIGNUM_LENGTH (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.
   */
  assert(carry == 0);
  return;
#undef product_high
}

static void
bignum_destructive_add(bignum_type bignum, bignum_digit_type n)
{
  bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  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
bignum_divide_unsigned_large_denominator(bignum_type numerator,
                                         bignum_type denominator,
                                         bignum_type * quotient,
                                         bignum_type * remainder,
                                         int q_negative_p,
                                         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;
  assert(length_d > 1);
  {
    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);
      if (remainder != ((bignum_type *) 0))
        (*remainder) = (bignum_trim (u));
      else
        BIGNUM_DEALLOCATE (u);
    }
  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))
        (*remainder) = bignum_destructive_unnormalization (u, shift);
      else
        BIGNUM_DEALLOCATE(u);
    }
  if (quotient != ((bignum_type *) 0))
    (*quotient) = (bignum_trim (q));
  return;
}

static void
bignum_divide_unsigned_normalized(bignum_type u, bignum_type v, 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;
  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
}

static bignum_digit_type
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);
}

static void
bignum_divide_unsigned_medium_denominator(bignum_type numerator,
                                          bignum_digit_type denominator,
                                          bignum_type * quotient,
                                          bignum_type * remainder,
                                          int q_negative_p,
                                          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);
    }
  {
    bignum_digit_type r = 0;
    bignum_digit_type * start = (BIGNUM_START_PTR (q));
    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
bignum_destructive_normalization(bignum_type source, bignum_type 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 = ((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
    assert(carry == 0);
  return;
}

/* This will also trim the number if necessary */
static bignum_type
bignum_destructive_unnormalization(bignum_type bignum, int shift_right)
{
  bignum_length_type length = BIGNUM_LENGTH(bignum);
  bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  bignum_digit_type * scan = (start + length);
  bignum_digit_type digit;
  bignum_digit_type carry = 0;
  int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
  bignum_digit_type mask = ((1L << shift_right) - 1);
  
  while (!(*--scan)) {
    if (start == scan) { /* Don't bother. */
      BIGNUM_SET_HEADER (bignum, 0, 0);
      BIGNUM_REDUCE_LENGTH (bignum, bignum, 0);
      return bignum;
    }
    --length;
  }

  digit = (*scan);
  (*scan) = (digit >> shift_right);
  length -= (*scan == 0); /* Add 1 or 0 */
  carry = ((digit & mask) << shift_left);
  
  while (start < scan)
    {
      digit = (*--scan);
      (*scan) = ((digit >> shift_right) | carry);
      carry = ((digit & mask) << shift_left);
    }
  assert(carry == 0);
  assert(BIGNUM_LENGTH(bignum) != length);
  assert(length != 1 || BIGNUM_REF(bignum, 0) != 0);
  BIGNUM_SET_HEADER
    (bignum, length, (BIGNUM_NEGATIVE_P (bignum)));
  BIGNUM_REDUCE_LENGTH (bignum, bignum, length);

  return bignum;
}

/* 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])));         \
}

static bignum_digit_type
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;
  bignum_digit_type v2;
  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;                                                        \
    }                                                                   \
}

static bignum_digit_type
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

static void
bignum_divide_unsigned_small_denominator(bignum_type numerator,
                                         bignum_digit_type denominator,
                                         bignum_type * quotient,
                                         bignum_type * remainder,
                                         int q_negative_p,  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
bignum_destructive_scale_down(bignum_type bignum, 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 (bignum));
  bignum_digit_type * scan = (start + (BIGNUM_LENGTH (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
bignum_remainder_unsigned_small_denominator(bignum_type 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;
  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 C_word
big_comp_big(C_word x, C_word y)
{
  bignum_type bigx = big_of(x), bigy = big_of(y);
  return C_fix((BIGNUM_NEGATIVE_P (bigx))
               ? ((BIGNUM_NEGATIVE_P (bigy))
                  ? (bignum_compare_unsigned (bigy, bigx))
                  : (bignum_comparison_less))
               : ((BIGNUM_NEGATIVE_P (bigy))
                  ? (bignum_comparison_greater)
                  : (bignum_compare_unsigned (bigx, bigy))));
}

static enum bignum_comparison
bignum_compare_unsigned(bignum_type x, bignum_type y)
{
  bignum_length_type x_length;
  bignum_length_type y_length;
  if (x == y) /* Objects are the same? */
    return (bignum_comparison_equal);

  x_length = (BIGNUM_LENGTH (x));
  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);
}

static void
big_abs(C_word c, C_word self, C_word k, C_word big)
{
  if (!BIGNUM_NEGATIVE_P(big_of(big)))
    C_kontinue(k, big);
  else
    C_return_bignum(k, bignum_new_sign(big_of(big), 0));
}

static void
big_quotient_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx = big_of(x);
  y = C_unfix(y);

  if (y == 1)
    C_kontinue(k, x);
  else if (y == -1)
    C_return_bignum(k, bignum_new_sign(bigx, !(BIGNUM_NEGATIVE_P(bigx))));

  /* Too bad, we really need to do some work... */
  {
    int neg_p = (y < 0) ? !(BIGNUM_NEGATIVE_P(bigx)) : BIGNUM_NEGATIVE_P(bigx);
    bignum_digit_type abs_y = (y < 0) ? -y : y;
    bignum_type quotient;
    
    if (y == C_MOST_NEGATIVE_FIXNUM) {
      if (!BIGNUM_NEGATIVE_P(bigx) && BIGNUM_LENGTH(bigx) == 1
          && BIGNUM_REF(bigx, 1) == 1 && BIGNUM_REF(bigx, 0) == 0) {
        /*
         * Very very special case:
         * quotient(MOST_NEGATIVE_FIXNUM, -(MOST_NEGATIVE_FIXNUM)) => -1
         */
        C_kontinue(k, C_fix(-1));
      } else {
        /* This is the only case we need to go allocate a bignum for */
        bignum_type bigy =
          bignum_allocate_from_fixnum(C_fix(C_MOST_NEGATIVE_FIXNUM));

        bignum_divide_unsigned_large_denominator
          (bigx, bigy, (&quotient), ((bignum_type *) 0), neg_p, 0);
        BIGNUM_DEALLOCATE(bigy);
        C_return_bignum(k, quotient);
      }
    } else if (abs_y < BIGNUM_RADIX_ROOT) {
      bignum_divide_unsigned_small_denominator
        (bigx, abs_y, (&quotient), ((bignum_type *) 0), neg_p, 0);
      C_return_bignum(k, quotient);
    } else {
      bignum_divide_unsigned_medium_denominator
        (bigx, abs_y, (&quotient), ((bignum_type *) 0), neg_p, 0);
      C_return_bignum(k, quotient);
    }
  }
}

static void
big_quotient_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type numerator = big_of(x);
  bignum_type denominator = big_of(y);
  int neg_p = ((BIGNUM_NEGATIVE_P (denominator))
               ? (! (BIGNUM_NEGATIVE_P (numerator)))
               : (BIGNUM_NEGATIVE_P (numerator)));

  switch (bignum_compare_unsigned (numerator, denominator))
    {
    case bignum_comparison_equal:
      C_kontinue(k, neg_p ? C_fix(-1) : C_fix(1));
    case bignum_comparison_less:
      C_kontinue(k, C_fix(0));
    case bignum_comparison_greater:
    default:                                  /* to appease gcc -Wall */
      {
        bignum_type quotient;
        assert(BIGNUM_LENGTH(denominator) > 1);
        bignum_divide_unsigned_large_denominator
          (numerator, denominator, (&quotient), ((bignum_type *) 0), neg_p, 0);
        C_return_bignum(k, quotient);
      }
    }
}

static void
big_remainder_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
{
   bignum_type remainder, bigy, bigx;
   int x_neg, y_neg;
   C_word abs_y;

   switch (y) {
   /* case 0:  SHOULD NOT HAPPEN (checked in Scheme)
     C_kontinue(k, C_SCHEME_FALSE); */
   case C_fix(1):
   case C_fix(-1):
     C_kontinue(k, C_fix(0));
   case C_fix(C_MOST_NEGATIVE_FIXNUM):
     bigx = big_of(x);
     if (BIGNUM_LENGTH(bigx) == 2 && BIGNUM_REF(bigx, 0) == 0 && BIGNUM_REF(bigx, 0))
       C_kontinue(k, C_fix(0));
     /* Don't handle 0 <= length(bigx) <= 1 since then it should be a fixnum */
     assert(BIGNUM_LENGTH(bigx) >= 2);
     
     bigy = bignum_allocate_from_fixnum(y);
     bignum_divide_unsigned_large_denominator
      (bigx, bigy, (bignum_type *)0, &remainder, 0, BIGNUM_NEGATIVE_P(bigx));
     BIGNUM_DEALLOCATE(bigy);
     
     C_return_bignum(k, remainder);
   default:
     bigx = big_of(x);
     y = C_unfix(y);
     y_neg = (y < 0);
     abs_y = y_neg ? -y : y;
     x_neg = BIGNUM_NEGATIVE_P(bigx);
     
     if (abs_y < BIGNUM_RADIX_ROOT)
       remainder =
         bignum_remainder_unsigned_small_denominator(bigx, abs_y, x_neg);
     else
       bignum_divide_unsigned_medium_denominator(bigx, abs_y,
                                                 (bignum_type *) 0, &remainder,
                                                 x_neg, x_neg);
     C_return_bignum(k, remainder);
   }
}

static void
big_remainder_big(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type numerator = big_of(x), denominator = big_of(y);
  
  switch (bignum_compare_unsigned (numerator, denominator))
    {
      case bignum_comparison_equal:
        C_kontinue(k, C_fix(0));
      case bignum_comparison_less:
        C_kontinue(k, x);
      case bignum_comparison_greater:
      default:                                  /* to appease gcc -Wall */
        {
          bignum_type remainder;
          bignum_divide_unsigned_large_denominator
            (numerator, denominator,
             ((bignum_type *) 0), (&remainder),
             0, BIGNUM_NEGATIVE_P(numerator));
          C_return_bignum(k, remainder);
        }
    }
}

static void
flo_to_int(C_word c, C_word self, C_word k, C_word x)
{
  C_return_bignum(k, double_to_bignum(C_flonum_magnitude(x)));
}

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

static bignum_type
double_to_bignum(double x)
{
  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_type result = (bignum_allocate (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 (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

static void 
int_bitwise_int(C_word c, C_word self, C_word k, C_word op, C_word x, C_word y)
{
  bignum_type bigx, bigy, result;
  
  if((x & C_FIXNUM_BIT) != 0)
    bigx = bignum_allocate_from_fixnum(x);
  else
    bigx = big_of(x);

  if((y & C_FIXNUM_BIT) != 0)
    bigy = bignum_allocate_from_fixnum(y);
  else
    bigy = big_of(y);

  result = ((BIGNUM_NEGATIVE_P (bigx))
             ? (BIGNUM_NEGATIVE_P (bigy))
               ? bignum_negneg_bitwise_op(C_unfix(op), bigx, bigy)
               : bignum_posneg_bitwise_op(C_unfix(op), bigy, bigx)
             : (BIGNUM_NEGATIVE_P (bigy))
               ? bignum_posneg_bitwise_op(C_unfix(op), bigx, bigy)
               : bignum_pospos_bitwise_op(C_unfix(op), bigx, bigy));
        
  if((x & C_FIXNUM_BIT) != 0)
    BIGNUM_DEALLOCATE(bigx);
  if((y & C_FIXNUM_BIT) != 0)
    BIGNUM_DEALLOCATE(bigy);

  C_return_bignum(k, result);
}

static void 
int_not(C_word c, C_word self, C_word k, C_word x)
{
  bignum_type bigx, result;
  
  if((x & C_FIXNUM_BIT) != 0)
    bigx = bignum_allocate_from_fixnum(x);
  else
    bigx = big_of(x);

  result = bignum_bitwise_not(bigx);

  if((x & C_FIXNUM_BIT) != 0)
    BIGNUM_DEALLOCATE(bigx);

  C_return_bignum(k, result);
}

bignum_type
bignum_bitwise_not(bignum_type x)
{
  static C_word invbits[] = { BIGNUM_RADIX | 1, 1 }; /* bignum representing -1 */
  return bignum_subtract((bignum_type)invbits, x);
}

/*
 * From Hacker's Delight by Frank Warren
 * based on a modified nlz() from section 5-3 (fig. 5-7)
 */
static C_word
ilen(C_uword x) {
  C_uword y;
  C_word n = 0;

#ifdef C_SIXTY_FOUR
  y = x >> 32; if (y != 0) { n += 32; x = y; }
#endif
  y = x >> 16; if (y != 0) { n += 16; x = y; }
  y = x >>  8; if (y != 0) { n +=  8; x = y; }
  y = x >>  4; if (y != 0) { n +=  4; x = y; }
  y = x >>  2; if (y != 0) { n +=  2; x = y; }
  y = x >>  1; if (y != 0) return n + 2;
  return n + x;
}

static void
int_length(C_word c, C_word self, C_word k, C_word x)
{
  if (x & C_FIXNUM_BIT) {
    x = C_unfix(x);
    C_kontinue(k, C_fix(ilen((x < 0) ? ~x : x)));
  } else { /* bignum */
    bignum_type bigx;
    bignum_digit_type bd;
    int len_1;

    bigx = big_of(x);
    if (BIGNUM_NEGATIVE_P(bigx)) {
      /* TODO: This can probably be done faster, by juggling with borrow */
      bignum_type notx = bignum_bitwise_not(bigx);
      len_1 = BIGNUM_LENGTH(notx) - 1;
      bd = *(BIGNUM_START_PTR(notx) + len_1); /* Most significant digit */
      BIGNUM_DEALLOCATE(notx);
    } else {
      len_1 = BIGNUM_LENGTH(bigx) - 1;
      bd = *(BIGNUM_START_PTR(bigx) + len_1); /* Most significant digit */
    }
    C_kontinue(k, C_fix(ilen(bd) + len_1 * BIGNUM_DIGIT_LENGTH));
  }
}

static void 
int_shift_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  bignum_type bigx, result;

  if (y == C_fix(0)) C_kontinue(k, x); /* Done too (no shift) */
  
  /* Ensure x is a bignum */
  if((x & C_FIXNUM_BIT) != 0) {
    if (x == C_fix(0)) /* Skip everything else */
      C_kontinue(k, x);
    
    bigx = bignum_allocate_from_fixnum(x);
  } else {
    bigx = big_of(x);
  }
   
  result = bignum_arithmetic_shift(bigx, C_unfix(y));
  if ((x & C_FIXNUM_BIT) != 0)
    BIGNUM_DEALLOCATE(bigx);
  C_return_bignum(k, result);
}

static bignum_type
bignum_arithmetic_shift(bignum_type arg1, C_word n)
{
  bignum_type tmp1, tmp2, result;
  if (BIGNUM_NEGATIVE_P(arg1) && n < 0) {
    tmp1 = bignum_bitwise_not(arg1);
    tmp2 = bignum_magnitude_ash(tmp1, n);
    BIGNUM_DEALLOCATE(tmp1);
    result = bignum_bitwise_not(tmp2);
    BIGNUM_DEALLOCATE(tmp2);
    return result;
  } else {
    return bignum_magnitude_ash(arg1, n);
  }
}

static bignum_type
bignum_magnitude_ash(bignum_type arg1, C_word n)
{
  bignum_type result;
  bignum_digit_type *scan1;
  bignum_digit_type *scanr;
  bignum_digit_type *end;

  C_word digit_offset,bit_offset;
  assert(n != 0); /* int_shift_fix checks this */

  if (n > 0) {
    digit_offset = n / BIGNUM_DIGIT_LENGTH;
    bit_offset =   n % BIGNUM_DIGIT_LENGTH;
    
    result = bignum_allocate_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 = bignum_allocate_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 ;
  }
  
  return (bignum_trim (result));
}

static bignum_type
bignum_pospos_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
{
  bignum_type result;
  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);

  result = bignum_allocate(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;
    /*
    fprintf(stderr, "[pospos op = %d, i = %ld, d1 = %lx, d2 = %lx]\n",
            op, endr - scanr, digit1, digit2);
            */
    *scanr++ = (op == bignum_and_op) ? digit1 & digit2 :
               (op == bignum_ior_op) ? digit1 | digit2 :
                                       digit1 ^ digit2;
  }
  return bignum_trim(result);
}

static bignum_type
bignum_posneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
{
  bignum_type result;
  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 == bignum_ior_op || op == bignum_xor_op;

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

  result = bignum_allocate(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 == bignum_and_op) ? digit1 & digit2 :
               (op == bignum_ior_op) ? digit1 | digit2 :
                                       digit1 ^ digit2;
  }
  
  if (neg_p)
    bignum_negate_magnitude(result);

  return bignum_trim(result);
}

static bignum_type
bignum_negneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
{
  bignum_type result;
  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 == bignum_and_op || op == bignum_ior_op;

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

  result = bignum_allocate(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 == bignum_and_op) ? digit1 & digit2 :
               (op == bignum_ior_op) ? digit1 | digit2 :
                                       digit1 ^ digit2;
  }

  if (neg_p)
    bignum_negate_magnitude(result);

  return bignum_trim(result);
}

static void
bignum_negate_magnitude(bignum_type 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;
  }
}

#define BIGNUM_STR_BLOCK_SIZE 8

/*
 * Contains s48_bignum_to_digit_stream.  We may want to separate that out
 * again in case of huge numbers; it may be more efficient to write
 * those straight to files without going through the string conversion.
 */
static void
big_to_string(C_word c, C_word self, C_word k, C_word value, C_word radix)
{
  char *buf, *index;
  C_word len, counter=0;
  char *characters = "0123456789abcdef";
  C_word *tmp, ret;
  bignum_type bignum;

  bignum = big_of(value);
  len = BIGNUM_STR_BLOCK_SIZE;
  buf = C_malloc(len);
  if (buf == NULL) {
    fprintf(stderr, "out of memory - can not allocate string");
    exit(EXIT_FAILURE);
  }
  index = buf + len - 1;

  if (BIGNUM_NEGATIVE_P(bignum)) {
    *index-- = '-';
    counter++;
  }

  radix = C_unfix(radix);
  
  assert((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  if (! (BIGNUM_ZERO_P (bignum)))
    {
      bignum_type working_copy = (bignum_copy (bignum));
      bignum_digit_type * start = (BIGNUM_START_PTR (working_copy));
      bignum_digit_type * scan = (start + (BIGNUM_LENGTH (working_copy)));
      bignum_digit_type digit;
      while (start < scan)
        {
          if ((scan[-1]) == 0) {
            scan -= 1;
          } else {
            digit = bignum_destructive_scale_down (working_copy, radix);
            *index-- = characters[digit];
            if (++counter == len) {
              char *newbuf;
              newbuf = C_malloc(len + BIGNUM_STR_BLOCK_SIZE);
              if (newbuf == NULL) {
                fprintf(stderr, "out of memory - can not allocate string");
                exit(EXIT_FAILURE);
              }
              C_memcpy(newbuf + BIGNUM_STR_BLOCK_SIZE, buf, len);
              C_free(buf);
              buf = newbuf;
              index = newbuf + BIGNUM_STR_BLOCK_SIZE - 1;
              len += BIGNUM_STR_BLOCK_SIZE;
            }
          }
        }
      BIGNUM_DEALLOCATE (working_copy);
    }

  if (BIGNUM_NEGATIVE_P(bignum))
    *index-- = '-';

  tmp = C_alloc(C_SIZEOF_STRING(counter));
  ret = C_string(&tmp, counter, index+1);
  C_free(buf);
  C_kontinue(k, ret);
}

static void
digits_to_big(C_word c, C_word self, C_word k, C_word n,
              C_word start, C_word end, C_word radix, C_word negp)
{
  char *str;
  size_t n_digits;
  int negative_p = (negp != C_SCHEME_FALSE);
  int digit;
  int hash = 0;

  str = C_c_string(n) + C_unfix(start);
  n_digits = C_unfix(end)-C_unfix(start);
  radix = C_unfix(radix);
  hash = /* abs(radix / 2) */ 0;
  
#define DIGIT_TO_INT(x)         \
  (((x) == '#') ? hash :        \
   (((x) >= (int)'a') ?((x) - (int)'a' + 10) : ((x) - (int)'0')))

  assert((radix > 1) && (radix < BIGNUM_RADIX_ROOT));
  if (n_digits == 0)
    C_kontinue(k, C_SCHEME_FALSE);
  if (n_digits == 1)
    {
      digit = DIGIT_TO_INT(C_tolower((int)*str));
      if (digit >= radix || digit < 0)
        C_kontinue(k, C_SCHEME_FALSE);
      else
        C_return_bignum(k, bignum_digit_to_bignum(digit, negative_p));
    }
  {
    bignum_length_type length;
    {
      unsigned int radix_copy = radix;
      unsigned int log_radix = 0;
      while (radix_copy > 0)
        {
          radix_copy >>= 1;
          log_radix += 1;
        }
      /* This length will be at least as large as needed. */
      length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix));
    }
    {
      bignum_type result = (bignum_allocate_zeroed (length, negative_p));
      while ((n_digits--) > 0)
        {
          digit = DIGIT_TO_INT(C_tolower((int)*str));
          str++;
          if (digit >= radix || digit < 0) {
            BIGNUM_DEALLOCATE(result);
            C_kontinue(k, C_SCHEME_FALSE);
          }
        
          bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
          bignum_destructive_add (result, digit);
        }
      C_return_bignum (k, bignum_trim (result));
    }
  }
}

/**
 * Below you will find the functions that have been refactored to
 * match the "core" style.
 *
 * Naming convention idea: Maybe name these fixnum ops differently,
 * _p_ for "promoting"?  OTOH, it may be safer and cleaner to rename
 * the old fixnum ops to have _fx_ to indicate they run in fixnum mode.
 */
static void bignum_negate_2(C_word c, C_word self, C_word new_big);
static void allocate_bignum_2(C_word c, C_word self, C_word bigvec);

/* Eventually this will probably need to be integrated into C_2_plus. */
void C_ccall
C_u_2_fixnum_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
{
  C_word z, ab[C_SIZEOF_BIGNUM(2)], *a = ab;

  /* Exceptional situation: this will cause a real overflow */
  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
    C_kontinue(k, C_bignum2(&a, 1, 0, 2));
  } else {
    z = C_unfix(x) + C_unfix(y);

    /* This code "knows" that both fixnums and bignums have 2 reserved bits */
    if(!C_fitsinfixnump(z)) {
      C_kontinue(k, C_bignum2(&a, (z < 0), labs(z) & (C_uword)BIGNUM_DIGIT_MASK, 1));
    } else {
      C_kontinue(k, C_fix(z));
    }
  }
}

/* TODO: This should probably be renamed C_fixnum_negate to replace
 * what's in core.  Unfortunately, that one is allocating inline.
 * That one may be renamed to C_u_i_fixnum_negate() or some such.
 * TODO: Convert this to be an inline function and move to header?
 */
void C_ccall
C_u_fixnum_neg(C_word c, C_word self, C_word k, C_word x)
{
  /* Exceptional situation: this will cause an overflow to itself */
  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM)) { /* C_fitsinfixnump(x) */
    C_word ab[C_SIZEOF_BIGNUM(2)], *a = ab;
    C_kontinue(k, C_bignum2(&a, 0, 0, 1));
  } else {
    C_kontinue(k, C_fix(-C_unfix(x)));
  }
}

void C_ccall
C_fixnum_gcd(C_word c, C_word self, C_word k, C_word x, C_word y)
{
   C_word r;
   
   x = C_unfix(x);
   y = C_unfix(y);
   
   if (x < 0) x = -x;
   if (y < 0) y = -y;
   
   while(y != 0) {
     r = x % y;
     x = y;
     y = r;
   }
   C_kontinue(k, C_fix(x));
}

void C_ccall
C_u_bignum_negate(C_word c, C_word self, C_word k, C_word x)
{
  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, negp, size;

  negp = C_i_not(C_u_i_bignum_negativep(x));
  k2 = C_closure(&ka, 3, (C_word)bignum_negate_2, k, x);
  
  size = C_u_i_bignum_size(x);
  C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
}

static void
bignum_negate_2(C_word c, C_word self, C_word new_big)
{
  C_word k = C_block_item(self, 1),
         old_big = C_block_item(self, 2);
  C_memcpy(C_bignum_digits(new_big), C_bignum_digits(old_big),
           /* TODO: This is currently in bytes.  If we change the
            * representation that needs to change!
            * We subtract the size of the header, too.
            */
           C_header_size(C_internal_bignum(old_big))-C_wordstobytes(1));
  C_kontinue(k, new_big);
}

void C_ccall
C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp)
{
  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, init;
  k2 = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp);

  init = C_and(initp, C_make_character('\0'));
  C_allocate_vector(6, (C_word)NULL, k2,
                    C_bytes(C_fixnum_plus(size, C_fix(1))), /* Add header */
                    /* Byte vec, initialization, align at 8 bytes */
                    C_SCHEME_TRUE, init, C_SCHEME_FALSE);
}

static void
allocate_bignum_2(C_word c, C_word self, C_word bigvec)
{
  C_word ab[C_SIZEOF_STRUCTURE(2)], *a = ab, bignum,
         k = C_block_item(self, 1),
         negp = C_truep(C_block_item(self, 2)),
         size = C_bytestowords(C_header_size(bigvec))-1;

  C_word tagvec = CHICKEN_gc_root_ref(tags);

  C_set_block_item(bigvec, 0, negp ? C_BIGNUM_HEADER_SIGN_BIT | size : size);

  bignum = C_structure(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
  C_kontinue(k, bignum);
}