素性测试和质因数分解

Miller-Rabin 素性测试

对于质数 \(p>2\)，方程 \(x^2\equiv 1\pmod p\) 有且仅有两个不同的解 \(x\equiv 1,-1\pmod p\)。证明考虑移项 \((x+1)(x-1)\equiv 0\pmod p\)。然而，如果模数不是质数，那么解可能多于两个，例如 \(3^2\equiv 1\pmod 8\)。

由费马小定理可以得知，若 \(p\) 是质数，则对任意的正整数 \(a\in [1,p-1]\) 有 \(a^{p-1}\equiv 1\pmod p\)。然而费马小定理的逆定理并不成立，最小的反例为 \(p=561\)。

考虑结合以上两者，分解 \(p-1=r\times 2^t\) 其中 \(r\) 是奇数。考虑数列 \(a^r,a^{2r},a^{4r},a^{8r},\cdots,a^{p-1}\)，由于序列最后一项等于 \(1\)，因此序列形如一个后缀等于 \(1\)，前一个数为 \(-1\pmod p\)。当然也可以全是 \(1\)。Miller-Rabin 指出，只需要对 \(a=2,3,5,7,11,13,17,37\) 总共 \(8\) 个数检验序列是否形如后缀一段 \(1\)，前一个为 \(-1\)，或者全是 \(1\)，就可以覆盖 long long 范围内的所有合数。

这样做时间复杂度为 \(O(8\log p)\)。

代码

bool is_prime(ll x) {
    if(x <= 1) return false;
    static ll base[] = {2, 3, 5, 7, 11, 13, 17, 37};
    ll t = 0, r = x - 1; while(~r & 1) t++, r >>= 1;
    for(ll a : base) {
        if(a == x) return true;
        a = qpow(a, r, x);
        if(a == 0) return false;
        if(a == 1) continue;
        for(int i = 1; i < t; i++) {
            if(a == x - 1) break;
            a = mul(a, a, x);
        }
        if(a != x - 1) return false;
    }
    return true;
}

Pollard-Rho 质因数分解

Pollard-Rho 可以在期望 \(O(n^{1/4})\) 的复杂度内给出一个 \(n\) 的非平凡（非 \(1\) 非 \(n\)）因子。

考虑函数 \(f(x)=x^2+C\) 其中 \(C\) 是一个常数。此时假如我们要分解一个大数 \(n\)（\(n\) 是合数），那么考虑数列 \(a,f(a),f(f(a)),\cdots\pmod n\)。设 \(d\) 是 \(n\) 的一个非平凡因子，不妨设 \(d\le \sqrt n\)。我们会发现数列 \(a,f(a),f(f(a)),\cdots\pmod d\) 会在期望 \(O(\sqrt d)\) 的时间内出现第一个重复数字。由于 \(\forall x_1\equiv x_2\pmod d,~f(x_1)\equiv f(x_2)\pmod d\)，因此接下来的数字将会重复之前的序列。

如果我们找到了这个位置，即找到了 \(p_1,p_2\) 满足 \(f^{(p_1)}(a)\equiv f^{(p_2)}(a)\pmod d\) 但 \(f^{(p_1)}(a)\not\equiv f^{(p_2)}(a)\pmod n\)，那么取 \(t=f^{(p_1)}(a)-f^{(p_2)}(a)\)，于是计算 \(gcd(t,n)\) 即可得到一个因子。

我们可以搞两个指针，在数列上一个一次跳一步，一个一次跳两步，只需要环长的步数就可以找到满足上述条件的 \(t\)。代码中我们可以把每次的差值都乘起来对 \(n\) 取模，然后每隔 \(\log\) 步求一次 \(\gcd\)。

这里有些东西需要特判，比如如果进入了模 \(n\) 的环，那就退出再随一个 \(a\) 和 \(C\)。如果累乘的结果突然变成 \(0\)，就说明当前差值和之前的累乘结果中各至少有 \(n\) 的一个因子，直接算 \(\gcd\) 即可。

代码

inline ll mul(ll a, ll b, ll mod) { return ((ull)a * b - (ull)((ld)a / mod * b) * mod + mod) % mod; }
void get_fac(ll x) {
    static mt19937_64 rng((random_device){}());
    if(x < LIM) {
        while(x != 1) { fac.push_back(npri[x]), x /= npri[x]; }
        return;
    }
    if(is_prime(x)) return fac.push_back(x);
    for(int cnt = 1; cnt <= 20; cnt++) {
        ll a = rng() % (x - 1) + 1, c = rng() % x, b = (mul(a, a, x) + c) % x, sum = 1;
        for(int i = 1; i < 64 || (ll)i * i * i * i / 4 <= x; i++) {
            ll delta = (b - a + x) % x; if(delta == 0) break;
            sum = mul(sum, delta, x); if(sum == 0) { ll d = gcd(delta, x); if(d > 1) return get_fac(d), get_fac(x / d), void(); }
            if((i & 63) == 0) { ll d = gcd(sum, x); if(d > 1) return get_fac(d), get_fac(x / d), void(); }
            a = (mul(a, a, x) + c) % x, b = (mul(b, b, x) + c) % x, b = (mul(b, b, x) + c) % x;
        }
    }
    throw runtime_error("x is not a prime but we cannot get any factor\n");
}

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

按三次根号分治，小于 \(1000\) 的质数只有 \(168\) 个，直接暴力前缀和；剩下的数最多有两个质因子，直接 PR 分解出来跑莫队即可。

代码

#include<iostream>
#include<algorithm>
#include<random>
#define ll long long
using namespace std;
const int N = 1e5 + 10;
const int M = 1010;
const int MOD = 19260817;

struct hash_table {
    static const int MOD = 500041;
    struct Node {
        int key, v, next;
    } pool[3 * N];
    int ne, head[MOD];
    inline hash_table() { ne = 0; for(int i = 0; i < MOD; i++) head[i] = 0; }
    inline int &operator[](int key) {
        for(int i = head[key % MOD]; i; i = pool[i].next) {
            if(pool[i].key == key) return pool[i].v;
        }
        pool[++ne] = {key, 0, head[key % MOD]};
        head[key % MOD] = ne;
        return pool[ne].v;
    }
};

int blen = 330, blo[N];
struct Qr {
    int l, r, id;
    inline bool operator<(const Qr &b) const {
        if(blo[l] != blo[b.l]) return blo[l] < blo[b.l];
        return (blo[l] & 1) ? r < b.r : r > b.r;
    }
} q[N];

inline int qpow(int a, int b, int p) {
    int res = 1;
    while(b) {
        if(b & 1) res = (ll)res * a % p;
        a = (ll)a * a % p;
        b >>= 1;
    }
    return res;
}

int gcd(int a, int b) {
    if(b == 0) return a;
    return gcd(b, a % b);
}

bool ispri(int n) {
    if(n <= 3) return true;
    static const int aa[] = {2, 3, 5, 7, 11, 13, 17, 37};
    int t = 0, r = (n - 1);
    while(!(r & 1)) ++t, r >>= 1;
    for(int i = 0; i < 8; i++) {
        if(n == aa[i]) return true;
        int a = qpow(aa[i], r, n);
        if(a == 0) return false;
        if(a == 1) continue;
        for(int j = 1; j < t; j++) {
            if(a == n - 1) break;
            a = (ll)a * a % n;
        }
        if(a != n - 1) return false;
    }
    return true;
}

int get_factor(int n) {
    if(ispri(n)) return 1;
    static mt19937 rng((random_device){}());
    for(int cnt = 1; cnt <= 20; cnt++) {
        ll a0 = rng() % (n - 1) + 1, C = rng() % n;
        ll x = a0, y = (a0 * a0 + C) % n, s = (x - y + n) % n;
        if(n % a0 == 0) return a0;
        if(s == 0) continue;
        for(int cnt = 1; cnt <= 1e5; cnt++) {
            x = (x * x + C) % n; y = (y * y + C) % n; y = (y * y + C) % n;
            if((x - y + n) % n == 0) break;
            s = s * (x - y + n) % n;
            if(s == 0) return gcd(n, (x - y + n) % n);
            if(!(cnt & 63)) {
                int d = gcd(n, s);
                if(d != 1) return d;
            }
        }
    }
    return 1;
}

int n, m;
int a[N], a1[N]; ll ans1[N], ans2[N];
int s[180][N];

int npri[M], pri[180], pcnt;
int inv[MOD];

int cl, cr; ll c_ans = 1;
hash_table mp;

void add_v(int v, int w = 1) {
    c_ans = c_ans * inv[mp[v] + 1] % MOD;
    mp[v] += w;
    c_ans = c_ans * (mp[v] + 1) % MOD;
}

void add_p(int p) {
    if(a[p] == 1) return;
    if(a1[p] != 1) add_v(a1[p]), add_v(a[p] / a1[p]);
    else add_v(a[p]);
}

void del_p(int p) {
    if(a[p] == 1) return;
    if(a1[p] != 1) add_v(a1[p], -1), add_v(a[p] / a1[p], -1);
    else add_v(a[p], -1);
}

int main() {

    for(int i = 2; i < M; i++) {
        if(!npri[i]) pri[++pcnt] = i;
        for(int j = 1; j <= pcnt; j++) {
            if(i * pri[j] >= M) break;
            npri[i * pri[j]] = 1;
            if(i % pri[j] == 0) break;
        }
    }
    inv[1] = 1;
    for(int i = 2; i < MOD; i++) inv[i] = MOD - (ll)inv[MOD % i] * (MOD / i) % MOD;

    cin >> n >> m;
    for(int i = 1; i <= n; i++) cin >> a[i];
    for(int i = 1; i <= m; i++) cin >> q[i].l >> q[i].r, q[i].id = i;

    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= pcnt; j++) {
            while(a[i] % pri[j] == 0) {
                a[i] /= pri[j];
                ++s[j][i];
            }
        }
    }
    for(int i = 1; i <= pcnt; i++) {
        for(int j = 1; j <= n; j++) s[i][j] += s[i][j - 1];
    }
    for(int i = 1; i <= m; i++) {
        ans1[i] = 1;
        for(int j = 1; j <= pcnt; j++) ans1[i] = ans1[i] * (s[j][q[i].r] - s[j][q[i].l - 1] + 1) % MOD;
    }

    for(int i = 1; i <= n; i++) {
        a1[i] = get_factor(a[i]);
    }

    for(int i = 1; i <= n; i++) blo[i] = (i + blen - 1) / blen;
    sort(q + 1, q + 1 + m);
    cl = 1, cr = 0;
    for(int i = 1; i <= m; i++) {
        int l = q[i].l, r = q[i].r, id = q[i].id;
        while(l < cl) add_p(--cl);
        while(r > cr) add_p(++cr);
        while(l > cl) del_p(cl++);
        while(r < cr) del_p(cr--);
        ans2[id] = c_ans;
    }

    for(int i = 1; i <= m; i++) cout << ans1[i] * ans2[i] % MOD << '\n';

    return 0;
}