patch bignums

[Ondrej Bilka <neleai@seznam.cz>]

I am so far with implementing faster bignums:

First refactoring are num_op where first argument can be
fixnum or bignum. It simplifies writing expresions.
pow refactored using num_ops.

Bugfix and implented bignum odd? and even? 

Now I use karatsuba multiplication and squaring. Array allocated was 1
bigger because result fits into x->len+y->len (Maximum is (2**x->len-1)*(2**y->len-1)<=(2**x+y-1))

Conversions: 
First optimalization is grouping charsperdig chars into one big digit.
When base is power of 2 then conversion is trivial.
For large strings is str2big algorithm fast as multiplication. (Unless fft is used when its log n times slower)
For large strings I have algorithm where big2str is fast as division.
I comented it out because division at O(n**2) makes it slower than trivial alg.

libtommath
Some volunteer who rewrites it that it uses all bits in its digit istead
using highest for carry in addition?
Then it would have ruby bignum compatible format.

BTW. How do you debug ruby. At README.EXT is this missing. 
I dont know easy way to inspect ruby objects inside gdb. I used printf+guessing what obj contain.
But now its "obvious bug" free.


Index: numeric.c
===================================================================
RCS file: /src/ruby/numeric.c,v
retrieving revision 1.148
diff -u -r1.148 numeric.c
--- numeric.c	21 Sep 2006 06:09:26 -0000	1.148
+++ numeric.c	21 Sep 2006 14:46:19 -0000
@@ -2,8 +2,8 @@
 
   numeric.c -
 
-  $Author: matz $
-  $Date: 2006/09/21 06:09:26 $
+  $Author: nobu $
+  $Date: 2006/09/17 01:42:28 $
   created at: Fri Aug 13 18:33:09 JST 1993
 
   Copyright (C) 1993-2003 Yukihiro Matsumoto
@@ -1724,12 +1724,14 @@
 static VALUE
 int_even_p(VALUE num)
 {
-    if (rb_funcall(num, '%', 1, INT2FIX(2)) == INT2FIX(0)) {
+    if (int_odd_p(num)==Qfalse) {
 	return Qtrue;
     }
     return Qfalse;
 }
 
+
+
 /*
  *  call-seq:
  *     int.next    => integer
@@ -1930,7 +1932,7 @@
  * result.
  */
 
-static VALUE
+VALUE
 fix_plus(VALUE x, VALUE y)
 {
     if (FIXNUM_P(y)) {
@@ -1963,7 +1965,7 @@
  * result.
  */
 
-static VALUE
+VALUE
 fix_minus(VALUE x, VALUE y)
 {
     if (FIXNUM_P(y)) {
@@ -1997,7 +1999,7 @@
  * result.
  */
 
-static VALUE
+VALUE
 fix_mul(VALUE x, VALUE y)
 {
     if (FIXNUM_P(y)) {
@@ -2236,30 +2238,11 @@
  *    2 ** 0.5    #=> 1.4142135623731
  */
 
+VALUE rb_num_pow(VALUE x,VALUE y);
 static VALUE
 fix_pow(VALUE x, VALUE y)
 {
-    if (FIXNUM_P(y)) {
-	long a, b;
-
-	b = FIX2LONG(y);
-	if (b == 0) return INT2FIX(1);
-	if (b == 1) return x;
-	a = FIX2LONG(x);
-	if (b > 0) {
-	    return rb_big_pow(rb_int2big(a), y);
-	}
-	return rb_float_new(pow((double)a, (double)b));
-    }
-    switch (TYPE(y)) {
-      case T_BIGNUM:
-	x = rb_int2big(FIX2LONG(x));
-	return rb_big_pow(x, y);
-      case T_FLOAT:
-	return rb_float_new(pow((double)FIX2LONG(x), RFLOAT(y)->value));
-      default:
-	return rb_num_coerce_bin(x, y);
-    }
+	return rb_num_pow(x,y);
 }
 
 /*
@@ -2831,13 +2814,14 @@
     return Qfalse;
 }
 
+
 /*
  *  call-seq:
  *     fix.odd? -> true or false
  *  
  *  Returns <code>true</code> if <i>fix</i> is an odd number.
  */
-
+ 
 static VALUE
 fix_odd_p(VALUE num)
 {
@@ -2846,7 +2830,7 @@
     }
     return Qfalse;
 }
-
+ 
 /*
  *  call-seq:
  *     fix.even? -> true or false
@@ -2863,6 +2847,7 @@
     return Qtrue;
 }
 
+
 void
 Init_Numeric(void)
 {
@@ -2976,8 +2961,9 @@
     rb_define_method(rb_cFixnum, "to_f", fix_to_f, 0);
     rb_define_method(rb_cFixnum, "size", fix_size, 0);
     rb_define_method(rb_cFixnum, "zero?", fix_zero_p, 0);
-    rb_define_method(rb_cInteger, "odd?", fix_odd_p, 0);
-    rb_define_method(rb_cInteger, "even?", fix_even_p, 0);
+    rb_define_method(rb_cFixnum, "odd?", fix_odd_p, 0);
+    rb_define_method(rb_cFixnum, "even?", fix_even_p, 0);
+
 
     rb_cFloat  = rb_define_class("Float", rb_cNumeric);
 
Index: bignum.c
===================================================================
RCS file: /src/ruby/bignum.c,v
retrieving revision 1.135
diff -u -r1.135 bignum.c
--- bignum.c	4 Sep 2006 20:10:45 -0000	1.135
+++ bignum.c	21 Sep 2006 14:46:25 -0000
@@ -38,6 +38,13 @@
 
 #define BIGZEROP(x) (RBIGNUM(x)->len == 0 || (RBIGNUM(x)->len == 1 && BDIGITS(x)[0] == 0))
 
+VALUE rb_num_plus(VALUE x,VALUE y);
+VALUE rb_num_minus(VALUE x,VALUE y);
+VALUE rb_num_mul(VALUE x,VALUE y);
+VALUE rb_num_sqr(VALUE x);
+VALUE rb_num_pow(VALUE x,VALUE y);
+
+
 static VALUE
 bignew_1(VALUE klass, long len, int sign)
 {
@@ -46,7 +53,6 @@
     big->sign = sign?1:0;
     big->len = len;
     big->digits = ALLOC_N(BDIGIT, len);
-
     return (VALUE)big;
 }
 
@@ -118,6 +124,36 @@
     return x;
 }
 
+/*
+ *  call-seq:
+ *     big.odd? -> true or false
+ *  
+ *  Returns <code>true</code> if <i>big</i> is an odd number.
+ */
+
+static VALUE
+big_odd_p(VALUE num)
+{   
+	return ( ((BDIGIT*)RBIGNUM(num)->digits)[0]%2!=0) ? Qtrue : Qfalse;
+}
+
+/*
+ *  call-seq:
+ *     int.even? -> true or false
+ *  
+ *  Returns <code>true</code> if <i>int</i> is an even number.
+ */
+
+static VALUE
+big_even_p(VALUE num)
+{
+    if (big_odd_p(num)==Qfalse) {
+	return Qtrue;
+    }
+    return Qfalse;
+}
+
+
 VALUE
 rb_big_norm(VALUE x)
 {
@@ -296,6 +332,33 @@
 
 #endif
 
+/*how much chars at base fits into BDIGIT*/
+void charperdig(int base,BDIGIT *hbase,long *charsperdig) {
+	*hbase = base;
+	*charsperdig=1;
+	/*if (base==10){//TODO implement as BDIGMAX aware table
+		*hbase=1000000000;
+		*charsperdig=9;
+		return;
+	}*/
+	while (*hbase<BDIGMAX/base){
+	  *hbase*=base;
+	  (*charsperdig)++;
+	}
+}
+
+static void bigdivmod(VALUE x, VALUE y, VALUE *divp, VALUE *modp);
+VALUE rb_big_sqr(VALUE x);
+static long big_log2(BDIGIT x){/*TODO table*/
+	int y=1,r=0;
+	while (x>y){
+		y*=2;
+		r++;
+	}
+	return r;
+}
+
+#define CSTR_TO_BIG_TRESHOLD 32
 VALUE
 rb_cstr_to_inum(const char *str, int base, int badcheck)
 {
@@ -305,7 +368,7 @@
     int c;
     BDIGIT_DBL num;
     long len, blen = 1;
-    long i;
+    long i,j,k;
     VALUE z;
     BDIGIT *zds;
 
@@ -354,30 +417,20 @@
     }
     switch (base) {
       case 2:
-	len = 1;
 	if (str[0] == '0' && (str[1] == 'b'||str[1] == 'B')) {
 	    str += 2;
 	}
 	break;
-      case 3:
-	len = 2;
-	break;
       case 8:
 	if (str[0] == '0' && (str[1] == 'o'||str[1] == 'O')) {
 	    str += 2;
 	}
-      case 4: case 5: case 6: case 7:
-	len = 3;
-	break;
       case 10:
 	if (str[0] == '0' && (str[1] == 'd'||str[1] == 'D')) {
 	    str += 2;
 	}
-      case 9: case 11: case 12: case 13: case 14: case 15:
-	len = 4;
 	break;
       case 16:
-	len = 4;
 	if (str[0] == '0' && (str[1] == 'x'||str[1] == 'X')) {
 	    str += 2;
 	}
@@ -386,22 +439,21 @@
 	if (base < 2 || 36 < base) {
 	    rb_raise(rb_eArgError, "illegal radix %d", base);
 	}
-	if (base <= 32) {
-	    len = 5;
-	}
-	else {
-	    len = 6;
-	}
 	break;
     }
     if (*str == '0') {		/* squeeze preceeding 0s */
 	while (*++str == '0');
 	--str;
     }
-    len *= strlen(str)*sizeof(char);
+	len=strlen(str);
+	BDIGIT hbase;
+	long charsperdig;
+	charperdig(base,&hbase,&charsperdig);
+	long si=len/charsperdig+1; 
+
 
-    if (len <= (sizeof(VALUE)*CHAR_BIT)) {
-	unsigned long val = strtoul(str, &end, base);
+    if (si<=sizeof(long)/sizeof(BDIGIT) ) {
+	long val = strtol(s, &end, base);
 
 	if (str < end && *end == '_') goto bigparse;
 	if (badcheck) {
@@ -409,67 +461,63 @@
 	    while (*end && ISSPACE(*end)) end++;
 	    if (*end) goto bad;	      /* trailing garbage */
 	}
-
-	if (POSFIXABLE(val)) {
-	    if (sign) return LONG2FIX(val);
-	    else {
-		long result = -(long)val;
-		return LONG2FIX(result);
-	    }
-	}
-	else {
-	    VALUE big = rb_uint2big(val);
-	    RBIGNUM(big)->sign = sign;
-	    return bignorm(big);
-	}
+	return LONG2NUM(val);
+	
     }
   bigparse:
-    len = (len/BITSPERDIG)+1;
     if (badcheck && *str == '_') goto bad;
-
-    z = bignew(len, sign);
+    z = bignew(si, sign);
     zds = BDIGITS(z);
-    for (i=len;i--;) zds[i]=0;
-    while (c = *str++) {
-	if (c == '_') {
-	    if (badcheck) {
-		if (nondigit) goto bad;
-		nondigit = c;
-	    }
-	    continue;
-	}
-	else if (!ISASCII(c)) {
-	    break;
-	}
-	else if (isdigit(c)) {
-	    c -= '0';
-	}
-	else if (islower(c)) {
-	    c -= 'a' - 10;
-	}
-	else if (isupper(c)) {
-	    c -= 'A' - 10;
-	}
-	else {
-	    break;
-	}
-	if (c >= base) break;
-	nondigit = 0;
-	i = 0;
-	num = c;
-	for (;;) {
-	    while (i<blen) {
-		num += (BDIGIT_DBL)zds[i]*base;
-		zds[i++] = BIGLO(num);
-		num = BIGDN(num);
-	    }
-	    if (num) {
-		blen++;
-		continue;
-	    }
-	    break;
+    for (i=si;i--;) zds[i]=0;
+	for(;;){
+		num=0;
+		j=charsperdig;
+		while (j--){			
+retry: 			c = *str++; 
+				if (c == '_') {
+					if (badcheck) {
+						if (nondigit) goto bad;
+						nondigit = c;
+					}
+					goto retry;
+				}
+				/*else if (!ISASCII(c)) {
+					goto end;
+				}*/
+				else if (isdigit(c)) {
+					c -= '0';
+				}
+				else if (islower(c)) {
+					c -= 'a' - 10;
+				}
+				else if (isupper(c)) {
+					c -= 'A' - 10;
+				}
+				else {
+					goto end;
+				}
+			if (c >= base) goto end;
+			nondigit = 0;
+			num=base*num+c;
+		}
+			
+		
+		if ((base&(base-1)/*base isn't power of 2*/)&& si<CSTR_TO_BIG_TRESHOLD){
+			i=0;
+			while (i<blen) {
+				num += (BDIGIT_DBL)zds[i]*hbase;
+				zds[i++] = BIGLO(num);
+				num = BIGDN(num);
+			}
+			if (num){
+				zds[blen++]=num;
+			}
+		}else{
+			zds[blen++]=num;
+		}
 	}
-    }
+
+			
     if (badcheck) {
 	str--;
 	if (s+1 < str && str[-1] == '_') goto bad;
@@ -479,6 +527,69 @@
 	    rb_invalid_str(s, "Integer");
 	}
     }
+end:
+	if  (!(base&(base-1))){/*its power of 2*/
+		BDIGIT_DBL num2=0;
+		int k=0,bits=0,usedbits;
+		VALUE x;
+		BDIGIT *xds;
+		usedbits=big_log2(hbase);
+		x = bignew(blen+1, sign);
+		xds = BDIGITS(x);
+		for (i=blen;i--;){		
+			num2|=((BDIGIT_DBL)zds[i])<<bits;	
+			bits+=usedbits;
+			if (bits>=BITSPERDIG){
+				bits-=BITSPERDIG;
+				xds[k++]=BIGLO(num2);
+				num2 = BIGDN(num2);
+			}
+		}
+		xds[k]=num2;
+		xds[k+1]=0;
+		RBIGNUM(x)->len=k+2;
+		z=x;
+		zds=xds;
+	}else if (si>=CSTR_TO_BIG_TRESHOLD){
+		blen--;
+		VALUE *temps=ALLOCA_N(VALUE,blen );
+		VALUE radix=ULONG2NUM(hbase);
+		for (i=0;i<blen;i++)
+			temps[i]=ULONG2NUM(zds[i+1]);
+		while (blen>1){	/*temps represent number written at hbase radix*/
+			if (blen%2){  /*we transform it into hbase**2 radix until it becomes 1 number*/
+				for (i=1;i<=blen/2;i++){
+					temps[i]=rb_num_plus(rb_num_mul(temps[2*i-1],radix),
+							temps[2*i]);	}				
+				blen=blen/2+1;
+			}else{
+				for (i=0;i<blen/2;i++){
+					temps[i]=rb_num_plus(rb_num_mul(temps[2*i],radix),
+							temps[2*i+1]);		}
+				blen/=2;
+			}
+			radix=rb_num_sqr(radix);		
+		}
+		z=temps[0];
+		blen=RBIGNUM(z)->len;
+		REALLOC_N(RBIGNUM(z)->digits, BDIGIT, ++RBIGNUM(z)->len);
+		BDIGITS(z)[RBIGNUM(z)->len-1]=0;
+		zds = BDIGITS(z);
+		RBIGNUM(z)->sign=sign;
+	}
+
+	hbase/=base;
+	while (j--)hbase/=base;
+		i=0;
+		while (i<blen) {
+			num += (BDIGIT_DBL)zds[i]*hbase;
+			zds[i++] = BIGLO(num);
+			num = BIGDN(num);
+		}
+		if (num){
+			zds[blen]=num;
+			blen++;
+		}
 
     return bignorm(z);
 }
@@ -577,13 +688,15 @@
     return rb_str_to_inum(str, base, base==0);
 }
 
+#define BIG_TO_STR_TRESHOLD LONG_MAX
 const char ruby_digitmap[] = "0123456789abcdefghijklmnopqrstuvwxyz";
 VALUE
 rb_big2str(VALUE x, int base)
 {
     volatile VALUE t;
     BDIGIT *ds;
-    long i, j, hbase;
+    long i, j; 
+	BDIGIT hbase;	long charsperdig;
     VALUE ss;
     char *s, c;
 
@@ -622,36 +735,78 @@
 	break;
     }
     j += 2;
-
-    hbase = base * base;
-#if SIZEOF_BDIGITS > 2
-    hbase *= hbase;
-#endif
-
+	charperdig( base,&hbase, &charsperdig);
+   
     t = rb_big_clone(x);
     ds = BDIGITS(t);
     ss = rb_str_new(0, j);
     s = RSTRING_PTR(ss);
 
     s[0] = RBIGNUM(x)->sign ? '+' : '-';
-    while (i && j) {
-	long k = i;
-	BDIGIT_DBL num = 0;
-
-	while (k--) {
-	    num = BIGUP(num) + ds[k];
-	    ds[k] = (BDIGIT)(num / hbase);
-	    num %= hbase;
-	}
-	if (ds[i-1] == 0) i--;
-	k = SIZEOF_BDIGITS;
-	while (k--) {
-	    c = (char)(num % base);
-	    s[--j] = ruby_digitmap[(int)c];
-	    num /= base;
-	    if (i == 0 && num == 0) break;
+	if (base&(base-1)/*not power of 2*/){
+		/*if (RBIGNUM(x)->len<BIG_TO_STR_TRESHOLD){*/
+			while (i && j) {
+				long k = i;
+				BDIGIT_DBL num = 0;
+
+				while (k--) {
+					num = BIGUP(num) + ds[k];
+					ds[k] = (BDIGIT)(num / hbase);
+					num %= hbase;
+				}
+				if (ds[i-1] == 0) i--;
+				k = charsperdig;
+				while (k--) {
+					c = (char)(num % base);
+					s[--j] = ruby_digitmap[(int)c];
+					num /= base;
+					if (i == 0 && num == 0) break;
+				}
+			}
+		/*}else{
+			VALUE powers[40];
+			int len=RBIGNUM(x)->len,po=0;
+			VALUE pow=ULONG2NUM(hbase);
+			while (rb_big_cmp(pow,x)!=INT2FIX(-1)){
+				powers[po++]=pow;
+				pow=rb_big_sqr(pow);
+			}
+			VALUE *split=ALLOCA_N(VALUE,4*len),*splitn=ALLOCA_N(VALUE,4*len),*tmp;
+			split[0]=x;
+			len=1;
+			while (po--)
+			{
+				for (i=0;i<len;i++)
+					divmod(split[i],powers(po),splitn+2*i+1,splitn+2*i);
+				len=2*len;
+				if BIGZEROP(splitn[len-1])	len--;
+				tmp=split;split=splitn;splitn=tmp;
+			}
+			for (i=0;i<len;i++)
+			{
+				k = charsperdig;
+				num=split[i];
+				while (k--) {
+					c = (char)(num % base);
+					s[--j] = ruby_digitmap[(int)c];
+					num /= base;
+				}
+			}
+		}*/
+	}else{
+		BDIGIT_DBL num2=0;
+		int bits=0;
+		int bitsperchar=big_log2(base);
+		for (i=0;i<RBIGNUM(x)->len;i++){		
+			num2|=((BDIGIT_DBL)BDIGITS(x)[i])<<bits;	
+			bits+=BITSPERDIG;
+			while (bits>=bitsperchar){
+				bits-=bitsperchar;
+			    s[--j] = ruby_digitmap[(int)(num2&(base-1))];
+				num2 >>=bitsperchar ;
+			}
+		}
 	}
-    }
     while (s[j] == '0') j++;
     i = RSTRING_LEN(ss)-(RBIGNUM(x)->sign?j:j-1);
     memmove(RBIGNUM(x)->sign?s:s+1, s+j, i);
@@ -1010,6 +1165,33 @@
     return bignorm(z);
 }
 
+
+static void
+big_minus_bang(VALUE x,VALUE y)
+{  
+	BDIGIT *xx,*yy;
+	if (FIXNUM_P(y)) y=rb_int2big(FIX2LONG(y));
+	xx=BDIGITS(x);
+	yy=BDIGITS(y);
+    BDIGIT_DBL_SIGNED num;
+    long i;
+
+
+    for (i = 0, num = 0; i < RBIGNUM(y)->len; i++) { 
+	num += (BDIGIT_DBL_SIGNED) BDIGITS(x)[i] - BDIGITS(y)[i];
+	BDIGITS(x)[i] = BIGLO(num);
+	num = BIGDN(num);
+    } 
+    while (num ) {
+	num += BDIGITS(x)[i];
+	BDIGITS(x)[i++] = BIGLO(num);
+	num = BIGDN(num);
+    }
+	if (i>RBIGNUM(x)->len){
+		rb_raise(rb_eFatal, "substracted greater number");
+	}
+}
+
 static VALUE
 bigsub(VALUE x, VALUE y)
 {
@@ -1056,6 +1238,30 @@
     return z;
 }
 
+static void 
+big_plus_bang(VALUE x, VALUE y, int shift)
+{
+    BDIGIT_DBL num;
+    long i, len;
+	if (FIXNUM_P(y)) y=rb_int2big(FIX2LONG(y));
+
+
+	len = RBIGNUM(y)->len;
+    for (i = 0, num = 0; i < len; i++) {
+	num += (BDIGIT_DBL)BDIGITS(x)[i+shift] + BDIGITS(y)[i];
+	BDIGITS(x)[i+shift] = BIGLO(num);
+	num = BIGDN(num);
+    }
+    while (num ) {
+	num += BDIGITS(x)[i+shift];
+	BDIGITS(x)[i+shift] = BIGLO(num);
+	i++;
+	num = BIGDN(num);
+    }	
+	if (i>RBIGNUM(x)->len)
+		rb_raise(rb_eFatal, "added beyond capacity");
+}
+
 static VALUE
 bigadd(VALUE x, VALUE y, int sign)
 {
@@ -1108,7 +1314,7 @@
 
 VALUE
 rb_big_plus(VALUE x, VALUE y)
-{
+{	
     switch (TYPE(y)) {
       case T_FIXNUM:
 	y = rb_int2big(FIX2LONG(y));
@@ -1149,6 +1355,19 @@
     }
 }
 
+static void big_set(VALUE x,VALUE y,int shift)
+{
+		if FIXNUM_P(y)y=rb_int2big(FIX2LONG(y));
+		MEMCPY(BDIGITS(x)+shift, BDIGITS(y), BDIGIT, RBIGNUM(y)->len);
+}
+
+static VALUE big_part(VALUE x,int beg,int end)
+{
+	VALUE z=bignew(end-beg,1);
+	MEMCPY(BDIGITS(z), BDIGITS(x)+beg, BDIGIT, end-beg);
+	return z;
+}
+
 static VALUE
 rb_big_mul0(VALUE x, VALUE y)
 {
@@ -1171,17 +1390,50 @@
       default:
 	return rb_num_coerce_bin(x, y);
     }
-
-    j = RBIGNUM(x)->len + RBIGNUM(y)->len + 1;
+	
+	if (RBIGNUM(x)->len < RBIGNUM(y)->len)	{z=x;x=y;y=z;}
+    j = RBIGNUM(x)->len + RBIGNUM(y)->len;
     z = bignew(j, RBIGNUM(x)->sign==RBIGNUM(y)->sign);
+	if (RBIGNUM(y)->len>=64){/*TODO find what at what value is best*/		
+		/*Karatsuba multiplication
+		 * Its based at folowing
+		 * x*y=(xhi<<shift+xlo)*(yhi<<shift+ylo)=
+		 *    =xhi*yhi<<(2*shift)+(xhi*ylo+yhi*xlo)<<shift+xlo*ylo
+		 *    =xhi*yhi<<(2*shift)+(zmid)<<shift+xlo*ylo
+		 * where zmid=xhi*ylo+yhi*ylo= (xhi+xlo)*(yhi+ylo)-xhi*yhi-xlo*ylo
+		 * */
+		memset(RBIGNUM(z)->digits,0,sizeof(BDIGIT)*RBIGNUM(z)->len);
+		int shift=RBIGNUM(x)->len/2;
+		VALUE xhi=big_part(x,shift,RBIGNUM(x)->len),xlo=big_part(x,0,shift);
+		VALUE zhi,zlo;
+		if (RBIGNUM(y)->len<=shift) {
+			zlo=rb_big_mul(xlo,y);
+			zhi=rb_big_mul(xhi,y);
+			big_set(z,zlo,0);
+			big_plus_bang(z,zhi,shift);
+			return z;
+		}
+		VALUE yhi=big_part(y,shift,RBIGNUM(y)->len),ylo=big_part(y,0,shift);
+		zhi=rb_big_mul(xhi,yhi),zlo=rb_big_mul(xlo,ylo);	
+		VALUE zmid=rb_num_mul(rb_big_plus(xhi,xlo),rb_big_plus(yhi,ylo));
+		if (FIXNUM_P(zmid)) zmid=rb_int2big(FIX2LONG(zmid));
+		big_minus_bang(zmid,zhi);
+		big_minus_bang(zmid,zlo);
+		/*zmid=(xhi+xlo)*(yhi+ylo)-xhi*yhi-xlo*ylo=xhi*ylo+yhi*ylo */		
+
+		big_set(z,zlo,0);
+		big_set(z,zhi,2*shift);/*zlo and zhi dont overlap*/
+		big_plus_bang(z,zmid,shift);
+		return z;
+	}
     zds = BDIGITS(z);
     while (j--) zds[j] = 0;
-    for (i = 0; i < RBIGNUM(x)->len; i++) {
-	BDIGIT_DBL dd = BDIGITS(x)[i]; 
+    for (i = 0; i < RBIGNUM(y)->len; i++) {
+	BDIGIT_DBL dd = BDIGITS(y)[i]; 
 	if (dd == 0) continue;
 	n = 0;
-	for (j = 0; j < RBIGNUM(y)->len; j++) {
-	    BDIGIT_DBL ee = n + (BDIGIT_DBL)dd * BDIGITS(y)[j];
+	for (j = 0; j < RBIGNUM(x)->len; j++) {
+	    BDIGIT_DBL ee = n + (BDIGIT_DBL)dd * BDIGITS(x)[j];
 	    n = zds[i + j] + ee;
 	    if (ee) zds[i + j] = BIGLO(n);
 	    n = BIGDN(n);
@@ -1193,6 +1445,26 @@
     return z;
 }
 
+VALUE rb_big_sqr(VALUE x)/*karatsuba sqr is faster than sf x*x */
+{
+    if (RBIGNUM(x)->len>16){
+	  int shift=RBIGNUM(x)->len/2;
+	  VALUE  z = bignew(RBIGNUM(x)->len*2, 1);
+ 	  memset(RBIGNUM(z)->digits,0,sizeof(BDIGIT)*RBIGNUM(z)->len);
+	
+	  VALUE xhi=big_part(x,shift,RBIGNUM(x)->len),xlo=big_part(x,0,shift);
+	  /* z= xhi**2 <<2*shift+2*xhi*xlo<<shift+xlo**2 */
+      VALUE zlo=rb_big_sqr(xlo),zhi=rb_big_sqr(xhi),zmid=rb_big_mul(xlo,xhi);
+	  big_set(z,zlo,0);/*zlo and zhi dont overlap*/
+	  big_set(z,zhi,2*shift);
+	  big_plus_bang(z,zmid,shift);
+  	  big_plus_bang(z,zmid,shift);
+      return z;
+	}
+    return rb_big_mul(x,x);
+}
+
+
 /*
  *  call-seq:
  *     big * other  => Numeric
@@ -1216,7 +1488,8 @@
     BDIGIT_DBL t2;
     BDIGIT_DBL_SIGNED num;
     BDIGIT dd, q;
-
+	
+   	if (FIXNUM_P(x)) x=rb_int2big(FIX2LONG(x));
     if (BIGZEROP(y)) rb_num_zerodiv();
     yds = BDIGITS(y);
     if (nx < ny || (nx == ny && BDIGITS(x)[nx - 1] < BDIGITS(y)[ny - 1])) {
@@ -1339,6 +1612,7 @@
 {
     VALUE mod;
 
+   	if (FIXNUM_P(x)) x=rb_int2big(FIX2LONG(x));
     bigdivrem(x, y, divp, &mod);
     if (RBIGNUM(x)->sign != RBIGNUM(y)->sign && !BIGZEROP(mod)) {
 	if (divp) *divp = bigadd(*divp, rb_int2big(1), 0);
@@ -1507,6 +1781,7 @@
     return rb_float_new(dx / dy);
 }
 
+
 /*
  *  call-seq:
  *     big ** exponent   #=> numeric
@@ -1519,54 +1794,9 @@
  *    123456789 ** 1.2    #=> 5126464716.09932
  *    123456789 ** -2     #=> 6.5610001194102e-17
  */
-
-VALUE
-rb_big_pow(VALUE x, VALUE y)
-{
-    double d;
-    long yy;
-    
-    if (y == INT2FIX(0)) return INT2FIX(1);
-    switch (TYPE(y)) {
-      case T_FLOAT:
-	d = RFLOAT(y)->value;
-	break;
-
-      case T_BIGNUM:
-	rb_warn("in a**b, b may be too big");
-	d = rb_big2dbl(y);
-	break;
-
-      case T_FIXNUM:
-	yy = FIX2LONG(y);
-	if (yy > 0) {
-	    VALUE z = x;
-
-	    if (RBIGNUM(x)->len * SIZEOF_BDIGITS * yy > 1024*1024) {
-		rb_warn("in a**b, b may be too big");
-		d = (double)yy;
-		break;
-	    }
-	    for (;;) {
-		yy -= 1;
-		if (yy == 0) break;
-		while (yy % 2 == 0) {
-		    yy /= 2;
-		    x = rb_big_mul0(x, x);
-		    if (!BDIGITS(x)[RBIGNUM(x)->len-1]) RBIGNUM(x)->len--;
-		}
-		z = rb_big_mul0(z, x);
-		if (!BDIGITS(z)[RBIGNUM(z)->len-1]) RBIGNUM(z)->len--;
-	    }
-	    return bignorm(z);
-	}
-	d = (double)yy;
-	break;
-
-      default:
-	return rb_num_coerce_bin(x, y);
-    }
-    return rb_float_new(pow(rb_big2dbl(x), d));
+VALUE 
+rb_big_pow(VALUE x,VALUE y){
+	return rb_num_pow(x,y);
 }
 
 /*
@@ -1584,6 +1814,7 @@
     long i, l1, l2;
     char sign;
 
+   	if (FIXNUM_P(xx)) xx=rb_int2big(FIX2LONG(xx));
     x = xx;
     y = rb_to_int(yy);
     if (FIXNUM_P(y)) {
@@ -1638,7 +1869,7 @@
     BDIGIT *ds1, *ds2, *zds;
     long i, l1, l2;
     char sign;
-
+   	if (FIXNUM_P(xx)) xx=rb_int2big(FIX2LONG(xx));
     x = xx;
     y = rb_to_int(yy);
     if (FIXNUM_P(y)) {
@@ -1691,7 +1922,8 @@
 VALUE
 rb_big_xor(VALUE xx, VALUE yy)
 {
-    volatile VALUE x, y;
+   	if (FIXNUM_P(xx)) xx=rb_int2big(FIX2LONG(xx));
+volatile VALUE x, y;
     VALUE z;
     BDIGIT *ds1, *ds2, *zds;
     long i, l1, l2;
@@ -1753,6 +1985,7 @@
 VALUE
 rb_big_lshift(VALUE x, VALUE y)
 {
+
     BDIGIT *xds, *zds;
     int shift = NUM2INT(y);
     int s1 = shift/BITSPERDIG;
@@ -1761,6 +1994,7 @@
     BDIGIT_DBL num = 0;
     long len, i;
 
+	if FIXNUM_P(x)x=rb_int2big(FIX2LONG(x));
     if (shift < 0) return rb_big_rshift(x, INT2FIX(-shift));
     len = RBIGNUM(x)->len;
     z = bignew(len+s1+1, RBIGNUM(x)->sign);
@@ -1796,7 +2030,7 @@
     BDIGIT_DBL num = 0;
     long i, j;
     volatile VALUE save_x;
-
+	if FIXNUM_P(x)x=rb_int2big(FIX2LONG(x));
     if (shift < 0) return rb_big_lshift(x, INT2FIX(-shift));
 
     if (s1 > RBIGNUM(x)->len) {
@@ -1827,6 +2061,92 @@
     return bignorm(z);
 }
 
+VALUE 
+rb_num_plus(VALUE x,VALUE y)
+{
+	if FIXNUM_P(x) 
+		return fix_plus(x,y);
+	else 
+		return rb_big_plus(x,y);
+}
+
+VALUE 
+rb_num_minus(VALUE x,VALUE y)
+{
+	if FIXNUM_P(x) 
+		return fix_minus(x,y);
+	else 
+		return rb_big_minus(x,y);
+}
+
+VALUE 
+rb_num_mul(VALUE x,VALUE y)
+{
+	if FIXNUM_P(x) 
+		return fix_mul(x,y);
+	else 
+		return rb_big_mul(x,y);
+}
+
+VALUE 
+rb_num_sqr(VALUE x)
+{
+	if FIXNUM_P(x) 
+		return fix_mul(x,x);
+	else
+		return rb_big_sqr(x);
+}
+
+VALUE 
+rb_num_pow(VALUE x,VALUE y)
+{
+    double d;
+    long yy;
+    
+    if (y == INT2FIX(0)) return INT2FIX(1);
+    switch (TYPE(y)) {
+      case T_FLOAT:
+	d = RFLOAT(y)->value;
+	break;
+
+      case T_BIGNUM:
+	rb_warn("in a**b, b may be too big");
+	d = rb_big2dbl(y);
+	break;
+
+      case T_FIXNUM:
+	yy = FIX2LONG(y);
+	if (yy > 0) {
+	    VALUE z = x;
+
+	    if (!FIXNUM_P(x)&&RBIGNUM(x)->len * SIZEOF_BDIGITS * yy > 1024*1024) {
+		rb_warn("in a**b, b may be too big");
+		d = (double)yy;
+		break;
+	    }
+	    for (;;) {
+		yy -= 1;
+		if (yy == 0) break;
+		while (yy % 2 == 0) {
+		    yy /= 2;
+		    x = rb_num_sqr(x);
+		}
+		z = rb_num_mul(z, x);
+	    }
+	    return z;
+	}
+	d = (double)yy;
+	break;
+
+      default:
+	return rb_num_coerce_bin(x, y);
+    }
+    return rb_float_new(pow(rb_big2dbl(x), d));
+}
+
+
+
+
 /*
  *  call-seq:
  *     big[n] -> 0, 1
@@ -2007,4 +2327,8 @@
     rb_define_method(rb_cBignum, "to_f", rb_big_to_f, 0);
     rb_define_method(rb_cBignum, "abs", rb_big_abs, 0);
     rb_define_method(rb_cBignum, "size", rb_big_size, 0);
+    rb_define_method(rb_cBignum, "odd?", big_odd_p, 0);
+    rb_define_method(rb_cBignum, "even?", big_even_p, 0);
+
+
 }
Index: sample/test.rb
===================================================================
RCS file: /src/ruby/sample/test.rb,v
retrieving revision 1.101
diff -u -r1.101 test.rb
--- sample/test.rb	10 Jul 2006 08:36:43 -0000	1.101
+++ sample/test.rb	21 Sep 2006 14:46:32 -0000
@@ -1572,7 +1572,11 @@
 test_ok(1299022.to_s(36) == "ruby")
 test_ok(-1045307475.to_s(36) == "-hacker")
 test_ok("Just_another_Ruby_hacker".to_i(36) == 265419172580680477752431643787347)
+s="thisisaverylargestringthisisaverylargestringthisisaverylargestringthisisaverylargestringthisisaverylargestringthisisaverylargestringsisaverylargestringthisisaverylargestring"
+test_ok(s.to_i(36).to_s(36)==s)
 test_ok(-265419172580680477752431643787347.to_s(36) == "-justanotherrubyhacker")
+s="123456789abcdef"
+test_ok(s==s.to_i(16).to_s(16))
 
 a = []
 (0..255).each {|n|