Miller-Rabin判质数和Pollared-Rho因数分解
朴素判质数：$ 在[2..\sqrt{n}]$范围内枚举逐一判断是不是$ n$的因数
时间复杂度:$ O(\sqrt{n})$
当n达到$ 10^{18}$级别时，显然效率过低
Miller-Rabin算法
这种算法本质上是一种基于概率的素数判断方法，因为复杂度小以及有极大的正确率而常被应用
费马小定理
对于一个素数$ p$和任意正整数$ x$有$ x^{p-1} \equiv 1 (mod\ p)$
那对于这个命题的逆命题呢？有高概率成立
因而很容易想到选用几个不同的底数逐一验证，只要有一个不满足则非素数
但由于有强伪素数的存在（即对于任何底数x这种伪素数都有素数性质）使得仅如此不管用多少底数都有误差
二次探测定理
若$ p$为素数且$ x^2 \equiv 1(mod\ p)$
则有$ x \equiv 1(mod\ p)\ or\ x \equiv p-1(mod\ p)$
证明：
由$ x^2 \equiv 1(mod\ p)$可得$ (x+1)(x-1) \equiv 0(mod\ p)$
即$ p|(x+1)(x-1)$
又∵$ p$是素数因而$ p|(x+1)\ or\ p|(x-1)$
得证
应用
以测试数$ 341$，底数$ 2$为例，这个伪素数满足$ 2^{340} \equiv 1 (mod\ 341)$
那么如果$ 341$是素数，根据二次探测定理显然有$ 2^{170} \equiv 1 (mod\ 341)$或$ 2^{170} \equiv 340 (mod\ 341)$
计算发现确实如此，由于$ 170$仍为偶数，递归验证$ 2^{85}\mod\ 341$
发现结果是$ 22$，并非$ 1$或者$ 340$
如此就把$ 341$这个伪素数揪出来了
找规律发现
对于一个素数$ p$，这样递归得到的对应结果序列应为$ 1\ \ 1\ \ 1 ...p-1$，即经过$ p-1$（$ p-1$之后会无序）或全为$ 1$（即不会经过$ p-1$)对于一个非素数$ p$，递归得到的对应结果序列一般不会以$ p-1$结尾，绝大多数情况下甚至不以1开头
那么我们直接从最大的$ p-1$的因数的奇数开始（这时候模后结果由于可能在$ p-1$后因而无序）
不断自身平方直到对应结果为$ p-1$或已经跳到$ p-1$
如果经过了对应结果$ p-1$的时候则有极高概率为素数（虽然并不完全正确）
当然也有可能这个奇数的对应结果直接为$ 1$，这时候也满足素数数性质（即整个结果序列均为$ 1$）(虽然也不完全保证是素数)
这样对于一个底数的测试复杂度仅有一个$ log$
如果只使用$ 2,3,7,61$四个因数至少可以保证$ 10^9$之内不存在任何反例
如果选用$ 2,3,7,61,24251$那么在$ 10^{16}$内只存在$ 46856248255981$这一个反例
正确率高到基本能满足所有$ long\ long$范围内的素数测试
代码：（此份是int范围内的）
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define rt register int
#define ll long long
using namespace std;
ll read()
{
ll x = 0; int zf = 1; char ch;
while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar();
if (ch == '-') zf = -1, ch = getchar();
while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar(); return x * zf;
}
int ksm(ll x,int y,const int p)
{
if(!y)return 1;int ew=1;
while(y>1)
{
if(y&1)y--,ew=x*ew%p;
x=x*x%p;y>>=1;
}
return x*ew%p;
}
int pri[4]={2,3,7,61};
bool check(const int x)//用2,3,7,61检验x
{
if(x==1)return 0;
if(x==2||x==3||x==5||x==7||x==11||x==61)return 1;
if(x%2==0||x%3==0||x%5==0||x%7==0||x%11==0)return 0;//小剪枝
for(rt i=0;i<=3;i++)
{
int d=x-1;
while(!(d&1))d>>=1;//找到最大的奇数因数
int s=ksm(pri[i],d,x);
while(s!=1&&s!=x-1&&d!=x-1)d<<=1,s=(ll)s*s%x;
//只要不是1或x-1或者整段序列已经遍历完毕就一直网上跳
if(s!=x-1&&!(d&1))return 0;//如果不存在x-1且序列最后一个不是1则为合数
}
return 1;
}
int main()
{
int n=read(),m=read();
while(m--)puts(check(read())?"Yes":"No");
return 0;
}
$ Pollared-Rho$算法
Pollared-Rho是用于质因数分解的随机算法
模板: 求$ \phi(n)$， $ n<=10^{18}$ here
朴素算法：在$ [2..\sqrt{n}]$范围内枚举质数判断是否为n的因数，然后直接计算
如果$ n$是$ 10^9$规模的两个质数相乘，复杂度显然过大
随机算法:
在$ [2..n]$范围内每次随机一个数，判断是否是x的因数
这种算法的复杂度均摊是$ O(n)$规模的，甚至劣于暴力
引出$ Pollard-Rho$算法
主要思想：对于一个大整数$ n$，我们取任意一个数$ x$使得$ x$是$ n$的因数的几率很小，但如果取两个数作差使得差为n的因数的几率就提高了，如果判断的是$ gcd(|x1-x2|,n)>1$那么概率会更高，这就是$ Pollard-Rho$的主要思想.
如何取随机数？
生成一个伪随机数列，这里定义伪随机函数为对于一个值$ x$和种子$ seed$会返回固定的新值
一般令当前$ seed$的$ x$的下一项为$ Nex(x)=x*x+seed (mod \ n)$
每次选取数列中相邻两个作差判断
可能会形成循环？
对！要的就是循环！
复杂度证明：
根据生日悖论，有在整数域$ [1..n]$范围内随机约$ \sqrt{n}$个数，就有约二分之一的概率存在两个相同的数.我们伪随机k个数，标号$ a1..ak$，令$ n$中最小的质因子为$ n1$,则有$n1<= \sqrt{n}$
将$ a$数组对$ n1$取模后的数组称为$ b$数组
则当$ b$数组出现循环的时候(约在$ n^{\frac{1}{4}}$处)$ a$数组有极大概率不进入循环
因为$ a$数组出现循环的期望位置约在$ n^{\frac{1}{2}}$处
此时若a尚未进入循环，则必存在某$ i ≠j$满足$ b[i]==b[j]$且$ a[i]!=a[j]$, 此时必有$ n1|abs(a[i]-a[j])$
因而在约进行$ n^{\frac{1}{4}}$次运算后跑出一个因数
复杂度：$ O(n^{\frac{1}{4}}log_2n)$ //不能忽略$ gcd$的复杂度
如果$ a$数组和$ b$数组同时进入循环？
暴力一点，更换一个$ seed$再次重复上述内容
如何判环？
可以采用$ floyd$判环
即对于每个当前随机数$ x，a$每次走一步，$ b$每次走两步，每次计算$ gcd(|a-b|,n)$
即每次$ a=Nex(a),b=Nex(Nex(b))$
而不是像之前一样$ b$在$ a$前一步然后每次$ a，b$都走一步
如果$ a,b$在同一环上显然存在某一时刻$ b$追上$ a$，即$ b$的值等于$ a$的值
这时候就直接跳出换一个新的$ seed$即可
期望换$ seed$的次数很小，因而复杂度可以满意
得到一个因数之后$ Miller-Rabin$判断是否是质数，如果不是递归继续，把得到的质数用$ map$存下来计算即可
代码：
#include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #include<map> #define rt register int #define l putchar('\n') #define ll long long #define r read() using namespace std; ll read() { ll x = 0; int zf = 1; char ch; while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar(); if (ch == '-') zf = -1, ch = getchar(); while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar(); return x * zf; } ll lowpow(ll x,ll y,ll p) { //计算x*y %p 防止爆long long; ll ew=0; while(y>1) { if(y&1)y--,ew=(ew+x)%p; else x=(x<<1)%p,y>>=1; } return (x+ew)%p; } ll ksm(ll x,ll y,const ll p)//快速幂 { if(!y)return 1;ll ew=1; while(y>1) { if(y&1)y--,ew=lowpow(x,ew,p); x=lowpow(x,x,p);y>>=1; } return lowpow(x,ew,p); } int pri[4]={2,3,7,61}; bool check(const ll x)//用2,3,7,61检验素数 { if(x==2||x==3||x==5||x==7||x==11||x==61)return 1; if(x%2==0||x%3==0||x%5==0||x%7==0||x%11==0)return 0; if(x==1)return 0; for(rt i=0;i<=3;i++) { ll d=x-1; while(!(d&1))d>>=1; ll s=ksm(pri[i],d,x); while(s!=1&&s!=x-1&&d!=x-1)d<<=1,s=(ll)lowpow(s,s,x); if(s!=x-1&&!(d&1))return 0; } return 1; } ll seed;//随机种子 ll Nex(ll x,ll p) { return (lowpow(x,x,p)+seed)%p;//生成伪随机序列中x的下一项 } ll gcd(ll x,ll y) { return (!y)?x:gcd(y,x%y); } map<ll,int>s;ll ans=1; void Pollared_Rho(ll n) { if(check(n))//如果n已经是素数 { s[n]++; if(s[n]==1)ans=ans*(n-1);else ans=ans*n;//计算φ return; } while(1)//如果a,b在同一环中就不断重复 { seed=rand()%(n-1)+1;//新建随机种子 ll a=rand()%(n-1)+1,b=Nex(a,n);//随机a,b初值 while(a!=b) { ll num=gcd(n,abs(a-b)); if(num>1) { Pollared_Rho(num); Pollared_Rho(n/num); return; } a=Nex(a,n);b=Nex(Nex(b,n),n);//b每次比a多走一步 } } } int main() { ll u=r; if(u==1) { cout<<1; return 0; } Pollared_Rho(u); cout<<ans; return 0; }
