|
数学関数 |
|
立方根 |
| 最終更新日：2006/07/01 書き直し |
●概要
除算を逆数法で行っても、あいかわらず、除算は乗算に比べて数倍以上遅い。従って、開根(平方、立方など)を行うニュートン法の漸化式に現われる多倍長除算(以下参照)を排除すれば、数倍の高速化が期待できる。
通常、Y = X3 - a とし、漸化式を
Xn+1 = Xn - 1 / 3 * (X3 - a) / X
となり、漸化式に逆数あるいは除算が現われ、速度的に不利となる。
●原理
根の逆数を求めることにすれば、関数は、
f(X) = 1 / Xm - a
となる。立方根なので m は 3 となる。これを、漸化式にすると、
Xn+1 = Xn - (1 / 3) * Xn * (a * Xn3 - 1)
となる。ここには、多倍長同士の除算は含まれていない。多倍長/Integer なる除算はあるが、これは速い。また、δ = X * (a * X3 - 1) とすれば、収束条件にもなる。
●方法
求めた、X の逆数を取れば、aの立方根となる。逆数はニュートン法で、しかも一回のみなので、良しとする。指数関数にて、a の1/3乗を求めるより圧倒的に速く求まる。
例のごとく、素直に適用すると初期値設定が困難あるいは不能になるので、以下のようにする。
a = M * 10P
= M * 10(3*Q+R)
と変形し、QとRは、
Q = P \ 3
R = P Mod 3 (|R|は、0,1,2）
である。こうすれば、
a = M * 10R * 103*Q
= m * 103*Q
となる。従って、1/ Cbrt(a) は、結局、
1/ Cbrt(a) = 1 /(m * 103*Q)1/3
= 1/Cbrt(m) * 10-Q
∴Cbrt(a) = Cbrt(m) * 10Q
で、求められる。即ち、1/Cbrt(m) を求めることに帰着される。10Q は、指数部の調整のみとなる。Rは、-2、-1、0、1、2 の値を取るので、m は、仮数部をもとにRに従って生成する。初期値としては、標準の数学関数にて、
X0 = (Pow(1 / m, 1 / 3))＜15桁をUltraLong化＞
と設定する。これにて、初期値の精度が、常に、15桁前後と安定する。a が、負の場合は、絶対値として演算し、結果の符号を調整する。
●高速化
ニュートン法にてまとめて解説しているので、ここでは結果のみ掲げる。
○単一初期値予測
初期値による立方根(逆数)の演算時間を専用ツールにてシミュレーションすると、下図のようになった。
3列目は初期値15桁の場合の基準演算時間、その右側に、列見出しの値を初期値にした演算時間である。薄緑色の時間は、その目標精度での最小値である。これによれば、概ね、目標精度の1/4か1/8の精度の初期値で求めた時間が最小となり、基準時間の1/2から1/3に改善されることが分かる。当初、初期値精度を512桁と固定していたが、本法の方が改善できる。
○多段初期値予測
実測は時間が掛かるので、同じく専用ツールでシミュレーションして見た。 結果は下図。
列見出しの15 は、基準時間、単一(1/8) は、目標精度の1/8 の初期値精度で求めた時間、1/4単位は、目標精度の1/4m 系列の初期値、1/8単位は、目標精度の1/8m 系列の初期値で求めた時間である。