円周率の近似計算
円周率の近似計算には、ラマヌジャンの公式を使うと良いそうなので、試しに自分でもコードを書いてみた。
ラマヌジャンについて
円周率の公式(ラマヌジャン)
奇妙な数式に見える。
π以外は全て有理数なので、プログラムしやすい。このラマヌジャンの公式を使って、円周率πを近似計算するコードをSchemeで書いた。次の2つのサイトを参考にした。
(3回くらいミスと修正を繰り返した。たぶん、正しい答えが出るコードだと思う)
;;power: aのn乗を求める関数。 ;;ただしnは非負の整数とする。 (define (power a n) (define (power-p a n tmp) (if (= n 0) tmp (power-p a (- n 1) (* tmp a)))) (if (and (number? a) (integer? n) (>= n 0)) (power-p a n 1) (error "power:(power a n). a must be a number! n must be a nonnegative integer!"))) ;;sum-ramanujan-formula: Ramanujanの公式より 882 * 4 / pi を近似する有限和を求める関数。 ;;注. 分子と分母の対を整数で返す。 ;;注. この有限和の際、通分は簡単であった。約分は一切していない。 (define (sum-ramanujan-formula m) (define (sum-ramanujan-formula-p last-numerator last-denominator flg tmp4nfacto n m) (if (> n m) (cons last-numerator last-denominator) (let ((step-d (* (power (* 4 n) 4) (power 882 2))) (next4n (* (* 4 n) (- (* 4 n) 1) (- (* 4 n) 2) (- (* 4 n) 3)))) (sum-ramanujan-formula-p (+ (* (* -1 flg) next4n tmp4nfacto (+ 1123 (* 21460 n))) (* step-d last-numerator)) (* step-d last-denominator) (* -1 flg) (* next4n tmp4nfacto) (+ n 1) m)))) (sum-ramanujan-formula-p 1123 1 1 1 1 m)) ;;pi-ramanujan: Ramanujanの公式でpiの近似値を計算する関数。 ;;mは有限和の取り方。第m項までの有限和を計算する。 ;;lenは表示する桁数。例えばlen=3ならば3.14 (define (pi-ramanujan m len) ;; 準備 /disp 分数を小数表示する関数。 (define (/disp x y len) (define (/disp-p x y len flg) (cond ((< len 0) (error "wrong len")) ((= len 0) "done") (else (let ((tmp (floor (/ x y)))) (begin (cond ((= flg 0) (display tmp) (display ".\n")) ((= flg 1) (display tmp))) (/disp-p (* (- x (* tmp y)) 10) y (- len 1) 1)))))) (/disp-p x y len 0)) ;;ここから本体 (let ((x (sum-ramanujan-formula m))) (/disp (* 882 4 (cdr x)) (car x) len)))
以下が実行結果。
50桁表示で実行してみた。
;;有限和 5項目まで (pi-ramanujan 5 50) gosh> 3. 1415926535897932384626433832795028976444922297274"done" ;;有限和 6項目まで (pi-ramanujan 6 50) gosh> 3. 1415926535897932384626433832795028841971532961927"done"
収束は速い。私の感覚だが、ありえないくらい速い。
そして、数式の意味は全然わからない。