Pollard's Rho

Quadratic residue

二次剩余

定义：

当存在某个$x$，式子$x^2\equiv d \pmod p$成立，称"$x$是$d$是模$p$下的二次剩余"。

性质：

• 对于一个数$0 \sim p - 1$的数，因为$x^2 \equiv (p-x)^2$，所以关于$p$的二次剩余的不超过$\frac{p}{2}+1$个。
• 二次剩余成二次剩余也是二次剩余

Miller-Rabin

考虑一种判断质数的方法，根据费马小定理，$a^{p-1} \equiv 1 \pmod p$。为了检测$p$是否为质数，可以带入几个数，判断费马小定理是否成立，如果都成立，那么可以认为这个数是质数。

但是有一些数没办法判断出来，如$561=3\times11\times17$，这一类数叫做Carmichael数（卡迈尔克数）。这一类数虽然不多，但是在$1e8$内还是有$255$个，会严重影响算法的正确率，Miller-Rabin这个算法便诞生了。

设要检测的数为$n$，$n$为奇数。

$$n=r\times 2^t + 1$$

根据$\text{Chinese Remainder Theorem}$，如果$n$是质数，$1$的在模$n$下二次剩余只有$1$和$n-1$。例如$6^2\equiv -1\pmod {37}$，$6^3 \equiv 1 \pmod {37}$，反之对于非质数会有更多其他的数，例如$11 \not \equiv -1$，$11^2 \equiv 1 \pmod {15}$，说明$15$是非质数。

对于不被$n$整除的一个数$A$，如果$n$为质数，根据费马小，有$A^{n-1} \equiv 1\to A^{r\times 2^t} \equiv 1\pmod n$。

考虑这个$1$是怎么来的，那么根据二次剩余，下面两个条件之一必然成立。

• $$A^r \equiv1 \pmod n$$
• $$\exists_{k=0}^{t-1}A^{r2^k} \equiv -1 \pmod n$$

于是我们检测到第一个$A^{r2^{k+1}}\equiv1$，那么$A^{r2^{k}} \equiv -1 \pmod n$。如果不是，则说明$n$不是质数。所以可以随机一个底数$A$，或者是取前几个小质数作为$A$。

对于一个奇合数$x$，它的证据个数有$\frac{x-1}{2}$个，对于一个证据，检测出的概率是$\frac{1}{2}$。如果选择了$n$个底数，那么正确的概率约为$2^{-n}$。

对于$1e18$的范围的数，取前12个质数就可以正确了。

int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 21, 23, 29, 31, 37};
IL bool Miller_Rabin(ll n) {
    if (n <= 1) return false;
    for (int i = 0; i < 9; ++i) if (n % pri[i] == 0) return n == pri[i];
    ll r, x, y; int t;
    for (r = n - 1, t = 0; ~r & 1; r >>= 1, t++);
    for (int i = 0; i < 9; ++i) {
        x = fpw(pri[i], r, n); 
        for (int j = 0; j < t && x > 1; ++j) {
            y = mul(x, x, n);
            if (y == 1 && x != n - 1) return false; 
            x = y;
        }
        if (x != 1) return false;
    }
    return true;
}

Birthday problem

生日悖论，23个人中两个人有相同生日的概率大于$50\%$。

对于值域为$P$，出现至少两个相同的数的概率为$1-\frac{P-1}{P}^{\binom{n}{2}}$

Pollard’s Rho

对于一个非质数$n$，存在$a_1$和$a_2$使得$n=a_1\times a_2$，那么递归分解$a_1$和$a_2$就可以实现质因数分解。

假设$a_1$是$n$的因子，暴力地可以随机一个数然后判断是否为$a_1$，这样找到$a1$的概率为$\frac{1}{n}$实在太小。如果得到了两个数$x$和$y$，其中$x\equiv y \pmod {a_1}$且$x\not \equiv y \pmod n$。如果$gcd(x−y,n)$不为$1$或$n$，那么可以得到 $n$ 的一个非平凡因子。

往一个纸带上填数，值域为$P$，出现有两个数就停下（显然步数最大为$P+1$)，期望的长度是多少？一共存在$\binom{len}{len-1}$对数，而每队相等的概率为$\frac{1}{p}$

$$len\times (len-1)\times \frac{1}{P}=1\\len\approx\sqrt{P}$$

我们可以算出$len$的期望长度为$\sqrt P$，也就是在模$n$意义下的期望长度为$\sqrt n$。

对于一组$(x,y)$，$y\equiv x\pmod a$的概率为$O(\frac{1}{n})$，根据生日悖论，如果选出大于$\sqrt n$组数，其中有两个相等的概率大于$\frac{1}{2}$，而后会更加快速地递增。

如果我们每次合理的随机选择$x,y$，使得$x-y$不相等，在大小为$n$的集合中，选出合理的$(x,y)$概率为$\sqrt n$，因为集合大小为$\sqrt n$，那么找到相等的概率为$\sqrt {\sqrt n}$，也就说$O(n^{\frac{1}{4}})$。

我们构造一个随机的函数$f(x)$使得$f(x)$在$\text{mod p}$意义下几乎随机，不断将$x$变为$f(x)$，$x$就会变为一个随机的数，经过测试，当$f(x)=x^2+c$的时候很合适。经过若干步以后一定会走到一个定长的环中。最开始$x,y$设为随机数，在每一步的操作中把$x$变为$f(x)$，$y$变为$f(f(y))$，检查$gcd(y-x,n)$。这样$y$每次比$x$多走一步，经过环长步以后$x,y$会到达相同位置。

因为对于每一步都要求一遍$gcd$，那么复杂度为$O(n^\frac{1}{4}\log n)$。

IL ll func(ll x, ll mod, ll a) { return (mul(x, x, mod) + a) % mod; }
IL ll Find(ll n) {
    static const int Test_Limit = 150000;
    ll a = rand(), x, y; int tim = 0;
    for (int i = 0; i < 12; ++i) if (n % pri[i] == 0) return pri[i];
    x = func(rand(), n, a), y = x;
    do {
        ull g = gcd((x - y + n) % n, n);
        if (g != 1 && g != n) return g; 
        x = func(x, n, a), y = func(func(y, n, a), n, a); ++tim; 
    } while (y != x && tim <= Test_Limit);
    return -1;
}
void Pollard_s_Rho(ll n, ll &a) {
    ll d; if (n <= a) return;
    while (!(n & 1)) n >>= 1, a = 2;
    if (n == 1 || Miller_Rabin(n)) { a = max(a, n); return; }
    for (d = Find(n); d == -1; d = Find(n));
    if (d > n / d) d = n / d; Pollard_s_Rho(d, a); Pollard_s_Rho(n / d, a);
}

Pollard’s Rho Update

但这并不是$\text{Pollard's Prho}$能达到的最优复杂度，可以优化为$O(n^\frac{1}{4})$。

首先$\gcd$的写法可以加上一点常数优化。

IL ll gcd(ll x, ll y) {
    int Shift = ctzll(x | y);
    y >>= ctzll(y);
    while (x) {
        x >>= ctzll(x); 
        if (x < y) swap(x, y);
        x -= y;
    }
    return y << Shift;
}

根据奇偶性，每隔两次$x$至少要除$2$，而换成非递归会快一点。

我们可以注意到$\gcd(t,n) | gcd(tq,n)$， 设定一个步长$M$，固定$x$，$y$一共跳$M$步，把每一步的$y-x$相等然后取$gcd$。虽然还是遍历了$\sqrt n$个元素，但是每隔$M$次取一次gcd可以降低复杂度，使得它为$O(\frac{n^{\frac{1}{4}}}{M}\times \log n)$。

如上图所示，图中的$Y$将从$X$跳到$2^{l+1}$的位置，因为之前已经算过蓝色部分并没有得到结果，那么因为图中的红色部分和蓝色部分等长，重复算红色部分没有意义，$Y$可以直接跳到现在这个位置；$Y$现在从图中$Y$的位置开始，每次跳$M$步，把每一步的$y-x$乘起来然后和$n$取$\gcd$检验。同样的，如果不是$1$或者是$n$说明找到了直接return；如果是$1$，视为没有找到；如果是$n$，可能这个$n$是由前面的几个数相乘得到的，也可能是走完了一个环，此时$x=y,gcd(y-x,n)=n$，那么我们需要再走一遍$M$这么长，判断中间是否出现结果，找到答案然后返回。

具体实现看代码（unsigned long long和__uint128_t真快)。

#define ctzll __builtin_ctzll
namespace Rho {
    IL ll gcd(ll x, ll y) {
        int Shift = ctzll(x | y);
        y >>= ctzll(y);
        while (x) {
            x >>= ctzll(x);
            if (x < y) swap(x, y);
            x -= y;
        }
        return y << Shift;
    }
    IL ll func(ll x, ll Mod, ll c) { return (mul(x, x, Mod) + c) % Mod; }
    IL ll Find(ll n) {
        static int Step = 1 << 9; int c = rand();
        ll x, y, temp; x = y = rand() % n; 
        for (int l = 1; ; l <<= 1) {
            x = y; 
            for (int i = 0; i < l; ++i) y = func(y, n, c);
            for (int k = 0; k < l; k += Step) {
                ll Prod_g = 1;
                temp = y;
                int rem = min(l - k, Step);
                for (int i = 0; i < rem; ++i) y = func(y, n, c), Prod_g = mul(Prod_g, (y + n - x) % n, n);
                Prod_g = gcd(Prod_g, n);
                if (Prod_g == 1) continue;
                if (Prod_g == n) for (y = temp, Prod_g = 1; Prod_g == 1; ) y = func(y, n, c), Prod_g = gcd(y + n - x, n); 
                return Prod_g;
            } 
        }           
    }   
    IL void Pollard_s_Rho(ll n) {} // 同上
}