はじめに
この度、自作のEasy-ISlispのbignum（多倍長整数）のデータ構造をすっかり書き直し効率を向上させました。そのあたりのことについて書き残しておこうと思います。
車輪の再発明
多倍長整数計算については既にGMPがあります。SBCLなど有名処理系はこれを採用しています。これはとても良くできていて乗算にカラツバ、FFTといった効率的なアルゴリズムが採用されています。十分に高品質なライブラリがあるならそれを利用するのが職業プラグラマの常道です。しかし、そこはアマチュアのLisperたるものこのくらい自分で作らないとね。
データ構造
10^10進法によっています。32bit intの配列を用意してこれをbignumに充てました。各配列の要素に０～10^10^1の非負整数が格納されます。ポインターを２つ用意しました。big_pt0は一時的利用の領域を管理します。big_pt1は恒久的利用の領域を管理します。通常の計算は全社の領域で行います。この時、除算ではその計算過程で使い捨ての乗算、減算、加算が使われます。この計算をするときには一時的にbig_pt0を作業領域に切り替えて計算をします。ここは使い捨ての領域であり再利用されます。big_pt1の管理する領域は例えば大域変数に束縛されたbignumです。defglobalを使う場合、lambdaが生成させるときにbignumの自由変数があるときにこの領域へコピーして恒久化します。
bigcell[int] heap[cell structure]
|(temporarly)| | |
|LSB ^ | |-------<|BIGX(car) |
| | | | | |
|MSB |<------ | |
|(working) | big_pt0 | |
| | | |
|(parmanent) | big_pt1 | |
bignumが生成される場合にはheap領域にbignumであることを示すセルが生成されます。このセルにはその実体をかくのうしたbigcellの配列添え字を保持します。他にその桁数、正負の区分、bignumであることを示すクラスが記憶されます。
GCとの関係
ガベージコレクションが起きた時には束縛されていないheap領域のセルのみを回収します。それが指示していた実体は回収しません。bigcellは2000万要素を用意してあります。千桁のbignumを表すのにおよそ１１１桁を要します。bigcellはおよそ18万個の千桁の数を保持できます。通常のLisp利用においてこれを使い切ることはないだろうと思います。敢えてGCで回収をしていません。計算が終了してREPLに戻ったときにポインタをゼロクリアするようにしています。
クヌース先生のDアルゴリズム
加算、減算、乗算は普通の筆算と同様のアルゴリズムとしています。乗算についてはカラツバ、FFTの応用が将来的に考えられます。除算については筆算と同様なアルゴリズムですが計算効率化のためにクヌース先生のDアルゴリズムを採用しています。序数が一定の条件を満たす場合には割り算の見込み違いによる計算回復が最悪2回で済むというものです。詳しくはThe art of computer programming 準数値算法/算術演算 第4巻p89を参照してください。
筆算の10進法程度ならば失敗による回復計算はせいぜい9回ですが、10^10進法だと見込み違いによる加算のしなおしは最悪10^10-2になります。当初、これを入れないでいてリソースエラーをおこしていました。
isqrt 整数平方根
brent-salaminのアルゴリズムによ多桁円周率の計算ではisqrtの性能がものをいいます。ニュートン法によって計算しています。ニュートン法は収束が速いのですが、初期値を適切に与える必要があります。bignumでの計算は多量にセルを消費します。できるだけ近い値を初期値として与える必要があります。次の方法によりました。
bignumの桁数が奇数の場合
最上位の32bit整数の平方根をCの関数を利用して計算します。最上位を除いた桁数は偶数です。そこでこの1/2の桁数分のシフト計算をします。これにより近似したbignumが得られます。ニュートン法では大きな数から接近、収束させるためcで計算したsqrtに１を加算しています。
令和4年7月21日追加
上記の方法だと問題があることがわかりました。セルが奇数で最上位セルの値が１の場合です。例えば (* (expt 10 99)(expt 10 99))はそうなります。これに上述の１を加算すると初期値が真の値の2倍にもなります。桁数が大きくなると収束に時間とメモリを食います。最上位の次のセルの値も含めてsqrtをとる必要があります。
bignumの桁数が偶数の場合
最上位と2番目の32bit整数を合成して64bit整数とします。これのsqrtをC関数で計算します。残りの桁数は偶数ですのでこれを1/2分、シフト計算します。
これによりおおよそ5回程度の計算で収束するようになりました。
性能比較
brent-salaminの漸化式を用いた次のコードで多桁の円周率計算をしました。
;; for BIGNUM test
;; calculate PI with bignum using brent-salamin method
;; e.g. (pi-brent-salamin 1000)
(defun square-root (x)
(isqrt x))
(defun square (x)
(* x x))
#|
;;For Common-Lisp
(defun div (x y)
(floor x y))
|#
; Compute pi using the 'brent-salamin' method.
(defun pi-brent-salamin (nb-digits)
(let ((one (expt 10 nb-digits)))
(pi-brent-salamin1 one
one
(square-root (div (square one) 2))
(div one 4)
1)))
(defun pi-brent-salamin1 (one a b tt x)
(if (= a b)
(div (square (+ a b)) (* 4 tt))
(let ((new-a (div (+ a b) 2)))
(pi-brent-salamin1 one
new-a
(square-root (* a b))
(- tt (div (* x (square (- new-a a))) one))
(* 2 x)))))
最大で10000桁程度までは計算可能でした。1000桁程度なら実用時間で計算可能です。しかし10000桁ともなると5秒程度かかります。GMPを採用しているSBCLでは1万桁であっても瞬時であり1秒もかかりません。
T
> (pi-brent-salamin 1000)
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201876
>
> (time (pi-brent-salamin 1000))
Elapsed Time(second)=0.034954
<undef>
> (time (pi-brent-salamin 8000))
Elapsed Time(second)=2.422727
<undef>
> (time (pi-brent-salamin 10000))
Elapsed Time(second)=4.409718
<undef>
>
桁数が増加すると指数関数的に時間が増加しています。これは乗算が素朴な筆算アルゴリズムであるためO(n^2)の時間を要するからと思われます。
将来展望
乗算の効率化により改善が期待できます。通常ですとカラツバの算法を用いるところですが、敢えてこれを使いません。2022年2月～のロシア人の蛮行に対する私のささやかな抵抗です。FFTなど他の高速アルゴリズムを試してみたいと思っています。
実装
Githunにおいてあります。よかったらお試しください。
https://github.com/sasagawa888/eisl