自作LispにKaratsuba乗算を組み込んだ話

改心

昨今のロシアによるウクライナ軍事進行に対する抗議の気持ちでKaratsuba乗算を取り入れないつもりでいました。しかし、先日あるロシア人の人から自作Prolog処理系に改良のプルリクエストをいただきました。macOSでも動くようにしていただきました。国としてのロシアには大いに問題がありますが、それを個々のロシア人にまで拡張解釈するのは良くないと反省したのでした。ここにKaratsubaさんの業績に敬意を表しつつ自作LISPにKaratsuba乗算を組み込みました。

アルゴリズム

アルゴリズムについては他所に詳細な記事がありますので割愛させていただきます。ロシア人は数学のような抽象的概念に優れた人の割合が多いようにに思います。数学者も多数排出しています。念のためKaratsubaさんのざっくりとしたアイディアは次のようになっています。

* KARATSUBA-Multiply 
*    X = x1*b + x0
*    Y = y1*b + y0
*    z2 = x1 * y1
*    z0 = x0 * y0
*    z1 = z2 + z0 − (x1 − x0)*(y1 − y0)
*    Z = z2*b^2 + z1*b + z0  
*

コード

次のようにアルゴリズムをコードに書き下しました。

int 
bigx_first_half(int x)
{
  int pointer,len,res;

  pointer = get_pointer(x);
  len = get_length(x);

  len = len / 2;
  res = gen_big();
  set_pointer(res,pointer);
  set_length(res,len);
  set_sign(res,get_sign(x));
  return(res);
}

int 
bigx_second_half(int x)
{
  int pointer,len,res;

  pointer = get_pointer(x);
  len = get_length(x);

  len = len / 2;
  pointer = pointer - len;
  res = gen_big();
  set_pointer(res,pointer);
  set_length(res,len);
  set_sign(res,get_sign(x));
  return(res);
}


int 
bigx_karatsuba_mult1(int x, int y)
{
  int len,x0,y0,x1,y1,z0,z1,z2,z;

  len = get_length(x);
  if(len < 10)
    return(bigx_mult(x,y));
  else{
    x1 = bigx_first_half(x);
    x0 = bigx_second_half(x);
    y1 = bigx_first_half(y);
    y0 = bigx_second_half(y);
    z2 = bigx_karatsuba_mult1(x1,y1);
    z0 = bigx_karatsuba_mult1(x0,y0);
    //z1 := z2 + z0 - (x1 - x0)*(y1 - y0) 
    z1 = minus(plus(z2,z0),mult(minus(x1,x0),minus(y1,y0)));
    //Z = z2*b^2 + z1*b + z0  
    z = plus(plus(bigx_shift_right(z2,len),bigx_shift_right(z1,len/2)),z0);
    return(z);
  }
}

int
bigx_karatsuba_mult(int x, int y)
{
  int lenx, leny, max_len, n, x1, y1, res;

  lenx = get_length (x);
  leny = get_length (y);
  if (lenx >= leny)
    {
      max_len = lenx;
    }
  else
    {
      max_len = leny;
    }

  n = 1;
  while (max_len > n)
    {
      n = 2 * n;
    }
  
  x1 = bigx_zero_supress(x,n-lenx);
  y1 = bigx_zero_supress(y,n-leny);
  
  // karatuba recursion
  res = bigx_karatsuba_mult1(x1,y1);
  return (res);
}

テスト

eisl-test という関数名にKaratsuba乗算を組み込んで試してみました。

> (eisl-test (expt 9 100) (expt 8 99))
67633156394733110901380138312618943356149760456325646953903999913949464703379968840706729314566411649098297553062918226864848164763933222398163799762670647441457568323756117490686492672

老舗のsbclと答え合わせをします。OKですね。

* (* (expt 9 100) (expt 8 99))
67633156394733110901380138312618943356149760456325646953903999913949464703379968840706729314566411649098297553062918226864848164763933222398163799762670647441457568323756117490686492672

通常の筆算アルゴリズムと比較してどのくらい高速になっているのでしょう？ 自作Lispには数論的FFTがくみこんでありますので、これをOFFにして比較してみました。

(defglobal a (expt 8 70000))
(defglobal b (expt 9 70000))


(defun bigtest0 (n)
    (if (= n 0)
        t
        (progn (eisl-test a b) (bigtest0 (- n 1)))))


(defun bigtest1 (n)
    (if (= n 0)
        t
        (progn (* a b) (bigtest1 (- n 1)))))

Easy-ISLisp Ver2.53
> (time (bigtest0 10))
Elapsed Time(second)=2.406790
<undef>
> (time (bigtest1 10))
Elapsed Time(second)=3.798342
<undef>
> 

計算時間は63％ほどに短縮されていました。

試験的組み込み

乗算のルーチンにKaratsuba乗算も組み込みました。乗数、被乗数のセル数の和が１００（およそ１０００桁）を超えたらNTT乗算を起動、２０（およそ２００桁）超えたらkaratsuba乗算を起動するようにしています。私の実装ではKaratsubaは高速になる代わりにメモリを食います。数十桁程度の乗算なら筆算アルゴリズムの方が良さそうです。大きくなるとNTT乗算が現在のところ最高速にようです。それにしてもかなり初期の段階でこの算法を見出したkaratsubaさんは才能に恵まれたお方だったのでしょう。

他にもコルモロゴフさんとか、ジューコフスキーさんとか優れた業績を残したロシア人がいました。彼らの業績に泥をかけないためにもロシア政府は直ちに軍事進行、力による政治をやめるべきですね。少数の愚か者のために民族全体が愚者と見られるのは耐えられないことでしょう。

ソースコード

自作Lisp処理系　Easy-ISlispはOSSとして下記にて公開しております。