でざわです

bignum に sqrt 作ってみました。

で、テストしていて気づいたのですが bignum が 32bit を切ると
Fixnum に変わってしまうんですね。
c=200000000000000; c.sqrt は通るが (c/100000000000).sqrt は
	 undefined local variable or method `sqrt' for 2000(Fixnum)
はて、
すべての Numneric に作る必要がでてきた、、。
と言うことで、Fixnum と Float にも作りました。

Known BUG：
	負の Fixnum を sqrt しても エラーにならない。

次の warning が出ます。

bignum.c:1260: warning: assignment makes pointer from integer without a cast
bignum.c:1263: warning: assignment makes pointer from integer without a cast

bignum な x を x>>1 するにはどうするのが正しい作法？

	    x = big_rshift(x, INT2FIX(1) );

--- bignum.c.orig	Sat Jan 25 09:47:03 1997
+++ bignum.c	Sat Jan 25 18:11:01 1997
@@ -1,3 +1,4 @@
+#include <stdio.h>
 /************************************************
 
   bignum.c -
@@ -1215,6 +1216,55 @@
     return INT2FIX(big->len*2);
 }
 
+static VALUE
+big_sqrt(x)
+    struct RBignum *x;
+{
+    struct RBignum *z;
+    USHORT *zds, *xds;
+    unsigned long y, a;
+    UINT len, sqrt_len, iz, ix , j ;
+    UINT mask;
+  
+  /* 符号確認*/
+    if (x->sign == 0)ArgError("square root for negative number");
+
+  /* <= 2^32 */
+    /* あり得ない */
+  /* 正規化 2^30 <=  < 2^32 にする*/
+    xds = BDIGITS(x);
+    iz = sqrt_len;  ix = x->len;
+    y = ((long)xds[x->len-1])<<16 ;
+    j=16;
+    mask = ~0177777 ; /* 0011 1111 1111 1111  or 0111 1111 1111 1111 */
+    while ( y & mask ){
+      j -= 2; y >>=2;      /* j bits 左shift した */
+    }
+    y = (y << 16) + (((long)xds[x->len-2])<<j)
+             + (( (xds[x->len-3])<<j & 0377 )>>16 );
+
+    /* 初期値 */
+    a = 143257 - 202860544/(1582+( y >>22));
+    a = a << ( 16 - (j>>1)); 
+    sqrt_len = (( x->len+1) >>1) ;
+    z = (struct RBignum*)bignew( sqrt_len, 1);
+    zds = BDIGITS(z);
+    zds[sqrt_len-1] = a>>16; zds[sqrt_len - 2 ] = a & 0177777;
+   
+    /* 精度計算 */
+    iz = (x->len *16 - j)>>1 ;
+    
+    /* Newton (until j * 2^i > iz ) */
+    /* X(i+1) = (X(i)^2 + Y)/2X(i) --> = X(i)>>1 + (Y>>1)/X(i)  */
+    /* x = x>>1 */
+    x = big_rshift(x, INT2FIX(1) );
+    for(j=15; j < iz; j = j<<1){
+      /* z = ( z >> 1 )             + ( x >> 1 )/z */
+      z = big_plus(big_rshift(z, INT2FIX(1)), big_div( x,z) );
+    }
+    return bignorm(z);
+}
+
 void
 Init_Bignum()
 {
@@ -1246,4 +1296,6 @@
     rb_define_method(cBignum, "to_f", big_to_f, 0);
     rb_define_method(cBignum, "abs", big_abs, 0);
     rb_define_method(cBignum, "size", big_size, 0);
+    rb_define_method(cBignum, "sqrt", big_sqrt, 0);
+
 }
*** numeric.c.orig      Sat Jan 25 09:44:25 1997
--- numeric.c   Sat Jan 25 18:41:49 1997
***************
*** 1040,1045 ****
--- 1040,1065 ----
      }
      return num;
  }
+ static VALUE
+ flo_sqrt(flt)
+     struct RFloat *flt;
+ {
+     double val = flt->value;
+     if (val < 0.0) ArgError("square root for negative number");
+     return float_new(sqrt(val));
+ }
+ static VALUE
+ fix_sqrt(fix)
+     VALUE fix;
+ {
+   double dbl;
+   fix = FIX2INT(fix) ;
+   if ( fix < 0 ) ArgError("square root for negative number");
+   dbl = fix;
+   fix = sqrt(dbl)+0.5 ;
+   fix =INT2FIX(fix) ;
+   return fix ;
+ }
  
  extern VALUE mComparable;
  extern VALUE eException;
***************
*** 1138,1141 ****
--- 1158,1163 ----
      rb_define_method(cFloat, "to_i", flo_to_i, 0);
      rb_define_method(cFloat, "to_f", flo_to_f, 0);
      rb_define_method(cFloat, "abs", flo_abs, 0);
+     rb_define_method(cFloat, "sqrt", flo_sqrt, 0);
+     rb_define_method(cFixnum, "sqrt", fix_sqrt, 0);
  }