Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

素因数分解

この例では、Symbolic Math Toolbox™ を使って sym オブジェクトに対して初等関数を使用する方法を示します。

MATLAB® の組み込み整数型は、2^64 未満の整数に適しています。しかし、ここでは、より大きな整数の素因数分解について統計調査を行うことにします。これを行うには、サイズに制限の無い、シンボリック整数を使用します。N0+1N0+100 の間にある整数を調べます。ここで、N0=3*1023 です。組み込みデータ型はこのような値を正確に格納できません。そのため、最も内部の数値を sym でラップし、計算ではシンボリック表現を使用します。これで丸めやオーバーフローによる誤差を回避します。

N0 = 3*sym(10)^23;
disp(['Roundng error for doubles: ' char(3*10^23 - N0)]);
Roundng error for doubles: -25165824
disp(['Overflow error for integers: ' char(3*uint64(10)^23 - N0)]);
Overflow error for integers: -299981553255926290448385

算術演算では、シンボリック数は double との組み合せが可能であり、演算は変換後に実行されます。したがって、次の定義で丸め誤差が発生することはありません。

A = N0 + (1:100);

factor を使用して A の要素の素因数分解を計算します。素因数の数はさまざまです。配列には異なる長さのベクトルを格納できませんが、cell 配列には格納できます。メモリの再割り当てを回避するには、cell 配列を初期化してから、ループで因数分解を計算します。

Bcell = cell(1, 100);
for i=1:100
   Bcell{i} = factor(A(i)); 
end

より効率的な方法は、arrayfun を使用することです。UniformOutputfalse に設定すると、結果が cell 配列として返されます。

Bcell = arrayfun(@factor, A, 'UniformOutput', false);

たとえば、主な素因数分解には以下があります。

Bcell{1:5}
ans = (13432332303316007278478583)
ans = (2171739911223244939171805943)
ans = (311471392531549797184491917)
ans = (22330131295383776910994983)
ans = (5627126502668233610146697)

max を使用して最大の素因数を取得します。出力が sym オブジェクトで構成される場合、出力が一様であっても、オプション UniformOutput を常に false に設定する必要があります。

Mcell = cellfun(@max, Bcell, 'UniformOutput', false);

たとえば、主な最大素因数には以下があります。

Mcell{1:5}
ans = 2303316007278478583
ans = 171805943
ans = 549797184491917
ans = 76910994983
ans = 3610146697

cell 配列をシンボリック ベクトルに変換し、最大の素因数とその数値の長さとの比を全体的に調べます。double の場合と同様に、シンボリック ベクトルに対し、算術演算を要素ごとに適用できます。ほとんどの統計関数では、引数を倍精度の数値にする必要があることに注意してください。

M = [Mcell{:}];
histogram(double(log(M)./log(A)), 20);
title('Ratio of lengths of the largest prime factor and the number');

Figure contains an axes object. The axes object with title Ratio of lengths of the largest prime factor and the number contains an object of type histogram.

同様に、素因数の数の分布を調べます。ここで、結果には一様な数値データが含まれます。したがって、UniformOutputfalse に設定する必要はありません。

omega = cellfun(@numel, Bcell);
histogram(omega);
title('Number of prime factors');

Figure contains an axes object. The axes object with title Number of prime factors contains an object of type histogram.

調査対象の区間には 2 つの素数が含まれます。

A(omega == 1)
ans = (300000000000000000000037300000000000000000000049)

最大素因数は、4 を法とする剰余類 1 および 3 にほぼ同じ頻度で属することを確認します。ここで、sym オブジェクトの方程式はシンボリック オブジェクトそのものであり、論理値ではありません。シンボリック オブジェクトを変換してからその和を計算する必要があります。

sum(logical(mod(M, 4) == 1))
ans = 49
sum(logical(mod(M, 4) == 3))
ans = 51