基于AKS素性检测的素数生成器

1.      解题思路

AKS算法整体包括六个步骤，它本身相当于一个过滤器，只有经过各种条件筛选以后的数才是素数。步骤如图 1所示。

1.  Step1就是判断n是否是一个数的幂次方形式，如果是的话就表明n是合数，程序结束。
2.  Step2需要找出比(log₂n)²大的最小的欧拉函数r。
3.   Step3是找出gcd(a,n)，如果存在a<=r的情况，则输出合数
4.  Step4表明如果n<=r，则输出素数。
5.  Step5是最麻烦的一部，其中涉及到了多项式取模运算，对于多项式的取模运算我们这里使用了一个叫做NTL的C++类库。对一个多项式取模结果进行判断，如果对于某个小范围的a，只要存在多项式，模运算结果不等的情况，则输出宿舍。
6.    Step6表示经过了前面多个step的过滤，依然存在的数就是素数。

程序源代码如下：（用到了NTL类库，使用方法参考前面的几篇博客）

#include<math.h>
#include<time.h>
#include<stdio.h>
#include<stdlib.h>
#include<iostream>
#include <NTL/ZZ.h>
#include <NTL/ZZ_p.h>
#include <NTL/ZZ_pX.h>
NTL_CLIENT
using namespace std;
//根据课本中的伪代码，使用Euclid算法求最大公约数。
int Euclid(int a,int b)//a>=b>=0
{
    if(a>=b)
    {
        if(b==0)
            return a;
        return Euclid(b,a%b);
    }
    else
    {
        if(a==0)
            return b;
        return Euclid(a,b%a);
    }
}

/*
摘自维基百科：
In number theory, given an integer a and a positive integer n with gcd(a,n) = 1, 
the multiplicative order of a modulo n is the smallest positive integer k with a^k ≡ 1 (modulo n). 
所以k>=1，且(a^k-1)%n==0.
For example, to determine the multiplicative order of 4 modulo 7, 
we compute 4^2 = 16 ≡ 2 (mod 7) and 4^3 = 64 ≡ 1 (mod 7), so ord7(4) = 3.
a^b%r=(a%r)^b%r
A(k)=a^k%r=(a*a^(k-1))%r=(a%r*a^(k-1)%r)%r=(A(k-1)*a%r)%r
*/
int Multiplicative_Order2(int n,int r)//对于大数n此方法无法使用。
{
    int k=1;
    while(true)
    {
        int temp=pow((double)n,k);//double pow(double,int),因此这里必须进行强制类型转换，不然编译无法通过。
        if((temp-1)%r==0)
            break;
        //if(((pow(n,k)-1)%r)==0)
        //    break;
        else
            k++;
    }
    return k;
}
int Multiplicative_Order(int n,int r)
{
    if(Euclid(n,r)!=1)//如果n,r不互质，不存在Multiplicative_Order
        return -1;
    else
    {
        int k=1,temp=1;
        while(true)
        {
//根据A(k)=a^k%r=(a*a^(k-1))%r=(a%r*a^(k-1)%r)%r=(A(k-1)*a%r)%r
            temp=(temp*(n%r))%r;
            if(temp==1)
                break;
            else
                k++;
        }
        return k;
    }
    
}

long MyRand(int n)//产生随机数算法。最高时9位。超过9位时会出现负数。
{
    static int inited=0;
    long X,k=1;
    int i;
    if(!inited)
    {
        srand((unsigned)time(0));
        inited=1;
    }
    if(n==1)
        return rand()%10;
    else if(n==2)
        return (rand()%9+1)*10+rand()%10;
    for(i=0;i<n/2+1;i++)
        k*=10;
    for(i=n/2+1;i<n-1;i++)
        k*=10;
    X=(rand()%9+1)*k+((long)(rand()%k)*(rand()%k))%k;
    if(X%2==0)//如果产生的随机数是偶数，那么+1使之成为奇数。偶数肯定不是素数。
        X=X+1;
    return X;
}
//step1：
//If n = a^b for integers a > 0 and b > 1, output composite.
//a,b为整数，且b>1，那么b>=2，那么a<=sqrt(n)
bool IsPower(int num)
{
    if(num==1)//如果num=1,那么1=1^x成立，其中x为任意大于1的整数。
        return true;
    else
    {
        for(int a=2;a<=floor(sqrt((double)num));++a)//判断n = a^b是否成立只需要让n不断除以a,若最后结果为1则表明成立。
        {
            int temp=num;
            while(1)
            {
                if(temp%a!=0)//出口1:判断num能够整除a，如果不能整除，n = a^b肯定不成立
                    break;
                else//如果能整除则除以a以后再去判断能够整除a。
                    temp=temp/a;
                if(temp==1)//出口2:temp=1
                    return true;
            }
        }
    }
    return false;
}
//step2：
//Find the smallest r such that or(n) > log2(n).
int get_Smallest_r(int n)
{
    int r=2;
    while(true)
    {
        if(Euclid(n,r)==1)
        {
            int multi_Order= Multiplicative_Order(n,r);
            //cout<<n<<" "<<r<<" "<<multi_Order<<endl;
            if(multi_Order>pow(log((double)n),2))
                return r;
                //break;
        }
        r++;
    }
    //return r;
}

//step3：
//If 1 < gcd(a,n) < n for some a ≤ r, output composite.
bool a_r(int n)//返回true表示合数，返回false表示素数。
{
    bool flag=false;
    int r=get_Smallest_r(n);
    for(int a=1;a<=r;a++)
    {
        if(Euclid(a,n)>1&&Euclid(a,n)<n)
        {
            flag=true;
            break;
            //return true;//output composite.
        }
    }
    return flag;
}
/*欧拉函数
摘自维基百科：在数论中，对正整数n，欧拉函数是小于或等于n的正整数中与n互质的数的数目。
例如Euler_Totient(8)=4，因为1,3,5,7均和8互质。
*/
int Euler_Totient(int r)//r>=1
{
    int count=0;
    int i;
    for(i=1;i<=r;i++)
        if(Euclid(i,r)==1)//如果i和r互质
        {
        //    cout<<i<<endl;
            count++;
        }
    return count;
}

int Congruence ( long a , ZZ n , ZZ r ) //step5
{
    ZZ_p::init(n); //mod n
    ZZ_pX b = ZZ_pX(to_long(r),1)- 1;// b = xˆr−1;
    ZZ_pX c = ZZ_pX( 1 , 1 ) - a;// c = x−a ;
    ZZ_pX f = PowerMod( c , n , b) ;// f=(x−a ) ˆn mod c , n which is the RHS
    ZZ_pX e = ZZ_pX( 1 , 1 ) ;
    ZZ_pX g = PowerMod( e , n , b) ;// xˆn mod b , n
    g = g-a ;// g1 = xˆn−a mod c , n .
    if ( f==g )
    {
        //cout<<"prime"<<endl;
        return  1 ; 
    }
    else 
    {
        //cout<<"not prime"<<endl;
        return  0  ; 
    }
}

//返回true表示素数，返回false表示合数。
bool AKS_IsPrimality(int num)
{
        if(IsPower(num))//step1
            return false;
        int r=get_Smallest_r(num);//step2
        if(a_r(num))//step3,存在a ≤ r，使得1 < gcd(a,n) < n 成立。返回合数,
            return false;
        if(num<=r)//step4
                return true;
        //step5
        long a;
        long b=floor(sqrt((double)Euler_Totient(r))*log((double)num));
        for(a=1;a<=b;a++)//对于区间的所有数，如果都不满足条件，则返回合数
        {
            int f= Congruence (a, to_ZZ(num),to_ZZ(r) ) ;
            if(f==0)
            {
            //    cout<<"end in step5"<<endl;
                return false;
            }        
        }
        return true;//step6
}
int main()
{
    //主程序测试
    int n,i;
    srand((unsigned)time(0));
    cout<<"随机数位数:";
    cin>>n;
    for(i=0;i<20;i++)
    {
        long mr=MyRand(n);
        bool flag=AKS_IsPrimality(mr);
        if(flag)
            cout<<mr<<" is prime"<<endl;
        else
            cout<<mr<<" is not prime"<<endl;
    }
cin>>n;
    
    /*测试：int Congruence ( long a , ZZ n , ZZ r )
    //88903,5
    ZZ a, b, c;
    cout<<"a=";cin >> a;//n
    cout<<"b=";cin >> b;//r
    for(c=1;c<b;c++)
        int f= Congruence ( to_long(c) , a, b ) ;
    */

/*
    int n=1342342;
    bool flag=AKS_IsPrimality(n);
    if(flag)
        cout<<n<<" is prime"<<endl;
    else
        cout<<n<<" is not prime"<<endl;
*/
    /*测试：int Euclid(int a,int b)
    int a=7001,b=54643;
    int c=Euclid(a,b);
    cout<<c;
    */
    /*测试：int Multiplicative_Order(int n,int r)
    int a=9,b=7;
    int c= Multiplicative_Order(a,b);
    cout<<c;
    */
    /*测试：long MyRand(int n)
    int n,i;
    srand((unsigned)time(0));
    scanf("%d",&n);
    for(i=0;i<20;i++)
    printf("%ld/n",MyRand(n));
    */
    /*测试：bool IsPower(int num)        
    int n=10001;
    if(IsPower(n))
        cout<<n<<" is power"<<endl;
    else
        cout<<n<<" is not power"<<endl;
*/
    /*测试：int get_Smallest_r(int n)
    int n=7000;
    cout<<"Smallest_r:"<<get_Smallest_r(n)<<endl;
    */
    /*测试：int Euler_Totient(int r)//r>=1
    int r=8;
    cout<< Euler_Totient(r)<<endl;
    */
    return 0;
}