AKS.mws

AKS素数判定法のアルゴリズム

AKS素数判定法のアルゴリズムをMapleで実装しました.C言語,UBASICでも書いて
みたのですがいろいろと制限がでてきたのでMapleで書き直しました.
多項式の剰余の計算はFFTを用いずオーダーはO~(log^19 n)になってしまいましたが
どっちにしろ時間がかかってしまうので…16ビットの素数で20日ほどかかります.(PIII 1GHz)

そこで,PRIMES is in Pの6章にある予想をプログラムしてみました.
n=10^10までは検証されているそうです.

AKS(Agrawal-Kayal-Saxena)素数判定法のアルゴリズムは、2002年8月にインドの
カンプール工科大学コンピュータサイエンス学部、M Agrawal教授、N.Kayal氏、
N.Saxena氏らにより多項式時間で決定的な素数判定を行うアルゴリズムです。
http://www.tcs.hut.fi/~helger/crypto/link/primality_tests/aks.html
http://www.cse.iitk.ac.in/news/primality.html

Mapleでの実装は、http://www.luschny.de/math/primes/aks.txtや
http://www.cybernet.co.jp/maple/hiroba/prime/prime.htmlにも紹介
されておりますが、このプログラムはアルゴリズムの実装と検証を目的として、
Cで実装したものをそのままMapleに移したので、
配列の多項式の計算スペースを rのオーダーに抑えました。

変数nに判定したい値を入力すると結果を返します。

>    restart;
  #多項式の剰余の計算
  PolyMult := proc(x,y,z,n,r)
      local i,j,k,tp;
      tp := array(0..r);
      for i from 0 to r-1 do
          tp[i] := 0;
     end do;
 
     for i from 0 to r-1 do
         for j from 0 to r-1 do
             k := (i+j) mod r;
             tp[k] := (((x[i]*y[j]) mod n)+tp[k]) mod n;
         end do;
     end do;

     z:=tp;

 end proc:
 
 
 #素数判定したい数を入力.
 n:=10^20+39;
 check_p:="PRIME":
 
 #AKS予想の条件を満たす r を求める
 for r from 2 to n do
     if gcd(n,r)<>1 then  check_p:="COMPOSITE"; break;    end if;
     if ((n mod r)*(n mod r) mod r) <> 1 then break;    end if;
 end do:

 if check_p="COMPOSITE" then next; end if:
 
 #配列の生成と初期化
 coef := array(0..r):
 tmp1 := array(0..r):
 tmp2 := array(0..r):
 pow  := array(0..r):
 ctmp := array(0..r):
 
 for i from 0 to r-1 do
     coef[i]:=0;
     tmp1[i]:=0;
     tmp2[i]:=0;
     pow[i] :=0;
     ctmp[i]:=0;
 end do:
 
 ctmp[0] := 1:
 tmp1[0] := n-1:
 tmp1[1] := 1:
 tmp2[0] := n-1:
 tmp2[1] := 1:
 pow[0]  := n-1:
 pow[1]  := 1:
 
 #多項式の剰余を計算する.
 binn := n:
     while binn > 0 do
       if binn mod 2 = 1 then
             PolyMult(pow,ctmp,coef,n,r);
             for j from 0 to r-1 do
                 ctmp[j] := coef[j];
             end do;     
         end if;
         for j from 0 to r-1 do
             tmp1[j] := pow[j];
             tmp2[j] := pow[j];
         end do;    
         PolyMult(tmp1,tmp2,pow,n,r);
         binn := trunc(binn / 2);
         if binn = 0 then break;    end if;
         for j from 0 to r-2 do
             coef[j];
         end do:
     end do:
    
     coef[0]:= (coef[0]+1)mod n:
     coef[n mod r]:= (coef[n mod r]-1)mod n:
 
     for i from 0 to r-1 do
         if coef[i] <> 0 then check_p:="composite"; break;    end if;
     end do:
 
 #もしnが素数なら"PRIME"と,合成数なら"composite"と表示
     if check_p="PRIME" then
     print("PRIME");
     end if:
 
     if check_p="composite" then
     print("composite");
     end if:

n := 100000000000000000039

PRIME

>