高级数论

· · 个人记录

Miller-Rabin

如何判定一个数 p 是否是质数？

我们考虑费马小定理：设 p 是任意素数，则对于所有 a 均有：

a^{p-1}\equiv1\pmod p

那么它的逆定理是否成立呢？答案是错误的，最小反例 p=561=3\times11\times17。

然后我们考虑另一个定理：二次探测定理，设 p 是任意素数，则方程

a^2\equiv 1\pmod p

的解为 a\equiv\pm1\pmod p。

我们用这个定理的逆否定理：设 p 是任意正整数，且方程

a^2\equiv 1\pmod p

有不满足 a\equiv\pm1\pmod p 的解，则 p 不是素数。

那么我们就可以设计出这样的算法：

首先把 p-1 分解为 2^kt 的形式，然后随机选出一个数 a，计算出 a^t，然后让其不断自乘，即：a\leftarrow a^2，那么重复 k 次后计算出的就是 a^{p-1}，同时结合二次探测定理判断，如果 a\not\equiv\pm1\pmod p，但 a^2\equiv1\pmod p，则 p 不是素数。

正确性：

若 p 通过了一次测试，则 p 不是质数的概率至多为 \dfrac14。

Pollard-Rho

Lucas & exLucas

求 \binom nm\bmod p 的值。

按照惯例，我们先考虑 p 为质数的情况。

如果 n,m<p，那么我们显然可以用逆元求。

然后我们又有一个定理：

\binom nm\equiv\binom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\times\binom{n\bmod p}{m\bmod p}\pmod p

证明不会。

那么我们考虑 p 是合数怎么办。

先把 p 质因数分解：p=\prod p_i^{a_i}，那么求出 \binom nm\bmod {p_i^{a_i}} 之后 CRT 合并即可。

考虑如何求出 \binom nm\pmod{p_i^{a_i}}，考虑阶乘式：

\binom nm\equiv\frac{n!}{m!(n-m)!}\pmod{p^a}

显然无法直接求，因为 m! 和 (n-m)! 不一定有逆元。

我们考虑求出 m! 和 (n-m)! 的因子中有多少个 p，设 n! 的因子中有 f(n) 个 p，则：

\binom nm\equiv\frac{\frac{n!}{p^{f(n)}}}{\frac{m!}{p^{f(m)}}\frac{(n-m)!}{p^{f(n-m)}}}p^{f(n)-f(m)-f(n-m)}\pmod{p^a}

那么现在问题就是求出：

\frac{n!}{p^{f(n)}}\bmod{p^a}

考虑把 n! 拆开：

\begin{aligned}n!&=1\times2\times\cdots\times n\\&=\left(p\times2p\times\cdots\times\left\lfloor\frac np\right\rfloor p\right)(1\times2\times\cdots\times(p-1)\times(p+1)\times\cdots)\\&=p^{\left\lfloor\frac np\right\rfloor}\left\lfloor\frac np\right\rfloor!\prod_{1\le i\le n,i\bmod p\not=0}i\end{aligned}

我们发现因为 \bmod\ p^a 了所以后面的很明显有循环节。

所以：

\begin{aligned}n!&=p^{\left\lfloor\frac np\right\rfloor}\left\lfloor\frac np\right\rfloor!\left(\prod_{1\le i\le p^k,i\bmod p\not=0}i\right)^{\left\lfloor\frac n{p^k}\right\rfloor}\left(\prod_{\lfloor\frac n{p^k}\rfloor\le i\le n,i\bmod p\not=0}i\right)\end{aligned}

然后我们发现，\left\lfloor\dfrac np\right\rfloor! 貌似需要继续递归，前面的 p^{\left\lfloor\frac n{p^k}\right\rfloor} 会被直接除掉，那么我们直接设一个函数 g(n)=\dfrac{n!}{p^{f(n)}}，则：

g(n)=g\left(\left\lfloor\frac np\right\rfloor\right)\left(\prod_{1\le i\le p^k,i\bmod p\not=0}i\right)^{\left\lfloor\frac n{p^k}\right\rfloor}\left(\prod_{\lfloor\frac n{p^k}\rfloor p^k\le i\le n,i\bmod p\not=0}i\right)

边界条件 g(0)=1，那么现在就有：

\binom nm\equiv\frac{g(n)}{g(m)g(n-m)}p^{f(n)-f(m)-f(n-m)}\pmod{p^k}

现在的问题就是怎么求 f(n)，还是像 g(n) 一样分组，我们就可以得到：

f(n)=f\left(\left\lfloor\dfrac np\right\rfloor\right)+\left\lfloor\dfrac np\right\rfloor

边界条件 f(n)=0(n<p)。

CRT & exCRT

求解以下同余方程组：

\begin{cases}x\equiv a_1\pmod{p_1}\\x\equiv a_2\pmod{p_2}\\\cdots\cdots\\x\equiv a_n\pmod{p_n}\end{cases}

首先考虑 p_i 都互质的情况，即 CRT 算法能处理的情况。

我只说最普遍的解法：设 M=\prod p_i,M_i=\dfrac M{p_i},t_i=M_i^{-1}\pmod{p_i}，那么我们就可以构造出 x=\sum a_iM_it_i，容易验证这是对的。

现在考虑 p_i 不互质的情况，考虑合并两个方程：

\begin{cases}x\equiv a_1\pmod{p_1}\\x\equiv a_2\pmod{p_2}\end{cases}

变形：

\begin{cases}x=a_1+up_1\\x=a_2+vp_2\end{cases} up_1-vp_2=a_2-a_1

然后用 exgcd 求出 u,v，然后带入 a_1+up_1 求出一个特解 x_1。

然后显然就有：

x\equiv x_1\pmod{\operatorname{lcm}(p_1,p_2)}

然后把 n 个方程都合并一下就行了。

BSGS & exBSGS

BSGS 用于求解以下方程：

a^x\equiv b\pmod p

其中 p 为质数。

考虑根号分治，设 t=\left\lfloor\sqrt p\right\rfloor，令 x=it-j，其中 0\le j<t，由此可得：

a^{it-j}\equiv b\pmod p (a^t)^i\equiv ba^j\pmod p

那么考虑求出所有 ba^j\bmod p 后，枚举 i，看 (a^t)^i 是否出现在了 ba^j 中，如果出现了，那么就得到了一个解，否则就没有解。

具体实现可以把 ba^j\bmod p 扔进 map 或者哈希表里。

原根与 N 次剩余

原根是个非常有用的东西。

首先我们定义阶：由欧拉定理可知如果 \gcd(a,m)=1，那么 a^{\varphi(m)}\equiv1\pmod{m}，因此满足 a^n\equiv1\pmod m 的最小正整数 n 存在，此时我们就称 n 是 a 模 m 的阶，记作 \delta_m(a)。

定义原根：我们称一个数 g 是 m 的原根当且仅当 \delta_m(g)=\varphi(m)。

一些定理：

  3. 若已知 m 的一个原根 g，那么对于所有 \gcd(i,\varphi(m))=1，都有 g^i 是原根。

  4. 原根存在定理：若 m 可以被表示为以下形式：

    ① m=2 或 m=4。

    ② 存在任意奇素数 p 和一正整数 a 使得 m=p^a 或 m=2p^a。

    则 m 有原根。

  5. 原根判定定理：若 g^{\varphi(m)}\equiv1\pmod m 且对于 \varphi(m) 的每个质因子 p，g 都满足 g^{\frac{\varphi(m)}p} \not\equiv1\pmod m，则 g 是 m 的一个原根。

现在我们要求 N 次剩余，即：找出 x^n\equiv a\pmod p 的所有解。

我们先讲一下 n=2，即二次剩余的情况，先变换一下记号：解 x^2\equiv n\pmod p，并称 n 是模 p 的二次剩余当且仅当上述方程有解。

1. p 为任意奇素数

先判断是否有解，由费马小定理可知 n^{p-1}\equiv1\pmod p，对其开根并由二次探测定理得 n^{\frac{p-1}2}\equiv\pm1\pmod p。

假设 n 是一个二次剩余，则存在 x 使得 n^{\frac{p-1}2}\equiv(x^2)^{\frac{p-1}2}\equiv x^{p-1}\equiv1\pmod p。

假设 n 满足 n^{\frac{p-1}2}\equiv1\pmod p，则令 n\equiv g^k\pmod p，其中 g 是 p 的原根，那么有 g^{k\frac{p-1}2}\equiv1\pmod p，由于 g 是原根，所以必然有 \varphi(p)|k\frac{p-1}2，又因为 p 是素数，所以 \varphi(p)=p-1，所以 k 显然是偶数，所以 n 开根之后得到的 g^{\frac k2} 一定存在。

综上所述，我们证明了 n^{\frac{p-1}2}\equiv1\pmod p 是 n 是一个二次剩余的充要条件（即，能相互推导），那么我们也就同时说明了 n 是一个非二次剩余等价于 n^{\frac{p-1}2}\equiv-1\pmod p。

ps：其实上面这个 n^{\frac{p-1}2}\bmod p 叫做勒让德符号，记作 \left(\dfrac np\right)。

还有个定理：模 p 下的二次剩余共有 \dfrac{p+1}2 个。

然后就是神仙构造了。

我们先随机找到一个 a 使得 a^2-n 不是一个二次剩余，然后把 \sqrt{a^2-n} 当成虚数单位 i，那么有 (a+i)^{p+1}\equiv1\pmod p。

证明：

引理 1

i^p\equiv-i\pmod p

证明：i^p\equiv i(i^2)^{\frac{p-1}2}\equiv i(a^2-n)^{\frac{p-1}2}=-i\pmod p。

引理 2

(a+b)^p\equiv a^p+b^p\pmod p

考虑二项式定理展开，因为有个阶乘，所以除了 \dbinom0p,\dbinom pp 之外的所有系数 \bmod\ p 之后都是 0，所以原式等于 a^p+b^p。

现在证明原定理：

(a+i)^{p+1}\equiv(a+i)(a+i)^p\equiv(a+i)(a^p+i^p)\equiv(a+i)(a-i)\equiv a^2-i^2\equiv n\pmod p

那么 n 的二次剩余即为 (a+i)^{\frac{p+1}2}\bmod p，其实部为一个解，实部的相反数为另一个解。

算一个简单的期望可以知道期望 2 次找到 a，所以时间复杂度 O(\log p)。

然后考虑 p 不是奇素数的情况：

2.p=p_0^a，其中 p_0 为任意奇素数

设 n=p_0^cm，且 m\bmod p_0\not=0。

若 c\ge k，则显然答案为 0。

如果 c<k，则有解当且仅当 c\bmod2=0。

代码

// Miller-Rabin
#include <cstdio>
using namespace std;
typedef __int128 ll;
ll pow(ll a, ll b, ll p) {
    ll ans = 1;
    while (b != 0) {
        if ((b & 1) == 1) {
            ans = ans * a % p;
        }
        a = a * a % p;
        b = b / 2;
    }
    return ans;
}
const int pcnt = 12;
const ll prime[pcnt] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
ll check(ll a, ll p) {
    if (p % a == 0) {
        return false;
    }
    if (pow(a, p - 1, p) != 1) {
        return false;
    }
    ll d = p - 1;
    while ((d & 1) == 0) {
        d /= 2;
        ll temp = pow(a, d, p);
        if (temp == p - 1) {
            return true;
        } else if (temp != 1) {
            return false;
        }
    }
    return true;
}
bool miller_rabin(ll n) {
    if (n <= 40) {
        for (int i = 0; i < pcnt; i++) {
            if (n == prime[i]) {
                return true;
            }
        }
        return false;
    } else {
        for (int i = 0; i < pcnt; i++) {
            if (n == prime[i]) {
                return true;
            }
            if (n % prime[i] == 0) {
                return false;
            }
            if (!check(prime[i], n)) {
                return false;
            }
        }
        return true;
    }
}
int main() {
    long long n;
    while (scanf("%lld", &n) != EOF) {
        if (miller_rabin(n)) {
            printf("Y\n");
        } else {
            printf("N\n");
        }
    }
    return 0;
}
// Lucas
#include <cstdio>
using namespace std;
typedef long long ll;
const int P = 100005;
ll fac[P], invfac[P];
ll pow(ll a, ll b, ll p) {
    ll ans = 1;
    while (b != 0) {
        if (b & 1 == 1) {
            ans = (ans * a) % p;
        }
        a = a * a % p;
        b = b / 2;
    }
    return ans;
}
ll C(ll n, ll m, ll p) {
    if (m > n) {
        return 0;
    } else {
        return (fac[n] * invfac[m]) % p * invfac[n - m] % p;
    }
}
ll Lucas(ll n, ll m, ll p) {
    if (m == 0) {
        return 1;
    } else {
        return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
    }
}
void inline solve() {
    ll n, m, p;
    scanf("%lld %lld %lld", &n, &m, &p);
    fac[0] = 1;
    for (int i = 1; i < p; i++) {
        fac[i] = (fac[i - 1] * i) % p;
    }
    invfac[p - 1] = pow(fac[p - 1], p - 2, p);
    for (int i = p - 2; i >= 0; i--) {
        invfac[i] = invfac[i + 1] * (i + 1) % p;
    }
    printf("%lld\n", Lucas(n + m, n, p));
}
int main() {
    int t;
    scanf("%d", &t);
    for (int i = 1; i <= t; i++) {
        solve();
    }
    return 0;
}
// BSGS
#include <cstdio>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
ll BSGS(ll a, ll b, ll p) {
    ll t = sqrt(p);
    map < ll , ll > hash;
    for (ll bpowerai = b, i = 0; i < t; i++, bpowerai = bpowerai * a % p) {
        hash[bpowerai] = i;
    }
    ll powerat = 1;
    for (ll i = 0; i < t; i++) {
        powerat = powerat * a % p;
    }
    for (ll powerati = 1, i = 0; i <= t; i++) {
        ll j = hash.find(powerati) == hash.end() ? -1 : hash[powerati];
        if (j != -1) {
            return i * t - j;
        }
        powerati = powerati * powerat % p;
    }
    return -1;
}
int main() {
    ll p, a, b;
    scanf("%lld %lld %lld", &p, &a, &b);
    ll ans = BSGS(a, b, p);
    if (ans == -1) {
        printf("no solution");
    } else {
        printf("%lld", ans);
    }
    return 0;
}
// CRT
#include <cstdio>
using namespace std;
typedef long long ll;
const int N = 100005;
int n;
ll a[N], b[N];
ll ans;
void exgcd(ll a, ll b, ll & x0, ll & y0) {
    if (b == 0) {
        x0 = 1;
        y0 = 0;
    } else {
        exgcd(b, a % b, y0, x0);
        y0 -= a / b * x0;
    }
}
ll inv(ll a, ll b) {
    ll x0, y0;
    exgcd(a, b, x0, y0);
    x0 = (x0 + b) % b;
    return x0;
}
void CRT() {
    ll M = 1;
    for (int i = 1; i <= n; i++) {
        M *= b[i];
    }
    for (int i = 1; i <= n; i++) {
        ll Mi = M / b[i];
        ll ti = inv(Mi, b[i]);
        ll xi = Mi * ti;
        ans = (ans + a[i] * xi) % M;
    }
}
signed main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        scanf("%lld %lld", b + i, a + i); // 题目里是 x = a (mod b)
        a[i] %= b[i];
    }
    CRT();
    printf("%lld", ans);
    return 0;
}
// exCRT
#include <cstdio>
using namespace std;
typedef __int128 ll;
const int N = 100005;
ll mult(ll n, ll k, ll mod){
    ll ans = 0;
    while (k != 0) {
        if ((k & 1) == 1) {
            ans = (ans + n) % mod;
        }
        k >>= 1;
        n = (n + n) % mod;
    }
    return ans;
}
ll inline fix(ll & a, ll mod) {
    return a = ((a % mod) + mod) % mod;
}
ll gcd(ll a, ll b) {
    if (b == 0) return a;
    else return gcd(b, a % b);
}
ll lcm(ll a, ll b) {
    return a / gcd(a, b) * b;
}
// 解方程 ax + by = gcd(a, b) 
void exgcd(ll a, ll b, ll & x0, ll & y0) {
    if (b == 0) {
        x0 = 1;
        y0 = 0;
    } else {
        exgcd(b, a % b, y0, x0);
        y0 -= a / b * x0;
    }
}
void merge(ll & a0, ll & b0, ll a1, ll b1) {
    ll k0, k1;
    exgcd(b0, b1, k0, k1);
    k0 *= a1 - a0;
    k0 /= gcd(b0, b1); // 只需要用一个 
    a0 = a0 + b0 * k0; // 一个特解 
    b0 = lcm(b0, b1);
    fix(a0, b0);
}
int n;
ll a[N], b[N];
ll ansa, ansb;
ll ans;
void excrt() {
    // 解方程组 
    // x = a1 (mod b1) 
    // x = a2 (mod b2) 
    // ...... 
    // x = an (mod bn) 
    ansa = a[0];
    ansb = b[0];
    for (int i = 1; i < n; i++) {
        merge(ansa, ansb, a[i], b[i]);
    }
    ans = ansa;
}
signed main() {
    scanf("%d", &n);
    for (int i = 0; i < n; i++) {
        scanf("%lld %lld", b + i, a + i); // 题目里是 x = a (mod b) 
        a[i] %= b[i];
    }
    excrt();
    printf("%lld", ans);
    return 0;
}
// 原根
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const ll N = 1000005;
bool isnp[N], have_g[N];
ll prime[N], pcnt;
ll phi[N];
void init() {
    isnp[1] = false;
    phi[1] = 1;
    for (ll i = 2; i < N; i++) {
        if (!isnp[i]) {
            prime[++pcnt] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= pcnt && i * prime[j] < N; j++) {
            isnp[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
    have_g[2] = have_g[4] = true;
    for (ll i = 2; i <= pcnt; i++) {
        for (ll j = 1; j * prime[i] < N; j *= prime[i]) {
            have_g[j * prime[i]] = true;
        }
        for (ll j = 2; j * prime[i] < N; j *= prime[i]) {
            have_g[j * prime[i]] = true;
        }
    }
}
ll factor[N], fcnt;
void get_factor(int n) {
    fcnt = 0;
    for (ll i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            factor[++fcnt] = i;
            while (n % i == 0) {
                n /= i;
            }
        }
    }
    if (n > 1) {
        factor[++fcnt] = n;
    }
}
ll pow(ll a, ll b, ll p) {
    ll ans = 1;
    while (b != 0) {
        if ((b & 1) == 1) {
            ans = ans * a % p;
        }
        a = a * a % p;
        b = b / 2;
    }
    return ans;
}
bool check(ll g, ll n) {
    if (pow(g, phi[n], n) != 1) {
        return false;
    }
    for (ll i = 1; i <= fcnt; i++) {
        if (pow(g, phi[n] / factor[i], n) == 1) {
            return false;
        }
    }
    return true;
}
ll get_min_g(ll n) {
    get_factor(phi[n]);
    for (ll i = 1; i < n; i++) {
        if (check(i, n)) {
            return i;
        }
    }
    return 0;
}
ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a % b);
}
ll all_g[N], gcnt;
void get_all_g(ll n) {
    gcnt = 1;
    ll prod = get_min_g(n);
    all_g[gcnt] = prod;
    for (int i = 2; i < phi[n]; i++) {
        prod = prod * all_g[1] % n;
        if (gcd(i, phi[n]) == 1) {
            all_g[++gcnt] = prod;
        }
    }
    sort(all_g + 1, all_g + gcnt + 1);
}
int main() {
    init();
    int t;
    scanf("%d", &t);
    for (ll n, d; t != 0; t--) {
        scanf("%lld %lld", &n, &d);
        if (have_g[n]) {
            get_all_g(n);
            printf("%lld\n", gcnt);
            for (int i = 1; i <= gcnt / d; i++) {
                printf("%lld ", all_g[i * d]);
            }
            printf("\n");
        } else {
            printf("0\n\n");
        }
    }
    return 0;
}
// Cipolla 算法
#include <cstdio>
#include <cstdlib>
using namespace std;
typedef long long ll;
ll n, p;
ll w;
struct Complex {
    // a + b sqrt(w) 其中 w 为模 p 意义下的非二次剩余 
    ll a, b;
    Complex() {
        a = b = 0;
    }
};
Complex operator * (const Complex & a, const Complex & b) {
    Complex c;
    c.a = (a.a * b.a % p + w * a.b % p * b.b % p) % p;
    c.b = (a.b * b.a % p + a.a * b.b % p) % p;
    return c;
}
Complex pow(Complex a, ll b) {
    Complex ans;
    ans.a = 1;
    while (b != 0) {
        if ((b & 1) == 1) {
            ans = ans * a;
        }
        b /= 2;
        a = a * a;
    }
    return ans;
}
ll pow(ll a, ll b, ll p) {
    int ans = 1;
    while (b != 0) {
        if ((b & 1) == 1) {
            ans = ans * a % p;
        }
        b /= 2;
        a = a * a % p;
    }
    return ans;
}
ll Legendre(ll n, ll p) {
    return pow(n, (p - 1) / 2, p);
}
ll min(ll a, ll b) {
    return a < b ? a : b;
}
ll max(ll a, ll b) {
    return a > b ? a : b;
}
void inline solve() {
    scanf("%lld %lld", &n, &p);
    if (n == 0) {
        printf("0\n");
    } else if (Legendre(n, p) == p - 1) {
        printf("Hola!\n");
    } else {
        ll a;
        while (true) {
            a = rand() % p;
            w = (a * a % p - n + p) % p;
            if (Legendre(w, p) == p - 1) {
                break;
            }
        }
        Complex res;
        res.a = a;
        res.b = 1;
        res = pow(res, (p + 1) / 2);
        ll ans = res.a;
        printf("%lld %lld\n", min(ans, p - ans), max(ans, p - ans));
    }
}
int main() {
    int t;
    scanf("%d", &t);
    while (t--) {
        solve();
    }
    return 0;
}