2011年8月16日火曜日

自然対数の底 e の計算 in Elisp

eを計算する方法はたくさんある。ここでは二つの方法を書いておく。

1. 級数

exの x=1 におけるテイラー展開を求めてやると、次のような級数表現が得られる。

e=\sum^{\inf}_{n=0}\frac{1}{n!}=1+\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\cdots

この級数表現をプログラムにすると、次のようになる。

(defun e (n)
  "級数表現にもとづいてeを計算"
  (if (< n 0) 0.0
    (+ (/ 1.0 (fact n)) (e (- n 1)))))
=> e

(defun fact (n)
  "階乗を計算"
  (if (<= n 0) 1
    (* n (fact (- n 1)))))
=> fact

(e 5)
=> 2.7166666666666663

n=1 から n=20 までの実行結果。

(2.0
2.5
2.6666666666666665
2.708333333333333
2.7166666666666663
2.7180555555555554
2.7182539682539684
2.71827876984127
2.7182815255731922
2.7182818011463845
2.718281826198493
2.718281808918177
2.7182818042763013
2.7182818091495133
2.718281802164987
2.7182817951863503
2.718281799212947
2.7182818049171673
2.7182818140377822
2.7182818360880554)

プロットしてみた。

2. 一般連分数

具体的な式は、「ネイピア数の表現 - Wikipedia」の「II. 一般連分数による表現」に書いてある。これを漸化式の形で書くと、

\begin{eqnarray}\begin{array}{rcl}f_n&=&2+\frac{1}{f_{n-1}}\\f_{n-1}&=&1+\frac{1}{f_{n-2}}\\f_{n-2}&=&2+\frac{2}{f_{n-3}}\\&\vdots&\\f_0&=&n\end{array}\end{eqnarray}

1, 2, ... と増えていく部分を i で表すと、

\begin{eqnarray}f_{n-i}=\left\{\begin{array}{l}2+\frac{1}{f_{n-1}} \ \ \ \ (i=0)\\i+\frac{i}{f_{n-i-1}}\ \ \ \ (1\leq i\leq n-1)\\n \ \ \ \ (i=n)\end{array}\end{eqnarray}

さらに、 j = n - i つまり i = n - j と変換すると、

\begin{eqnarray}f_{j}=\left\{\begin{array}{l}2+\frac{1}{f_{j-1}} \ \ \ \ (j=n)\\n-j + \frac{n-j}{f_{j-1}} \ \ \ \ (1\leq j\leq n-1)\\n \ \ \ \ (j=0)\end{array}\end{eqnarray}

これをプログラムで表現すると、

(defun e (n)
  "一般連分数による表現にもとづいてeを計算"
  (labels ((f (j)
              (cond
               ((= n j) (+ 2 (/ 1.0 (f (1- j)))))
               ((zerop j) n)
               (t (+ (- n j) (/ (float (- n j)) (f (1- j))))))))
    (f n)))
=> e

(e 5)
=> 2.7184466019417477

※浮動小数点数で計算するために、一箇所だけfloat関数を使っている。

n=1 から n=20 までの実行結果。

(3.0
2.6666666666666665
2.7272727272727275
2.7169811320754715
2.7184466019417477
2.718263331760264
2.71828369389345
2.7182816576664037
2.7182818427778273
2.7182818273518743
2.718281828538486
2.718281828453728
2.7182818284593786
2.7182818284590256
2.7182818284590464
2.718281828459045
2.7182818284590455
2.7182818284590455
2.7182818284590455
2.7182818284590455)

1の方法で求めた図に重ねてプロットしてみた（赤い点）。

0 件のコメント:

コメントを投稿