1.解题思路
AKS算法整体包括六个步骤,它本身相当于一个过滤器,只有经过各种条件筛选以后的数才是素数。步骤如图1所示。
- Step1就是判断n是否是一个数的幂次方形式,如果是的话就表明n是合数,程序结束。
- Step2需要找出比(log2n)2大的最小的欧拉函数r。
- Step3是找出gcd(a,n),如果存在a<=r的情况,则输出合数
- Step4表明如果n<=r,则输出素数。
- Step5是最麻烦的一部,其中涉及到了多项式取模运算,对于多项式的取模运算我们这里使用了一个叫做NTL的C++类库。对一个多项式取模结果进行判断,如果对于某个小范围的a,只要存在多项式,模运算结果不等的情况,则输出宿舍。
- Step6表示经过了前面多个step的过滤,依然存在的数就是素数。
程序源代码如下:(用到了NTL类库,使用方法参考前面的几篇博客)
View Code
#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;
}
原文链接: https://www.cnblogs.com/xwdreamer/archive/2011/03/28/2297035.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/22987
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!