C++Pollard_rho分解质因数及其例题————

前言：

在观看此博客之前请学习miller_rabin。

我们在分解质因子时也许只会用试根法（也就是暴力）。

而在此我们将学习一个玄学的算法——Pollard_rho。

概念：

Pollard_rho是一种基于随机的算法，它的思路是先用miller_rabin来判断当前数是否已经是素数了，如果是的话记录并返回。如果不是，我们设要分解的数为n,那么我们考虑去找一个当前数的因数p，找到之后再分别对p和n/p分解质因数，有一点类似分治但是不是分治。

算法如图：

那么，我们要如何快速地找到当前数的因数呢？

这里，我们采用以下方法：先随机生成两个1~n的数，利用这两个数差的绝对值和n比较，来判断这两个数差的绝对值是不是n的一个因数。

可这又是为什么？

那么我们就要讲到一个叫生日悖论的东西了。

这里又出现了一个问题：我们为什么要使用两个随机数差的绝对值来判断是否为n的因数之一，而不是直接生成一个随机数来进行判断呢？

这里我们来讲一下生日悖论：

生日悖论，指如果一个房间里有23个或23个以上的人，那么至少有两个人的生日相同的概率要大于50%。这就意味着在一个典型的标准小学班级(30人)中，存在两人生日相同的可能性更高。对于60或者更多的人，这种概率要大于99%。

如何计算：

第一个人的生日是 365选365

第二个人的生日是 365选364

第三个人的生日是 365选363

第n个人的生日是 365选365-(n-1)

那么，所有人的生日都不同的概率为:

所有人中，至少有两个人生日相同的概率为：

当n=23时，这个概率为0.507，n=100时，概率为0.999999692751072。

当然，选取两个数也是有道理的。

一开始，我们先随机生成两个数x和c，每次都使x=x*x+c(mod n) ，另外再设置一个数y，y的初始值为x，每进行2^k次操作，就让y记为现在的x。在计算的过程中，我们可以发现x=x*x+c(mod n)这个过程是会出现循环的，具体图像会像希腊字母ρ（rho）一旦出现循环我们还没有找到可以分解的因数，就意味着这次随机失败了，需要我们再次进行随机。我们为了记录什么时候出现了循环，我们每个2^k次计算就把y记录成现在的x，如果在接下来的2^k次操作中，出现了y=x，就意味着出现了循环。所以，y的作用一是为了当随机的第二个数，二是为了判断循环。

代码如图：

ll Pollard_rho(ll n1,ll c)//ll 为long long
{ll i=1,k=2;ll x=rand()%n1;ll y=x;while(1){i++;x=(qkch(x,x,n1)+c)%n1;ll d=gcd(y-x,n1);if(d>1&&d<n)return d;if(x==y)return n1;if(i==k){y=x;k=k*2;}}
}

题目描述：

Given a big integer number, you are required to find out whether it's a prime number.

输入：

The first line contains the number of test cases T (1 <= T <= 20 ), then the following T lines each contains an integer number N (2 <= N < 2 54).

输出：

For each test case, if N is a prime number, output a line containing the word "Prime", otherwise, output a line containing the smallest prime factor of N.

样例输入：

2
5
10

样例输出：

Prime
2

此题为模板题就不需多讲。

代码实现：

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<ctime>
#define INF 1ll << 60
using namespace std;
#define ll long long
ll n,ans,t;
ll qkch(ll x,ll y,ll mod)
{ll res=0;while(y){if(y&1)res=(res+x)%mod;x=(x*2)%mod;y/=2;}return res;
}
ll gcd(ll a,ll b)
{if(a==0)return 1;if(a<0)return gcd(-a,b);while(b>0){ll t=a%b;a=b;b=t;}return a;
}
ll qkpow(ll x,ll y,ll mod)
{ll ans=1;while(y){if(y&1)ans=qkch(ans,x,mod);x=qkch(x,x,mod);y/=2;}return ans;
}
bool miller_robin(ll n1)
{if(n1==2)return 1;if(n1<2||!(n1&1))return 0;ll x,u,pre;ll k=0;u=n1-1;while(!(u&1)){k++;u>>=1;}for(int i=1;i<=3;i++){x=rand()%(n1-1)+1;x=qkpow(x,u,n1);pre=x;for(int j=1;j<=k;j++){x=qkch(x,x,n1);if(x==1&&pre!=n1-1&&pre!=1)return 0;pre=x;}if(x!=1)return 0;}return 1;
}
ll Pollard_rho(ll n1,ll c)
{ll i=1,k=2;ll x=rand()%n1;ll y=x;while(1){i++;x=(qkch(x,x,n1)+c)%n1;ll d=gcd(y-x,n1);if(d>1&&d<n)return d;if(x==y)return n1;if(i==k){y=x;k=k*2;}}
}
void findit(ll p)
{if(miller_robin(p)==1){if (ans>p)ans=p;return;}ll p1=p;while(p1>=p){p1=Pollard_rho(p1,rand()%(p-1)+1);//printf("%lld %lld\n",p,p1);}findit(p1);findit(p/p1);
}
int main()
{srand(time(0));scanf("%lld",&t);while(t--){ans=INF;scanf("%lld",&n);if(miller_robin(n))printf("Prime\n");else{findit(n);printf("%lld\n",ans);}}
}