34桁目のせめぎあい - きどたかのブログ

なかなか34桁目まで合わすというのは苦労が伴う。

とりあえずMathContextに69桁くらい設定したら、なんとか計算出来たっぽい。
平方根の計算で検算しながらやってみた。
9.99999999999999999999999999999999999997E6144の平方根の難易度が高い。
検算(2乗)すると先頭から68桁目にわずかな差が出るようだ。
68桁目まで計算しておけと････。

$x^y=e^{y\times \ln (x)}=10^{y\times \log (x)}=10^{y\times \frac{\ln (x)}{\ln (10)}}$
出来るだけ10進数で計算したいので展開してる。
割り算が厄介な誤差を生むんだろうなぁ。

$e^x=\sum^{\infty}_{n=0}{\frac{1}{n!}\times x^n$
これから指数関数をJavaで実装してる。
当然割り算が出てくるから誤差は出る。

$\ln (x)=\sum^{\infty}_{n=1}{\frac{(-1)^{n+1}}{n}{(x-1)}^n}$
これから対数関数を実装してる。
ただし、これはxの範囲に制限があるので、 $\ln (10)$ の計算は出来ない。

$\ln (x)= 2 \left ( \frac{x-1}{x+1} + \frac{1}{3}{\left(\frac{x-1}{x+1}\right)}^3 + \frac{1}{5}{\left(\frac{x-1}{x+1}\right)}^5 + \cdots \right ) = 2 \sum_{n=0}^\infty \frac{1}{2n+1} {\left ( \frac{x-1}{x+1} \right ) }^{2n+1}$
$\ln (10)$ を求めるのにはこちらの式を使って求めてる。

$\ln (x)=\ln(m \times 10^n)=\ln (m) + n \times \ln (10)$
この式も大事だ。xの小数点位置を調整することで、mを小さくして、xの範囲がある方の計算式でまとめられるようにする。
また、当然のことながら $\ln (10)$ を先に計算しておくことで、多少の演算時間の短縮が狙える？？
いや、どちらかというと、 $\ln (10)$ を事前に準備しておくという特性を考えて、
$\ln (10)$ だけめちゃくちゃ高い精度で計算を済ませておくことで、精度の向上に繋がるだろう。
いちおう有効桁100桁で計算することに成功しているので、
その文字列をストックしておくことで、JVMがあがる度に計算し直す必要はなくなるだろう。
この計算式はmを1より小さい値まで小数点を調整するように考えているので、 $\ln (m)$ の計算結果は必ずマイナス値になる。
ふむ、だからどうしたってんだろうけど、34桁目に直結する要素を探してるだけだ・・・。

指数関数での計算をやめる有効桁数が大事なんだと思う。
$x^y=10^{y\times \frac{\ln (x)}{\ln (10)}}$
この式において、指数部分は、整数部と小数部に分解する。
10の整数乗はスケールで解決出来るので、有効桁には直結しない。
つまりは10の累乗における指数部分が1未満の計算が大事。
実はここは底eの指数関数に戻して計算している。
$10^x=e^{(x \times \ln (10))}$
$x=\log (y)=\frac{\ln (y)}{\ln (10)}$
$\ln (y) =x \times \ln(10)$

ここでeの指数部分にくる数字の桁数次第で、何桁目まで計算してるのかが左右される。
それは指数関数の中身を見ればわかる。
自分の場合は最低でも34プラスアルファの桁で与えていた。
2乗で68桁、3乗で102桁、4乗で136桁、5乗で170桁、6乗で204桁。
ただ、先頭34桁目に影響を与えない累乗の数が6乗で十分かというとそうではない。

0.9で考えてみる。
n= 0,x^n/n!=1
n= 1,x^n/n!=0.9
n= 2,x^n/n!=0.405
n= 3,x^n/n!=0.1215
n= 4,x^n/n!=0.0273375
n= 5,x^n/n!=0.00492075
n= 6,x^n/n!=0.0007381125
n= 7,x^n/n!=0.000094900178571428571428571428571429 (9.4900178571428571428571428571429e-5)
n= 8,x^n/n!=0.000010676270089285714285714285714286 (1.0676270089285714285714285714286e-5)
n= 9,x^n/n!=0.0000010676270089285714285714285714286 (1.0676270089285714285714285714286e-6)
n=10,x^n/n!=0.000000096086430803571428571428571428571 (9.6086430803571428571428571428571e-8)
n=11,x^n/n!=0.0000000078616170657467532467532467532468 (7.8616170657467532467532467532468e-9)

こいつらを足していき、かつ目標の有効桁での丸めに差が出ないところまで計算するのか。
途方も無い。
割り切れなくて誤差が出てる部分もあるんだよな、きっと。

$\exp (x)=\lim_{n \to \infty}\left( 1 + \frac{x}{n}\right)^n$
かと言って、こちらの公式で求めるのは非効率すぎるか？
さきほどの0.9の実験からも言えるのだが、nが10を超えると、必ず最低1桁は減っていく。
とすると、n=34+10と、n=34+11あたりを計算して、丸めてみるというアリか。

あと、経験則ではあるが、有効桁34桁を算出する際には3桁余分な精度を与えると、
大半の計算は正しい答えを導き出せる。完璧にではない。
これは、java.math.BigDecimal#pow(int,MathContext)の説明の中にある、
mc.precision()+elength+1の内容と一致しているなと思える。
ただ、あれは整数乗しかしないから、大半の確率がだいぶ違うのだろう。
あのメソッドは2ulp以内の誤差を述べているけど、
いまのところ僕のは1ulpの誤差がたまに見つかる程度だ、当然めちゃくちゃ計算して遅いのだが。

まだまだやらないといけないことはあるようだ。
xの値が1未満ならexpm1も実装しておいた方がよさそうだ。
どうやるんだろう。

こうかな。
$e^x-1=\sum^{\infty}_{n=1}\frac{1}{n!}x^n$

まあ、やってみよう。