#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") #if 0 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; } #endif 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_sign_shift) - ysign; } } else { s1 = y; n1 = ylimbs; s2 = x; n2 = xlimbs; result_sign = (1 << bignum_sign_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; }