おはようございます。
学校の体育の授業でサッカーができなさすぎて悲しいまがりかどです。
そんなことより今回は、
C++でネイピア数を求めてみよう
の回です。
数学で円周率の次ぐらいにとても有名な無理な定数ですね。
\(e^x\)は微分しても変わりませんっていうのはよく知られている性質ですよね。
で、今回はそんなネイピア数をC++のプログラミングによって求めていきたいと思います。
どうやってC++でネイピア数を求めるか
eの定義式は
\[\displaystyle\lim_{n\rightarrow\infty}\left(1 + \frac{1}{n}\right)^n\]
となっています。
\(e^x\)のマクローリン展開
まあ\(n\)にできるだけ大きい値を入れたらまあ大体求まるんですけどそれではあまり楽しくないので今回は\(e^x\)をマクローリン展開した式を使います。
https://mathtrain.jp/exptenkai
マクローリン展開の方法は上のページに分かりやすく書かれていたのでどうぞご参照を。
で、マクローリン展開の結果、
\[e^x = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + ・・・\]
となるらしいです。
この式、きれいきれいですね。
そして、\(x = 1\)のとき、
\[e = 1 + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + ・・・\]
となるのでこの式を使います。
ループでネイピア数eを求めていく
ということで、上のマクローリン展開の式を一項ずつ求めていって足していくのをループして求めていきます。
無限回やっても仕方がないのでこちらからnをコンピュータに入力させてn回ループさせます。
方針が決まりましたので、求めます。
C++でネイピア数を求めるプログラム
#include<bits/stdc++.h>
using namespace std;
const double e = 2.718281828459;
int main(){
int n;
cin >> n;
double ans = 0;
double div = 1;
for(int i = 0; i < n; i++){
if(i > 0) div *= 1 / (double)i;
ans += div;
}
cout << setprecision(13) << ans << endl;
if(fabs(ans - e) <= 1e-13) puts("Success");
else puts("Failure");
}一行目にあるのは求めたansの答え合わせ用の答えです。
で、上で書いたように普通にn回ループしながら求めていっているわけなんですけど、
という疑問が生まれるかもしれません。
そこについて解説します。
階乗によるオーバーフロー
マクローリン展開の式では分母に階乗が出てきていますよね。
これを何かしらの変数に入れると最大で\(n!\)が出てきます。
例えば\(n = 100\)のとき大体オーダーでいうと\(10^{157}\)が出てきてしまい、到底intでもlong long でもとにかく何でも入りません。
なので単純に階乗を求めていてはだめです。
ここで、n項目に注目します。(0項目から数えてます。）
\[\frac{1}{n!}\]
となっていますが、これって
\[\frac{1}{1}\times\frac{1}{2}\times・・・\times\frac{1}{n-1}\times\frac{1}{n}\]
というふうに変形できますよね。
これを一項ずつ求めていき掛けていくと一回一回逆数を取る形になるので膨大な数にならずに済みます。
この方法でやると計算量は\(O(n^2)\)になりますかね。
え、けど二重リープなんかないじゃん。
どうなってるの？
となるかもしれませんが、安心してください、説明します。
divに値を保存しておく
ゼロ項目で求めた\(\frac{1}{1!}\)って一項目以降でも使いますよね。
同様に一項目で求めた\(\frac{1}{2!}\)も二項目以降再び使うことになりますよね。
同じものを何回も求めるのは面倒なので省略しましょう。
ここでi項目の値を\(val_i\)とすると、
\[val_i = val_{i-1} + \frac{1}{i!}\]
となりますよね。
i回目のループの時もう\(val_{i-1}\)は使わないのでdivで更新していきましょう。
こんなかんじですね。
正誤判定
そして求めたネイピア数eが正しいかどうかを判定します。
13桁くらい目で確認すればいいのですがまあお勉強のためにも正誤判定のプログラムを作りましょう。
これ普通にans == eで判定すればいいのかと思ったんですけどこれでやるとtrueの場合も結構falseになってしまいました。
doubleの誤差
今回は13桁までをsetprecisionで指定して表示してしますが実際にはansにはもう少し桁があります。
なので14桁以降が0でない限りfalseとなってしまいます。
ということで今回はansとeの誤差が\(10^{-13}\)未満かどうかで正誤判定をしました。
求めたいのは\(10^{-12}\)の桁までだからです。
ネイピア数を求めるC++プログラムを実行
まずn = 100で実行してみると
100
2.718281828459
Success余裕ですね。
ここでちょっとプログラムをいじってnを20まで回してみます。
#include<bits/stdc++.h>
using namespace std;
const double e = 2.718281828459;
int main(){
int n;
for(int n = 1; n <= 20; n++){
double ans = 0;
double div = 1;
for(int i = 0; i < n; i++){
if(i > 0) div *= 1 / (double)i;
ans += div;
}
cout << "n = "<< n << " ans = " << setprecision(13) << ans << ' ';
if(fabs(ans - e) <= 1e-13) puts("Success");
else puts("Failure");
}
}こうすると、
n = 1 ans = 1 Failure
n = 2 ans = 2 Failure
n = 3 ans = 2.5 Failure
n = 4 ans = 2.666666666667 Failure
n = 5 ans = 2.708333333333 Failure
n = 6 ans = 2.716666666667 Failure
n = 7 ans = 2.718055555556 Failure
n = 8 ans = 2.718253968254 Failure
n = 9 ans = 2.718278769841 Failure
n = 10 ans = 2.718281525573 Failure
n = 11 ans = 2.718281801146 Failure
n = 12 ans = 2.718281826198 Failure
n = 13 ans = 2.718281828286 Failure
n = 14 ans = 2.718281828447 Failure
n = 15 ans = 2.718281828458 Failure
n = 16 ans = 2.718281828459 Success
n = 17 ans = 2.718281828459 Success
n = 18 ans = 2.718281828459 Success
n = 19 ans = 2.718281828459 Success
n = 20 ans = 2.718281828459 Successとなりました。
n = 16からSuccessですね。
これだとlong longを使えば階乗求めてもいけますね。
まあ最初はnがどれくらいか見当がつかないしなんかこっちのほうがきれいなプログラムなのでまあいいんじゃないですかね。
ありがとうございました。
さっきの式に階乗があったと思うけどどこにあるの？
階乗っぽいのなくない？