前回の記事でネイピア数の計算について少し書きました。計算の仕方の見通しが付くと、次は速さを求めてみたくなってきますよね。ということで、今回は本気でソフトを作りこんでいきます。
式変形
まずは式をできるだけコンピューターで計算しやすい形に変形して無駄な計算を減らします。おさらいですが、ネイピア数は次の範囲にあります。
\[1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}+\frac{1}{(n+1)!}<e<1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}+\frac{e}{(n+1)!}\]まずは左半分を整理します。
\[\begin{eqnarray*}1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}+\frac{1}{(n+1)!}&<&e\\\frac{{}_{n+1}P_{n+1}+{}_{n+1}P_{n}+_{n+1}P_{n-1}+\cdots+{}_{n+1}P_{1}+1}{(n+1)!}&<&e\end{eqnarray*}\]ただし${}_nP_r$は順列です。高校の数学でよく出てくるやつです。
つづいて右半分を整理します。
\[\begin{eqnarray*}e&<&1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}+\frac{e}{(n+1)!}\\\left\{1-\frac{1}{(n+1)!}\right\}e&<&1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}\\e&<&\left(1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}\right)\frac{(n+1)!}{(n+1)!-1}\\e&<&\frac{{}_{n}P_{n}+{}_{n}P_{n-1}+_{n}P_{n-2}+\cdots+{}_{n}P_{1}+1}{n!}\frac{(n+1)!}{(n+1)!-1}\\e&<&\frac{\left({}_{n}P_{n}+{}_{n}P_{n-1}+_{n}P_{n-2}+\cdots+{}_{n}P_{1}+1\right)(n+1)}{(n+1)!-1}\\e&<&\frac{{}_{n+1}P_{n+1}+{}_{n+1}P_{n}+_{n+1}P_{n-1}+\cdots+{}_{n+1}P_{1}}{(n+1)!-1}\end{eqnarray*}\]勘の良いひとは気付くかもしれませんが、$S$と$F$を
\[\begin{eqnarray*}S&=&{}_{n+1}P_{n+1}+{}_{n+1}P_{n}+_{n+1}P_{n-1}+\cdots+{}_{n+1}P_{1}\\F&=&(n+1)!\end{eqnarray*}\]と置くと、最初の不等式は次のように表せます。
\[\frac{S+1}{F}<e<\frac{S}{F-1}\]ああ、なんと美しいのでしょう。
ちなみに、${}_{n}P_{k}=(n-k+1)\cdot{}_{n}P_{k-1}$ですので、$S$は右側の項から順に計算していけば、一項につき掛け算を1回するだけで求まっていきます。そして$n!={}_nP_n$ですので、Sの各項を求めると最後はFが求まることになります。事前に式変形で計算量を減らせるのはこのあたりが限界でしょう。
実装
さて、上の式変形をもとにプログラムを書いてみましょう。
static (Fraction lower, Fraction higher) NapierRange(int order)
{
var permutation = BigInteger.One;
var totalNumerator = BigInteger.Zero;
for(int i = order + 1; i > 0; i--) {
permutation *= i;
totalNumerator += permutation;
}
var lower = new Fraction(totalNumerator + BigInteger.One, permutation);
var higher = new Fraction(totalNumerator, permutation - BigInteger.One);
return (lower, higher);
}
こんな感じでしょうか。要するにtotalNumeratorは$S$、permutationは$F$を意味しています。
なお、前回の記事では以下のようなプログラムを書きました。
static (Fraction lower, Fraction higher) NapierRange_F(int order)
{
Fraction napier = Fraction.One;
Fraction factorial = Fraction.One;
Fraction factorialFrom = Fraction.Zero;
for(int i = 0; i < order; i++) {
factorialFrom += Fraction.One;
factorial *= factorialFrom;
napier += Fraction.One / factorial;
}
factorialFrom += Fraction.One;
factorial *= factorialFrom;
Fraction lower = napier + Fraction.One / factorial;
Fraction higher = napier * factorial / (factorial - Fraction.One);
return (lower, higher);
}
どちらのメソッドでも同じ値を返してきますが、例えば、第1万項までを計算してみたところ、NapierRangeは0.12秒だったのに対し、NapierRange_Fは48秒でした。実に400倍も実行速度が違います。それだけ、通分や約分にはとてつもない計算量が必要ということですね。
並列化
さて、前回の記事より格段と速くなりましたが、もう一声いきたいですね。
現在の計算式はシングルスレッドです。ですが、近年のパソコンはCPUコアがたくさんついているので、並列化をできるならしたほうが計算速度の向上が見込めます。しかし、見ての通り、NapierRangeのforループは前回までの計算結果を参照しますので、各回が独立していません。すなわち、このままでは並列化できません。
ここで、順列の以下の性質を思い出してください。
\[{}_nP_k\cdot{}_{n-k}P_l={}_nP_{k+l}\]これをうまく使えば、
\[{}_{n}P_{1}+{}_{n}P_{2}+\cdots+{}_{n}P_{n}=\left({}_{n}P_{1}+{}_{n}P_{2}+\cdots+{}_{n}P_{\lfloor n/2 \rfloor}\right)\left(1+{}_{\lceil n/2 \rceil}P_{1}+{}_{\lceil n/2 \rceil}P_{2}+\cdots+{}_{\lceil n/2 \rceil}P_{\lceil n/2 \rceil}\right)\]と式変形できます。ただし$\lfloor \rfloor$は床関数、$\lceil \rceil$は天井関数です。こうすることで、左右のカッコ内の計算はそれぞれ独立で、計算量も均等です（任意のビット数を全く同じ速度で計算できるならばですが…実際はそんなことないのですが…）。今回は2分割ですが、うまくやれば好きな数で分割することができます。そうやって実装したのがこちらになります。
static (Fraction lower, Fraction higher) NapierRange(int order, int thread)
{
var opt = new ParallelOptions() { MaxDegreeOfParallelism = thread };
// 事前に通分すれば、分子は順列の総和になる
// 分子をthreadの個数で分割し、並列で計算する
var parallel = new (BigInteger permutation, BigInteger total)[thread];
Parallel.For(0, thread, opt, t => {
int from = (int)((long)order * (thread - t) / thread);
int to = (int)((long)order * (thread - t - 1) / thread);
BigInteger permutation = BigInteger.One;
BigInteger totalNumerator = BigInteger.Zero;
for(int i = from; i > to; i--) {
permutation *= i;
totalNumerator += permutation;
}
parallel[t].permutation = permutation;
parallel[t].total = totalNumerator;
});
// 分割した分を統合する
BigInteger permutation = BigInteger.One;
BigInteger totalNumerator = BigInteger.One;
for(int i = 0; i < thread; i++) {
totalNumerator += parallel[i].total * permutation;
permutation *= parallel[i].permutation;
}
// ここの時点でtotalNumeratorにはorder次までの分子の総和、permutationにはorder!が入っている
// 分母分子を(order+1)倍するとこの後の計算が超綺麗になる
permutation *= (order + 1);
totalNumerator *= (order + 1);
// 上限と下限は同時に計算する（独立なので）
Fraction lower = default, higher = default;
Parallel.Invoke(opt,
() => lower = new Fraction(totalNumerator + BigInteger.One, permutation),
() => higher = new Fraction(totalNumerator, permutation - BigInteger.One));
return (lower, higher);
}
区間を分割して並列で計算しますが、そのあと、それぞれを掛け合わせる再は結局はループ計算になってしまいます。でも多少は速くなるようです。ついでに、時間がべらぼうにかかる約分も上限側と下限側で別々にするようにしました。
試しに第10万項まで計算してみると、シングルスレッド版で5.5秒、マルチスレッド版で12秒でした。倍くらい速くなりましたね。項数が増えてくると差が顕著になると思われます。
分数の十進小数化
さて、ネイピア数の範囲が分数で表せるようになったところで、今度はそれを小数にしてやらねばなりません。これはむしろネイピア数を求めるよりも重たい処理になります。
前回は以下のようなコードを書きました。
static IEnumerable<char> ToDecimalString(Fraction lower, Fraction higher)
{
var lower_dec = lower.ToBigInteger();
var higher_dec = higher.ToBigInteger();
if(lower_dec == higher_dec) {
foreach(char c in lower_dec.ToString())
yield return c;
yield return '.';
var lower_rem = lower.Numerator;
var higher_rem = higher.Numerator;
while(true) {
lower_rem = lower_rem % lower.Denominator * 10;
higher_rem = higher_rem % higher.Denominator * 10;
char lower_digit = (lower_rem / lower.Denominator).ToString()[0];
char higher_digit = (higher_rem / higher.Denominator).ToString()[0];
if(lower_digit != higher_digit)
break;
else
yield return lower_digit;
}
}
}
大まかな流れはこれでOKです。いろいろ考えて試してみましたが、このパターンが一番速いようです。ですが、少しだけ追加で方針を加えます。
- DivRemメソッドを使用する
- 並列化をする
- 1桁ずつじゃなくて、n桁(nは1000くらい)ずつ計算する
まずDivRemメソッドですが、除算をしたときに商と剰余を一緒に出してくるメソッドになります。C#には除算と剰余算で別々の演算子を持っていますが、これら2つは片方の計算をしようとするともう片方も必然的に出てくるものです。別々に計算していてはほぼ同じ計算を2度することになってしまいます。
並列化は字のままです。独立している計算、すなわち下限と上限それぞれを並列で計算しましょう。
最後ですが、10で割った商を次の桁、その余りを10倍してループの次へとやっていくと、1桁を求めるのに1回ずつ除算と掛け算をしなければならなくなってしまいます。しかし、それらはどちらも計算コストの高い計算です。そこで、100で割った商を次の2桁とすれば、除算や乗算をする回数を半分にできます。1000で割れば…と、ここの数を増やすことで格段と計算の回数を減らせるのです。その代わり商は大きくなってしまって、商を十進数にするのに少し時間がかかるようになってしまいますが、DivRemのほうが遥かに重い処理であるようです。そのため、1000桁、すなわち10^1000くらいで割って、一気に1000桁分を求めるほうが遥かに計算速度は速くなります。
そうして実装したコードが次となります。これは1000桁ごとです。digit_unitをいじることで何桁ごとか変えることはできます。
static IEnumerable<char> ToDecimalString(Fraction lower, Fraction higher, int thread)
{
var opt = new ParallelOptions() { MaxDegreeOfParallelism = thread };
const int digit_unit = 1000;
BigInteger digitPow = BigInteger.Pow(10, digit_unit);
// まずは整数部分を比較
var lower_dec = lower.ToBigInteger();
var higher_dec = higher.ToBigInteger();
if(lower_dec == higher_dec) {
foreach(char c in lower_dec.ToString())
yield return c;
yield return '.';
var lower_deno = lower.Denominator;
var higher_deno = higher.Denominator;
var lower_rem = lower.Numerator % lower_deno;
var higher_rem = higher.Numerator % higher_deno;
while(true) {
string lower_digit_str = "", higher_digit_str = "";
// DivRemの処理がほかの処理の500倍くらい時間かかるので、これ以上の高速化は難しい。
Parallel.Invoke(opt, () => {
(var lower_digit, lower_rem) = BigInteger.DivRem(lower_rem * digitPow, lower_deno);
lower_digit_str = lower_digit.ToString().PadLeft(digit_unit, '0');
}, () => {
(var higher_digit, higher_rem) = BigInteger.DivRem(higher_rem * digitPow, higher_deno);
higher_digit_str = higher_digit.ToString().PadLeft(digit_unit, '0');
});
foreach(var (lower_c, higher_c) in lower_digit_str.Zip(higher_digit_str)) {
if(lower_c == higher_c)
yield return lower_c;
else
yield break;
}
}
}
}
意外とこれで速くなって、従来コードだと第10万項まで（＝45.6万桁）の計算で380秒かかりましたが、今回のコードで11秒まで縮まりました。38倍とは…かなりですね。
1000万項（約6560万桁）への挑戦
ここまで計算が速くなったら、大きな桁数のネイピア数に挑戦してみたくなりますよね。ただし計算量は桁数が増えるとともに爆発的に増えていきます。
ということで、1000万項までの計算をしてみました。それで求まったネイピア数は65,657,065桁、計算には丸々4日かかりました。
こちらからダウンロードしてください。
65.6 million digits.zip (33.1MB)
1億桁くらい計算してみたい気もしますが、結構時間がかかりそうですのでまあ何かの機会にでも。
0 件のコメント:
コメントを投稿