● 自然対数の底 e の近似値を求めます。

　ずいぶんと、いろいろに表現をされる e です。

・ e の x 乗はテイラー展開式によって近似値を求める。
　e^x = 1 + x + x^2/2! + ... + x^n/n!

・ Σ[n=0 - ∞](1/n!)

・ その値は、2.718281828459045235360287471352662497757247093699959574966967627724...
と続く無理数である。これは、e^x をマクローリン展開し、x = 1 を代入すると求められる；
e = 1/0! + 1/1! + 1/2! + ... ＝ Σ_n=0^∞ { 1 / n! }

また、以下の式で定義することもある。
e = lim _n→∞ ( 1 + 1 / n )ⁿ

　しかし、ここでは数値計算します。 とりあえず mpf_t 型を使います。

　GMP 浮動小数点データ型は mpf_t といいます。 仮数部の精度は mpf_set_default_prec() で指定します。 指数部は固定長で、多くの場合 1 machine word です、limb のビット数と同じサイズです。
32-bit システムでの指数部はおよそ 2^-68719476768 to 2^68719476736 の範囲です。
　但し、mpf_get_str は mp_exp_t に収まる範囲の指数部のみを返します。 mpf_set_str は long に収まりきらない指数部を受け付けません。

　まず double 型の変数を使って計算手順を確認してから mpf_t 型に移行しました。
デフォールトの仮数部 64bit ではちょっとさびしいので、128bit に増やしました。

　　ソースプログラム

　実行結果

2: f 2, ta 0.5, a 0.50000000000000000000　＜* この部分は double 型での練習
3: f 6, ta 0.166667, a 0.66666666666666662966
4: f 24, ta 0.0416667, a 0.70833333333333325932
5: f 120, ta 0.00833333, a 0.71666666666666656305
6: f 720, ta 0.00138889, a 0.71805555555555544700
7: f 5040, ta 0.000198413, a 0.71825396825396814471
8: f 40320, ta 2.48016e-05, a 0.71827876984126970417
9: f 362880, ta 2.75573e-06, a 0.71828152557319213667
10: f 3628800, ta 2.75573e-07, a 0.71828180114638440212
11: f 39916800, ta 2.50521e-08, a 0.71828182619849278989
12: f 479001600, ta 2.08768e-09, a 0.71828182828616848887
13: f 6227020800, ta 1.6059e-10, a 0.71828182844675891872
14: f 87178291200, ta 1.14707e-11, a 0.71828182845822963198
15: f 1307674368000, ta 7.64716e-13, a 0.71828182845899435360
16: f 20922789888000, ta 4.77948e-14, a 0.71828182845904209319
17: f 355687428096000, ta 2.81146e-15, a 0.71828182845904486875
18: f 6402373705728000, ta 1.56192e-16, a 0.71828182845904497977
19: f 121645100408832000, ta 8.22064e-18, a 0.71828182845904497977
20: f 2432902008176640000, ta 4.11032e-19, a 0.71828182845904497977
a 2 + 0.71828182845904497977
(2.71828182845904509080 (double)M_E)
(2.7182818284590452354 (char)M_E)

mp_bits_per_limb 32　　　　　　　　　　　　＜* ここから mpf_t 変数を使います
old mpf_get_default_prec 64 bit
new mpf_get_default_prec 128 bit
31: f 8.22284e+33, de 1.21613e-34
0.7182818284590452353602874713526624938382
32: f 2.63131e+35, de 3.80039e-36
0.7182818284590452353602874713526624976386
33: f 8.68332e+36, de 1.15163e-37
0.7182818284590452353602874713526624977538
34: f 2.95233e+38, de 3.38716e-39
0.7182818284590452353602874713526624977571
35: f 1.03331e+40, de 9.67759e-41
0.7182818284590452353602874713526624977572
36: f 3.71993e+41, de 2.68822e-42
0.7182818284590452353602874713526624977572
37: f 1.37638e+43, de 7.26546e-44
0.7182818284590452353602874713526624977572

2.7182 8182 8459 0452 3536 0287 4713 5266 2497 7572 4709 3699 ...　＜* これは数表の値

★ 前半は double 型を使った練習です。 "2." を省いて、小数部のみを計算しています。
f は n の階乗、ta は 1/f、a は e の値の計算途中経過です。
a 2 + 0.71828182845904497977　　　　　＜* これが計算結果です。
(2.71828182845904509080 (double)M_E)　＜* これは (double)M_E を単に printf したもの
(2.7182818284590452354 (char)M_E)　　 ＜* これは math.h 定義の M_E を printf したもの

★ 後半が mpf_t 変数を使って実行した結果です。　"2." を省いて、小数部のみを計算しています。
mp_bits_per_limb 32　　　　　　　　＜* limb は 32bit でした
old mpf_get_default_prec 64 bit　　＜* mpf_get_default_prec の値は 64bit でした
new mpf_get_default_prec 128 bit　 ＜* mpf_set_default_prec(128) で 128bit に変えて確認
31: f 8.22284e+33, de 1.21613e-34　＜* n < 30 以下は報告省略、f = n!, de = 1/f です
0.7182818284590452353602874713526624938382
　　...
37: f 1.37638e+43, de 7.26546e-44　　　　　　 ＜* de は 1/n! の値です
0.7182818284590452353602874713526624977572　　＜* これが計算結果

2.7182 8182 8459 0452 3536 0287 4713 5266 2497 7572 4709 3699 ...　＜* これは数表の値

★ 「"2." を省いて、小数部のみを計算しています。」 という意味は、1/0! + 1/1! を計算に含めなかった、ということです。

　とりあえず、40 桁くらいの e の値は計算できました。 (本格的な誤差の評価は行っておりません。)