学习 Pollard-Rho

前言

• woc这月赛是啥a
• msjing在写 Ynoi 需要用这个于是来学了
• pb好强强强TTC好强强强

\(\mathcal{Miller-Rabin}\) 算法

• 这是一个判素数的算法，想要学 \(\mathcal{Pollard-Rho}\) 就要先学这个
• 不是知更鸟（\(\mathcal{robin}\)）哦
• 学这个要先学费马小，这个百度一下就行了
• \(\mathcal{Rabin}\) 算法本质是一种随机化算法，不过我们要先介绍一个错误判素数的方法，就是上面说的费马小
• 我们知道费马小长这样：

\[a^{p-1} \equiv 1 \>\>\> (\!\!\!\!\!\mod p) \]

• 其中 \(p\) 为素数
• 我们发现貌似使这个成立的 \(p\) 都是素数
• 但是是假的，因为一类数被构造了

\(\mathcal{Carmichael}\) 数

• 这个数的定义是：对于一个合数 \(n\)，所有与 \(n\) 互质的正整数 \(a\)，都满足费马小，那称 \(n\) 为 \(\mathcal{Carmichael}\) 数

• 然后费马小就假了

• 所以我们要换一个了

• 但是其实是没啥好办法了，于是请出随机化

• 不过我们要先证一个东西

• 我们设 \(x\) 在 \(p\) 的剩余系中，当 \(p\) 为素数时，\(x\) 只有两个解 \(1\) 或 \(p-1\)，这个是充要的，即：

\[x^2 \equiv 1 \>\>\>(\!\!\!\!\! \mod p) \Longleftrightarrow x = 1 \>\>\> or \>\>\> p-1 \]

• 这个叫二次探测定理
• 我们证明一下（TTC好强谢谢你喵）

证明

• 先证必要性
• \(x = 1\) 是显然成立的
• 我们把 \(x = p-1\) 带入，有：

\[(p-1)^2 \]

• 完全平方展开有：

\[p^2 - 2p + 1 \]

• 提出 \(p\)，有：

\[p(p - 2) + 1 \]

• \(\!\!\!\mod p\)，余 \(1\)，证毕
• 再证充分性
• 我们把柿子换个形式：

\[x^2 = k \times p + 1 \]

• 其中，\(k\) 是一个整数
• 移项，有：

\[x^2 - 1 = k \times p \]

• 即：

\[x^2 - 1 \mid p \]

• 平方差展开，有：

\[(x + 1)(x - 1) \mid p \]

• 即：

\[(x + 1) \mid p \>\>\> or \>\>\> (x - 1) \mid p \]

• 前面的是 \(x = p + 1\)，后面的是 \(x = p - 1\)，由于定义 \(x\) 在 \(p\) 剩余系内，所以前面解为 \(x = 1\)，证毕

• 有了这个，我们就可以判素数了
• 具体做法很神奇，我们直接说步骤
• 枚举 \(k\) 个数满足 \(1 < a_i < p\)，带入验证
• 如何带入验证
• 我们令 \(p - 1 = 2^k \times b\)，\(k\) 尽可能大，\(b\) 为整数
• 我们可以拆出 \(p - 1\) 的约数个数为 \(k\)，将一个 \(a\) 带入，并与拆剩下的数做快速幂运算，结果设为 \(s\)，如果 \(s = 1\)，成立，\(p\) 为素数
• 不成立，枚举 \(k\) 次，每次让 \(s\) 自乘，就是把幂再累加上去，每次判 \(s\) 是否等于 \(p - 1\)，成立，\(p\) 为素数，直接出，枚举完 \(k\) 后还不成立，\(p\) 不是素数，注意先判再自乘
• \(a\) 的个数 \(k\) 大概在 \(8\)，直接枚举素数就行，记得特判一下 \(p\) 是否是 \(a\) 中的数
• 这个算法就没了
• 欸，但是这样不会出错吗
• e，事实上，在 long long 内的数这个出错概率基本为 \(0\)，实际出错率为每次最多 $ 1 \over 4$
• 复杂度 \(O(k \log n)\)

板子SP288 PON - Prime or Not

点击查看代码

• 可能会爆，所以用龟速乘

#include <bits/stdc++.h>
#define int long long
using namespace std;
constexpr int maxn=2e6+10;
int read()
{
    int x=0,f=1;
    char ch=getchar();
    while (ch<'0' || ch>'9')
    {
        if (ch == '-') f=-1;
        ch=getchar();
    }
    while (ch>='0' && ch<='9')
    {
        x=(x<<1)+(x<<3)+ch-'0';
        ch=getchar();
    }
    return x*f;
}
int a[9]={0,2,3,5,7,11,13,17,19};
int mul(int x,int y,int p)
{
    int res=0;
    while (y)
    {
        if (y&1) (res+=x)%=p;
        (x+=x)%=p;
        y>>=1;
    }
    return res;
}
int power(int x,int y,int p)
{
    int res=1;
    while (y)
    {
        if (y&1) res=mul(res,x,p);
        x=mul(x,x,p);
        y>>=1;
    }
    return res;
}
int chk(int p,int a)
{
    int k=p-1,cnt=0;
    while (!(k&1))
        cnt++,k>>=1;
    int x=power(a,k,p);
    if (x == 1) return 1;
    for (int i=1;i<=cnt;i++)
    {
        if (x == p-1) return 1;
        x=mul(x,x,p);
    }
    return 0;
}
int prime(int x)
{
    for (int i=1;i<=8;i++)
    {
        if (x == a[i]) return 1;
        if (!(x%a[i])) return 0;
        if (!chk(x,a[i])) return 0;
    }
    return 1;
}
signed main()
{
    int T=read();
    while (T--)
    {
        int n=read();
        prime(n) ? puts("YES") : puts("NO");
    }
    return 0;
}

• \(\mathcal{Rabin}\) 就没了，接下来是 \(PR\)

• 不过我们在学这个之前，要先学 \(\mathcal{Floyd}\) 判圈算法

\(\mathcal{Floyd}\) 判圈算法

• \(\mathcal{Floyd}\) 判圈算法，又称龟兔赛跑算法，由美国计算机科学家 \(\mathcal{Robert\>\>\>W.Floyd}\) 提出，那个最短路也是他提出的
• \(\mathcal{Floyd}\) 判圈算法一般用于判链表是否存在环
• 那 \(\mathcal{Floyd}\) 判圈算法是什么原理呢？
• 我们开两个指针，一个快一个慢，快的每次走两步，好比兔子，慢的每次走一步，好比乌龟
• 当没有环时，慢指针一定追不上快指针，而有环时，快指针会在环内追上慢指针
• 实现还好，msjing不打算写了
• 复杂度线性

\(\mathcal{Pollard-Rho}\)

• 这个算法可以以 \(O(n^{1 \over 4})\) 的期望时间算出一个 \(n\) 的非平凡因子

非平凡因子：若 \(x\) 能整除 \(n\) 且 \(1<x<n\) ，则称 \(x\) 是 \(n\) 的非平凡因子

• 这个算法本质是通过一个迭代的函数（比如 \(x_i = f(x_{i-1})\)）在取模（模数为 \(n\) 的最小质因子，设为 \(p\)）意义下进行迭代，由于有取模限制，所以 \(x\) 的取值是有限的，因此它会进一个环，期望进环时间为 \(O(\sqrt p)\)
• 当有两个数满足：\(x_i \equiv x_j \>\>\>(\!\!\!\! \mod p)\) 时，如果 \(x_i \neq x_j\)，那就是找到环了，同时，有 \(p \mid \gcd(|x_i ​− x_j​|,n)\)，如果 \(\gcd(|x_i ​− x_j​|,n) > 1\)，\(p\) 就是 \(n\) 的一个非平凡因子
• \(\mathcal{Pollard-Rho}\) 算法通过一个序列 \(a_i = a_{i-1}^2 + c \>\>\>(\!\!\!\! \mod n),a_0 = x\) 来实现
• 我们通过随机参数，发现点的分布有周期性（desmos驯服失败，图没搓出来）
• 那对于一些值，会形成有个尾巴的环，和 \(\rho\) 很像
• 期望进环的时间证明需要生日悖论，介绍一下

生日悖论

• 在不少于 \(23\) 人的群体中出现相同生日的概率超过 \(50%\)，\(60\) 人时概率可达 \(99%\)
• 其函数为：

\[f\left(x\right)\ =\ 1-\ e^{-\frac{x\left(x-1\right)}{2 \times N}} \>\>(N = 365) \]

• 可以用 desmos 画一下

• 比例没设好

• woc我真要证这个吗
• 对于这个问题，用生日悖论可以证明期望进环时间，但msjing不会证
• 那你直接随机参数带入验证，迭代计算，就可以得到结论
• 但是你发现有可能会由于运气脸黑而导致找不到环，所以我们用上面的 \(\mathcal{Floyd}\) 判环
• 复杂度是 \(O(\sqrt p)\) 的
• 不过你要把 \(4\) 判掉，要不然跑不出来
• 但是算法复杂度还是有点高，我们发现 \(\gcd\) 有点慢
• 发现由于算法是与 \(n\) 取 \(\gcd\)，所以我们发现如果将 \(|x_i - x_j|\) 累加起来再做 \(\gcd\) 是不影响的，求出来的数仍然是 \(n\) 的非平凡因子
• 累加阈值可以调，一般用 \(128\)
• 整个算法就没有了，我们去把板子写了

板子P4718 【模板】Pollard-Rho

• 找个好看的代码贺一下

点击查看代码

• \(\mathcal{Rabin}\) 算法的伪随机 \(a\) 个数设为 \(8\) 会掉一个，设为 \(9\) 是对的，msjing设为 \(8\) 掉了
• __int128 好像没法用 \(abs\)，所以手写了一个

#include <bits/stdc++.h>
#define int __int128
using namespace std;
constexpr int maxn=2e6+10;
int read()
{
    int x=0,f=1;
    char ch=getchar();
    while (ch<'0' || ch>'9')
    {
        if (ch == '-') f=-1;
        ch=getchar();
    }
    while (ch>='0' && ch<='9')
    {
        x=(x<<1)+(x<<3)+ch-'0';
        ch=getchar();
    }
    return x*f;
}
namespace MR
{
    int a[11]={0,2,3,5,7,11,13,17,19,23,29};
    int power(int x,int y,int p)
    {
        int res=1;
        while (y)
        {
            if (y&1) (res*=x)%=p;
            (x*=x)%=p;
            y>>=1;
        }
        return res;
    }
    int chk(int p,int a)
    {
        int k=p-1,cnt=0;
        while (!(k&1))
            cnt++,k>>=1;
        int x=power(a,k,p);
        if (x == 1) return 1;
        for (int i=1;i<=cnt;i++)
        {
            if (x == p-1) return 1;
            (x*=x)%=p;
        }
        return 0;
    }
    int prime(int x)
    {
        for (int i=1;i<=10;i++)
        {
            if (x == a[i]) return 1;
            if (!(x%a[i])) return 0;
            if (!chk(x,a[i])) return 0;
        }
        return 1;
    }
}using namespace MR;
namespace PR
{
    mt19937 rd(time(0));
    int f(int x,int c,int p) {return (x*x%p+c)%p;}
    int _abs(int x,int y)
    {return x-y<0 ? y-x : x-y;}
    int rho(int n)
    {
        int mo=n;
        int x=rd()%mo,c=rd()%(mo-1)+1,p=1;
        for (int i=2,j=2,d=x;i;i++)
        {
            x=f(x,c,mo);
            (p*=_abs(x,d))%=mo;
            if (!(i%127) && __gcd(p,n)!=1)
                return __gcd(p,n);
            if (i == j)
            {
                j<<=1,d=x;
                if (__gcd(p,n)!=1)
                    return __gcd(p,n);
            }
        }
    }
    int pr(int n)
    {
        if (prime(n)) return n;
        int p=n;
        while (p == n) p=rho(n);
        return max(pr(p),pr(n/p));
    }
}using namespace PR;
signed main()
{
    int T=read();
    while (T--)
    {
        int n=read();
        int ans=pr(n);
        if (ans == n) puts("Prime");
        else printf("%lld\n",(long long)ans);
    }
    return 0;
}

例题P5071 [Ynoi Easy Round 2015] 此时此刻的光辉

• 就是因为这个我才学的
• 这里看
• 这真是个黑科技我瞎搞做法能过

后记

• 没了，写 Ynoi 去了
• upd on 8.22 21:30：写完了，来补例题