Mathematics(2):二次剩余与Pollard's Rho 2019-05-10

Quadratic residue

二次剩余

定义:

当存在某个\(x\),式子\(x^2\equiv d \pmod p\)成立,称"\(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)\)

E2IZWR.png 如上图所示,图中的\(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) {} // 同上
}