For Gauche 0.9.14Search (procedure/syntax/module):

Next: , Previous: , Up: ライブラリモジュール - ユーティリティ   [Contents][Index]

12.37 math.prime - 素数

Module: math.prime

このモジュールは、素数を扱うユーティリティ関数を提供します。

素数のシーケンス

Variable: *primes*

{math.prime} 素数の無限遅延シーケンスです。

;; show 10 prime numbers from 100-th one.
(take (drop *primes* 100) 10)
 ⇒ (547 557 563 569 571 577 587 593 599 601)
Function: reset-primes

{math.prime} *primes*から大きな素数を取り出すと、それ以前の素数も全てメモリに 残り続けます。*primes*変数がシーケンスの頭を抱えているからです。 もう素数を必要としないことがわかっている場合、それらのメモリがガベージコレクト されることが望ましいかもしれません。reset-primes手続きは *primes*をまだ現実化されていない遅延シーケンスに再束縛し、 次のGCで計算済みの素数シーケンスが回収されるようにします。

Function: primes

{math.prime} 新たな素数の遅延シーケンスを返します。 その時だけ素数を使いたい、という時に便利です。返されたシーケンスへの 参照が無くなれば、シーケンスはガベージコレクトされます。 (ある素数の計算には素数のシーケンスが最初から必要なので、 たとえシーケンスの頭ではなく途中への参照だけを持っていたとしても、 遅延シーケンスの中のサンクにはシーケンスの頭への参照が保持されています。 シーケンスがGCされるためには、いかなる部分への参照も残さないようにしなければなりません)。

primesが返す各シーケンスは独立しているので、素数の計算もそれぞれで (重複して)行われることになります。

単純なルールとして、プログラム中で何度も素数を使う必要があるのなら 変数*primes*を利用するのが良いでしょう。各素数の計算は一度しか 行われず、余分な計算を省くことができます。しかしその場だけ素数が欲しいなら、 primesを呼んで、仕事が済んだらシーケンスを捨ててしまえば、 不要なシーケンスがメモリに残りつづけることを心配しなくても済みます。

素数かどうかを調べる

Function: small-prime? n

{math.prime} 比較的小さな正整数 (*small-prime-bound*以下の正整数) に対して、 それが素数であるかどうかを判定し、素数なら#tを返します。 nがそれ以上である場合は常に#fを返します。

この手続きは確実に素数であるとわかるものを素早く判別する時に便利です。 #tが返れば確実に素数であるとわかるからです (入力が大きな素数の時に #fを返すことはありえますが)。 これに対し、下に述べるMiller-Rabin法では、合成数は確実に判別できますが、 素数であるかどうかは確実には言えません。

Variable: *small-prime-bound*

{math.prime} これより小さな数に対しては、small-prime?は決定的に 素数かどうかを判別します。現在の実装ではこの数は3.4e14よりちょっと大きな数です。

Function: miller-rabin-prime? n :key num-tests random-integer

{math.prime} 2以上の正確な整数nが素数かどうかを、確率的なMiller-Rabin法を使って判定します。 この手続きが#fを返したなら、nは確実に合成数です。 この手続きが#tを返した場合、nはおそらく素数ですが、 疑陽性である確率もわずかにあります。

ただし、nがある数(*small-prime-bound*) より小さければ、アルゴリズムは決定的で、#tが返るnは確実に素数です。

n*small-prime-bound*以上の場合は 確率的テストを用います。デフォルトでは7回、ランダムにベース整数値を選んで Miller-Rabinテストを適用します。試行回数はnum-testsキーワード引数で 変更可能です。合成数に対して誤って#tを返してしまう確率は たかだか(expt 4 (- num-tests))です。

確率的テストでは、miller-rabin-prime?はデフォルトで この手続き固有の、固定したランダムシードを使います。固定値なのは再現性を確保するためです。 異なる乱数系列を使いたければ、ランダムな整数生成手続きを random-integerキーワード引数に与えてください。 手続きは正整数kを取り、0からk-1までのランダムな整数値を 返すものでなければなりません。

Function: bpsw-prime? n

{math.prime} nが素数かどうかをBaillie-PSW法を用いて判定します (http://www.trnicely.net/misc/bpsw.html)。 このアルゴリズムは2^64 (約1.8e19) 以下の入力に対しては決定的であり、 正しい答えを返します。入力がそれ以上の場合、合成数に対して#tが返る可能性が あります (具体的な数はまだ見つかっていませんが)。素数に対して#fが返ることは 決してありません。

Miller-Rabin法より遅いですがカジュアルに使う分には十分に速いので、 上記の入力範囲で確実な答えを得たい場合は便利でしょう。

素因数分解

Function: naive-factorize n :optional divisor-limit

{math.prime} 正整数nを、(sqrt n)までの素数で順に割ってみることで 素因数分解します。戻り値は小さい順に並べられた素因数(2以上の整数)のリストです。

(naive-factorize 142857)
  ⇒ (3 3 3 11 13 37)

(naive-factorize 1)()を返します。

この方法は極めてナイーブなものですが、目安としてどの素因数も1e7程度以下であれば それなりに使えます。例えば次の例は2.4GHz Core2マシンで0.4秒で答えが返ります (ただし、初回の実行は遅延素数シーケンスの実現化があるので1.3秒ほどかかりますが)。

(naive-factorize 3644357367494986671013))
  ⇒ (10670053 10670053 32010157)

もちろんnがより大きなオーダーの素因数を含んでいると、性能は急激に 悪化します。安全に使うにはnを1e14程度のオーダーに止めておくのが良いでしょう。

オプショナル引数divisor-limitを与えると、試行する素数の上限を指定 できます。この引数がある場合、naive-factorizeは因数fdivisor-limit以下の素数で割りきれなければ、そこで諦めてfを 結果に含めます。この場合、結果の最後の要素は合成数であるかもしれないわけです。 これは、より高度な素因数分解アルゴリズムを適用する前にありきたりの素因数を 除外するのに便利です。

(naive-factorize 825877877739 1000)
  ⇒ (3 43 6402154091)

;; whereas
(naive-factorize 825877877739)
  ⇒ (3 43 4591 1394501)

この手続きは高速化のために小さなnに対する結果はメモ化しています。

Function: mc-factorize n

{math.prime} 正整数nをモンテカルロ素因数分解法 (R. P. Brent, An improved Monte Carlo factorization algorithm, BIT 20 (1980), 176-184. http://maths-people.anu.edu.au/~brent/pub/pub051.html)により 素因数分解します。

この手続きはnaive-factorizeよりも大きな数に使えます (目安としては1e20程度まで)。

アルゴリズムは確率的なので、同じnに対しても実行時間はばらつきますが、 nの素因数が全て2^64より小さければ、かならず確定的な答えを返します。

今のところ、nが2^64以上の素因数を含んでいる場合、この手続きは 永遠にそれを分割しようとしてループしてしまいます。現実的なアプリケーションは 何らかの方法で一定の時間でルーチンを中断して諦めるメカニズムが必要でしょう。 全ての入力に大して確定的な素数判定が実装されれば、この欠陥も修正されます。

その他の関数

Function: jacobi a n

{math.prime} Jacobi symbol (a/n) を計算します (http://en.wikipedia.org/wiki/Jacobi_symbol)。

Function: totient n

{math.prime} オイラーのトーシェント関数です。nは非負整数です。

現在の実装は上のmc-factorizeを使っており、 nが大きな素因数を持っている場合は非常に長い時間がかかります。


Next: , Previous: , Up: ライブラリモジュール - ユーティリティ   [Contents][Index]


For Gauche 0.9.14Search (procedure/syntax/module):