素数判定&质因数分解(数论)(Miller Rabin)(Pollard Rho)

太玄学了！

我真的被概率的魅力折服了。此前我认为1便是1，0.9999999999…便是0.9999999999…。
但实际上它们有着千丝万缕的关系。
试想，如果一件事发生的概率是0.9999999999，你觉得它是不是必然事件？
数学逻辑严谨的人说：“emmmm，这得回到‘必然事件’的定义上去。所以嘛，其实还是差了一丢丢的。”
然而，抛开束缚人的数学定义，从直观上讲，我们完全可以相信，这件事全然是一件必然事件。
那人还说：“我信仰！我就是要在十个铭文位中选一个祸源！”
那好吧，或许我说服不了你，但你且看这个算法，我相信你会改变你的想法。

（介于作者水平，我就只讲做法。至于原理，你感兴趣的话自己下去证吧。不会把不会吧，不会真有人感兴趣吧[奸笑] ）

今天有个好朋友将从头到尾一直陪伴着我们，它是rand()。

素数判定:Miller Rabin算法

初阶:费马小定理

玄学的素数判定法：Miller_Rabin算法
先搬出费马小定理：
$$x^{n-1}\equiv 1(\mathrm{mod}n)$$
这个式子成立的前提是$\gcd (x,n)=1$。那么我们反过来用，如果我们让x先n-1次幂，如果对n取模结果不为1，那么$\gcd (x,n)$便是$n$的一个因子了(注意, 并不是说$x\mid n$,而是$gcd(x,n)\mid n$,这就是高级之处所在) 。而这个x当然不会让它从$2$到$\sqrt{n}$这样for一遍，这样效率是没有变化的。出其不意地，我们会让x=rand()%(n-1)+1，即在$[1,n-1]$从随机选一个数。

非常神奇，这样选出来的数，按以往的方法，去尝试$n\%x==0?$命中的概率十分低，然而去尝试$x^{n-1}\%n!=1?$，概率却能大大UP！你不服可以自己去试一下，据说至少是$\frac{3}{4}$的概率能成功，而实测更高，我还没找到需要随机两次的。。。即使真是保底$\frac{3}{4}$的概率成功，我们做它10次，$1-(\frac{1}{4}{)}^{10}=0.99999904632568359375$，成功率上升到了这个数，你觉得稳了吗？

据此，注意细节，写出代码：

bool miller_rabbin(__int64 n)
{
    __int64 i,a;
    for(i=0;i<50;i++)//自己调节次数
    {
        a=rand()%(n-2)+2;
        //printf("test%d:%I64d\n",i,a);
        if(mod_exp(a,n-1,n)!=1)
            return false;
    }
    return true;
}

进阶:二次勘测定理

这还只是一个低配版的素数判断。我们再摆出一条式子，二次探测定理：
$$x^2\equiv 1(\mathrm{mod}p)\Rightarrow x1=1,x2=n-1$$
用人话说，如果p是素数，如果x的平方的满足左边的式子，那么x一定是1或n-1（在对p取模下）。我的理解是，如果p为素数，数A是平方得到的数，且A满足$A^2\%n=1$，那么平方前的数必为1或n-1。

这是另一个判断方法。

然后我们考虑把上面两种方法撮合在一起。
随机一个数x，我们想做费马小定理检验，则需要$x^{p-1}$，对p-1进行变化，看成$p-1=2^k\ast u$, 整体可以看成$x^{p-1}=x^{2^k\cdot u}=(x^u{)}^{2\cdot \mathrm{2...}\cdot 2}$.

我们对$x=x^u$进行二次探测定理检验，不符合定理前提或者检验通过，就把$x=x^2$，反复进行. 直到x变到了$x^{p-1}$.

这个时候就是费马小定理的检验时刻, 直接判断x是否为1. 不是的话, 说明不是素数, 这其实就回到了我们上面给出的第一种方案.

第二种方案, 是把第一种方案的快速幂过程拆解了一下, 在中途插入二次探测定理的检验.

bool miller_rabin(ll n)
{
	if(n==2) return true;
	if(n<2||n%2==0) return false;
	ll u=n-1,pre,x;
	int k=0;
	while(!(u&1)) u>>=1,k++;
	for(int T=1;T<=10;T++)
	{
		x=rand()%(n-1)+1;
		x=qkpow(x,u,n);
		for(int i=1;i<=k;i++)
		{
			pre=x;
			x=qkmul(x,x,n);
			if(x==1 && pre!=1 && pre!=n-1) return false;
		}
		if(x!=1) return false;
	}
	return true;
}

质因数分解:Pollard Rho算法

这个其实没怎么搞懂…
拿题目poj1811来讲,要求一个数的最小质因数. 思路如图, 先找到一个因数d, 那么n/d是它的另一个因数, 分别去求d和n/d的因数,知道为不可再分的质数.

如何找因数是一个玄学问题. 由于我也不能理解, 便只是将步骤记录如下:

1. 随机两个数x, c. 为的是一个玄学的函数$f(x)=x^2+c$
2. 定义y,i,k. (i=1) i,k都是用来计数的,每进行$2^k$次运算就要重新将y赋为当前的x.所以一开头就要了(y=x)
3. 开始循环, x=f(x).
4. 求$\gcd (|x-y|,n)$, 我们的主角并非x,也并非y, 是x与y的产物$|x-y|$
5. gcd()!=1? 如果gcd的结果不是1！！！那我们就成功了, 返回gcd的结果即可.
6. 否则判断x==y? 是否进入了$\rho$的头部, 进入了, 接下来的工作没有意义
7. i++,开始下一个循环 同时注意k是否等于i
   

好了,希望大家能明白…
代码如下(提醒用C++交,G++会RE, poj1811)

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
typedef long long ll;
const ll INF=0x3f3f3f3f;
using namespace std;

ll qkmul(ll a,ll b,ll mod)
{
	ll re=0;
	while(b)
	{
		if(b&1) re=(re+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return re;
}
ll qkpow(ll a,ll b,ll mod)
{
	ll re=1;//debug re=0;
	while(b)
	{
		if(b&1) re=qkmul(re,a,mod);
		a=qkmul(a,a,mod);
		b>>=1;
	}
	return re;
}
ll Gcd(ll x,ll y)
{
	if(x==0) return y;
	if(x<0) return Gcd(-x,y);
	return y==0?x:Gcd(y,x%y);
}

bool miller_rabin(ll n)
{
	if(n==2) return true;
	if(n<2||n%2==0) return false;
	ll u=n-1,pre,x;
	int k=0;
	while(!(u&1)) u>>=1,k++;
	for(int T=1;T<=10;T++)
	{
		x=rand()%(n-1)+1;
		x=qkpow(x,u,n);
		for(int i=1;i<=k;i++)
		{
			pre=x;
			x=qkmul(x,x,n);
			if(x==1 && pre!=1 && pre!=n-1) return false;
		}
		if(x!=1) return false;
	}
	return true;
}

ll pollard_rho(ll n,ll c)
{
	int i=1,k=2;
	ll x=rand()%n,y=x;
	while(1)
	{
		x=(qkmul(x,x,n)+c)%n;
		ll gcd=Gcd(y-x,n);
		if(gcd>1&&gcd<n) return gcd;
		if(x==y) return n;
		if(i==k) y=x,k<<=1;
		i++;
	}
}

ll minn;
void Find(ll n)
{
	if(miller_rabin(n))
	{
		minn=min(minn,n);
		return ;
	}
	ll p=n;
	while(p==n)
		p=pollard_rho(p,rand()%(n-1)+1);
	Find(p);
	Find(n/p);
}

int main()
{
	srand(time(0));
	int BRI;scanf("%d",&BRI);
	while(BRI--)
	{
		ll n;scanf("%lld",&n);
		if(miller_rabin(n))
		{
			puts("Prime");continue;
		}
		minn=INF;
		Find(n);
		printf("%lld\n",minn);
	}
	return 0;
}