高精度な算術を使用したほとんど整数の検索

PSLQ アルゴリズムのような整数関係を求めるアルゴリズム [1] では、正確な結果を得るために高精度な算術が必要となります。アルゴリズムの入力の精度が高いほど、単なる数値アーティファクトではない整数関係をアルゴリズムが見つけられる可能性が高くなります。PSLQ アルゴリズムでは、ほとんど整数などが原因で誤検知が発生する場合があります。

この例では、Symbolic Math Toolbox™ の可変精度の算術を使用して、ほとんど整数 (整数に非常に近い数) を見つける方法を説明します。この例では、整数 $\mathit{n}=1,...,200$ について、$\exp \left(\pi \sqrt{\mathit{n}}\right)$ または $\exp \left(\pi \mathit{n}\right)$ の形式をもつほとんど整数 (整数に近い数) を検索します。

まず、ほとんど整数 [2] のよく知られた例、すなわち実数 $\exp \left(\pi \sqrt{163}\right)$ について考えます。この実数を正確なシンボリック数として作成します。

r = exp(pi*sqrt(sym(163)))

r = $${\mathrm{e}}^{\pi \hspace{0.17em}\sqrt{163}}$$

vpa を使用して、この数を可変精度の算術で評価します。既定では、vpa は有効桁数 32 桁まで値を計算します。

vpa(r)

ans = $$262537412640768743.99999999999925$$

digits を使用すると、有効桁数を変更して精度を高めることができます。同じ数を有効桁数 40 桁で評価します。

digits(40)
vpa(r)

ans = $$262537412640768743.9999999999992500725972$$

この数は整数に非常に近くなっています。この実数および最も近い整数との差を調べます。vpa を使用して、この結果を有効桁数 40 桁で評価します。

dr = vpa(round(r)-r)

dr = $$0.0000000000007499274028018143151532171817442410122968$$

次に、整数 $\mathit{n}=1,...,200$ について、$\exp \left(\pi \sqrt{\mathit{n}}\right)$ の形式をもつほとんど整数を検索します。これらの数を正確なシンボリック数として作成します。

A = exp(pi*sqrt(sym(1:200)));

整数部の桁数が $\exp \left(\pi \sqrt{200}\right)$ 桁、小数部の桁数が 20 桁となるように、有効桁数を設定します。

d = log10(A(end));
digits(ceil(d) + 20)

この数字列および最も近い整数との差を評価します。丸め誤差が 0.0001 未満となるほとんど整数を見つけます。結果を厳密なシンボリック型で表示します。

B = vpa(round(A)-A);
A_nearint = A(abs(B)<0.0001)'

A_nearint = 
$$\left(\begin{array}{c}{\mathrm{e}}^{\pi \hspace{0.17em}\sqrt{37}}\\ {\mathrm{e}}^{\pi \hspace{0.17em}\sqrt{58}}\\ {\mathrm{e}}^{\pi \hspace{0.17em}\sqrt{67}}\\ {\mathrm{e}}^{\pi \hspace{0.17em}\sqrt{163}}\end{array}\right)$$

小数点以下 20 桁以上の精度でほとんど整数を評価します。

A_nearint = vpa(A_nearint)

A_nearint = 
$$\left(\begin{array}{c}199148647.9999780465518567665009238753359\\ 24591257751.99999982221324146957619235527\\ 147197952743.9999986624542245068292613126\\ 262537412640768743.9999999999992500725972\end{array}\right)$$

差のヒストグラムをプロットし、その分布を示します。この分布を見ると、0 に近い差の度数が多くなっています。これは、$\exp \left(\pi \sqrt{\mathit{n}}\right)$ の形式がほとんど整数である場所です。

histogram(double(B),40)

次に、整数 $\mathit{n}=1,...,200$ について、$\exp \left(\pi \mathit{n}\right)$ の形式をもつほとんど整数を検索します。これらの数を正確なシンボリック数として作成します。

A = exp(sym(pi)*1:200);

整数部の桁数が $\exp \left(\pi \cdot 200\right)$ 桁、小数部の桁数が 20 桁となるように、有効桁数を設定します。

d = log10(A(end));
digits(ceil(d) + 20)

この数字列および最も近い整数との差を評価します。丸め誤差が 0.0001 未満となるほとんど整数を見つけます。結果は空の配列 sym となります。これは、この条件を満たすほとんど整数が存在しないことを意味します。

B = vpa(round(A)-A);
A_nearint = A(abs(B)<0.0001)

 
A_nearint =
 
Empty sym: 1-by-0
 

差のヒストグラムをプロットします。このヒストグラムは比較的均等に分布しており、$\exp \left(\pi \mathit{n}\right)$ の形式がほとんど整数となる場合があまりないことを示しています。この例では、丸め誤差が 0.0001 未満であるほとんど整数が存在しません。

histogram(double(B),40)

最後に、既定の精度である有効桁数 32 桁に戻して計算を続けます。

digits(32)

参考

sym | vpa | digits

関連するトピック

• 数値演算またはシンボリック演算の選択
• 数値のシンボリックへの変換
• 数値計算精度の引き上げ