Miller-Rabin & Pollard-Rho

Miller - Rabin 的 用途 : 快速判断一个 1e18 级的数是不是质数

如果试除的话复杂度较高, 于是就有了 Miller - Rabin

首先考虑 费马小定理   , 如果 2 ^ (p-1) 次方 模 p 为 1 可不可以说明 p 是质数呢

当然不行, 因为有反例 ---- 伪素数  

这样的话出错率很高, 于是我们考虑继续以3为底数, 5为底数测试

但是还是有一些极端的伪素数, 将小于它的底数试完过后还是未能将它筛出, 例如 561

以下来自 https://www.cnblogs.com/Norlan/p/5350243.html

Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。

新的测试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1, 那么要么x=1，要么x=p-1。

这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。

由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时 x=p-1)。
    

2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，面具后面真实的嘴脸显现了出来.
    

这就 是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。

那么我们需要计算的东西就 变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。

如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足a^(d * 2^i) mod n = n-1

或者最后指数中的2用完了得到的a^d mod n=1或n-1。

也就是说, 尽可能提取因子2， 把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1

或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) （注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）

---

代码 : 选用 2, 3, 7, 61, 24251, 1e16 范围内只有一个46856248255981不满足

bool test(ll n, ll p, ll d){
	while(!(d&1)) d>>=1;
	ll t = fpow(p, d, n);
	while(d != n-1 && t != n-1 && t != 1){
		t = mul(t, t, n); d<<=1;
	} return t == n-1 || (d&1);
}
bool Miller(ll x){
	if(x == 46856248255981ll || x<2) return false;
    if(x==2 || x==3 || x==7 || x==61 || x==24251) return true;
    if(!(x%2) || !(x%3) || !(x%7) || !(x%61) || !(x%24251)) return false;
	if(!test(x, 2, x-1)) return false;
	if(!test(x, 61, x-1)) return false;
	return true;
}

---

Pollard-Rho : P4718 【模板】Pollard-Rho算法

我们可以每次猜一个数, 如果它跟 n 的gcd 不为1, 那么它就是 n 的一个因子

但猜中的概率可能有点小, 但我们每次随记两个数作差, 再运用一些随机的玄学算法, 就可以比较客观地猜出一个因子

来自luogu题解

用一个悖论来提高成功率

随机算法的优化 

 

Pollard Rho 和他的随机函数 

每当我们猜到一个质因子 t时, 我们递归求解 t 与 n/t, 知道 Miller-Rabin 测出来是质数 

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll read(){
	ll cnt = 0; char ch = 0;
	while(!isdigit(ch)) ch = getchar();
	while(isdigit(ch)) cnt = (cnt<<1) + (cnt<<3) + (ch-'0'), ch = getchar();
	return cnt;
}
int T; ll n, Ans, tot;
ll gcd(ll a, ll b){ return !b ? a : gcd(b, a%b);}
ll mul(ll a, ll b, ll p){
    ll d = ((long double) a / p * b + 1e-8);
    ll r = a * b - d * p;
    return r<0 ? r + p : r;
} 
ll fpow(ll a, ll b, ll p){
	ll ans = 1; for(;b;b>>=1){
		if(b&1) ans = mul(ans, a, p);
		a = mul(a, a, p);
	} return ans;
}
bool test(ll n, ll p, ll d){
	if(!(n&1)) return false;
	while(!(d&1)) d>>=1;
	ll t = fpow(p, d, n);
	while(d != n-1 && t != n-1 && t != 1){
		t = mul(t, t, n); d<<=1;
	} return t == n-1 || (d&1);
}
bool Miller(ll x){
	if(x == 46856248255981ll || x<2) return false;
    if(x==2 || x==3 || x==7 || x==61 || x==24251) return true;
    if(!(x%2) || !(x%3) || !(x%7) || !(x%61) || !(x%24251)) return false;
	if(!test(x, 2, x-1)) return false;
	if(!test(x, 61, x-1)) return false;
	return true;
}
ll calc(ll x){
	ll c = rand() % (x-1) + 1, n = 0, m = 0, g = 1, q = 1; 
	for(ll k=2; ;k<<=1, m = n, q = 1){
		for(ll i = 1; i <= k; i++){
			n = (mul(n, n, x) + c) % x;
			q = mul(q, abs(m-n), x);
		} g = gcd(q, x); if(g > 1) return g;
	}
}
void Pollard_Rho(ll n){
	if(n == 1) return;
	if(Miller(n)){ Ans = max(Ans, n); tot++; return;}
	ll t = n; while(t == n) t = calc(n);
	Pollard_Rho(t); Pollard_Rho(n/t);
}
int main(){
	srand(time(0)); 
	T = read();
	while(T--){
		n = read(); tot = Ans = 0; Pollard_Rho(n);
		if(tot == 1) printf("Prime\n");
		else{ printf("%lld\n", Ans);}
	} return 0;
}