[模版] 数论基础模版

欧几里得算法

• gcd，(O(log n))

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}

• lcm，(O(log n))

// 注意数据范围，过大的数要用__int128
ll lcm(ll a, ll b) {return a * b / gcd(a, b);}

• 扩展欧几里得，(log n)

ll exgcd(ll a, ll b, ll &x, ll &y) { // ax + by = gcd(a, b)
    if(b == 0) {
        y = 0, x = 1;
        return a;
    }  
    ll res = exgcd(b, a%b, x, y);
    ll tx = x,ty = y;
    x = ty;
    y = tx - (a / b) * ty;
    return res;
}
bool solve(ll a, ll b, ll c, ll &x, ll &y) { // ax + by = c, gcd(a, b)|c, x>=0
    if(!c) {
        x = y = 0;
        return true;
    }
    ll d = exgcd(a, b, x, y);
    if(c % d) return false;
    a /= d, b /= d, c /= d;
    ll ret = c % b;
    x = x * ret;
    x = x + (b - x) / b * b; // 调整此处控制x下界，x可以加上容易倍数的b
    y = (c - a * x) / b;
    return true;
}

欧拉定理&费马小定理

[a^bequiv egin{cases} a^{b mod varphi(p)},\,&gcd(a,\,p)=1\ a^b,&gcd(a,\,p) e1,\,blt varphi(p)\ a^{b mod varphi(p)+varphi(p)},&gcd(a,\,p) e1,\,bgevarphi(p) end{cases} pmod p ]

注意模数和底数不互质时的处理。

ll getphi(ll n) {
    ll p = n;
    for(int i = 0; 1ll * pri[i] * pri[i] < p; i++) {
        if(n % pri[i] == 0) {
            p = p / pri[i] * (pri[i] - 1);
            while(n % pri[i] == 0) n /= pri[i];
        } 
    }
    if(n > 1) p = p / n * (n - 1);
    return p;
}

乘法逆元

• 线性求逆元

// i无逆元时，inv[i]是未定义的
inv[1] = 1;
for(int i = 2; i <= n; i++) {
    inv[i] = (ll)(p - p / i) * inv[p % i] % p;
}

• 扩展欧几里得

ll inv(ll a, ll p) {
    ll x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p; 
}

中国剩余定理(CRT)

// a, r下标均从1开始
// crt
ll crt(int k, ll a[], ll r[]) { // xi mod ai = ri, ai模数必须互质
    ll n = 1, ans = 0;
    for(int i = 1; i <= k; i++) n = n * a[i];
    for(int i = 1; i <= k; i++) {
        ll m = n / a[i], b, y;
        exgcd(m, a[i], b, y); // b * m mod a[i] = 1
        ans = (ans + r[i] * m * b % n) % n;
    }
    return (ans % n + n) % n;
}
// 任意模数crt
ll excrt(int n, ll a[], ll r[]) { //xi mod ai = ri, 无解时返回-1
    ll res, m, x, y;
    res = r[1];
    m = a[1];
    for(int i = 2; i <= n; i++) {
        ll d = ((r[i] - res) % a[i] + a[i]) % a[i];
        ll g = exgcd(m, a[i], x, y);
        if(d % g) return -1;
        ll k = d / g;
        ll rm = lcm(a[i], m);
        res = ((__int128)m * x % rm * k + res) % rm;
        m = rm;
    }
    return (res % m + m) % m;
}

二次剩余

• cipolla，要求模数必须为奇素数，期望时间复杂度(O(log n))

namespace cipolla {
    ll w;
    struct num {  //建立一个复数域
      ll x, y;
    };
    num mul(num a, num b, ll p) {  //复数乘法
      num ans = {0, 0};
      ans.x = ((a.x * b.x % p + a.y * b.y % p * w % p) % p + p) % p;
      ans.y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
      return ans;
    }
    ll qpow_img(num a, ll b, ll p) {  //虚部快速幂
      num ans = {1, 0};
      while (b) {
        if (b & 1) ans = mul(ans, a, p);
        a = mul(a, a, p);
        b >>= 1;
      }
      return ans.x % p;
    }
    ll solve(ll n, ll p) { // x^2 mod p = n
          n %= p;
          if (p == 2) return n;
          if (qpow(n, (p - 1) / 2, p) == p - 1) return -1;
          ll a;
          while (1) {  //生成随机数再检验找到满足非二次剩余的a
                a = rand() % p;
                w = ((a * a % p - n) % p + p) % p;
                if (qpow(w, (p - 1) / 2, p) == p - 1) break;
          }
          num x = {a, 1};
          return qpow_img(x, (p + 1) / 2, p);
    }
}

• BSGS，求任意次开根，即求(x^aequiv b mod p)，(p)为奇素数。时间复杂度(O(sqrt{n}))

ll getroot(ll p) { // p是素数，求p的原根
    ll phi = p - 1, tp = phi;
    vector<ll> factors;
    for(int i = 2; 1ll * i * i <= phi; i++) {
        if(tp == 1) break;
        if(tp % i == 0) {
            while(tp % i == 0) tp /= i;
            factors.push_back(i);
        }
    }
    ll g = 2;
    while(g < p) {
        bool ok = true;
        for(ll d : factors) {
            if(qpow(g, phi / d, p) == 1) {
                ok = false;
                break;
            }
        }
        if(ok) return g;
        g++;
    } 
    return -1;
}

vector<ll> ans;
void solve(ll a, ll b, ll p) { // 获得x^a % p = b所有解，保存在数组ans中
    if(b == 0) {
        ans.push_back(0);
        return ;
    }
    ll g = getroot(p), phi = p - 1;
    ll d = phi / gcd(a, phi);
    ll c = bsgs(qpow(g, a, p), b, p);
    if(c == -1) return ;
    for(ll i = c % d; i < phi; i += d) {
        ans.push_back(qpow(g, i, p));
    }
    sort(ans.begin(), ans.end());
}

离散对数

• exbsgs，求解(a^xequiv bmod p)，其中(a)不必和(p)互质，时间复杂度(O(sqrt{n}))。

ll inv(ll a, ll p) {
    ll x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p; 
}

map<ll, int> ex;
ll bsgs(ll a, ll b, ll p) { // a^x % p = b; a, p必须互质
    ex.clear();
    int sqp = ceil(sqrt(p));
    ll pab = b % p, pw = qpow(a, sqp, p), pa = pw;
    for(int i = 0; i <= sqp; i++) {
        ex[pab] = i;
        pab = pab * a % p;
    }
    for(int i = 1; i <= sqp; i++) {
        if(ex.count(pa)) {
            ll res = 1ll * i * sqp - ex[pa];
            return res;
        }
        pa = pa * pw % p;
    }
    return -1;
}

ll exbsgs(ll a, ll b, ll p) {
    b %= p;
    a %= p;
    ll d = 1, tp = p, tb = b, pwa = 1;
    int k = 0;
    while(tb != 1 && (d = gcd(a, tp)) != 1) {
        if(b % d != 0) return -1;
        if(tb % d != 0) break;
        tp /= d;
        tb /= d;
        k++; 
    }
    for(int i = 0; i <= k; i++) {
        if(pwa == b) return i;
        pwa = pwa * a % p;
    }
    if(gcd(a, tp) != 1) return -1;
    ll res = bsgs(a % tp, b * inv(qpow(a, k, tp), tp) % tp, tp);
    if(res == -1) return -1;
    return res + k;
}

卢卡斯定理

• 递推式求组合数

for(int i = 1; i <= n; i++) C[i][0] = C[i][i] = 1;
for(int i = 1; i <= n; i++) {
    for(int j = 1; j < i; j++) {
        C[i][j] = (C[i-1][j-1] + C[i-1][j]) % M; // M = mod
    }
}

• 公式法求组合数，注意若模数为(p)。当(n< p)可以用公式法，(nge p)时请使用其他方法，否则可能会出现模数为0导致结果出错。

[C(n,m)=frac{n!}{(n-m)!m!} ]

• 卢卡斯定理：求(C(n, m)mod p)​，(p)​为素数，时间复杂度(O(f(p) + g(p)log n))​。其中(f(p))​代表预处理值域为前(p)​组合数的复杂度，(g(p))​为计算单个组合数的复杂度。若直接暴力计算，时间复杂度(O(plog p))​，

ll lucas(ll n, ll m, ll p) {
    if(m == 0) return 1;
    return C(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}

• 扩展卢卡斯定理：求(C(n, m)mod p)，(p)为任意模数，时间复杂度(O(plog p))。预处理+p固定+若干询问：(O(p+Tlog p))，(T)为询问次数。

ll calc(ll n, ll x, ll P) {
    if(!n) return 1;
    ll s = 1;
    for(ll i = 1; i <= P; i++) // 此处可以预处理
        if(i % x) s = s * i % P;
    s = qpow(s, n / P, P);
    for(ll i = 1; i <= n % P; i++) // 此处可以预处理
        if(i % x) s = i * s % P;
    return s * calc(n / x, x, P) % P;
}

ll multilucas(ll n, ll m, ll x, ll P) {
    int cnt = 0;
    for(ll i = n; i; i /= x) cnt += i / x;
    for(ll i = m; i; i /= x) cnt -= i / x;
    for(ll i = n - m; i; i /= x) cnt -= i / x;
    // inv使用exgcd的
    return qpow(x, cnt, P) % P * calc(n, x, P) % P * inv(calc(m, x, P), P) % P * inv(calc(n - m, x, P), P) % P;
}

ll exlucas(ll n, ll m, ll P) { // n >= m, C(n, m)
    if(m > n) return 0; 
    int cnt = 0;
    ll p[20], a[20];
    for(ll i = 2; i * i <= P; i++) {
        if(P % i == 0) {
            p[++cnt] = 1;
            while (P % i == 0) p[cnt] = p[cnt] * i, P /= i;
            a[cnt] = multilucas(n, m, i, p[cnt]);
        }
    }
    if(P > 1) p[++cnt] = P, a[cnt] = multilucas(n, m, P, P);
    return excrt(cnt, p, a);
}