伪Miller Rabin—— JZOJ(GMOJ)4278. 【NOIP2015模拟10.29B组】质数

伪Miller Rabin——

JZOJ(GMOJ)4278. 【NOIP2015模拟10.29B组】质数

Description

定义质数为因数只含1和其本身的数，对于N组询问，试判断每个数是否为素数。

Input

第一行一个正整数N，表示有N组询问。
接下来N行，每行一个正整数M，表示询问M是否为质数。

Output

输出N行，每行一个字符串。
若是质数则输出‘Prime’，若不是质数则输出‘Not prime’。

Sample Input

5
2
10
89807289
9032482948
1000000007

Sample Output

Prime
Not prime
Not prime
Not prime
Prime
样例解释：
10=2*5
89807289=3 * 11 * 11 * 13 * 19031
9032482948=2 * 2 * 439 * 5143783

Data Constraint

20%的数据满足N≤100，1< M≤800,000。
50%的数据满足N≤1,000，1< M≤100,000,000。
100%的数据满足N≤1,000，1< M≤1,000,000,000,000.
显然，\(O(n\sqrt m)\)一定会T飞
1< M≤1,000,000,000,000.

Time Limits:2000 ms Memory Limits: 262144 KB

所以这个世界才出现了各种玄学算法

如 Miller Rabin

但我只会伪Miller Rabin

前置知识

1.费马小定理
\(x^{p-1}\equiv 1 (mod\; p)\)

若 有很多的x 都满足，那p就极有可能
是一个质数
反之 如果有 \(x \ne p\)且不满足，一定是合数

2.分块乘

我只会分块乘,快速乘的正确性据说不能保证?
\((ax+b)(cx+d)mod\; p=(acx^2+adx+bcx+bd)mod\;p\)
分块乘的极限是\(10^{13}\)左右

**豪同学狂言称霸10^18次方，
一到 某谷 某题上就萎了
如果用
2,3,7,43,61,24251,19491001
这组质数
即使不增加二次探测的加持
貌似\(10^{16}\)以内找不到伪素数

代码

//M<=10^12
#include<cstdio>
#include<cstring>
using namespace std;
#define Np 8
int wp[Np]={0,2,3,7,43,61,24251,19491001};
#define tmi 1048575ll
#define tty 1048576ll
#define tmr 1099511627776ll
inline long long supx(long long a,long long b,long long mod)
{
    long long a0=a&tmi,a1=a>>20;
    long long b0=b&tmi,b1=b>>20;
    long long ans=((a1%mod)%mod*(b1*tmr%mod)%mod+a0*b0%mod)%mod;
    ans=(ans+(a0*(b1%mod)%mod*(tty%mod)%mod)%mod)%mod;
    ans=(ans+(b0*(a1%mod)%mod*(tty%mod)%mod)%mod)%mod;
    return ans;
}
inline long long qpow(long long a,long long x,long long mod)
{
    long long ans=1,apow=a;
    while(x>0)
    {
        if(x&1) ans=supx(ans,apow,mod);
        apow=supx(apow,apow,mod);
    //      printf("  %lld %lld\n",apow,ans);
        x>>=1;
    }
    return ans;
}
int main()
{
    freopen("prime.in","r",stdin);
    freopen("prime.out","w",stdout);
    int n,i,j;long long wpow,m;
    bool ans;
    scanf("%d",&n);
    for(i=1;i<=n;i++)
    {
        scanf("%lld",&m);
        ans=1;
        for(j=1;j<Np;j++)
        {
            if(wp[j]!=m)
            {
                wpow=qpow(wp[j],m-1,m);
                if(wpow!=1){ans=0;break;}
            }
            else break;
        }
        if(ans) printf("Prime\n");
        else printf("Not prime\n");
    }
    return 0;
}

2020/08/11更新

如果需要进一步加大适用的范围
我们需要将分块乘替换掉,

龟速乘

long long mul(long long x,long long y,long long mod)
{
	long long ans=0;
	while(x>0)
	{
		if(x&1) ans=(ans+y)%mod;
		y=(y<<1)%mod;
		x>>=1;
	}
	return ans;
}

适用范围\(10^{18}\)以内
但复杂度要带个\(\log\)会稍微慢一些

二次探测定理

若有 \(x^2\equiv 1(\mod p)\)
则 \(x^2-1\equiv 0(\mod p)\)
\((x+1)(x-1)\equiv 0(\mod p)\)
若\(x\not=1\)且\(x\not=p-1\)
p不为质数

完整代码

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <ctime>
using namespace std;
long long mul(long long x,long long y,long long mod)
{
	long long ans=0;
	while(x>0)
	{
		if(x&1) ans=(ans+y)%mod;
		y=(y<<1)%mod;
		x>>=1;
	}
	return ans;
}
long long qpow(long long a,long long d,long long mod)
{
	long long s=a,ans=1;
	while(d>0)
	{
		if(d&1) ans=mul(ans,s,mod);
		d>>=1;s=mul(s,s,mod);
	}
	return ans;
}
int pd[7]={2,3,7,43,61,24251,19491001};
bool Miller_Rabin(long long p)
{
	bool bz=1;int tp=0,i,j;
	long long a,b,x=p-1;
	while(~x&1) x>>=1,tp++;
	if(tp==0&&p!=2) return 0;
	for(i=0;i<7;i++)
	{
		if(p==pd[i]) return 1;
		a=pd[i];
		a=qpow(a,x,p);
		for(j=1;j<=tp;j++)
		{
			b=mul(a,a,p);
			if(b==1&&(a!=1&&a!=p-1)) return 0;
			a=b;
		}
		if(a!=1) return 0;
	}
	return 1;
}
int main()
{
	freopen("prime.in","r",stdin);
	freopen("prime.out","w",stdout);
	srand((int) time(NULL));
	long long m;int n,i;
	scanf("%d",&n);
	for(i=1;i<=n;i++)
	{
		scanf("%lld",&m);
		if(Miller_Rabin(m)) printf("Prime\n");
		else printf("Not prime\n");
	}
	return 0;
}