素数の求め方。

このエントリーの最終目標
NScripterで、is_prime命令を実装する。
is_prime命令とは
「is_prime %0,%1」と記述した時、%1が素数であれば、%0に1を代入し、そうでなければ0を代入して返す命令。
おまけ目標
フラグ反転したis_not_prime(no_prime?)命令もついでに。

素数の性質。おさらい

  • 2以上の自然数に含まれる。
  • 1と、その数自身でしかきっちり割り切れない。

素数の判定方法

有名な物に「エラトステネスの篩」がある。これを逆に考えると、このように導きだされる。

ある自然数nが素数かどうか判定するには、2〜n-1までの数で、nを割ってみて、一度でも余りが0になれば素数ではないことがわかる。

数学者ならこれでもいいのかも知れないが、しかしこれをそのままNScriptで実装すると、多分死ぬほど実行時間がかかって使い物にならないと思われる。
そこで、計算量を減らすための工夫が必要になる。
(以下、現時点での実装)

numalias prime_try,%0:inc %0
numalias prime_flag,%0:inc %0
numalias prime_loop,%0:inc %0

defsub is_prime
defsub is_not_prime
defsub no_prime

======

*is_not_prime
*no_prime
getparam i%prime_flag,%prime_try
is_prime %%prime_flag,%prime_try
mov %%prime_flag,1-%%prime_flag
return

*is_prime
getparam i%prime_flag,%prime_try
mov %%prime_flag,0
if %prime_try<2 return ; 2以上の自然数だっつーとろーが!
mov %prime_loop,2
~
mov %%prime_flag,%prime_try mod %prime_loop
if %%prime_flag=0 return ; 割り切れちゃった。ここで終わり。
inc %prime_loop ; カウントアップ
notif %prime_loop=%prime_try jumpb ; ループ判定
mov %%prime_flag,1
return

ダイエットする。その1

「2〜n-1」までの数と言うのが怪しい。
これ、よく考えてみれば、「2〜nの平方根」までで足りることがわかる。証明はちょっとパスさせてもらうが。数学者じゃないし。
nが素数ではないと仮定すると、それは素数の積で表現できることになる。その場合、素因数の最大値は、nの平方根以下になる。
(以下、現時点での実装)

numalias prime_try,%0:inc %0
numalias prime_flag,%0:inc %0
numalias prime_loop,%0:inc %0

defsub is_prime
defsub is_not_prime
defsub no_prime

======

*is_not_prime
*no_prime
getparam i%prime_flag,%prime_try
is_prime %%prime_flag,%prime_try
mov %%prime_flag,1-%%prime_flag
return

*is_prime
getparam i%prime_flag,%prime_try
mov %%prime_flag,0
if %prime_try<2 return ; 2以上の自然数だっつーとろーが!
mov %prime_loop,2
~
mov %%prime_flag,%prime_try mod %prime_loop
if %%prime_flag=0 return ; 割り切れちゃった。ここで終わり。
inc %prime_loop ; カウントアップ
notif %prime_loop*%prime_loop>%prime_try jumpb ; ループ判定
mov %%prime_flag,1
return

ダイエットする。その2

また、2で割り切れるかどうかをチェックした後、4で割り切れるかどうかをチェックするのは時間の無駄である。つまり、「2〜nの平方根の範囲の素数」で順に割って行くのが一番いいと言うことになる。
そうなると、ある程度の素数表を先に作っておけばいいと言うことに気づく。
その素数表の最後の素数、その二乗の数値までは判定できるプログラムができる。
(以下、現時点での実装)

numalias prime_try,%0:inc %0
numalias prime_flag,%0:inc %0
numalias prime_loop,%0:inc %0
numalias prime_table,%0:inc %0
numalias prime_table_len,100 ; 仮にテーブルの長さを100とする。
dim ?prime_table[100]
mov ?prime_table[0],2
mov ?prime_table[1],3

defsub is_prime
defsub is_not_prime
defsub no_prime

; 素数テーブルの作成
mov %prime_try,3
for %1=2 to prime_table_len
	~
	add %prime_try,2
	is_prime %2,%prime_try
	if %2=0 jumpb
	mov ?prime_table[%1],%prime_try
next

======

*is_not_prime
*no_prime
getparam i%prime_flag,%prime_try
is_prime %%prime_flag,%prime_try
mov %%prime_flag,1-%%prime_flag
return

*is_prime
getparam i%prime_flag,%prime_try
mov %%prime_flag,1
if %prime_try<2 mov %%prime_flag,0:return ; 2以上の自然数だっつーとろーが!
if %prime_try=2 return ; 特別処理。
if %prime_try=3 return ; 特別処理。
mov %prime_table,0
~
notif %prime_table>prime_table_len mov %prime_loop,?prime_table[%prime_table]
if %prime_table>prime_table_len add %prime_loop,2
if %prime_loop*%prime_loop>%prime_try jumpf ; ループ判定
mov %%prime_flag,%prime_try mod %prime_loop
if %%prime_flag=0 return ; 割り切れちゃった。ここで終わり。
inc %prime_table ; カウントアップ
jumpb
~
mov %%prime_flag,1
return

最終的に。

素数テーブルを使った構造上、ある数値(素数テーブルの最大値の二乗)を超える自然数を与えた場合は、素数テーブルの最大値以降、判定対象の平方根の範囲の総ての奇数で割り切れるかどうかのチェックをするため、速度がさらに遅くなる仕様となった。
素数テーブルの長さを伸ばせば、この速度の遅くなる境界はさらに大きな数値になるが、素数テーブルの長さが長すぎればそれはそれでメモリーを圧迫する。
適宜、prime_table_lenの長さを調整したいところ。

追記

Q
NScripter素数が扱えるからって、だからどうなの?
A
……まあ、その。素因数分解をするのに必要だったから作ったんですよ。ええ。あ、あと、ブッチ神父のセリフを自動生成するとか!