[LOJ143]质数判定(Miller-rabin素数测试+O(1)快速乘讲解)

题面

题解

Miller-rabin素数测试

是一种随机化算法，能够在较短的时间内判断出一个数是合数，还是很可能为素数。其出错的概率极小；当用于测试的p取遍前10个素数时，在(3e18)的范围内不会出错。

原理：1、费马小定理

p为素数的必要不充分条件是(forall (a,p)=1,a^{p-1}{equiv}1{mod p})。证明不再赘述，可以参见百度百科或《数学奥林匹克命题人讲座：初等数论》2.3节。

2、二次探测

只使用费马小定理时，出错概率仍然较大，尤其是对于Carmichael数n（在(1e8)以下有255个，最著名的、最小的Carmichael数是561），对于所有的与n互质的正整数a，都可以满足(a^{n-1}{equiv}1{mod n})。因此，仅仅使用Fermat测试是不够的。

二次探测的原理是：如果(x^2{equiv}1{mod p})，p为素数，那么有(x{equiv}{pm}1{mod p})。证明显然，将1移项至左边，然后因式分解即可。

具体做法

设待测数为mod，若mod为偶数则可以直接判断。否则，我们将mod-1分解为(2^t*n)，其中(2{ mid}n)。然后，取与mod互质的数p，计算

[p^n、 p^{2n}、 p^{2^2n} … p^{2^tn}{\% mod} ]

从第二个数开始，每一个数可以通过前一个数平方得到。

如果最后一项不为1，那么根据原理1，p不是素数；
如果第一项就为1，那么此次测试失败，需要更换p。但是，这样的概率极小。
如果第一项不为1，且最后一项为1：我们寻找这个数列中的第一个1，其前一项如果不是n-1，那么根据原理2，p不是素数。

对于特定选取的p，经过一次test(p,mod)后，出错的概率在({frac{1}{4}})以下。证明见《算法导论》第31.8节。

在待测数(mod{leq}3e18)的情况下，选取p为2,3,5,7,11,13,17,19,23这前9个素数，可以保证不出错。

对于一个待测数mod，Miller-rabin素数测试的时间复杂度为(O(klogmod))。（其中k为测试的轮数）

防止相乘爆long long的小技巧（O(1)快速乘）

此方法适用范围：x,y,mod均在(1e18)级别。

#define ll long long
#define ld long double

inline void Adjust(ll &x,ll mod){
    x = (x % mod + mod) % mod;
}

inline void Tms(ll &x,ll y,ll mod){
    x = x * y - (ll)((ld)x * y / mod) * mod;
    Adjust(x,mod);
}

在Miller-rabin素数测试中，由于数据范围很大，会碰到求xy%mod（其中x,y,mod在1e18级别)的问题，此时如果直接乘会爆long long。为了解决这个问题，我们可以利用自然溢出：若两个long long型变量(x，y)，进行了某种运算({igodot})（包括加、减、乘）后超出了([-2^{63},2^{63}))的范围，那么返回值是使得

[x{igodot}y=2^{64}*t+r(-2^{63}{leq}r<2^{63}) ]

的r。在Tms()中，(x*y)和((ll)((ld)x * y / mod) * mod)均属于这种情况。二者的相减也是。

我们设（假设没有溢出）(x*y-(ll)((ld)x*y/mod)*mod)的值为r'，我们真实想要的数为r。由于long double的运算可能产生误差，所以实际上(r')可能等于r或者(r{pm}mod)。但是无论如何，(r')是([-2^{61},2^{61}])内，与r关于mod同余的一个数。

而考虑到溢出，真实情况下，(x*y-(ll)((ld)x*y/mod)*mod)的返回值应该是与(r')关于(2^{64})同余的、在([-2^{63},2^{63}))内的一个数(s)。但是，我们发现(r')不管是加上还是减去(2^{64})，都超出了([-2^{63},2^{63}))，因此，一定有(s=r')。

然后为了消除可能有的误差，再进行一次Adjust操作即可。

代码

#include<bits/stdc++.h>

using namespace std;

#define ll long long 
#define ld long double
#define rg register

namespace ModCalc{
	inline void Adjust(ll &x,ll mod){
		x = (x % mod + mod) % mod;
	}
	
	inline ll Check(ll x,ll mod){
		Adjust(x,mod);return x;
	}
	
	inline void Tms(ll &x,ll y,ll mod){
		x = Check(x * y - (ll)((ld)x*y/mod) * mod,mod); 
	}
	
	inline ll Mul(ll x,ll y,ll mod){
		Tms(x,y,mod);return x;
	}
}
using namespace ModCalc;

ll pri[9] = {2,3,5,7,11,13,17,19,23};

inline ll power(ll a,ll n,ll mod){
	ll x = a,s = 1;
	while(n){
		if(n & 1)Tms(s,x,mod);
		Tms(x,x,mod);
		n >>= 1;
	}
	return s;
}

inline bool test(ll mod,ll p){
	if(mod == p)return 1; //mod=p需要特判，因为p的倍数一定不符合费马小定理
	ll n = mod - 1;
	while(!(n&1))n >>= 1;
	ll r = power(p,n,mod);
	if(r == 1)return 1; //运气不好，本轮测试失败
	while(n < mod - 1){
		ll x = Mul(r,r,mod);
		if(x == 1)return r == mod - 1;
		r = x,n <<= 1;
	}
	return 0;
} 

inline bool MR(ll mod){
	if(mod < 2)return 0;     
	if(!(mod&1))return mod == 2;
	for(rg ll i = 0;i < 9;i++)if(!test(mod,pri[i]))return 0;
	return 1;
}

int main(){
	ll mod;
	while(~scanf("%lld",&mod)){
		puts(MR(mod) ? "Y" : "N");
	} 
	return 0;
}