円周率の近似計算

円周率の近似計算には、ラマヌジャンの公式を使うと良いそうなので、試しに自分でもコードを書いてみた。

ラマヌジャンについて

シュリニヴァーサ・ラマヌジャン - Wikipedia

円周率の公式(ラマヌジャン)

{\displaystyle
\frac{4}{\pi} = \frac{1}{882} \sum_{n=0}^{\infty} \frac{(-1)^n (4n)! (1123+21460n)}{(4^n n!)^4 882^{2n}}
}

奇妙な数式に見える。

π以外は全て有理数なので、プログラムしやすい。このラマヌジャンの公式を使って、円周率πを近似計算するコードをSchemeで書いた。次の2つのサイトを参考にした。

円周率 - Wikipedia

keisan.casio.jp


(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"

収束は速い。私の感覚だが、ありえないくらい速い。

そして、数式の意味は全然わからない。