其中c为 b的二进制数的后缀 (bi-1,bi-2,….b0)对应的十进制数,当c成倍增加时,算法保持条件 d=a^c mod n不变,直至 c=b.
下面为程序分析,函数modular_exp(long a,long b,long n)输入底数a,次幂b和模n后,通过反复平方法计算和返回a^bmod n的值.
long modular_exp(long a,long b,long n)//d≡a^b mod n
{
long d=1;
long t=a;
while(b>0)
{
if(b%2==1)
d=(d*t)%n;
b=b/2;
t=(t*t)%n;
}
return d;
}
二次探测定理 如果p是一个素数,且0<x<p,则方程x*x≡1(mod p)的解为x=1,p-1.
事实上, x*x≡1(mod p)等价于 x*x-1≡0(mod p).由此可知;
(x-1)(x+1) ≡1(mod p)
故p必须整除x-1或x+1.由p是素数且 0<x<p,推出x=1或x=p-1.
利用二次探测定理,我们可以在利用费马小定理计算 a^(n-1) mod n的过程中增加对于整数n的二次探测.一旦发现违背二次探测条件,即可得出n不是素数的结论.
下面的算法power用于计算a^p mod n,并在计算过程中实施对n的二次探测.
void power(unsigned long a,unsigned long p,unsigned long n,unsigned long &result,bool &composite)
//计算 a^p mod n,并实施对n的二次探测
{
unsigned long x;
if(p==0) result=1;
else
{
power(a,p/2,n,x,composite);//递归计算
result=(x*x)%n;//二次探测
if((result==1)&&(x!=1)&&(x!=n-1))
composite=true;
if((p%2)==1) //p是奇数
result=(result*a)%n;
}
}
在算法power的基础上,可设计Miller_Rabin素数测试的算法如下:
bool Miller_Rabin(unsigned long n)
{
RandomNumber rnd;
unsigned long a,result;
bool composite=false;
a=rnd.Random(n-3)+2;
power(a,n-1,n,result,composite);
if(composite||(result!=1)) return false;
else return true;
}
上述算法返回false时,整数n一定是一个合数,而当返回值为true时,整数n在高概率意义下是一个素数.仍然可能存在合数n,对于随机选取的基数a,算法返回true.但对于上述算法的深入分析表明,当n充分大时,这样的基数a不超过(n-9)/4个. Miller_Rabin算法的错误概率可通过多次重复调用而迅速降低.重复调用k次的Miller_Rabin算法可描述如下:
bool Miller_Rabin(unsigned long n,unsigned int k)
//重复k次调用
{
RandomNumber rnd;
unsigned long a,result;
bool composite=false;
for(int i=1;i<=k;i++)
{
a=rnd.Random(n-3)+2;
power(a,n-1,n,result,composite);
if(composite||(result!=1)) return false;
}
return true;
}
分析得上述算法的错误概率不超过
,这是一个很保守的估计,实际使用的效果要好得多.
附:上述程序中的类RandomNumber定义如下:
class RandomNumber
{
private ://当前种子
unsigned long randSeed;
public:
//构造函数,缺省值0表示由系统自动产生种子
RandomNumber( unsigned long s=0)
//产生0~n-1之间的随机整数
unsigned long Random( unsigned long n)
}
RandomNumber::RandomNumber(unsigned long s)
{
if(s==0)
{
randSeed= (unsigned long)time(0);
}
else
{
randSeed=s;
}
}
unsigned long RandomNumber:Random(unsigned long n)
{
randSeed=multiplier * randSeed + adder;
return randSeed % n;
}