Miller_Rabin素数测试算法模板对比

昨天在USACO做了一道判断素数的题，就想着学习一下Miller_Rabin素数测试算法，在网上找到两种模版，第一种十分简洁，运行速度也很快，但是会判错极少的几个非素数；第二种比较麻烦，运行速度很慢，所以我便想找到第一种模版不能判断的非素数特判一下，结果用了一天，电脑只找到10^8以下的，10^9内还有2个没找到，但正确的模版运行速度太慢，我的电脑又太渣，耗不起时间了，姑且先这样，等以后有深入理解有更好的方法再更新一下。

第一种：源自吉林大学ACM模版

刚开始用的是随机数测试，我想到以前了解过只要用2，7，61便可能判断unsigned内的数，便做了修改，交了USACO，由于特判了5，7，11，竟然水过了，然后就愉快地去完了。。。

后来做另一道素数题，终于暴露了它的缺点，但是感觉它判断的很快，而且只有极少的不能判断，便先找出它的特例特判一下，经过一天的无聊枚举终于找到了。

文本对比还是自己写来的快，网上下的用了结果内存爆了而且运行慢，自己写的10s出结果（50mb的文本）。

目前只针对2，3，5，找到的特例是：4097，1048577，16777217，25326001，（960946321是一个10^9内的特例，这个算法枚举一次10^9，我的渣电脑都要8min左右，别说更严密的那个了）

bool witness(long long a,long long n) {
    long long d=1,x;
    int i=ceil(log(n-1.0)/log(2.0))-1;//ceil()是向上取整
    for(;i>=0;i--) {
        x=d;
        d=(d*d)%n;
        if(d==1&&x!=1&&x!=n-1)
            return true;
        if(((n-1)&(1<<i))>0)
            d=(d*a)%n;
    }
    return d!=1;
}

bool miller_rabin(long long n) {
    int s[]={2,3,5};//只用2，3，5和下面的特例判断就能正确判断10^8以内的所有素数，枚举一遍用时40s左右
    if(n==2||n==3||n==5||n==7)
        return true;
    if(n==1||(n&1)==0||0==n%3||0==n%5||n==4097||n==1048577||n==16777217||n==25326001)
        return false;
    for(int i=0;i<3;i++)
        if(witness(s[i], n))
            return false;
    return true;
}
           

第二种模版：网上大多数用的模版（模版修改自：http://blog.csdn.net/z690933166/article/details/9860937）

就算只用2，7，61判断，也奇慢无比，枚举一次10^8要跑17min左右，但unsigned int内没有特例

typedef unsigned long long LL;

LL modular_multi(LL x,LL y,LL mo) {
	LL t;
	x%=mo;
	for(t=0;y;x=(x<<1)%mo,y>>=1)
		if (y&1)
			t=(t+x)%mo;
	return t;
}

LL modular_exp(LL num,LL t,LL mo) {
	LL ret=1,temp=num%mo;
	for(;t;t>>=1,temp=modular_multi(temp,temp,mo))
		if (t&1)
			ret=modular_multi(ret,temp,mo);
	return ret;
}

bool miller_rabin(LL n) {
	if (n==2||n==7||n==61)
        return true;
	if (n==1||(n&1)==0)
        return false;
	int t=0,num[3]={2,7,61};//2,7,61对unsigned int内的所有数够用了，最小不能判断的数为4 759 123 141；用2,3,7,61在 10^16 内唯一不能判断的数是 46 856 248 225 981
	LL a,x,y,u=n-1;
	while((u&1)==0)
        t++,u>>=1;
	for(int i=0;i<3;i++) {
		a=num[i];
		x=modular_exp(a,u,n);
		for(int j=0;j<t;j++) {
			y=modular_multi(x,x,n);
			if (y==1&&x!=1&&x!=n-1)
				return false;
            //其中用到定理，如果对模n存在1的非平凡平方根，则n是合数。
            //如果一个数x满足方程x^2≡1 (mod n),但x不等于对模n来说1的两个‘平凡’平方根：1或-1，则x是对模n来说1的非平凡平方根
			x=y;
		}
		if (x!=1)//根据费马小定理,若n是素数，有a^(n-1)≡1(mod n).因此n不可能是素数
			return false;
	}
	return true;
}