Main Content

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

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

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

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

r = exp(pi*sqrt(sym(163)))
r = eπ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

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

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

整数部の桁数が exp(π200) 桁、小数部の桁数が 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 = 

(eπ37eπ58eπ67eπ163)

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

A_nearint = vpa(A_nearint)
A_nearint = 

(199148647.999978046551856766500923875335924591257751.99999982221324146957619235527147197952743.9999986624542245068292613126262537412640768743.9999999999992500725972)

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

histogram(double(B),40)

Figure contains an axes object. The axes object contains an object of type histogram.

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

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

整数部の桁数が exp(π200) 桁、小数部の桁数が 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(π n) の形式がほとんど整数となる場合があまりないことを示しています。この例では、丸め誤差が 0.0001 未満であるほとんど整数が存在しません。

histogram(double(B),40)

Figure contains an axes object. The axes object contains an object of type histogram.

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

digits(32)

参考文献

[1] “Integer Relation Algorithm.” In Wikipedia, April 9, 2022. https://en.wikipedia.org/w/index.php?title=Integer_relation_algorithm&oldid=1081697113.

[2] “Almost Integer.” In Wikipedia, December 4, 2021. https://en.wikipedia.org/w/index.php?title=Almost_integer&oldid=1058543590.

参考

| |

関連するトピック