P4718 【模板】Pollard-Rho算法

发表于 | 更新于 | 评论数： | 阅读次数：
本文字数： 7.1k | 阅读时长 ≈ 6 分钟

P4718 【模板】Pollard-Rho算法

传送门

洛谷

题目描述

Miller Rabin算法是一种高效的质数判断方法。虽然是一种不确定的质数判断法，但是在选择多种底数的情况下，正确率是可以接受的。

Pollard Rho是一个非常玄学的方式，用于在$O(n^{1/4})$的期望时间复杂度内计算合数n的某个非平凡因子。事实上算法导论给出的是$O(\sqrt p)$，$p$是$n$的某个最小因子，满足$p$与$n/p$互质。但是这些都是期望，未必符合实际。但事实上Pollard Rho算法在实际环境中运行的相当不错。

这里我们要写一个程序，对于每个数字检验是否是质数，是质数就输出Prime；如果不是质数，输出它最大的质因子是哪个。

输入输出格式

输入格式：

第一行，$T$代表数据组数（不大于$350$）

以下$T$行，每行一个整数$n$，保证$1\le n\le 10^{18}$。

输出格式：

输出$T$行。

对于每组测试数据输出结果。

输入输出样例

输入样例#1：

6
2
13
134
8897
1234567654321
1000000000000

输出样例#1：

Prime
Prime
67
41
4649
5

分析

先讲讲算法。

Miller-Rabin算法

简介

Miller-Rabin算法是一种可以快速判定大质数的随机算法，有一定概率阵亡。
嗯…好像讲完了。

前置知识

首先，我们需要知道一个定理——费马小定理。

其次，我们还需要知道一个定理(不知道算不算定理，以下姑且将其称为定理)
对于一个奇质数$p$且$e \ge 1$，满足$a^2 \equiv 1 \pmod {p^e}$的根只有平凡平方根。
严格的证明不会，口胡一下还是会的。
首先我们化式子。

即，$p \mid (a+1)$或者$p \mid (a-1)$。
然后分情况讨论再合并，就是上面这个东西了。
于是Miller和Rabin就捣鼓出了传说中的Miller-Rabin素性测试算法。

实现

Miller-Rabin算法的主体思路很简单。
首先，随机一个数$a$。
想想费马小定理，我们想知道一个数$n$是不是素数，就验证一下$a^{n-1} \equiv 1 \pmod p$这个式子是否成立。
但这样是有几率误判的（因为费马小定理的逆定理本来就是错的）。
如取$a=2$，在前10亿个自然数中共有$50847534$个素数，而满足$2^{n-1} \equiv 1 \pmod p$的合数有$5597$个。这样算下来，算法出错的可能性约为$0.00011$。我们无法接受。
无法接受归无法接受，我们还是得用它。
所以想办法提升它的正确性。
还记得第二条定理吗？
我们将$a^{p-1}$拆开，拆成$a^{d \times 2^s}$，那么就像相当于将$a^d$一直平方，就可以套第二个式子了。
当然能通过上面两个定理的检验的假素数还是有的（藏在数群中带着素数的面具勾搭素MM）。
所以，我们还需要多次取数检验提高正确率。
通常认为，Miller-Rabin算法中随机选取$ k$个底数进行测试的失误率大概为$4^{-k}$。

Miller-Rabin算法中，一般情况下底数随机选取，但当待测数不太大时，选择测试底数就有一些技巧了。比如，如果被测数小于$4759123141$，那么只需要测试三个底数$2, 7$和$61$就足够了。如果你每次都用前$7$个素数($2, 3, 5, 7, 11, 13$和$17$)进行测试，所有不超过$341550071728320$的数都是正确的。如果选用$2, 3, 7, 61$和$24251$作为底数，那么$10^{16}$内唯一的强伪素数为$46856248255981$。这样的一些结论使得Miller-Rabin算法在OI中非常实用。
$\uparrow$ 这段底数选取的数据来自Norlan的素数与素性测试（Miller-Rabin测试）。

Pollard-Rho算法

Pollard-Rho算法的主要思想是找因子分解，这个因子不一定是质因子。
那么问题就在于如何找到它的因子。
首先我们知道Pollard-Rho算法是一个随机算法，那么我们找因子自然是随机。
但每次生成一个数显然撞上的可能性太小，那么根据生日悖论，我们每次多生成一些的话，且每次我们尝试的数为$x-y$（$x,y$为我们生成的随机数，下文一样），撞上的可能性就会增加。
但这样我们又面临这一个问题，我们的程序可以调用的内存是有限的，如果我们每次生成很多数，空间就爆炸了。
真是一个令人尴尬的问题。
好在已经被解决了，要不然要这算法有何用。
因为我们每次只调用两个随机数，提前把全部随机数生成好毫无意义。
那么我们就按照一定顺序生成数。
于是就有了一个伪随机数生成函数——

其中$a$是我们随机的一个数（个人觉得很像随机种子）。
那么我们就用了一系列随机数，还不用将其全部存放与内存中。
但Pollard-Rho算法的精髓还不止这一点。
我们发现尽管生成了一串随机数，但要求能整除$n$，这样的数被撞上的概率还是不足以让人满意。
于是Pollard又想出了一个办法。
概率不然人满意，是不是我们对生成的数要求太高了。
那换个要求，我们要求$gcd(x-y,n)>1$。
这样的数又$x,2x,…,(y-1)x,y,2y,…,(x-1)y$，共$x+y-2$个数。
于是成功率可以让人接受了。
但算法还是不怎么完美，因为我们的伪随机数生成函数是长这样的。

$\uparrow$ 图来自 努力努力再努力x 的 大数因式分解 Pollard_rho 算法详解

很像希腊字母$\rho$，这也是算法名字中Rho的来历。
于是新的问题再次被提了出来。
怎么判断环？？？
如果把生成的数都存下来判环，我们使用这个函数的意义就没有了。
好在前人已经帮我们解决了这个问题。
我们可以使用Floyd判环法。
想象两个人A,B在一个环形跑道上奔跑，我们让A的速度是B的两倍，当A第一次赶上B时，我们就知道B一定比A多跑了一圈。
但实际应用中我们采用的是brent判环法（因为更高效），具体的，就是倍增，让A每次比B多跑$2^i$步。
吐槽：我说我怎么看不懂他们在写什么，原来是因为我不知道我写的是brent判环法，还以为我们都是用的Floyd判环法。 连判环法都认不出来的蒟蒻一只，还有，觉得Floyd判环法好丑好麻烦。
来来来，我们分析一波复杂度。
根据生日悖论，每 $k = \sqrt n$个人中有 $50\%$的可能性产生生日冲突的话。
那么这个算法的期望复杂度应该是带一个根号的，粗略的说，我们试一个数，原来试到的概率为$\frac{1}{n}$，加上生日悖论的话，变成$\sqrt \frac {1}{n}$，现在有$x+y-2$个数可以成功，再乘一下，概率还是比较大。所以期望复杂度感觉上虽然带个根号，当应该不会很高。
事实上，算法导论上给出的期望复杂度为$O( \sqrt p )$，$p$为$n$的某个最小因子，满足$p$与$n/p$互质。证明见算法导论。
好像还可以当作$O( n ^ \frac{1}{4} )$来算。

参考资料

Miller-Rabin算法

努力努力再努力x的 大素数测试的Miller-Rabin算法
Norlan的 素数与素性测试（Miller-Rabin测试） 讲的非常好，语言也很生动形象，很多本文由于作者懒没有提到的的定义，证明，定理以及有趣的历史都可以在上面找到。

Pollard-Rho算法

book丶book丶的 Pollard_rho 因数分解 主要看判环法的代码。
Doggu的 大数质因解：浅谈Miller-Rabin和Pollard-Rho算法 以及由Doggu翻译的一篇文章 感觉两篇文章都非常好，很多本人没法证明的或者是解释不清楚的都可以在这两篇文章里找到。
努力努力再努力x的 大数因式分解 Pollard_rho 算法详解
cz_xuyixuan的 【学习笔记】Pollard’s rho算法 写的比较详细，不过本人理解起来较为困难，提供神奇的优化。
Great_Influence的洛谷题解 Pollard-Rho 变成了 Pollard−Pho ：) $46856248255981$好像可以被其他数字判掉？？？学到了。

代码

不要打龟速乘，不要打龟速乘，不要打龟速乘
打$O(1)$快速乘，亲测龟速乘跑$4s$，快速乘跑$0.5s$。
为什么我一个每个数$gcd$一下$Wa$了，优化算法复杂度的时候每$127$个数$gcd$一下，竟然$A$了，优化复杂度怎么优化到正确性上去了？？？
卡常卡死我了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<time.h>
using namespace std;
template <typename T>inline T read()
{
    register T sum=0;
    register char cc=getchar();
    int sym=1;
    while(cc!='-'&&(cc>'9'||cc<'0'))cc=getchar();
    if(cc=='-')sym=-1,cc=getchar();
    sum=sum*10+cc-'0';
    cc=getchar();
    while(cc>='0'&&cc<='9')sum=sum*10+cc-'0',cc=getchar();
    return sym*sum;
}
template <typename T>inline T read(T &a)
{
    a=read<T>();
    return a;
}
template <typename T,typename... Others> inline void read(T& a, Others&... b)
{
    a=read(a);
	read(b...);
}
int T;
long long ans;
const long long MR[]={2,3,7,61,24251};
//long long mul(long long x,long long y,long long p)
//{
//	long long ans=0;
//	while(y)
//	{
//		if(y&1)ans=(ans+x)%p;
//		x=(x+x)%p;
//		y>>=1;
//	}
//	return ans;
//}
inline long long mul(long long x,long long y,long long mod)
{
	long long tmp=(x*y-(long long)((long double)x/mod*y+1.0e-8)*mod);
	return tmp<0 ? tmp+mod : tmp;
}
inline long long add(long long x,long long y,long long p)
{
	long long tmp=x+y;
	if(tmp>p)tmp-=p;
	return tmp;
}
inline long long fpow(long long x,long long y,long long p)
{
	long long ans=1;
	while(y)
	{
		if(y&1)ans=mul(ans,x,p);
		x=mul(x,x,p);
		y>>=1;
	}
	return ans;
}
long long gcd(long long x,long long y)
{
	return y? gcd(y,x%y):x;
}
inline bool Miller_Rabin(long long n)
{
	if(n==2||n==3)return true;
	if(n%2==0||n==1||n==46856248255981)return false;
	long long d=n-1;
	int s=0;
	while(!(d&1))s+=1,d>>=1;
	for(int i=0;i<5;i++)
	{
		if(n==MR[i])return true;
		if(n%MR[i]==0)return false;
		long long x=fpow(MR[i],d,n),y=0;
		for(int j=1;j<=s;j++)
		{
			y=mul(x,x,n);
			if(y==1&&x!=1&&x!=n-1)return false;
			x=y;
		}
		if(y!=1)return false;
	}
	return true;
}
inline long long Irand(long long x)
{
    return 1ll*((rand()<<15^rand())<<30^(rand()<<15^rand()))%x;
}
inline long long Pollard_Rho(long long n,long long c)
{
	long long x=Irand(n-1)+1;
	long long y=x,k=2,val=1;
	for(int i=1;;i++)
	{
		x=add(mul(x,x,n),c,n);
		val=mul(val,add(y-x,n,n),n);
		if(i%127==0)
		{
			long long d=gcd(val,n);
        	if(d!=1&&d!=n)return d;
        	val=1;
		}
//		long long d=gcd(add(y-x,n,n),n);
//        if(d!=1&&d!=n)return d;
		if(x==y)return n;
		if(i==k)
		{
			long long d=gcd(val,n);
        	if(d!=1&&d!=n)return d;
        	val=1;
			y=x;
			k<<=1;
		}
	}
}
void Find(long long n,long long c)
{
	if(n==1||n<=ans)return ;
	if(Miller_Rabin(n))
	{
		ans=max(n,ans);
		return ;
	}
	long long p=n;
	while(p>=n)p=Pollard_Rho(p,c--);
	Find(p,c);
	Find(n/p,c);
}
int main()
{
	srand(19491001);
	T=read<int>();
	for(int tot=1;tot<=T;tot++)
	{
		long long x=read<long long>();
		if(!Miller_Rabin(x))
		{
			ans=0;
			Find(x,100);
			printf("%lld\n",ans);
		}
		else puts("Prime");
	}
    return 0;
}