数论算法小结

Miller Rabin

费马测试：

根据费马小定理 $a^{p-1} \equiv 1\mod p$，其中 $p$ 为质数。于是我们猜想若对于一个 $a$ 有 $a^{p-1} \equiv 1\mod p$ 成立则 $p$ 为质数。但这个猜想是错误的，我们只能知道若 $a^{p-1} \mod p$ 不为 $1$ 时，$p$ 一定是合数；若 $a^{p-1}$ 模 $p$ 为 $1$ 时，$p$ 只是有很大概率是素数。

因此我们可以想到可以随机地多选取一些 $a$ 对 $p$ 进行费马测试，以提高准确率。可惜的是存在一些合数比如 $516$，能够通过所有的费马测试，为了更高的准确率，在费马测试的基础上有了 Miller-Rabin 质数测试算法。

二次探测：

有一个定理是：如果 $p$ 是奇素数，则 $ x^2 ≡ 1\mod p$ 的解为 $x ≡ 1$ 或 $x ≡ p – 1$。那么如果有一个 $x$，它使得 $x^2$ 模 $p$ 不为 $1$ 或 $p-1$，那么 $p$ 就一定不为素数。

因为 $a^{p-1} ≡ 1\mod p$，所以将指数 $p-1$ 分解为 $2^k*x$ 的形式后，根据上面定理可以倒着推出 $a^{2^{k-1}*x} $ 在模 $p$ 下为 $1$ 或 $p-1$，如果该数为 $1$ 则可以继续往下推，直到一个数模 $p$ 为 $p-1$。反而言之，如果 $a^x$ 不为 $1$ 且从 $a^{x},a^{2x},a^{4x},…$ 开始，到 $a^{p-1}$ 都没有出现 $p-1$ 那么 $p$ 就是合数。

如果将上面两种测试合起来，多随机一些 $a$ 进行测试，这样我们就可以以很高的正确率来判断一个数是否是素数，而这个算法就是 Miller-Rabin 测试。

代码：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#include <vector>
#include <map>
#include <set>
#define rint register int
#define LL long long
#define I128 __uint128_t
#define LD long double
using namespace std;

const LL base[5]={2, 3, 7, 61, 24251};

int T;
LL n;

LL ksm(LL x, LL y, LL p)
{
    LL sum=1;
    while(y)
    {
        if(y&1) sum=(I128)x*sum%p;
        y>>=1, x=(I128)x*x%p;
    }
    return sum;
}

bool MR(LL p)
{
    LL x=p-1, k=0;
    while(!(x&1)) x>>=1, k++;
    for(rint i=0; i<5; ++i)
    {
        if(p==base[i]) return true;
        LL t=ksm(base[i], x, p);
        if(t==1 || t==p-1) continue;
        int tag=0;
        for(rint j=1; j<=k; ++j, t=(I128)t*t%p)
            if(t==p-1) {tag=1; break;}
        if(!tag) return false;
    }
    return true;
}

int main()
{
    scanf("%d", &T);
    while(T--)
    {
        scanf("%lld", &n);
        if(MR(n)) printf("Yes\n");
        else printf("No\n");
    }
    return 0;
}

Pollard Rho

首先，十分感谢这篇文章的作者，这篇文章给了我很大的帮助。

假设我们要分解的数 $n$ 能被写成 $p * q$ 的形式，那么有一种算法就是不断地取一个随机数，直到它为 $p$ 或 $q$，然后递归地下去再找 $p, q$ 的因数，如果这个因数为质数，那么就可以停止递归且这个数就是 $n$ 的一个质因数。如果 $n$ 只有两个因子的话，根据这个算法找到一个因子的概率就是 $\frac{2}{n}$，是完全不可接受的。于是我们接下来就考虑如何提高这个概率。

生日悖论：

如果一个班有 23 名学生，那么两个学生同一天生日的概率就超过了 $50\%$。类似的，如果在 $[2,n]$ 范围内有一个我们要找的数 $p$，我们大约随机 $\sqrt n$ 个数，那么其中有一对数的差等于 $p$ 的概率也超过了 $50\%$。这样，我们将这个概率由原来的 $\frac{1}{n}$ 级别提升至了 $\sqrt \frac{1}{n}$ 级别。 但是如果我们枚举随机出来的所有数对，并判断它们的差是否能被 $n$ 整除，这样的复杂度是 $O(\sqrt n^2) $ 即 $O(n)$ 的，依然有优化的余地。

考虑原来一对数 $(x_i-x_j) | n$ 的情况很少，如果我们改变判断的标准，判断 $gcd(x_i-x_j, n)$ 是否大于 $1$ ，这样的情况数就大大增加了，这个算法在期望意义下随机选取 $n^{1/4}$ 个数就会存在一对数 $gcd(x_i-x_j, n)$ 大于 $1$，是比较优秀的。

Pollard rho 算法:

而按照上面提到的方法，将随机出的 $n^{1/4}$ 个数储存起来，空间开销又会很大，Pollard-rho 算法就是解决了这个问题。考虑这样一个随机函数 $f(x)=(x^2+a) \mod n$，其中 $a$ 作为随机种子。我们从 $x_0$ 开始，利用 $x_{i+1}=f(x_i)$ 生成 $x_{i+1}$ 并判断数对 $x_{i}, x_0$ 是否合法，这样我们只需要储存两个数就好了。

那么问题又来了，按照这个随机方法一段时间以后随机出来的数就会产生环，如果出现了环，我们就需要修改随机种子，以免一直在环上走。

判环：

有一种 $Floyd$ 发明的判环方法，而这里讲一种似乎更为 玄学 高效的判环方法。设定一个上限 $k$，然后以 $x_0$ 为初值生成随机数，如果生成 $k$ 个数之后没有和 $x_0$ 相等的数，那么调整 $k$ 并重新设置 $x_0$；否则就已经在一个环上了。

代码：
指数形式分解质因数

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#include <vector>
#include <map>
#include <set>
#include <iostream>
#define MAXN 10005
#define rint register int
#define LL __uint128_t
#define LD long double
using namespace std;

inline void read(LL &x)
{
    char ch=getchar(); x=0;
    while(ch<'0' || ch>'9') ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0', ch=getchar();
}

const LL base[5]={2, 3, 7, 61, 24251};

int T, id, num[MAXN];
LL n;

map<LL, int> mp;
map<LL, int>::iterator it;

LL gcd(LL x, LL y)
{
    return y?gcd(y, x%y):x;
}

LL ksm(LL x, LL y, LL p)
{
    LL sum=1;
    while(y)
    {
        if(y&1) sum=x*sum%p;
        y>>=1, x=x*x%p;
    }
    return sum;
}

bool MR(LL p)
{
    LL x=p-1, k=0;
    while(!(x&1)) x>>=1, k++;
    for(rint i=0; i<5; ++i)
    {
        if(p==base[i]) return true;
        LL t=ksm(base[i], x, p);
        if(t==1 || t==p-1) continue;
        int tag=0;
        for(rint j=1; j<=k; ++j, t=t*t%p)
            if(t==p-1) {tag=1; break;}
        if(!tag) return false;
    }
    return true;
}

LL PR(LL n)
{
    LL x=1LL*rand()*rand()%(n-2)+1;
    LL y=x, c=1LL*rand()*rand()%n+1;
    for(rint i=1, k=2; ; ++i)
    {
        x=(x*x%n+c)%n;
        LL d=gcd(y-x, n);
        if(1<d && d<n) return d;
        if(x==y) return n;
        if(i==k) y=x, k<<=1;
    }
}

void find(LL x, int t)
{
    if(x==1) return;
    if(MR(x))
    {
        if(!mp[x]) mp[x]=++id;
        num[mp[x]]+=t;
        return;
    }
    LL y=x;
    int temp=0;
    while(x==y) y=PR(x);
    while(x%y==0) x/=y, temp++;
    find(x, t); find(y, temp*t);
}

int main()
{
    scanf("%d", &T);
    while(T--)
    {
        memset(num, 0, sizeof(num));
        mp.clear();
        read(n);
        find(n, 1);
        for(it=mp.begin(); it!=mp.end(); it++)
            printf("%lld %d\n", (long long)it->first, num[it->second]);
    }
    return 0;
}

BSGS

BSGS 算法是用来求解形如 $a^x \equiv y \mod p $ 的高次同余方程的算法。

首先根据欧拉定理的扩展，可以知道 $a^{x \% \varphi(p)} \equiv a^x$。所以如果方程有解，在 $[0, \varphi(p)]$ 范围内也一定会有解。考虑分块，设 $num=\sqrt p$ ，并将值域分为大小为 $num$ 块。预处理 $y*a^i \mod p$ $(i \in [1, num])$ 的值，并放入 map 中。然后可以将式子写成 $a^{k*num} \equiv y*a^i \mod p$ 的形式，接着从小到大枚举 $k$，在 $map$ 中找是否有对应 $a^{k*num}$ 的值 $i$，如果存在那么方程的最小非负解为 $x=k*num-i$。

这样的复杂度就是 $O(\sqrt p * \log\sqrt p)$，如果手写 hash 表代替 map 可以把复杂度的 $\log$ 去掉。

代码：
BZOJ2242 计算器 AC 代码。

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#include <vector>
#include <map>
#include <set>
#define MAXN
#define INF 0x3f3f3f3f
#define rint register int
#define LL long long
#define LD long double
using namespace std;

int T, k, x, y, z, p;

map<int, int> mp;

int ksm(int x, int y)
{
    int sum = 1;
    while (y)
    {
        if (y & 1)
            sum = 1LL * sum * x % p;
        x = 1LL * x * x % p;
        y >>= 1;
    }
    return sum;
}

int exgcd(int a, int b, int &x, int &y)
{
    if (!b)
    {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, x, y);
    int t = x;
    x = y, y = (t - a / b * y);
    return d;
}

int bsgs(int y, int z, int p)
{
    if (y % p == 0)
        return (z == 0) ? 1 : -1;
    int num = sqrt(p) + 0.5, base = 1;
    mp.clear();
    int t = z % p;
    mp[t] = 0;
    for (rint i = 1; i <= num; ++i)
    {
        t = 1LL * t * y % p;
        mp[t] = i;
        base = 1LL * base * y % p;
    }
    t = 1;
    for (rint i = 1; i <= num; ++i)
    {
        t = 1LL * t * base % p;
        if (mp.count(t))
        {
            x = i * num - mp[t];
            return (x % p + p) % p;
        }
    }
    return -1;
}

int main()
{
    scanf("%d%d", &T, &k);
    while (T--)
    {
        scanf("%d%d%d", &y, &z, &p);
        if (k == 1)
            printf("%d\n", ksm(y, z));
        if (k == 2)
        {
            int t, d = exgcd(y, p, x, t);
            if (z % d)
            {
                printf("Orz, I cannot find x!\n");
                continue;
            }
            x = (1LL * x * z / d % (p / d) + (p / d)) % (p / d);
            printf("%d\n", x);
        }
        if (k == 3)
        {
            x = bsgs(y, z, p);
            if (x == -1)
                printf("Orz, I cannot find x!\n");
            else
                printf("%d\n", x);
        }
    }
    return 0;
}