ikarus/bin/ikarus-numerics.c

#include "ikarus.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gmp.h>


#define most_positive_fixnum 0x1FFFFFFF
#define most_negative_fixnum 0x20000000

#define max_digits_per_limb 10

#ifdef NDEBUG
#define verify_bignum(x,caller) (x)
#else
static ikp
verify_bignum(ikp x, char* caller){
  if(tagof(x) != vector_tag){
    fprintf(stderr, "Error in (%s) invalid primary tag %p\n", caller, x);
    exit(-1);
  }
  ikp fst = ref(x, -vector_tag);
  int limb_count = ((unsigned int) fst) >> bignum_length_shift;
  if(limb_count <= 0){
    fprintf(stderr, 
        "Error in (%s) invalid limb count in fst=0x%08x\n", 
        caller, (int)fst);
    exit(-1);
  }
  int pos;
  if((int)fst & bignum_sign_mask){
    pos = 0;
  } else {
    pos = 1;
  }
  unsigned int last_limb = 
    (unsigned int) ref(x, off_bignum_data + (limb_count - 1) * wordsize);
  if(last_limb == 0){
    fprintf(stderr, 
        "Error in (%s) invalid last limb = 0x%08x", caller, last_limb);
    exit(-1);
  }
  if(limb_count == 1){
    if(pos){
      if(last_limb <= most_positive_fixnum){
        fprintf(stderr, 
                "Error in (%s) should be a positive fixnum: 0x%08x\n", 
                caller, last_limb);
        exit(-1);
      }
    } else {
      if(last_limb <= most_negative_fixnum){
        fprintf(stderr, 
                "Error in (%s) should be a negative fixnum: 0x%08x\n", 
                caller, last_limb);
        exit(-1);
      }
    }
  }
  /* ok */
  return x;
}
#endif

#define BN(x) verify_bignum(x,"BN")

ikp 
ikrt_isbignum(ikp x){
  if(tagof(x) == vector_tag){
    ikp fst = ref(x, -vector_tag);
    if (bignum_tag == (bignum_mask & (int)fst)){
      return true_object;
    }
  }
  return false_object;
}

ikp
ikrt_positive_bn(ikp x){
  ikp fst = ref(x, -vector_tag);
  if((int)fst & bignum_sign_mask){
    return false_object;
  } else {
    return true_object;
  }
}

ikp 
ikrt_even_bn(ikp x){
  int fst = (int)ref(x, wordsize-vector_tag);
  if(fst & 1){
    return false_object;
  } else {
    return true_object;
  }
}



ikp 
ikrt_fxfxplus(ikp x, ikp y, ikpcb* pcb){
  int n1 = unfix(x);
  int n2 = unfix(y);
  int r = n1 + n2;
  ikp q = fix(r);
  if(r == unfix(q)){
    return q;
  } 
  else {
    ikp bn = ik_alloc(pcb, align(disp_bignum_data + wordsize));
    if(r > 0){
      ref(bn, 0) = (ikp)(bignum_tag | (1 << bignum_length_shift)); 
      ref(bn, disp_bignum_data) = (ikp)r;
    }
    else {
      ref(bn, 0) = 
        (ikp)(bignum_tag | 
              (1 << bignum_length_shift) | 
              (1 << bignum_sign_shift));
      ref(bn, disp_bignum_data) = (ikp)-r;
    }
    return verify_bignum(bn+vector_tag, "fxfx+");
  }
}

ikp 
ikrt_fxbnplus(ikp x, ikp y, ikpcb* pcb){
  if(x == 0){ return y ; }
  ikp fst = ref(y, -vector_tag);
  int limb_count = ((unsigned int) fst) >> bignum_length_shift;
  int intx = unfix(x);
  if(intx > 0){
    if((bignum_sign_mask & (int)fst) == 0){
      /* positive fx + positive bn = even bigger positive */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+(limb_count+1)*wordsize));
      int carry = mpn_add_1((mp_limb_t*)(r+disp_bignum_data), 
                            (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                            limb_count,
                            intx);
      if(carry){
        ref(r, disp_bignum_data + limb_count*wordsize) = (ikp)1;
        ref(r, 0) = (ikp)
             (((limb_count + 1) << bignum_length_shift) |
              (0 << bignum_sign_shift) |
              bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn+1");
      } else {
        ref(r, 0) = (ikp)
          ((limb_count << bignum_length_shift) |
           (0 << bignum_sign_shift) |
           bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn+2");
      }
    }
    else {
      fprintf(stderr, "this case 0x%08x\n", intx);
      /* positive fx + negative bn = smaller negative bn */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+limb_count*wordsize));
      int borrow = mpn_sub_1((mp_limb_t*)(r+disp_bignum_data), 
                             (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                             limb_count,
                             intx);
      if(borrow){
        fprintf(stderr, "Error: BUG in borrow1 %d\n", borrow);
        exit(-1);
      }
      int result_size = 
        (ref(r, disp_bignum_data + (limb_count-1)*wordsize)) 
        ? limb_count
        : (limb_count - 1);
      if(result_size == 0){ 
        return 0;
      }
      if(result_size == 1){
        unsigned int last = 
          (unsigned int) ref(r, disp_bignum_data + (result_size-1)*wordsize);
        if(last <= most_negative_fixnum){
          return fix(-(int)last);
        }
      }
      ref(r, 0) = (ikp)
        ((result_size << bignum_length_shift) |
         (1 << bignum_sign_shift) |
         bignum_tag);
      return verify_bignum(r+vector_tag, "fxbn+3");
    }
  }
  else {
    if((bignum_sign_mask & (int)fst) == 0){
      /* negative fx + positive bn = smaller positive */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+limb_count*wordsize));
      int borrow = mpn_sub_1((mp_limb_t*)(r+disp_bignum_data), 
                             (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                             limb_count,
                             - intx);
      if(borrow){
        fprintf(stderr, "Error: BUG in borrow2\n");
        exit(-1);
      }
      int result_size = 
        (ref(r, disp_bignum_data + (limb_count-1)*wordsize) == 0) 
        ? (limb_count - 1)
        : limb_count;
      if(result_size == 0){
        return 0;
      }
      if(result_size == 1){
        unsigned int last = 
          (unsigned int) ref(r, disp_bignum_data + (result_size-1)*wordsize);
        if(last <= most_positive_fixnum){
          return fix((int)last);
        }
      }
      ref(r, 0) = (ikp)
        ((result_size << bignum_length_shift) |
         (0 << bignum_sign_shift) |
         bignum_tag);
      return verify_bignum(r+vector_tag, "fxbn+4");
    } else {
      /* negative fx + negative bn = larger negative */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+(limb_count+1)*wordsize));
      int carry = mpn_add_1((mp_limb_t*)(r+disp_bignum_data), 
                            (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                            limb_count,
                            -intx);
      if(carry){
        ref(r, disp_bignum_data + limb_count*wordsize) = (ikp)1;
        ref(r, 0) = (ikp)
             (((limb_count + 1) << bignum_length_shift) |
              (1 << bignum_sign_shift) |
              bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn+5");
      } else {
        ref(r, 0) = (ikp)
          ((limb_count << bignum_length_shift) |
           (1 << bignum_sign_shift) |
           bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn+5");
      }
    }
  }
}




ikp 
ikrt_bnbnplus(ikp x, ikp y, ikpcb* pcb){
  unsigned int xfst = (unsigned int)ref(x, -vector_tag);
  unsigned int yfst = (unsigned int)ref(y, -vector_tag);
  int xsign = xfst & bignum_sign_mask;
  int ysign = yfst & bignum_sign_mask;
  int xlimbs = xfst >> bignum_length_shift;
  int ylimbs = yfst >> bignum_length_shift;
  if(xsign == ysign){
    int n1,n2;
    ikp s1,s2;
    if(xlimbs > ylimbs){
      n1 = xlimbs; n2 = ylimbs; s1 = x; s2 = y;
    } else {
      n1 = ylimbs; n2 = xlimbs; s1 = y; s2 = x;
    }
    ikp res = ik_alloc(pcb, align(disp_bignum_data + (n1+1)*wordsize));
    mp_limb_t carry = mpn_add((mp_limb_t*) (res+disp_bignum_data),
                              (mp_limb_t*) (s1-vector_tag+disp_bignum_data),
                              n1,
                              (mp_limb_t*) (s2-vector_tag+disp_bignum_data),
                              n2);
    if(carry){
      ref(res, disp_vector_data + xlimbs*wordsize) = (ikp)1;
      ref(res, 0) = (ikp)
                    (((n1+1) << bignum_length_shift) |
                     xsign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn+1");
    } else {
      ref(res, 0) = (ikp)
                    ((n1 << bignum_length_shift) |
                     xsign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn+2");
    }
  }
  else {
    ikp s1=x, s2=y;
    int n1=xlimbs, n2=ylimbs;
    int result_sign = xsign;
    while((xlimbs == ylimbs) && 
          (ref(x, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize) == 
           ref(y, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize))){
      xlimbs -= 1;
      ylimbs -= 1;
      if(xlimbs == 0){ return 0; }
    }
    /* |x| != |y| */
    if(xlimbs <= ylimbs){
      if(xlimbs == ylimbs){
        if((ref(y, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize) > 
            ref(x, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize))){
          s1 = y; n1 = ylimbs;
          s2 = x; n2 = xlimbs;
          result_sign = ysign;
        }
      } else {
        s1 = y; n1 = ylimbs;
        s2 = x; n2 = xlimbs;
        result_sign = ysign;
      }
    }
    /* |s1| > |s2| */
    ikp res = ik_alloc(pcb, align(disp_bignum_data + n1 * wordsize));
    int burrow = mpn_sub((mp_limb_t*) (res + disp_bignum_data),
                         (mp_limb_t*) (s1 - vector_tag + disp_bignum_data),
                         n1,
                         (mp_limb_t*) (s2 - vector_tag + disp_bignum_data),
                         n2);
    if(burrow){
      fprintf(stderr, "BUG: Burrow error in bnbn+\n");
      exit(-1);
    }
    int len = n1;
    while(ref(res, disp_bignum_data + (len-1)*wordsize) == 0){
      len--;
      if(len == 0){
        return 0;
      }
    }
    if(result_sign == 0){
      /* positive result */
      if(len == 1){
        unsigned int fst_limb = (unsigned int) ref(res, disp_bignum_data);
        if(fst_limb <= most_positive_fixnum){
          return fix((int)fst_limb);
        }
      }
      ref(res, 0) = (ikp)
                    ((len << bignum_length_shift) |
                     result_sign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn+3");
    } else {
      /* negative result */
      if(len == 1){
        unsigned int fst_limb = (unsigned int) ref(res, disp_bignum_data);
        if(fst_limb <= most_negative_fixnum){
          return fix(-(int)fst_limb);
        }
      }
      ref(res, 0) = (ikp)
                    ((len << bignum_length_shift) |
                     result_sign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn+4");
    }
  }
}




ikp 
ikrt_fxfxminus(ikp x, ikp y, ikpcb* pcb){
  int n1 = unfix(x);
  int n2 = unfix(y);
  int r = n1 - n2;
  if(r >= 0){
    if(((unsigned int)r) <= most_positive_fixnum){
      return fix(r);
    } else {
      ikp bn = ik_alloc(pcb, align(disp_bignum_data + wordsize));
      ref(bn, 0) = (ikp) (bignum_tag | (1 << bignum_length_shift));
      ref(bn, disp_bignum_data) = (ikp)r;
      return verify_bignum(bn+vector_tag,"fxfx-1");
    }
  } else {
    ikp fxr = fix(r);
    if(unfix(fxr) == r){
      return fxr;
    } else {
      ikp bn = ik_alloc(pcb, align(disp_bignum_data + wordsize));
      ref(bn, 0) = (ikp) 
        (bignum_tag | 
         (1 << bignum_sign_shift) |
         (1 << bignum_length_shift));
      ref(bn, disp_bignum_data) = (ikp)(-r);
      return verify_bignum(bn+vector_tag, "fxfx-2");
    }
  }
}


ikp 
ikrt_bnnegate(ikp x, ikpcb* pcb){
  ikp fst = ref(x, -vector_tag);
  int limb_count = ((unsigned int) fst) >> bignum_length_shift;
  if(limb_count == 1){
    if((bignum_sign_mask & (int)fst) == 0){
      /* positive bignum */
      unsigned int limb = (unsigned int) ref(x, disp_bignum_data - vector_tag);
      if(limb == (most_positive_fixnum + 1)){
        return fix(-(int)limb);
      }
    }
  }
  ikp bn = ik_alloc(pcb, align(disp_bignum_data + limb_count * wordsize));
  memcpy(bn+disp_bignum_data,
         x-vector_tag+disp_bignum_data,
         limb_count*wordsize);
  ref(bn, 0) = (ikp)
    (bignum_tag |
     ((1 << bignum_sign_shift) - (bignum_sign_mask & (int)fst)) |
     (limb_count << bignum_length_shift));
  return verify_bignum(bn+vector_tag, "bnneg");
}

ikp 
ikrt_fxbnminus(ikp x, ikp y, ikpcb* pcb){
  if(x == 0){ return ikrt_bnnegate(y, pcb) ; }
  ikp fst = ref(y, -vector_tag);
  int limb_count = ((unsigned int) fst) >> bignum_length_shift;
  int intx = unfix(x);
  if(intx > 0){
    if(bignum_sign_mask & (int)fst){
      /* positive fx - negative bn = positive bn */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+(limb_count+1)*wordsize));
      int carry = mpn_add_1((mp_limb_t*)(r+disp_bignum_data), 
                            (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                            limb_count,
                            intx);
      if(carry){
        ref(r, disp_bignum_data + limb_count*wordsize) = (ikp)1;
        ref(r, 0) = (ikp)
             (((limb_count + 1) << bignum_length_shift) |
              (0 << bignum_sign_shift) |
              bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn-1");
      } else {
        ref(r, 0) = (ikp)
          ((limb_count << bignum_length_shift) |
           (0 << bignum_sign_shift) |
           bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn-2");
      }
    }
    else {
      /* positive fx - positive bn = smaller negative bn/fx */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+limb_count*wordsize));
      int borrow = mpn_sub_1((mp_limb_t*)(r+disp_bignum_data), 
                             (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                             limb_count,
                             intx);
      if(borrow){
        fprintf(stderr, "Error: BUG in borrow3\n");
        exit(-1);
      }
      int result_size = 
        (ref(r, disp_bignum_data + (limb_count-1)*wordsize)) 
        ? limb_count
        : (limb_count - 1);
      if(result_size == 0){ 
        return 0;
      }
      if(result_size == 1){
        unsigned int last = 
          (unsigned int) ref(r, disp_bignum_data + (result_size-1)*wordsize);
        if(last <= most_negative_fixnum){
          return fix(-(int)last);
        }
      }
      ref(r, 0) = (ikp)
        ((result_size << bignum_length_shift) |
         (1 << bignum_sign_shift) |
         bignum_tag);
      return verify_bignum(r+vector_tag, "fxbn-");
    }
  }
  else {
    if(bignum_sign_mask & (int)fst){
      /* negative fx - negative bn = smaller positive */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+limb_count*wordsize));
      int borrow = mpn_sub_1((mp_limb_t*)(r+disp_bignum_data), 
                             (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                             limb_count,
                             - intx);
      if(borrow){
        fprintf(stderr, "Error: BUG in borrow4\n");
        exit(-1);
      }
      int result_size = 
        (ref(r, disp_bignum_data + (limb_count-1)*wordsize) == 0) 
        ? (limb_count - 1)
        : limb_count;
      if(result_size == 0){
        return 0;
      }
      if(result_size == 1){
        unsigned int last = 
          (unsigned int) ref(r, disp_bignum_data + (result_size-1)*wordsize);
        if(last <= most_positive_fixnum){
          return fix((int)last);
        }
      }
      ref(r, 0) = (ikp)
        ((result_size << bignum_length_shift) |
         (0 << bignum_sign_shift) |
         bignum_tag);
      return verify_bignum(r+vector_tag,"fxbn-");
    } else {
      /* negative fx - positive bn = larger negative */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+(limb_count+1)*wordsize));
      int carry = mpn_add_1((mp_limb_t*)(r+disp_bignum_data), 
                            (mp_limb_t*)(y - vector_tag + disp_bignum_data),
                            limb_count,
                            -intx);
      if(carry){
        ref(r, disp_bignum_data + limb_count*wordsize) = (ikp)1;
        ref(r, 0) = (ikp)
             (((limb_count + 1) << bignum_length_shift) |
              (1 << bignum_sign_shift) |
              bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn-");
      } else {
        ref(r, 0) = (ikp)
          ((limb_count << bignum_length_shift) |
           (1 << bignum_sign_shift) |
           bignum_tag);
        return verify_bignum(r+vector_tag, "fxbn-");
      }
    }
  }
}

ikp 
ikrt_bnfxminus(ikp x, ikp y, ikpcb* pcb){
  if(y == 0){ return x; }
  ikp fst = ref(x, -vector_tag);
  int limb_count = ((unsigned int) fst) >> bignum_length_shift;
  int inty = unfix(y);
  if(inty < 0){
    if((bignum_sign_mask & (int)fst) == 0){
      /* - negative fx + positive bn = positive bn */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+(limb_count+1)*wordsize));
      int carry = mpn_add_1((mp_limb_t*)(r+disp_bignum_data), 
                            (mp_limb_t*)(x - vector_tag + disp_bignum_data),
                            limb_count,
                            -inty);
      if(carry){
        ref(r, disp_bignum_data + limb_count*wordsize) = (ikp)1;
        ref(r, 0) = (ikp)
             (((limb_count + 1) << bignum_length_shift) |
              (0 << bignum_sign_shift) |
              bignum_tag);
        return verify_bignum(r+vector_tag,"bnfx-");
      } else {
        ref(r, 0) = (ikp)
          ((limb_count << bignum_length_shift) |
           (0 << bignum_sign_shift) |
           bignum_tag);
        return verify_bignum(r+vector_tag,"bnfx-");
      }
    }
    else {
      /* - negative fx + negative bn = smaller negative bn/fx */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+limb_count*wordsize));
      int borrow = mpn_sub_1((mp_limb_t*)(r+disp_bignum_data), 
                             (mp_limb_t*)(x - vector_tag + disp_bignum_data),
                             limb_count,
                             -inty);
      if(borrow){
        fprintf(stderr, "Error: BUG in borrow5\n");
        exit(-1);
      }
      int result_size = 
        (ref(r, disp_bignum_data + (limb_count-1)*wordsize)) 
        ? limb_count
        : (limb_count - 1);
      if(result_size == 0){ 
        return 0;
      }
      if(result_size == 1){
        unsigned int last = 
          (unsigned int) ref(r, disp_bignum_data + (result_size-1)*wordsize);
        if(last <= most_negative_fixnum){
          return fix(-(int)last);
        }
      }
      ref(r, 0) = (ikp)
        ((result_size << bignum_length_shift) |
         (1 << bignum_sign_shift) |
         bignum_tag);
      return verify_bignum(r+vector_tag,"bnfx-");
    }
  }
  else {
    if((bignum_sign_mask & (int)fst) == 0){
      /* - positive fx + positive bn = smaller positive */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+limb_count*wordsize));
      int borrow = mpn_sub_1((mp_limb_t*)(r+disp_bignum_data), 
                             (mp_limb_t*)(x - vector_tag + disp_bignum_data),
                             limb_count,
                             inty);
      if(borrow){
        fprintf(stderr, "Error: BUG in borrow6\n");
        exit(-1);
      }
      int result_size = 
        (ref(r, disp_bignum_data + (limb_count-1)*wordsize) == 0) 
        ? (limb_count - 1)
        : limb_count;
      if(result_size == 0){
        return 0;
      }
      if(result_size == 1){
        unsigned int last = 
          (unsigned int) ref(r, disp_bignum_data + (result_size-1)*wordsize);
        if(last <= most_positive_fixnum){
          return fix((int)last);
        }
      }
      ref(r, 0) = (ikp)
        ((result_size << bignum_length_shift) |
         (0 << bignum_sign_shift) |
         bignum_tag);
      return verify_bignum(r+vector_tag, "bnfx-");
    } else {
      /* - positive fx + negative bn = larger negative */
      ikp r = ik_alloc(pcb, align(disp_bignum_data+(limb_count+1)*wordsize));
      int carry = mpn_add_1((mp_limb_t*)(r+disp_bignum_data), 
                            (mp_limb_t*)(x - vector_tag + disp_bignum_data),
                            limb_count,
                            inty);
      if(carry){
        ref(r, disp_bignum_data + limb_count*wordsize) = (ikp)1;
        ref(r, 0) = (ikp)
             (((limb_count + 1) << bignum_length_shift) |
              (1 << bignum_sign_shift) |
              bignum_tag);
        return verify_bignum(r+vector_tag, "bnfx-");
      } else {
        ref(r, 0) = (ikp)
          ((limb_count << bignum_length_shift) |
           (1 << bignum_sign_shift) |
           bignum_tag);
        return verify_bignum(r+vector_tag, "bnfx-");
      }
    }
  }
}



ikp 
ikrt_bnbnminus(ikp x, ikp y, ikpcb* pcb){
  if(x == y) { return 0; }
  unsigned int xfst = (unsigned int)ref(x, -vector_tag);
  unsigned int yfst = (unsigned int)ref(y, -vector_tag);
  int xsign = xfst & bignum_sign_mask;
  int ysign = yfst & bignum_sign_mask;
  int xlimbs = xfst >> bignum_length_shift;
  int ylimbs = yfst >> bignum_length_shift;
  if(xsign != ysign){
    int n1,n2;
    ikp s1,s2;
    if(xlimbs >= ylimbs){
      n1 = xlimbs; n2 = ylimbs; s1 = x; s2 = y;
    } else {
      n1 = ylimbs; n2 = xlimbs; s1 = y; s2 = x;
    }
    ikp res = ik_alloc(pcb, align(disp_bignum_data + (n1+1)*wordsize));
    mp_limb_t carry = mpn_add((mp_limb_t*) (res+disp_bignum_data),
                              (mp_limb_t*) (s1-vector_tag+disp_bignum_data),
                              n1,
                              (mp_limb_t*) (s2-vector_tag+disp_bignum_data),
                              n2);
    if(carry){
      ref(res, disp_vector_data + xlimbs*wordsize) = (ikp)1;
      ref(res, 0) = (ikp)
                    (((n1+1) << bignum_length_shift) |
                     xsign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn-");
    } else {
      ref(res, 0) = (ikp)
                    ((n1 << bignum_length_shift) |
                     xsign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn-");
    }
  }
  else {
    /* same sign */
    if(xlimbs == ylimbs){
      while((ref(x, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize) == 
             ref(y, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize))){
        xlimbs -= 1;
        if(xlimbs == 0){ return 0; }
      }
      ylimbs = xlimbs;
    }
    ikp s1=x, s2=y;
    int n1=xlimbs, n2=ylimbs;
    int result_sign = xsign;
    /* |x| != |y| */
    if(xlimbs <= ylimbs){
      if(xlimbs == ylimbs){
        if((ref(y, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize) > 
            ref(x, -vector_tag+disp_bignum_data+(xlimbs-1)*wordsize))){
          s1 = y; n1 = ylimbs;
          s2 = x; n2 = xlimbs;
          result_sign = (1 << bignum_length_shift) - ysign;
        }
      } else {
        s1 = y; n1 = ylimbs;
        s2 = x; n2 = xlimbs;
        result_sign = (1 << bignum_length_shift) - ysign;
      }
    }
    /* |s1| > |s2| */
    ikp res = ik_alloc(pcb, align(disp_bignum_data + n1 * wordsize));
    int burrow = mpn_sub((mp_limb_t*) (res + disp_bignum_data),
                         (mp_limb_t*) (s1 - vector_tag + disp_bignum_data),
                         n1,
                         (mp_limb_t*) (s2 - vector_tag + disp_bignum_data),
                         n2);
    if(burrow){
      fprintf(stderr, "BUG: Burrow error in bnbn-\n");
      exit(-1);
    }
    int len = n1;
    while(ref(res, disp_bignum_data + (len-1)*wordsize) == 0){
      len--;
      if(len == 0){
        return 0;
      }
    }
    if(result_sign == 0){
      /* positive result */
      if(len == 1){
        unsigned int fst_limb = (unsigned int) ref(res, disp_bignum_data);
        if(fst_limb <= most_positive_fixnum){
          return fix((int)fst_limb);
        }
      }
      ref(res, 0) = (ikp)
                    ((len << bignum_length_shift) |
                     result_sign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn-");
    } else {
      /* negative result */
      if(len == 1){
        unsigned int fst_limb = (unsigned int) ref(res, disp_bignum_data);
        if(fst_limb <= most_negative_fixnum){
          return fix(-(int)fst_limb);
        }
      }
      ref(res, 0) = (ikp)
                    ((len << bignum_length_shift) |
                     result_sign |
                     bignum_tag);
      return verify_bignum(res+vector_tag, "bnbn-");
    }
  }
}


ikp 
ikrt_fxfxmult(ikp x, ikp y, ikpcb* pcb){
  int n1 = unfix(x);
  int n2 = unfix(y);
  mp_limb_t lo = 0;
  mp_limb_t s1 = n1;
  mp_limb_t s2 = n2;
  int sign = 0;
  if(n1 < 0){
    s1 = -n1;
    sign = 1 - sign;
  } 
  if(n2 < 0){
    s2 = -n2;
    sign = 1 - sign;
  } 
  mp_limb_t hi = mpn_mul_1(&lo, &s1, 1, s2);
  if(hi == 0){
    if(sign){
      if(lo <= most_negative_fixnum){
        return fix(-((int)lo));
      } 
    } else {
      if(lo <= most_positive_fixnum){
        return fix((int)lo);
      }
    }
    ikp r = ik_alloc(pcb, disp_bignum_data + wordsize);
    ref(r, 0) = (ikp)
      (bignum_tag |
       (sign << bignum_sign_shift) |
       (1 << bignum_length_shift));
    ref(r, disp_bignum_data) = (ikp)lo;
    return BN(r+vector_tag);
  } else {
    ikp r = ik_alloc(pcb, align(disp_bignum_data + 2*wordsize));
    ref(r, 0) = (ikp)
      (bignum_tag | 
       (sign << bignum_sign_shift) |
       (2 << bignum_length_shift));
    ref(r, disp_bignum_data) = (ikp)lo;
    ref(r, disp_bignum_data+wordsize) = (ikp)hi;
    return BN(r+vector_tag);
  }
}

ikp
normalize_bignum(int limbs, int sign, ikp r){
  while(ref(r, disp_bignum_data + (limbs-1)*wordsize) == 0){
    limbs--;
    if(limbs == 0){ return 0;}
  }
  if(limbs == 1){
    unsigned int last = (unsigned int) ref(r, disp_bignum_data);
    if(sign == 0){
      if(last <= most_positive_fixnum){
        return fix((int)last);
      }
    } else {
      if(last <= most_negative_fixnum){
        return fix(-((int)last));
      }
    }
  }
  ref(r, 0) = (ikp)
    (bignum_tag |
     sign |
     (limbs << bignum_length_shift));
  return BN(r+vector_tag);
}


ikp 
ikrt_fxbnmult(ikp x, ikp y, ikpcb* pcb){
  int n2 = unfix(x);
  if(n2 == 0) { return 0; }
  mp_limb_t s2 = (n2>0) ? n2 : (- n2);
  ikp fst = ref(y, -vector_tag);
  int limb_count = ((unsigned int) fst) >> bignum_length_shift;
  ikp r = ik_alloc(pcb, align(disp_bignum_data + (limb_count+1)*wordsize));
  mp_limb_t hi = mpn_mul_1((mp_limb_t*)(r+disp_bignum_data),
                           (mp_limb_t*)(y-vector_tag+disp_bignum_data),
                           limb_count,
                           s2);
  ref(r, disp_bignum_data + limb_count * wordsize) = (ikp)hi;
  int sign = 
    ((n2 > 0) ?
     (bignum_sign_mask & (int)fst) : 
     ((1 << bignum_sign_shift) - (bignum_sign_mask&(int)fst)));
  return normalize_bignum(limb_count+1, sign, r);
}

ikp 
ikrt_bnbnmult(ikp x, ikp y, ikpcb* pcb){
  int f1 = (int)ref(x, -vector_tag);
  int f2 = (int)ref(y, -vector_tag);
  int n1 = ((unsigned int)f1) >> bignum_length_shift;
  int n2 = ((unsigned int)f2) >> bignum_length_shift;
  int nr = n1 + n2;
  ikp bn = ik_alloc(pcb, align(disp_bignum_data + nr*wordsize));
  mp_limb_t r;
  if(n1 >= n2){
    r = mpn_mul((mp_limb_t*)(bn+disp_bignum_data),
                (mp_limb_t*)(x-vector_tag+disp_bignum_data),
                n1,
                (mp_limb_t*)(y-vector_tag+disp_bignum_data),
                n2);
  } else {
    r = mpn_mul((mp_limb_t*)(bn+disp_bignum_data),
                (mp_limb_t*)(y-vector_tag+disp_bignum_data),
                n2,
                (mp_limb_t*)(x-vector_tag+disp_bignum_data),
                n1);
  } 
  int sign = 
    ((bignum_sign_mask & f1) ?
     ((1 << bignum_sign_shift) - (bignum_sign_mask & f2)) :
     (bignum_sign_mask & f2));
  return normalize_bignum(nr, sign, bn);
}




ikp
ikrt_bnbncomp(ikp bn1, ikp bn2){
  ikp f1 = ref(bn1, -vector_tag);
  ikp f2 = ref(bn2, -vector_tag);
  if(bignum_sign_mask & (int)f1){
    if(bignum_sign_mask & (int)f2){
      /* both negative */
      int n1 = ((unsigned int) f1) >> bignum_length_shift;
      int n2 = ((unsigned int) f2) >> bignum_length_shift;
      if(n1 < n2) {
        return fix(1);
      } else if(n1 > n2){
        return fix(-1);
      } else {
        int i;
        for(i=(n1-1); i>=0; i--){
          unsigned int t1 = 
            (unsigned int) ref(bn1,disp_bignum_data-vector_tag+i*wordsize);
          unsigned int t2 = 
            (unsigned int) ref(bn2,disp_bignum_data-vector_tag+i*wordsize);
          if(t1 < t2){
            return fix(1);
          } else if(t1 > t2){
            return fix(-1);
          }
        }
      }
      return 0;
    } else {
      /* n1 negative, n2 positive */
      return fix(-1);
    } 
  } else {
    if(bignum_sign_mask & (int)f2){
      /* n1 positive, n2 negative */
      return fix(1);
    } else {
      /* both positive */
      int n1 = ((unsigned int) f1) >> bignum_length_shift;
      int n2 = ((unsigned int) f2) >> bignum_length_shift;
      if(n1 < n2) {
        return fix(-1);
      } else if(n1 > n2){
        return fix(1);
      } else {
        int i;
        for(i=(n1-1); i>=0; i--){
          unsigned int t1 = 
           (unsigned int) ref(bn1,disp_bignum_data-vector_tag+i*wordsize);
          unsigned int t2 = 
            (unsigned int) ref(bn2,disp_bignum_data-vector_tag+i*wordsize);
          if(t1 < t2){
            return fix(-1);
          } else if(t1 > t2){
            return fix(1);
          }
        }
      }
      return 0;
    }
  }
}

ikp
ikrt_fxbnlogand(ikp x, ikp y, ikpcb* pcb){
  int n1 = unfix(x);
  ikp fst = ref(y, -vector_tag);
  if(n1 >= 0){
    if(bignum_sign_mask & (unsigned int) fst){
      /* y is negative */
      return fix(n1 & (1+~(int)ref(y, disp_vector_data-vector_tag))); 
    } else {
      /* y is positive */
      return fix(n1 & (int)ref(y, disp_vector_data-vector_tag)); 
    }
  } else {
    if(n1 == -1){ return y; }
    if(bignum_sign_mask & (unsigned int) fst){
      /* y is negative */
      int len = (((unsigned int) fst) >> bignum_length_shift);
      unsigned int nn = 
        (1+~((1+~(int)ref(y, disp_bignum_data - vector_tag)) & n1));
      if((len == 1) && (nn <= most_negative_fixnum)){
        return fix(-nn);
      }
      ikp r = ik_alloc(pcb, align(disp_bignum_data + len * wordsize));
      ref(r, 0) = fst;
      ref(r, disp_bignum_data) = (ikp) nn;
      int i;
      for(i=1; i<len; i++){
        ref(r, disp_bignum_data+i*wordsize) = 
          ref(y, disp_bignum_data-vector_tag+i*wordsize);
      }
      return BN(r+vector_tag);
    } else {
      /* y is positive */
      int len = (((unsigned int) fst) >> bignum_length_shift);
      ikp r = ik_alloc(pcb, align(disp_bignum_data + len * wordsize));
      ref(r, 0) = fst;
      ref(r, disp_bignum_data) = (ikp)
        (((int)ref(y, disp_bignum_data - vector_tag)) & n1);
      int i;
      for(i=1; i<len; i++){
        ref(r, disp_bignum_data+i*wordsize) = 
          ref(y, disp_bignum_data-vector_tag+i*wordsize);
      }
      return BN(r+vector_tag);
    }
  }
}


static inline int
count_leading_ffs(int n, unsigned int* x){
  int idx;
  for(idx=0; idx<n; idx++){
    if(x[idx] != -1){
      return idx;
    }
  }
  return n;
}

ikp
ikrt_bnbnlogand(ikp x, ikp y, ikpcb* pcb){
  ikp xfst = ref(x, -vector_tag);
  ikp yfst = ref(y, -vector_tag);
  int n1 = ((unsigned int) xfst) >> bignum_length_shift;
  int n2 = ((unsigned int) yfst) >> bignum_length_shift;
  if(bignum_sign_mask & (unsigned int) xfst){
    if(bignum_sign_mask & (unsigned int) yfst){
      fprintf(stderr, "not yet for bnbnlogand\n");
      exit(-1);
    } else {
      return ikrt_bnbnlogand(y,x,pcb);
    }
  } else {
    if(bignum_sign_mask & (unsigned int) yfst){
      /* x positive, y negative */
      /*  the result is at most n1 words long */
      unsigned int* s1 = ((unsigned int*)(x+disp_bignum_data-vector_tag));
      unsigned int* s2 = ((unsigned int*)(y+disp_bignum_data-vector_tag));
      if(n1 <= n2){
        int i = n1-1;
        while(i >= 0){
          unsigned int t = s1[i] & (1+~s2[i]);
          if(t != 0){
            if((i == 0) && (t <= most_positive_fixnum)){
              return fix(t);
            }
            ikp r = ik_alloc(pcb, align(disp_bignum_data+(i+1)*wordsize));
            ref(r, 0) = (ikp) (bignum_tag | ((i+1) << bignum_length_shift));
            unsigned int* s = (unsigned int*)(r+disp_bignum_data);
            s[i] = t;
            for(i--; i>=0; i--){
              s[i] = s1[i] & (1+~s2[i]);
            }
            return r+vector_tag;
          } else {
            i--;
          }
        }
        return 0;
      } else {
        fprintf(stderr, "not yet for bnbnlogand\n");
        exit(-1);
      }
    } else {
      /* both positive */
      int n = (n1<n2)?n1:n2;
      int i;
      for(i=n-1; i>=0; i--){
        int l1 = 
          (int) ref(x, disp_bignum_data-vector_tag+i*wordsize);
        int l2 = 
          (int) ref(y, disp_bignum_data-vector_tag+i*wordsize);
        int last = l1 & l2;
        if(last){
          if((i == 0) && (last < most_positive_fixnum)){
            return fix(last);
          }
          ikp r = ik_alloc(pcb, align(disp_bignum_data+(i+1)*wordsize));
          ref(r, 0) = (ikp) (bignum_tag | ((i+1)<<bignum_length_shift));
          ref(r, disp_bignum_data + i*wordsize) = (ikp)last;
          int j;
          for(j=0; j<i; j++){
            ref(r, disp_bignum_data + j*wordsize) = (ikp)
              (((int)ref(x, disp_bignum_data-vector_tag+j*wordsize))
               &
               ((int)ref(y, disp_bignum_data-vector_tag+j*wordsize)));
          }
          return r+vector_tag;
        }
      }
      return 0;
    }
  }
}

#if 0
From TFM:
void
mpn_tdiv_qr (
  mp limb t *qp,        /* quotient placed here */
  mp limb t *rp,        /* remainder placed here */
  mp size t qxn,        /* must be zero! */ 
  const mp limb t *np,  /* first number  */
  mp size t nn,         /* its length    */
  const mp limb t *dp,  /* second number */
  mp size t dn          /* its length    */
)

Divide {np, nn} by {dp, dn} and put the quotient at {qp,nn-dn+1}
and the remainder at {rp, dn}. The quotient is rounded towards 0.
No overlap is permitted between arguments. nn must be greater than
or equal to dn. The most significant limb of dp must be non-zero.
The qxn operand must be zero. 
#endif

ikp
ikrt_bnbndivrem(ikp x, ikp y, ikpcb* pcb){
  ikp xfst = ref(x, -vector_tag);
  ikp yfst = ref(y, -vector_tag);
  mp_size_t xn = ((unsigned int) xfst) >> bignum_length_shift;
  mp_size_t yn = ((unsigned int) yfst) >> bignum_length_shift;
  if(xn < yn){
    /* quotient is zero, remainder is x */
    ikp rv = ik_alloc(pcb, pair_size);
    ref(rv, disp_car) = 0;
    ref(rv, disp_cdr) = x;
    return rv+pair_tag;
  }
  mp_size_t qn = xn - yn + 1;
  mp_size_t rn = yn;
  ikp q = ik_alloc(pcb, align(disp_bignum_data + qn*wordsize));
  ikp r = ik_alloc(pcb, align(disp_bignum_data + rn*wordsize));
  mpn_tdiv_qr (
      (mp_limb_t*)(q+disp_bignum_data),
      (mp_limb_t*)(r+disp_bignum_data),
      0,
      (mp_limb_t*)(x+off_bignum_data),
      xn,
      (mp_limb_t*)(y+off_bignum_data),
      yn);

  if(((unsigned int) xfst) & bignum_sign_mask){
    /* x is negative => remainder is negative */
    r = normalize_bignum(rn, 1 << bignum_sign_shift, r);
  } else {
    r = normalize_bignum(rn, 0, r);
  }

  if(((unsigned int) yfst) & bignum_sign_mask){
    /* y is negative => quotient is opposite of x */
    int sign = 
      bignum_sign_mask - (((unsigned int)xfst) & bignum_sign_mask);
    q = normalize_bignum(qn, sign, q);
  } else {
    /* y is positive => quotient is same as x */
    int sign = (((unsigned int)xfst) & bignum_sign_mask);
    q = normalize_bignum(qn, sign, q);
  }
  ikp rv = ik_alloc(pcb, pair_size);
  ref(rv, disp_car) = q;
  ref(rv, disp_cdr) = r;
  return rv+pair_tag;
}


#if 0
[Function]

mp_limb_t
mpn_divrem_1 (
  mp limb t *r1p,
  mp size t qxn, 
  mp limb t *s2p,
  mp size t s2n,
  mp limb t s3limb 
) 

Divide {s2p, s2n} by s3limb, and write the quotient at r1p. Return the remainder. 
The integer quotient is written to {r1p+qxn, s2n} and in addition qxn fraction limbs are 
developed and written to {r1p, qxn}. Either or both s2n and qxn can be zero. For most 
usages, qxn will be zero. 
#endif

ikp
ikrt_bnfxdivrem(ikp x, ikp y, ikpcb* pcb){
  int yint = unfix(y);
  mp_limb_t* s2p = (mp_limb_t*)(x+off_bignum_data);
  ikp fst = ref(x, -vector_tag);
  mp_size_t s2n = ((unsigned int) fst) >> bignum_length_shift;
  ikp quot = ik_alloc(pcb,
      align(s2n*wordsize + disp_bignum_data));
  mp_limb_t rv = mpn_divrem_1(
      (mp_limb_t*)(quot+disp_bignum_data),
      0,
      s2p,
      s2n,
      abs(yint));

  ikp rem;

  if(yint < 0){
    /* y is negative => quotient is opposite of x */
    int sign = 
      bignum_sign_mask - (((unsigned int)fst) & bignum_sign_mask);
    quot = normalize_bignum(s2n, sign, quot);
  } else {
    /* y is positive => quotient is same as x */
    int sign = (((unsigned int)fst) & bignum_sign_mask);
    quot = normalize_bignum(s2n, sign, quot);
  }

  /* the remainder is always less than |y|, so it will
     always be a fixnum.  (if y == most_negative_fixnum, 
     then |remainder| will be at most most_positive_fixnum). */
  if(((unsigned int) fst) & bignum_sign_mask){
    /* x is negative => remainder is negative */
    rem = (ikp) -(rv << fx_shift);
  } else {
    rem = fix(rv);
  }
  ikp p = ik_alloc(pcb, pair_size);
  ref(p, disp_car) = quot;
  ref(p, disp_cdr) = rem;
  return p+pair_tag;
}






ikp
ikrt_bntostring(ikp x, ikpcb* pcb){
  /* FIXME: avoid calling malloc, instead, use the heap pointer itself
   * as a buffer to hold the temporary data after ensuring that it has enough
   * space */
  ikp fst = ref(x, -vector_tag);
  int limb_count = (((unsigned int)fst) >> bignum_length_shift);
  if(limb_count <= 0){
    fprintf(stderr, "BUG: nbtostring: invalid length %d\n", limb_count);
    exit(-1);
  }
  int sign_bit = bignum_sign_mask & (int) fst;
  int nbsize = limb_count * sizeof(mp_limb_t);
  int strsize = limb_count * max_digits_per_limb;
  int mem_req = nbsize + strsize + 1;
  unsigned char* mem = malloc(mem_req);
  if(! mem){
    fprintf(stderr, "Error allocating space for bignum\n");
    exit(-1);
  }
  memcpy(mem, x - vector_tag + disp_bignum_data, nbsize);
  mp_size_t bytes = 
    mpn_get_str(mem+nbsize,       /* output string */ 
                10,               /* base */
                (mp_limb_t*) mem, /* limb */
                limb_count        /* number of limbs */
        );
  unsigned char* string_start = mem + nbsize;
  while(*string_start == 0){
    string_start++;
    bytes--;
  }
  ikp str = ik_alloc(pcb, align(bytes + disp_string_data + (sign_bit?1:0)));
  ref(str, 0) = fix(bytes + (sign_bit?1:0));
  ikp dest = str + disp_string_data;
  if(sign_bit){
    *dest = '-';
    dest++;
  }
  { 
    int i = 0;
    while(i < bytes){
      dest[i] = string_start[i] + '0';
      i++;
    }
    dest[bytes] = 0;
  }
  free(mem);
  return str + string_tag;
}