[ruby-dev:11363] Re: bignum

田中です

>From: [e-mail] (Yukihiro Matsumoto)
>Subject: [ruby-dev:11360] Re: bignum

> Knuthの本ではなんだか再帰を使ってるようなんですが、なんといっ
> ても例題がMIX(要するにアセンブラ)で書いてあるんで、良く分か
> んなかったのでした。新版を英語で読んだ方が良いのかも。って新
> 版でもMIXなのかしら。

GMPを参考にしてCで書いたものをご参考までに添付します。
add() とかを作らないと動かないので、動作確認はしていません。
alloca を使うかどうかは要検討。GMPでは workspace を
最初に割り当てて、未使用領域を下請けに渡しているようです。

> ところでKaratsubaってなんとなく日本的な響きですがロシア(当時
> はソ連か)の人みたいですね。論文は1962年に出てます。

私もてっきり日本人かと思いました。

田中昌宏


#define cpy(a,b,n) memcpy(a,b,sizeof(BDIGIT)*(n))
#define alloc(n)   (BDIGIT*)alloca(sizeof(BDIGIT)*(n))
#define THRESH 8  /* or so */

int   gt(BDIGIT *x, int nx, BDIGIT *y, int ny);  /* x > y */
int  add(BDIGIT *x, int nx, BDIGIT *y, int ny);  /* x += y */
int  sub(BDIGIT *x, int nx, BDIGIT *y, int ny);  /* x -= y */
/* r[nx+ny] = x[nx] * y[ny] */
void mul_base(BDIGIT *r, BDIGIT *x, int nx, BDIGIT *y, int ny);


/* r[2*nm] = x[nm] * y[nm] */
void mul(BDIGIT *r, BDIGIT *x, BDIGIT *y, int nm)
{
    BDIGIT *u, *v, *w, *z;
    int n, m, n2, flag;

    if (m<THRESH)
        mul_base(r,x,nm,y,nm);
  
    /* nm == n + m,  n==m or n==m+1 */
    m = nm/2;
    n = nm-m;
    n2 = n*2;
    /*
    x1 = x[0..n-1], x2 = x[n..nm-1]
    y1 = y[0..n-1], y2 = y[n..nm-1]
    u = r[0..n2-1], v = r[n2..2nm-1]
    */
    u = r;
    v = r+n2;

    if ( gt( x,n, x+n,m ) ) {  /* x1 > x2 */
        cpy(u, x, n);
        sub(u, n, x+n, m);     /* x1-x2 */
	flag = 0;
    } else {
        cpy(u, x+n, m); u[m]=0;
        sub(u, m, x, n);       /* x2-x1 */
        flag = 1;
    }

    if ( gt( y,n, y+n,m ) ) {  /* y1 > y2 */
        cpy(v, y, n);
        sub(v, n, y+n, m);     /* y1-y2 */
    } else {
        cpy(v, y+n, m); v[m]=0;
        sub(v, m, y, n);       /* y2-y1 */
        flag ^= 1;
    }

    w = alloc(n2);
    mul(w, u, v, n);      /* w = +/- (x1-x2)*(y1-y2) */
    mul(u, x,   y,   n);  /* u = r[0..n2-1]   = x1*y1 */
    mul(v, x+n, y+n, m);  /* v = r[n2..2nm-1] = x2*y2 */

    z = alloc(n2+1);
    cpy(z, u, n2); z[n2]=0;
    add(z, n2+1, v, n2);      /* z = u+v */

    if (flag)
        add(z, n2+1, w, n2);  /* z = u+v+w */
    else
        sub(z, n2+1, w, n2);  /* z = u+v-w */

    add(r+n, n+2*m, z, n2+1); /* r[n..2nm-1] += u+v-w */
}