素数判定.mws


素数判定のアルゴリズム(PRIMES is in P)

2002年8月5日の週に素数判定の新しいアルゴリズムに関するニュースが世界を駆け巡りました。

インドの3人の科学者が発表した"PRIMES is in P"という論文に関するものです。

原著論文 "PRIMES in P"は下記から入手できます。
http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf

アルゴリズムは与えられた数が素数か否かを判定するものです。

Mapleではisprimeという関数があり素数かどうかの判定を行います。

素数は1あるいはその数以外の素因数(割り切れる正の整数)を持たない数です。

6は2と3で割り切れるので素数ではありません。

 

>    isprime(12011);

true

 

>    isprime(12001);

false

 

>    ifactor(12001);

``(11)*``(1091)

この一見普通のアルゴリズムがニュースになった理由は多項式時間と呼ばれるはっきり(多項式で)予測できる

時間内に判定ができるアルゴリズムであることがこの論文で証明されたからです。多くの数学者がどうやらこの

論文は正しいと公表しているようです。いままでのアルゴリズムはまだ確実に多項式時間で判定できるとは保証されていませんでした。

なぜこのことが注目をあびているかというと、素数は暗号解読(cryptography)において重要な役割を担うからです。

3人はインドのIndian Institute of Technology  (Kanpur)のManindra Agrawal(教授), Neeraj Kayal と Nitin Saxenaという人たちです。

このアルゴリズムをMapleのプログラムで作成された方(匿名希望)がプログラムを公開してくださいました。

論文の4ページ目の記述をそのままMAPLEで記述されたそうです。

以下作成者のコメントです。

「早速、Maple 8が役に立ちました。
8月8日頃にインドで発表された素数判定のアルゴリズムを
論文の通りに作ってみました。

実用的には、isprime関数の方がずっと早いのですが
よくわからない論文も数式処理で書くと、書けるようです。
論文の4ページ目の記述をそのままMAPLEで記述しました。

1番目のループは、判定対象が平方数かどうかのチェックです。

この方式の効果が出てくるのは12000付近の素数からで
while loopで rが1つ見つかります。
nが十分大きいなら、もっと効果が出るはずです。

2番目のfor loopは二項定理を使って工夫すれば
もっと速くなるはずですが、Mapleの数式処理機能で
そのまま処理できるところが、面白いですね。
ただし、メモリを異常に消費するので、512Mの
マシンでも13000くらいの素数の判定が限界でした。
(x-a)^n を n=13000でexpand処理するわけですから
できたこと自体が驚きです。

参考までに、ソースを添付します。
匿名扱いで転載自由(ただし、正しいかどうかは
保証できません (^ ^;))
利用する場合は、専門の方にチェックしてもらって下さい。」

 

>    restart;

 

>    n:=nextprime(100);

n := 101

 

>       isprime(n);

true

 

>     lnN:=(ln(n));

lnN := ln(101)

 

>    evalf(trunc(lnN/ln(2)));

6.

 

>    check_p:="next":

 

>     for b from 2 to trunc(evalf[100](lnN/ln(2))) do
 root_b_a:=(evalf[100](n^(1/b)));
 nx:=ceil(root_b_a)-floor(root_b_a);
 if round(nx*10e100)=0.0 then check_p:='Compsit'; break;
 fi;
od:

 

>      print(b,check_p);

7, next

 

>     for r from 2  while r < n do
 #print(r);
 if gcd(n,r)<>1 then  check_p:="Compsit(gcd)"; break;
 fi;
    if isprime(r) then
         en:=nops(ifactor(r-1) );
         Llist:=seq(op(1,(op(j,ifactor(r-1)))),j=1..en);
         en2:=nops([Llist]);
         op(1,(op(1,[Llist])));
        LLlist:=seq(op(1,(op(j,[Llist]))),j=1..en2);
         q:=max(LLlist);
         #print(r,q,evalf(4*sqrt(r)*log[2](n)),(n^((r-1)/q)) mod r );
      if q >= evalf(4*sqrt(r)*log[2](n))
         and ((n^((r-1)/q)) mod r <> 1) then  break;  
      fi;
    fi;
 od;

 

>     print(r,check_p);

101, next

 

>    for a from 1  to trunc(2*sqrt(r)*log[2](n)) do
 # print(a);
 s:=rem((x^n-a),(x^r-1),x) mod n; t:=rem((x-a)^n,(x^r-1),x) mod n;
 if s <> t then check_p:="Compsit"; break;
 fi;
od:

 

>     if check_p="next" then print('Prime!');
 fi;
  

Prime!

お送り頂いたプログラムでは最初の

n:=nextprime(100);

のところは

n:=nextprime(12110);

となっていましたが計算時間がかかるため100に変更しています。

ご研究にぜひMapleをお役立てください。