pollard_rho

复习整理一下PollardRho的相关知识把，刚好密码，acm都要用

Miller Rabin素性检验

Miller Rabin素性检验是一种素数判定的法则，由CMU的教授Miller首次提出，并由希大的Rabin教授作出修改，变成了今天竞赛人广泛使用的一种算法，故称Miller Rabin素性检验。
该算法本质上是一种随机化算法，能在$O(klo{g}^{3}n)$ 的时间复杂度下快速判断出一个数是否是素数，但具有一定的错误概率。不过在一定数据范围内，通过一些技巧可以使其不出错。

Fermat素性检验

由费马小定理我们得知：对任意素数$p$和与其互质的整数$a$ ，$a^{p-1}\equiv 1(modp)$ 。
随机选取一小于$p$的整数$a$，若$a^{p-1}\equiv 1(modp)$ ，则是否可以证明$p$是一素数？
答案是不行，因为存在一些合数满足这个等式，这些合数被称为费马伪素数，例如最小的费马伪素数为341。
通过多选取几个底数，我们很大程度上降低了错误的概率，比如341就成功被筛去了。但仍旧存在极少的一些合数，即便遍历$[2,p-1]$的每一个数字作为底数，也无法筛去。
这样的合数被称为卡迈克尔数，在一亿内有255个，最小的卡迈克尔数为561。
若$n$为卡迈克尔数，则$2^n-1$也是卡迈克尔数，故其个数是无穷的。

Miller Rabin素性检验

单单费马小定理行不通，所以我们考虑引入新的定理来提高检测的正确性。
二次探测定理 ：对于质数$p$，若$x^2\equiv 1(modp)$，则小于$p$的解只有两个，$x_1=1,x_2=p-1$。
由Fermat检测得到$a^{p-1}\equiv 1(modp)$，且 $p-1$为偶数（否则$p$为偶数直接特判筛掉），则 $a^{p-1}$就相当于$x^2$。
也就是说，我们可以将$p-1=u\times 2^t(u为奇数)$ ，对$a^u,a^{u\times 2},a^{u\times 2^2}\dots$ 这一系列数进行检验，他们的解要么全是1，要么出现$p-1$后全是1（之前不能有1），否则就不是素数。当然，要注意$p-1$ 不能出现在最后一个数，否则就连费马小定理都不满足了。
还要注意这过程中不能产生$p$的倍数。
总结一下，该检验的思路是：

1）先特判筛掉3以下的数与偶数。
2）将待检验数$n$，$n-1$化为$u\times 2^t$ 。
3）选取多个底数，分别对$a^u,a^{u\times 2},a^{u\times 2^2}\dots$ 进行检验，判断其解是否全为1，或在非最后一个数的情况下出现$p-1$。
4）如果都满足，我们就认为这个数是素数。

代码解决

那么思路的问题解决了，接下来考虑代码的实现。
首先是底数的选取：Miller Rabin检测依旧有错误的概率，虽然微乎其微，但我们显然不想接受。但通过选取适合的底数，可以避免这一情况的发生。
在int范围内，选取 { 2 , 7 , 61 } \{2,7,61\} {2,7,61} 三个数作为底数可以保证100%正确。
在long long范围内， { 2 , 325 , 9375 , 28178 , 450775 , 9780504 , 1795265022 } \{2,325,9375,28178,450775,9780504,1795265022\} {2,325,9375,28178,450775,9780504,1795265022} 七个数作为底数可保证100%正确。
故正常情况下时间复杂度最多为$O(lo{g}^{3}n)$ ，常数为7。
其次，由于这个过程中可能要计算long long型变量的平方，所以要考虑数据溢出的问题，解决方案是快速乘。
这里提供一个O(1) 的快速乘，原理是用溢出来解决溢出。

inline ll qmul(ll x,ll y,ll mod){
    return (x*y-(ll)((ld)x/mod*y)*mod+mod)%mod;     
}

总参考代码：

inline ll qmul(ll x, ll y, ll mod){
    return (x*y-(ll)((ld)x/mod*y)*mod+mod)%mod;     
}
ll qpow(ll a, ll n, ll mod)/{
    ll res=1;
    while(n){
        if(n&1) res = qmul(res,a,mod);
        a = qmul(a,a,mod);
        n >>= 1;
    }
    return res;
}

bool Miller_Rabin(ll n){
    if(n<3 || n&1==0) return n==2;
    ll u=n-1, t=0;
    while(u&1 == 0) u/=2, t++;
    ll ud[]={2,325,9375,28178,450775,9780504,1795265022};
    for(ll a:ud){
        ll v=qpow(a, u, n);
        if(v==1 || v==n-1 || v==0) continue;
        for(int j=1;j<=t;j++){
            v=qmul(v,v,n);
            if(v==n-1 && j!=t){v=1;break;}
            if(v==1) return 0;
        }
        if(v!=1) return 0;//Fermat检验
    }
    return 1;
}

Floyd判圈算法（龟兔赛跑算法）

算法简述

Floyd判圈算法（Floyd Cycle Detection Algorithm），又称龟兔赛跑算法（Tortoise and Hare Algorithm），是一个可以在有限状态机、迭代函数或者链表上判断是否存在环，以及判断环的起点与长度的算法。

基本思路

在某种关系下，顶点 i 到 k 拓扑有序，顶点 k 到 j 也是相同的顺序，那么 i 和 j 也存在这个顺序。要是某一个顶点出现了自己到自己的环，那么图中就有环，但是这种方法复杂度高一些，没有检测顶点出度或者DFS的方法快，但是非常简单。

（1）判断是否有环
龟兔解法的基本思想可以用我们跑步的例子来解释，如果两个人同时出发，如果赛道有环，那么快的一方总能追上慢的一方。进一步想，追上时快的一方肯定比慢的一方多跑了几圈，即多跑的路的长度是圈的长度的倍数。
基于上面的想法，Floyd用两个指针，一个慢指针（龟）每次前进一步，快指针（兔）指针每次前进两步（两步或多步效果时等价的，只要一个比另一个快就行）。如果两者在链表头以外的某一点相遇，那么说明链表有环，否则，如果（快指针）到达了链表的结尾，那么说明没环。
（2）求环的长度
相遇的时候，一定已经在环上了，然后两个人只要再次在环上接着跑，再次相遇的时候（也就是所谓的套圈），跑的快的那个人就比跑的慢的人整整多跑了一圈，所以环的长度也就出来了。
（3）如何确定环的起点
环的检测用上述原理，接下来我们来看一下如何确定环的起点，这也是Floyd解法的第二部分。方法是将其中一个指针移到链表起点，两者同时移动，每次移动一步，那么两者相遇的地方就是环的起点。

解析：

首先假设第一次相遇的时候慢指针走过的节点个数为i，设链表头部到环的起点的长度为m，环的长度为n，相遇的位置与起点与起点位置距离为k。
于是有：$$i=m+a\ast n+k$$,其中a为慢指针走的圈数。
因为快指针的速度是慢指针的2倍，于是又可以得到另一个式子：$$2\ast i=m+b\ast n+k$$,其中b为快指针走的圈数。
两式相减得：$$i=(b-a)\ast n$$,也就是说i是圈长的整数倍。
这是将其中一个节点放在起点，然后同时向前走m步时，此时从头部走的指针在m位置。而从相遇位置开始走的指针应该在距离起点i+m，i为圈长整数倍，则该指针也应该在距离起点为m的位置，即环的起点。

PollardRho 算法：

PollardRho 是一个很神奇的算法，用于在 $O(n^{1/4})$ 的期望时间复杂度内计算合数n的某个非平凡因子（除了1和它本身以外能整除它的数）。事书上给出的复杂度是 $O(\sqrt{p})$ ， p 是 n 的某个最小因子，满足 p 与 n/p 互质。虽然是随机的，但 PollardRho 算法在实际环境中运行的相当不错，不会被卡。

问题模型：

给一个数 n ，你需要快速求出它的一个非平凡因子。
对于一个比较小的数（$n\le 10^9$），我们直接暴力枚举质因子就行但太大了我们就必须考虑一下随机算法了。
对于一个非常大的合数 $n\le 10^{18}$ （如果可能是质数，我们可以用Miller Rabin判一下）我们要求 n 的某一个非平凡因子，如果 n 的质因子很多（就是约数很多）我们也可以暴力随机弄一弄，但如果是一些（像由两个差不多的而且很大的质数乘得的n）它的非平凡因子就只有两个而且 n 本身还很大，此时如果我们随机的话复杂度 $O(\sqrt{n})$ ，这个太难接受了，所以我们想办法优化一下。

1.求 gcd

直接随机求 n 的某一个约数复杂度很高，可以考虑求一下 gcd 。因为我们可以保证一个合数它绝对存在一个质因子小于 $\sqrt{n}$ ，所以在 n 就存在至少$\sqrt{n}$ 个数与 n 有大于一的公约数。于是我们随机的可能性就提高到了 $O(\sqrt{n}lgn)$ 。

2.生日悖论

生日悖论其实不是一个逻辑悖论，只是与很多人的第一感觉不符罢了。它可以表述为： 一个房间里有23个人，则他们中有两人生日相同的概率超过一半 （不考虑闰年）。其实这在数学上很好证明，即 $\frac{365}{365}\times \frac{364}{365}\times \frac{363}{365}\times \cdots \times \frac{365-22}{365}<\frac{1}{2}$ 。
生日悖论启示我们，如果我们不断在某个范围内生成随机整数，很快便会生成到重复的数，期望大约在根号级别。精确地说，对于一个$[1,N]$内整数的理想随机数生成器，生成序列中第一个重复数字前期望有$\sqrt{\frac{\pi N}{2}}$个数
应用到原题上，即是对于最坏情形$n=p^2$ ，如果我们不断在$[1,n-1]$ 间生成随机数，那么期望在生成大约$\sqrt{p}=n^{\frac{1}{4}}$ 个数后，可以出现两个在模$p$下相同的数（注意$[1,n-1]$ 间的随机数模$p$ 大致是$[0,p-1]$间的随机数）。那么这两个数的差的绝对值 $d$ ，就一定满足$d\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$ ，于是也满足$\gcd (d,n)>1$ 。
但这件事意义并没有那么大。正如虽然生日悖论是正确的，但你不一定能在班上遇到和自己生日相同的人，因为这个高概率是在两两比较下才成立的。对这$n^{\frac{1}{4}}$ 个数两两进行验证，复杂度立刻退回到 $O(\sqrt{n}\lg n)$ ，并没有什么进步。所以我们需要一些技巧。

3.伪随机数序列

Pollard使用一种特别的伪随机数生成器来生成 $[0,N-1]$ 间的伪随机数序列：设序列第一个数为 $x$ ，$f(x):=(x^2+c)\mathrm{m}\mathrm{o}\mathrm{d}N$ ，则$x,f(x),f(f(x)),\cdots$为一个伪随机数序列。
以下为$x=0,c=24,N=9409$时，生成的前100个数：

不过如果稍微换一下参数，把$N$换成$9400$的话……

其实这是显然的，因为每个数都是由前一个数决定的，可以生成的数又是有限的，那么迟早会进入 循环 。当然，这个循环很可能是 混循环 ，所以生成的序列常常形成这样的ρ形，这也是为什么Pollard把这个算法命名为rho算法：
[图片]
我们在这里使用 Floyd判环算法 （也叫 龟兔赛跑算法 ），设置两个变量$t,r$ ，每次判断是否有 $\gcd (|t-r|,N)>1$ ，如果没有，就令$t=f(t)$ ，$r=f(f(r))$ 。因为 r 跑得更快，如果没有找到答案，最终会与 t 在环上相遇，这时退出，换一个 c 重新生成伪随机数。
那么，这有什么好处呢？其实，这个伪随机数生成器生成的数具有一个性质。注意到，如果$|i-j|\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$ ，那么一定有$|f(i)-f(j)|=|i^2-j^2|=|i-j|\cdot |i+j|\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 。由此可得， 只要环上距离为 d 的两个数满足条件，那么所有距离为 d 的数都满足条件 。在Floyd判环的过程中，每次移动都相当于在检查一个新的距离 d ，这样就不需要进行两两比较了。
这个算法的复杂度依赖于这个伪随机数生成器的随机程度，还没有被严格证明。如果它是足够随机的，那么期望复杂度显然是$O(n^{\frac{1}{4}}\log n)$。

4.判环

这种随机生成方法虽然优秀，但也有一些需要注意之处，比如有可能会生成一个环，并不断在这个环上生成以前生成过一次的数，所以我们必须写点东西来判环：

1. 我们可以让y根据生成公式以比x快两倍的速度向后生成，这样当y再次与x相等时，x一定跑完了一个圈且没重复跑多少！
2. 我们可以用倍增的思想，让y记住x的位置，然后x再跑当前跑过次数的一倍的次数。

inline ll rho(ll n){
    ll x, y, c, g, i=0, j=1;
    x=y=rand(); c=rand();
    while(i++){
        x=(qmul(x,x,n)+c)%n;
        if(x==y)break;
        g=gcd(abs(y-x),n);
        if(g>1)return g;
        if(i==j)y=x,j<<=1;
    }
}

PollardRho 算法的二次优化：

我们发现其实 PollardRho 算法其实复杂度不止 $O(n^{1/4})$ ，它每一次操作都会求一次 gcd ，所以复杂度会变成 $O(n^{1/4}\times log)$ 我们发现这个 log 实在有些拖慢速度。于是一个很奇妙的优化诞生了！
二次优化 ：我们将若干个 (y−x) 乘到一起，如果我的若干个 (y−x) 中有一个与n有公约数，最后的结果定然也会含有这个公约数！所以我们完全可以多算几次 (y−x) 的乘积在来求 gcd （一般连续算127次再求一次gcd）
对原算法的一些影响： 任何一个优化都是要考虑其对原算法的影响的。这个优化也会有一些影响：首先，因为我们会等大概127次后再去gcd，然而我们有可能在生成时碰上一个环，我们有可能还没生成127次就跳出这个环了，这样就无法得出答案；其次，我们的 PollardRho 算法实在太玄学优秀了，可能我们跑127次之后，所有 (y−x) 的乘积就变成了n的倍数（模 n 意义下得到 0 ），所以我们不能完全就只呆板的等127次在取模。
我们可以用倍增，分别在生成（1次，2次，4次，8次，16次，32次…）之后进行 gcd ！这样就可以完全避免上述的两个影响

inline ll rho(ll p){
    ll x, y, z, c, g, i;
    while(1){
        y = x = rand()%p;
        z = 1;
        c = rand()%p;
        i=0, j=1;
        while(i++){
            x = ((x*x%p)+c)%p;
            z = z*abs(y-x)%p;
            if(x==y || !z)break;
            if(i%127==0 || i&(i-1)==1){//分别在i=2^n 和 i%127==0时求gcd
                g = __gcd(z,p);
                if(g > 1)return g;
                if(i == j)y=x, j<<=1;
            }
        }
    }
}

总参考代码

inline ll qmul(ll x,ll y,ll mod){
    return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;     
}

ll qpow(ll a,ll n,ll mod)/{
    ll res=1;
    while(n){
        if(n&1) res=qmul(res,a,mod);
        a=qmul(a,a,mod);
        n>>=1;
    }
    return res;
}

bool Miller_Rabin(ll n){
    if(n<3 || n&1==0) return n==2;
    ll u=n-1, t=0;
    while(u&1 == 0) u/=2, t++;
    ll ud[]={2,325,9375,28178,450775,9780504,1795265022};
    for(ll a:ud){
        ll v=qpow(a,u,n);
        if(v==1 || v==n-1 || v==0) continue;
        for(int j=1; j<=t; j++){
            v = qmul(v,v,n);
            if(v==n-1 && j!=t){ v=1; break; }
            if(v == 1) return 0;
        }
        if(v != 1) return 0;//Fermat检验
    }
    return 1;
}

inline ll rho(ll p){
    ll x, y, z, c, g, i, j;
    if(Miller_Rabin(p)) return p;
    while(1){
        y = x = rand()%p;
        c = rand()%p + 2;
        i = 0, j = 1, z = 1;
        while(i++){
            x = (qmul(x,x,p) + c) % p;
            z=qmul(z, abs(y-x), p);
            if(x==y || !z) break;
            if(i%127==0 || i&(i-1)==1){//分别在i=2^n 和 i%127==0时求gcd
                g = __gcd(z,p);
                if(g > 1) return g;
                if(i == j) y=x, j<<=1;
            }
        }
    }
}