pollard_rho算法

pollard_rho 算法

这是一个因数分解的算法，时间复杂度大概为O(n1/4)O(n^{1/4})O(n1/4)。
算法流程是这样的：
先对n进行素数判断，这一步采用miller_rabbin算法来做。如果是素数，则无需进行因数分解了。
否则，选取一个随机数x，使得x<n, 按如下方法产生一个序列：
{a1=x%nai=(ai−12+C)%n\begin{cases} a_1=x \%n \\ a_i=(a_{i-1}^2+C)\%n\end{cases}{a1=x%nai=(ai−12+C)%n
其中C是一个常数。
因为有模n的存在，所以，序列最后一定会陷入循环中。这个序列有点像希腊字母ρ\rhoρ的形状，这是本算法为什么会有“rho”的原因。
进入循环以后，我们就不用考虑后续的项了。
取序列的相邻两项，设为ai,ai−1a_i,a_{i-1}ai,ai−1,
则有ai−ai−1=(ai−12+C)%n−(ai−22+C)%n=(ai−1−ai−2)∗(ai−1+ai−2)%na_i-a_{i-1}=(a_{i-1}^2+C)\%n-(a_{i-2}^2+C)\%n=(a_{i-1}-a_{i-2})*(a_{i-1}+a_{i-2})\%nai−ai−1=(ai−12+C)%n−(ai−22+C)%n=(ai−1−ai−2)∗(ai−1+ai−2)%n
上面的式子告诉我们，相邻两数的差会包含前面所有相邻两数差的乘积。这些乘积中可能会包含n的因子，模n后n的因子仍然存在。
在序列中取ai,aja_i,a_jai,aj,假设i>ji>ji>j,
则有：(ai−aj)%n=(ai−12+C−aj−12+C)%n=(ai−1−aj−1)∗(ai−1+aj−1)%n(a_i-a_j)\%n=(a_{i-1}^2+C-a_{j-1}^2+C)\%n=(a_{i-1}-a_{j-1})*(a_{i-1}+a_{j-1})\%n(ai−aj)%n=(ai−12+C−aj−12+C)%n=(ai−1−aj−1)∗(ai−1+aj−1)%n
这就是说，序列中距离为k的两项的差，一定为前面任意距离为k的两项的差的倍数，同样的，模n后也不会导致n的因子丢失。

在序列中取ai,aja_i,a_jai,aj,假设i>ji>ji>j,
则有：(ai−aj)%n=(ai−ai−1+ai−1−⋯−aj)%n=b∗(aj−aj−1)%n(a_i-a_j)\%n=(a_i-a_{i-1}+a_{i-1}-\dots-a_j)\%n=b*(a_j-a_{j-1})\%n(ai−aj)%n=(ai−ai−1+ai−1−⋯−aj)%n=b∗(aj−aj−1)%n
可知：序列中任意两数的差，也一定可以转换为相邻两个数的差的倍数。

接下来我们了解一下“生日悖论”。
这个原理来源于“生日悖论”：23个人的班级，有50%以上的概率，该班存在两个人的生日相同。
设PnP_nPn表示任意选取nnn个人，不存在任意两人的生日相同的概率。
Pn=∏1≤j<n(1−j/365)P_n=\prod_{1\leq j<n}(1-j/365)Pn=1≤j<n∏(1−j/365)

由均值不等式可得：
(∏1≤j<i(1−j/365)n−1<1n−1∑j=1n−1(1−j/365)\sqrt[n-1]{(\prod_{1\leq j<i} (1-j/365)} <\frac{1}{n-1}\sum_{j=1}^{n-1}(1-j/365)n−1(1≤j<i∏(1−j/365)<n−11j=1∑n−1(1−j/365)
所以：Pi<(1n−1∗(n−1)(2−n/365)/2)n−1=(1−n/730)n−1Pi<(\frac{1}{n-1}*(n-1)(2-n/365)/2)^{n-1}=(1-n/730)^{n-1}Pi<(n−11∗(n−1)(2−n/365)/2)n−1=(1−n/730)n−1
根据泰勒定理的无穷小推论：e−x>(1−x)e^{-x}>(1-x)e−x>(1−x),可得：
Pi<e−n/730∗(n−1)P_i<e^{-n/730*(n-1)}Pi<e−n/730∗(n−1)
当n=23时，则Pi<0.5P_i<0.5Pi<0.5.
所以，存在两人生日相同的概率超过了50%。
即在n个数中选择k个数，存在两个数相等的概率，当k接近n∗log⁡e2\sqrt{n*\log_{e}2}n∗loge2,已经接近50%了。
两个数相等，可以看作两个数的差为0。
所以，稍加扩展，上述结论即为：
在n个数中选择k个数，存在两个数的差值为n的因数的概率，当k接近n∗log⁡e2\sqrt{n*\log_{e}2}n∗loge2,已经接近50%了。
然而，n可能很大，这样k=nk=\sqrt{n}k=n也很大，这k个数不好保存。

关键的来了。回忆我们开始产生的那个序列，假设序列的在未出现循环前的长度为lenlenlen,我们需要检测序列中的数的两两之差，其中大概率包含n的因子,但那是平方级别的时间复杂度。因为任意两个数都可以做差。我们不必那么做。
根据前面的分析，序列有一些特殊的性质，即下标距离为k的两个数之差，是前面的下标距离为k的某两个数之差的倍数。所以，对于距离为k的两个数的差，我们只需要检测最后那一对即可。这样，每一个距离都只需要检测一次。
考虑到序列是动态产生的，如何检测它陷入了循环呢。可以使用floyd判环。想象有两个人，两个人从起点出发，一个人每次走一步，另一个人每次走两步。当两个人又碰面了，可以肯定一定出现了环。
下面程序是使用floyd判环来做的。
但是，对于一个环，因为两人碰面了就结束当前序列的检测，实际上并没有检测完序列中所有可能的差值。所以需要不断的生成新的序列，而且必须让生成的序列的式子为：a^2+c,其中a和c都是随机数。
如果c=1，则可能检测不出n的因子，这样即有可能陷入死循环。
比如n=4，c=1出不了结果。
当a=0时，序列a为0,1,2,1,2,……
当a=1时，序列a为1,2,1,2，……
当a=2时，序列a为2,1,2,1,……
当a=3时，序列a为3,2,1,2,1,……
因为只有当a=3，差值在a1,a3a_1,a_3a1,a3之间产生，但是从第二个数开始，就是循环了，循环周期而2。所以，当我们判断出循环后，也没有计算a1a_1a1和a3a_3a3。只有寄希望于c为其他值，前两个数的差值就为2.

LL pollard_rho(LL x)
{while(1){LL a = 1ll*rand()%x;LL c=rand()%x;LL b = (ksc(a,a,x)+c)%x, d; while (1){//cerr<<"a,b"<<a<<','<<b<<endl; if(a==b)break;d = gcd(x,abs(b-a));if(d != 1 && d != x)return d;else if( d==x )break;b = (ksc(b,b,x)+c)%x;b = (ksc(b,b,x)+c)%x;a= (ksc(a,a,x)+c)%x;}}
}

也可以倍增检测区间长度，对该区间，左端点不变，右端点依次递增，检测左右端点值之差是否为n的因子。第一次检测区间[1,2],差为s2-s1;第二次检测区间[2,4],差有：s3-s2,s4-s2;第三次检测区间[4,8],差有s5-s4,s6-s4,s7-s4,s8-s4;……，直到当前检测的区间出现了差值为0的情况，即可以肯定该区间包含了序列的循环周期，停止检测该序列；并生成新的序列进行下一次处理。如果在过程中，发现了n的因子a，则可以递归处理a和(n/a)。当然，在之前都要进行素性测试。

LL pollard_rho(LL x)
{srand(time(0));while (1){LL a = rand() % (x - 1) ;LL c = rand() % (x-1);LL b = a, d;LL cnt = 0,k = 2;while (1){b = (ksc(b,b,x) + c) % x;d = gcd(x, abs(a - b));if (d != 1 && d != x)return d;else if (d == x)break;cnt++;if (k == cnt){a = b;k <<= 1;}}}
}

感觉这个算法对序列中的差值检测还是有些不严谨，最终是通过产生其他随机序列来弥补了漏洞。

例题1：
给出若干组整数，对每组整数，判断它是否刚好有四个约数，如果是，则由小到大输出其大于1的三个约数。
分析：如果一个数包含可以拆成两个不同质数的乘积，则它一定仅包含4个约数。但还有一种情况，如果一个数是一个质数的立方，则它也包含4个约数。
所以，我们用pollard_rho进行因数分解即可，同时配合素性测试检测因数是否为素数。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <ctime>
using namespace std;
#define LL long long int
LL ksc(LL a, LL b, LL p)
{LL res = 0;while (b){if (b & 1){res = res + a;if (res >= p)res -= p;}a = a + a;if (a >= p)a -= p;b >>= 1;}return res;
}
LL ksm(LL a, LL b, LL p)
{LL res = 1;while (b){if (b & 1)res = ksc(res, a, p);a = ksc(a, a, p);b >>= 1;}return res;
}
bool miller_rabbin(LL x)
{if (x == 2 || x == 3)return 1;if (x == 1)return 0;LL n = x - 1, res, lastres;int r = 0, r1, t = 10;while (!(n & 1)){n >>= 1;r++;}srand(time(0));while (t--){LL a = rand() % (x - 1) + 1;res = ksm(a, n, x);if (res == 1)continue;r1 = r;while (res != 1 && r1 > 0){lastres = res;res = ksc(res, res, x);r1--;if (res == 1 && lastres != x - 1)return 0;}if (res != 1)return 0;}return 1;
}
LL gcd(LL a, LL b)
{if (!b)return a;return gcd(b, a % b);
}
LL pollard_rho(LL x)
{srand(time(0));while (1){LL a = rand() % (x - 1) ;LL c = rand() % (x-1);LL b = a, d;LL cnt = 0,k = 2;while (1){b = (ksc(b,b,x) + c) % x;d = gcd(x, abs(a - b));if (d != 1 && d != x)return d;else if (d == x)break;cnt++;if (k == cnt){a = b;k <<= 1;}}}
}
int main()
{//freopen("1.in","w",stdout);LL num, fac1, fac2;bool flg;while (scanf("%I64d", &num) != -1){flg = 1;if (miller_rabbin(num) == 1||num == 1)flg = 0;else{fac1 = pollard_rho(num);fac2 = num / fac1;if (fac1 != fac2 && miller_rabbin(fac1) == 1 && miller_rabbin(fac2) == 1){flg = 1;if (fac1 > fac2)swap(fac1, fac2);}elseflg = 0;}if (flg == 1)printf("%I64d %I64d %I64d\n", fac1, fac2, num);else{printf("is not a D_num\n");}}return 0;
}