カラツバ法

CourseraのAlgorithm Part1で出てきてふと実装熱が再燃した。

カラツバ法はかなり面白い乗算のアルゴリズム(と個人的に思っている)で、最初見たときはいきなり意味不明なことしてるぞ、とか思ったりした。アルゴリズムの肝は以下のような感じ。

;; multiplicands, x y
(define x #x123456789)
(define y #x908765432)

(define B 16) ;; radix
(define n 9)  ;; n digits

#|
x = a*B^(n/2) + b
y = c*B^(n/2) + d
|#
(define a #x12345) ;; most significant part of x
(define b #x6789)
(define c #x90876) ;; most significant part of y
(define d #x5432)

#|
x*y = (a*B^(n/2) + b) * (c*B^(n/2) + d)
    = B^n*ac + B^(n/2)*(ad + bc) + bd
    = 16^n*ac + 16^(n/2)*(ad + bc) + bd
|#
(let ((ac (* a c))
      (bd (* b d))
      (xx (+ (* a d) (* b c))))
  (+ (* (expt B (* (div n 2) 2)) ac) 
     (* (expt B (div n 2)) xx)
     bd))

以前実装を諦めた理由を覚えていないのだが、なんとなく実装できたのでベンチマークをとってみた。ベンチマークには以下のコードとスクリプトを仕様。

(import (time) (sagittarius control))
(define ex 2048)
(define ey 2048)
(define x (expt 2 ex))
(define y (expt 2 ey))
(define count 50000)

(time (dotimes (i count) (* x y)))

#!/bin/sh
echo Non karatsuba
time sash multi.scm
echo Karatsuba
time ./build/sash multi.scm

何のことは単なるスクリプト。シェルでtimeコマンド使う必要は無かったかもしれない。っで、以下が結果。

$ ./time.sh
Non karatsuba

;;  (dotimes (i count) (* x y))
;;  2.238128 real    2.215000 user    0.000000 sys

real    0m2.771s
user    0m2.449s
sys     0m0.296s
Karatsuba

;;  (dotimes (i count) (* x y))
;;  1.882108 real    1.872000 user    0.000000 sys

real    0m2.387s
user    0m2.090s
sys     0m0.296s

400msほど高速に計算する。まぁ、5万回まわしてこの程度とも言える。

カラツバ法のアルゴリズムだからなのか実装が悪いのか分からないが、乗算の数値のサイズが違いすぎると逆に遅くなる。適当に計測した結果、Sagittariusの実装だとおよそ10ワードが境界となった。なのでそれ以上の差がある場合は従来の乗算を使う。(この辺例えばToom-3とか実装すると解決するのかも、とか思いながら腰は重い。)