指数関数を計算するアルゴリズム
7381 ワード
ことばを引く
前の随筆で自然対数を計算する高速アルゴリズムを紹介しました.指数関数を計算するアルゴリズムを見てみましょう.指数関数exはテーラー級数として展開できることを知っています.
この級数は全体実数xに対して収束し,xが零に近づくと収束が速い.
このアルゴリズムを実現するC〓〓手続き
前に述べたexのテーラー級数展開式によれば、以下のC〓手順を書き出してdecimalデータタイプにExp拡張方法を追加することができる.行目のexpmaxの値は、decimal.MaxValueの自然対数の近似値であり、Expメソッドがオーバーフローしているかどうかを検出するために使用される(22行目). 20から30行目のExp拡張方法は、指数関数を計算するために用いられます. 本の方法は、exy=exyという数式を用いて、パラメータxを整数部分nと小数部分x-nに分けて計算する. 整数部分nはまた、1、2、4、8、16、32、64の諸数のうちのいくつかの和に分解され、事前に計算された定数を用いて計算される. 第25行は、e 66.5421をe 67 e-0044579に分解することを防止するためであり、e 67を計算するとオーバーフローする.e 66 e 0.5241に分解します. 第32行から37行のExponential方法は、テイラー級数を用いてexを計算する.そのパラメータqはゼロに近いほど計算が速い. このアルゴリズムはまだ早いです.35行目のforサイクルの実行回数は22回を超えません. テストプログラム
以下はdecimalデータタイプのExp拡張方法のテストプログラムです. decimal.MaxValue=79,228,162,514,264,337,593,543,950,335, e 67=125,236,317,084,221,378,051,352,196,0704.3657675348274… e 67は既にdecimalの最大値を超えていることが分かる.
事前に計算した定数
Decimal Extens.csプログラムの9行目から18行目のexpsスタティック読み取り専用配列には、e 1、e 2、e 4、e 8、e 16、e 32、e 64の値が格納されている.これらの値はどうやって得られますか?これは簡単です.Linuxオペレーティングシステムに高精度計算機bcがあります.まず次のような内容のテキストファイルを編集できます.i.txt:の前の7行はeの各乗です. 最後の行はdecimal.MaxValueの自然対数です. 参考資料 Wikipedia:Exponential function Wikipedia:Taylor series Linux man pages:bc-An arbitrary precision callator lagge
前の随筆で自然対数を計算する高速アルゴリズムを紹介しました.指数関数を計算するアルゴリズムを見てみましょう.指数関数exはテーラー級数として展開できることを知っています.
この級数は全体実数xに対して収束し,xが零に近づくと収束が速い.
このアルゴリズムを実現するC〓〓手続き
前に述べたexのテーラー級数展開式によれば、以下のC〓手順を書き出してdecimalデータタイプにExp拡張方法を追加することができる.
1 using System;
2
3 namespace Skyiv.Extensions
4 {
5 static class DecimalExtensions
6 {
7 static readonly decimal expmax = 66.542129333754749704054283659m;
8 static readonly int[] mask = { 1, 2, 4, 8, 16, 32, 64 };
9 static readonly decimal[] exps =
10 {
11 2.71828182845904523536028747135m, // exp(1)
12 7.38905609893065022723042746058m, // exp(2)
13 54.5981500331442390781102612029m, // exp(4)
14 2980.95798704172827474359209945m, // exp(8)
15 8886110.52050787263676302374078m, // exp(16)
16 78962960182680.6951609780226351m, // exp(32)
17 6235149080811616882909238708.93m // exp(64)
18 };
19
20 public static decimal Exp(this decimal x)
21 {
22 if (x > expmax) throw new OverflowException("overflow");
23 if (x < -66) return 0;
24 var n = (int)decimal.Round(x);
25 if (n > 66) n--;
26 decimal z = 1, y = Exponential(x - n);
27 for (int m = (n < 0) ? -n : n, i = 0; i < mask.Length; i++)
28 if ((m & mask[i]) != 0) z *= exps[i];
29 return (n < 0) ? (y / z) : (y * z);
30 }
31
32 static decimal Exponential(decimal q)
33 { // q (almost) in [ -0.5, 0.5 ]
34 decimal y = 1, t = q;
35 for (var i = 1; t != 0; t *= q / ++i) y += t;
36 return y;
37 }
38 }
39 }
簡単な説明は以下の通りです.以下はdecimalデータタイプのExp拡張方法のテストプログラムです.
1 using System;
2 using Skyiv.Extensions;
3
4 class Tester
5 {
6 static void Main()
7 {
8 try
9 {
10 foreach (var x in new decimal[] {
11 -100, -66, -65, -1, 0, 1, 2.5m, 16, 66.5421m, 67 })
12 Console.WriteLine("{0,-30}: exp({1})", x.Exp(), x);
13 }
14 catch (Exception ex) { Console.WriteLine(ex.Message); }
15 }
16 }
運転結果は以下の通りです.work$ dmcs Tester.cs DecimalExtensions.cs
work$ mono Tester.exe
0 : exp(-100)
0.0000000000000000000000000000: exp(-66)
0.0000000000000000000000000001: exp(-65)
0.3678794411714423215955237702: exp(-1)
1 : exp(0)
2.7182818284590452353602874714: exp(1)
12.182493960703473438070175950: exp(2.5)
8886110.520507872636763023741 : exp(16)
79225838488862236701995526357 : exp(66.5421)
overflow
e 67の計算でオーバーフローが発見されたことが分かる.なぜなら、事前に計算した定数
Decimal Extens.csプログラムの9行目から18行目のexpsスタティック読み取り専用配列には、e 1、e 2、e 4、e 8、e 16、e 32、e 64の値が格納されている.これらの値はどうやって得られますか?これは簡単です.Linuxオペレーティングシステムに高精度計算機bcがあります.まず次のような内容のテキストファイルを編集できます.i.txt:
scale=30
e(1)
e(2)
e(4)
e(8)
e(16)
e(32)
e(64)
l(2^96-1)
quit
上のeはexpを代表して、lはlnを代表して、296-1はdecimal.MaxValueです.次のコマンドを実行します.work$ bc -l exps_in.txt > exps_out.txt
下記の内容の出力exps_を得ることができます.out.txt:2.718281828459045235360287471352
7.389056098930650227230427460575
54.598150033144239078110261202860
2980.957987041728274743592099452888
8886110.520507872636763023740781450350
78962960182680.695160978022635108224219956195
6235149080811616882909238708.928469744831391846235799914388
66.542129333754749704054283659972
ちょっと整理すれば上のC〓のプログラムに使えます.