高精度の数値計算

この例では、Symbolic Math Toolbox™ を使った可変精度の演算により、高精度の計算を得る方法を説明します。

ほとんど整数を表す式を求めます。典型的な例は次になります。$$\exp (\sqrt{163}\cdot \pi )$$ を 30 桁まで計算します。結果は、丸め誤差を伴って表示されますが、整数になるように見えます。

digits(30);
f = exp(sqrt(sym(163))*sym(pi));
vpa(f)

ans = $$262537412640768743.999999999999$$

同じ値を 40 桁まで計算します。これが整数にはならないことがわかります。

digits(40);
vpa(f)

ans = $$262537412640768743.9999999999992500725972$$

この現象をさらに調査します。以下では、最大 $$\exp (1000)$$ になる数値が出現します。そのため、調査には、小数点以下の正しい桁数が必要です。必要な作業精度を計算します。

d = log10(exp(vpa(1000)))

d = $$434.2944819032518276511289189166050822944$$

必要な精度を設定した後はじめて、それに左右される関数を呼び出します。他には、round、vpa および double がそのような関数です。

digits(ceil(d) + 50);

$$\exp \left(\sqrt{n}\pi \right)$$ という形式に類似した例を探します。もちろん、163 の平方根を乗算して、そのような数値 n を多数求めることができます。しかし、それ以外にも、この形式の数値の多くが何らかの整数に近い値を取ります。このことは、小数部のヒストグラム プロットからわかります。

A = exp(pi*sqrt(vpa(1:1000)));
B = A-round(A);
histogram(double(B), 50)

ほとんど整数となる式 $$\exp (n)$$ があるかどうかを計算します。

A = exp(vpa(1:1000));
B = A-round(A);
find(abs(B) < 1/1000)

ans =

  1x0 empty double row vector

今度は、A の要素の小数部がかなり均等に分布していることがわかります。

histogram(double(B), 50)