昨今の円周率計算ではここで語られてるChudnovskyのアルゴリズムが良く使われているようだ。
ただ、この論文で紹介されてるアルゴリズム(2例)だと、何故か計算が上手く行かないんだよな〜。
しかも、ここで使われてる再帰はあまり効率が良くない。
というわけで、例によって「数式を直接プログラム」していこう。
手続き型言語とは違って、関数型プログラミングだと数式を直接プログラムした方が早かったりする、ってのは何度も指摘している。
「考えなくて良い」のだ。人間コンパイラにならないで良い。
直感的に数式を見たままプログラミングすれば良いのでラクなのである。
そしてスピードだ何だ、ってのは極論、こっちの責任ではなくってインタプリタ、ないしはコンパイラのせいにしちまえば良いのだ。
他力本願。それが関数型プログラミングの真骨頂である。
多分。
とは言っても、まずは階乗関数factをプログラムする。
これ自体も再帰のネタによくなってるが、高階関数を使う段階に入ってるのなら、reduceやそれに類する関数で書くのがお洒落だろう。
ここではSRFI-1に含まれるfoldを用いる。
(require srfi/1)(define (fact n)
(fold (lambda (x y)
(if (zero? x)
y(* x y))) 1 (iota (+ n 1))))
これで準備はオッケーだ。
まずは論文の通り、定数A、B、C、を定義する。
(define A 13591409)
(define B 545140134)
(define C 640320)
次に、Chudnovsky級数の項をプログラムする。
次の部分だ。
これは当然次の一部分なのだが。
最終項が∞、となっている。
って事はまたもや考えるのが面倒臭いので無限長リストにした方が良さそうだ。
遅延評価も使うのが難しいんじゃなくって考えるのが面倒臭い時に使うシロモノである。
なんせ、遅延評価は英語でlazy evaluationだ。lazyとは「面倒臭い」と言う意味がある。だから惰弱な人向けの機能なのは間違いない。
考え方としてはこうだ。0から始まる無限長の整数列を用意する。そしてその各項に対して
を適用すれば良い。
また「0から始まる無限長の整数列」といった胡散臭いものはSRFI-41のstream-fromとして提供されている。
(define terms
(stream-map (lambda (n)
(/ (* (expt -1 n) (fact (* 6 n)) (+ A (* B n)))
(* (expt (fact n) 3) (fact (* 3 n))
(expt C (+ (* 3 n) 3/2)))))
(stream-from 0)))
これで終わりだ。
あとは円周率πを返す関数を書けば良い。
ここで、Chudnovskyの級数は円周率の逆数を返すように定義されてるが、Schemeだと"/"を使うだけで逆数を得ることが出来る。逆数の逆数だと元に戻るわけで、ここで我々は目出度く円周率πを得られる。
(define (π n)
(/ (* 12 (stream-ref (stream-scan +
(stream-car terms)
(stream-cdr terms)) n))))
関数stream-scanは無限長のStream用のfoldである。
SRFI-41には大まかに同機能で2種類の関数が用意されてる事が多く、一つは有限長用、もう一つは無限長用、となってる。stream-foldは有限長用だからそれを使うと上手く行かないのだ。
いずれにせよ、ここでstream-scanは初期値を(stream-car terms)としてtermsの和を取っていったStreamを返す。要するにここでΣのStreamを形成してる。
それからstream-refで任意番の和を返し、12をかけて逆数を取れば、円周率πの出来上がりである。
> (π 0)
3.1415926535897345
> (π 1)
3.1415926535897936
Chudnovsky級数は収束が早いのが特徴のようで、n = 1の時点で、小数点以下15桁までの精度が確約される。
とは言っても1桁1桁出力させるのはちと手間なカンジがする。
とまぁ、4回に渡って4種類の円周率計算のアルゴリズムを紹介してみたが如何だったろうか。
個人的には、やはり扱いやすさから言うとHaskellで書かれたJeremy Gibbonsのアルゴリズムが一番良いような気がした。
両者とも遅延評価を扱えるプログラミング言語なので、考え方も移植性もピカイチだと思う。
反面、Cで書かれることを前提としたアルゴリズムは、数学的にはともあれ、移植性には難がある、ってのが正直なトコだ。
