一番おかしな所は階乗の所です。for文の使い方というよりも、intで大きな階乗を扱うのが不味いです。階乗は大きくなるのがはやく、13!=6,227,020,800 であり、2**32=4,294,967,296 を越えます。現在のほとんどのパソコンではintは32bitですので、すぐに限界になってしまいます。
そのほか<complex.h>をincludeしているのに使ってないとか、ちょっと色々怪しいです。参考までに作ってみたので見比べてみてください。(C11オンリー、gccに"-std=c11"をつけないと無理かも知れない。)
C
1/**
2 * オイラーの等式の計算
3 * $e^{i\pi} = \sum_{n=0}^{\infty}\frac{\pi^n}{n!}i^n = -1$
4*/
5
6#include <complex.h>
7#include <limits.h>
8#include <stdio.h>
9#include <tgmath.h>
10
11void print_complex(long double complex c);
12
13int main(void)
14{
15 long double complex eipi = 1;
16 long double complex prev_eipi = 1;
17 long double pi = acos((long double)-1);
18 // $\frac{\pi^n}{n!}$
19 // 最初の \pi^0 が 0 になるように pi で割っておく
20 long double pin_n = 1 / pi;
21 for (long long n = 0; n < LLONG_MAX; n++) {
22 pin_n *= pi;
23 if (n != 0) {
24 pin_n /= n;
25 }
26 // 計算できないほど小さいなら終わりにする。
27 if (pin_n == 0) {
28 break;
29 }
30 prev_eipi = eipi;
31 switch (n % 4) {
32 case 0:
33 // $n^(4k+0) = (n^4)^k = 1$
34 eipi += pin_n;
35 break;
36 case 1:
37 // $n^(4k+1) = (n^4)^ki = i$
38 eipi += pin_n * I;
39 break;
40 case 2:
41 // $n^(4k+2) = (n^4)^ki^2 = -1$
42 eipi += -pin_n;
43 break;
44 case 3:
45 // $n^(4k+3) = (n^4)^ki^3 = -i$
46 eipi += -pin_n * I;
47 break;
48 }
49 printf("%lld: ", n);
50 print_complex(eipi);
51 printf("\n");
52 // long doubleの限界を超えて同じになるなら終わりにする。
53 if (prev_eipi == eipi) {
54 break;
55 }
56 }
57 printf("%s = ", "e^(i*pi)+1");
58 print_complex(eipi);
59 printf("\n");
60 return 0;
61}
62
63void print_complex(long double complex c)
64{
65 long double c_r = creal(c);
66 long double c_i = cimag(c);
67 if (c_i == 0) {
68 printf("%Lf", c_r);
69 } else if (c_r == 0) {
70 printf("%Lfi", c_i);
71 } else if (c_i > 0) {
72 printf("%Lg+%Lgi", c_r, c_i);
73 } else {
74 printf("%Lg%Lgi", c_r, c_i);
75 }
76}
Q/A
- $...$は? -> mathjaxを使うときのtexの記法。Qiitaあたりで使うと数式になるよ。https://stackedit.io/editorとかで確かめてみると良いかも。
<tgmath.h>
って何? -> _Genericを使ったmath。floatでもdoubleでlong doubleでもcomplexでも同じ関数名でできるようになる。
long double
って意味あるの? -> 4倍精度浮動小数点数対応のCPUとコンパイラならきっと…たぶん。
- 途中で止めてる? ->
long double
がdouble
と同じぐらいの環境では、34回目あたりで表示の限界になり、49回目で足しても値が変わらなくなる。これ以上は無意味なので、止めてます。
- 階乗するところは別の関数にしないの? -> 前に足した部分へ π/n かければ次になるのでそのままforで一緒に回している。関数にする場合はメモ化とかしないと遅くなる。
- M_PIは使わないの? -> M_PIはCの標準じゃ無いから、
acos(-1)
でわざわざ計算している。tgmath.hのせいであらかじめlong doubelにキャストしないとlong dobuleで計算してくれない。
- 最終的にはどうなるの? -> 手元では e^(i*pi)+1 = -2.27573e-19-1.43047e-19i となった。限りなく0に近くなるので間違ってはいないと思う。
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2016/07/12 21:29