【模板】数论

1 最大公约数 (gcd)

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

2 快速幂 (递归版本)

ll PowerMod(ll a, ll n, ll m = mod) {
    if (!n || a == 1) return 1ll;
    ll s = PowerMod(a, n >> 1, m);
    s = s * s % m;
    return n & 1 ? s * a % m : s;
}

3 快速幂 (非递归版本) 

ll PowerMod(ll a, int n, ll c = 1) 
{
    for (; n; n >>= 1, a = a * a % mod)
    if (n & 1) 
    c = c * a % mod;
    return c;
}

4 扩展 Euclid 算法 (exgcd)

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

5 基本预处理 

// factorials, fact-inverses, unbounded
void init() {
    int i;
    for (*fact = i = 1; i < N; ++i) 
        fact[i] = (ll)fact[i - 1] * i % mod;
    --i;
        finv[i] = PowerMod(fact[i], mod - 2);
    for (; i; --i) 
        finv[i - 1] = (ll)finv[i] * i % mod;
}

// factorials, fact-inverses, bounded
void init(int n) {
    int i;
    for (*fact = i = 1; i <= n; ++i)
        fact[i] = (ll)fact[i - 1] * i % mod;
    finv[n] = PowerMod(fact[n], mod - 2);
    for (i = n; i; --i) 
        finv[i - 1] = (ll)finv[i] * i % mod;
}

// inverses, factorials, fact-inverses, unbounded
void init() {
    int i;
    for (inv[1] = 1, i = 2; i < N; ++i)
        inv[i] = ll(mod - mod / i) * inv[mod % i] % mod;
    for (*finv = *fact = i = 1; i < N; ++i)
        fact[i] = (ll)fact[i - 1] * i % mod, finv[i] = (ll)finv[i - 1] * inv[i] % mod;
}

// inverses, factorials, fact-inverses, bounded
void init(int n) {
    int i;
    for (inv[1] = 1, i = 2; i <= n; ++i) 
        inv[i] = ll(mod - mod / i) * inv[mod % i] % mod;
    for (*finv = *fact = i = 1; i <= n; ++i) 
        fact[i] = (ll)fact[i - 1] * i % mod, 
        finv[i] = (ll)finv[i - 1] * inv[i] % mod;
}

6 线性筛素数 (筛最小素因子)

int pn = 0, c[N], p[78540];

void sieve(int n) {
    int i, j, v;
    memset(c, -1, sizeof c);
    for (i = 2; i <= n; ++i) {
    if (!~c[i]) p[pn] = i, c[i] = pn++;
    for (j = 0; (v = i * p[j]) <= n && j <= c[i]; ++j) 
        c[v] = j;
    }
}

7 Miller-Rabin 素数测试

// No pseudoprime below 2^64 to base 2, 3, 5, 7, 11, 13, 82, 373.
const int test[8] = {2, 3, 5, 7, 11, 13, 82, 373};
bool Miller_Rabin(ll n) {
    if (n < MAX) return !c[n] && n > 1;
    int c, i, j; ll s = n - 1, u, v;
    for (c = 0; !(s & 1); ++c, s >>= 1);
    for (i = 0; i < 8; ++i) {
       if (!(n % test[i])) return false;
       u = PowerMod(test[i], s, n);
       for (j = 0; j < c; ++j, u = v)
          if (v = MulMod(u, u, n), u != 1 && u != n - 1 && v == 1) return false;
          if (u != 1) return false;
    }
    return true;
}

8 分解质因数

// version 1
inline void push(ll prime, int alpha)
 {pr[cnt] = prime, al[cnt++] = alpha;}

// version 2
inline void push(ll prime, int alpha) {
    map::iterator it = result.find(prime);
    it == result.end() ? (void)result.emplace(prime, alpha) : (void)(it->second += alpha);
}

void factor(ll n) {
    int i, j; cnt = 0; // result.clear();
    for (i = 0; n > 1; ) {
       if (n >= N) {
       for (; (ll)p[i] * p[i] <= n && n % p[i]; ++i);
        if (n % p[i]) return push(n, 1);
        } else i = c[n];
        for (j = 0; !(n % p[i]); n /= p[i], ++j); 
                push(p[i], j);
    }
}

9 线性筛 Euler φ 函数

void sieve(int n) {
    int i, j, v; phi[1] = 1;
    memset(c, -1, sizeof c);
    for (i = 2; i <= n; ++i) {
    if (!~c[i]) 
        p[pn] = i, c[i] = pn++, phi[i] = i - 1;
    for (j = 0; (v = i * p[j]) <= n && j < c[i]; ++j) 
        c[v] = j, phi[v] = phi[i] * (p[j] - 1);
    if (v <= n) 
        c[v] = j, phi[v] = phi[i] * p[j];
    }
}

10 线性筛 Möbius μ 函数 

void sieve(int n) {
    int i, j, v; mu[1] = 1;
    memset(c, -1, sizeof c);
    for (i = 2; i <= n; ++i) {
    if (!~c[i]) 
        p[pn] = i, c[i] = pn++, mu[i] = -1;
    for (j = 0; (v = i * p[j]) <= n && j < c[i]; ++j) 
        c[v] = j, mu[v] = -mu[i];
    if (v <= n) c[v] = j;
    }
}

11 遍历 ⌊n / i⌋ 的所有可能值 (旧版整除分块)

#define next(n, i) (n / (n / i + 1))
for(i = n; i >= 1; i = j){
    j = next(n, i);
    // do something
}

#define _next(n, i) (n / (n / (i + 1)))
for(i = 0; i < n; i = j){
    j = _next(n, i);
    // do something
}

// 0 1 100 1 2 50 2 3 33 3 4 25 4 5 20 5 6 16 16 6 7 14 7 8 12 8 9 11 9 10 10 10 11 9 11 12 8 12 13 7 14 15 6 16 17 5 20 21 4 25 26 3 33 34 2 50 51 1 100 101 0 ?

12 新版整除分块

// int32 ver.
int n, cnt, srn;
int v[G];

inline int ID(int x) {return x <= srn ? x : cnt + 1 - n / x;}

{
    int i;
    srn = sqrt(n), cnt = srn * 2 - (srn * (srn + 1) > n);
    for (i = 1; i <= srn; ++i)
        v[i] = i, v[cnt - i + 1] = n / i;
}

// int64 ver.
ll n, v[G];
int cnt, srn;

inline int ID(ll x) {return x <= srn ? x : cnt + 1 - n / x;}

{
    int i;
    srn = sqrt(n), cnt = srn * 2 - (srn * (srn + 1) > n);
    for (i = 1; i <= srn; ++i) 
        v[i] = i, v[cnt - i + 1] = n / i;
}

13 杜氏求和法 / 杜教筛

void Dsum() {
    int i, j, k, l, n, sq, _;
    for (i = 1; i <= cnt && v[i] < N; ++i) 
        Mu[i] = smu[v[i]], Phi[i] = sphi[v[i]];
    for (; i <= cnt; ++i) {
        int &m = Mu[i];
            ll &p = Phi[i];
        n = v[i], sq = sqrt(n), m = 1, p = n * (n + 1ll) / 2, l = n / (sq + 1);
        for (j = sq; j; --j) {
        _ = (k = n / j) - l, m -= _ * smu[j], p -= _ * sphi[j], l = k;
        if (j != 1 && k > sq) m -= Mu[_ = ID(k)], p -= Phi[_];
        }
    }
}

14大小步离散对数

int Boy_Step_Girl_Step(){
    mii :: iterator it;
    ll boy = 1, girl = 1, i, t;
    for(i = 0; i < B; i++){
        m.insert(pr((int)boy, i)); // if there exists in map, then ignore it
        (boy *= a) %= p;
    }
    for(i = 0; i * B < p; i++){
        t = b * inv(girl) % p;
        it = m.find((int)t);
        if(it != m.end())
            return i * B + it->second;
        (girl *= boy) %= p;
    }
    return -1;
} 

// input p and then B = (int)sqrt(p - 1);
// inv(x) is the multiplicative inverse of x in Z[p].
// mii: map <int, int> and m is a variable of which.

15 Pollard-Rho 素因数分解

ll Pollard_Rho(ll n, int c) {
    ll i = 1, k = 2, x = rand() % (n - 1) + 1, y = x, p;
    for (; k <= 131072; ) 
       {
    x = (MulMod(x, x, n) + c) % n;
    p = gcd<ll>((ull)(y - x + n) % n, n);
    if (p != 1 && p != n) return p;
    if (++i == k) y = x, k <<= 1;
    }
    return 0;
}

inline void push(ll prime, int alpha) {
    map::iterator it = result.find(prime);
    it == result.end() ? (void)result.insert(pr(prime, alpha)) : (void)(it->second += alpha);
}

void factor(ll n) {
    int i, j, k; ll p;
    for (i = 0; n != 1; )
    if (n >= MAX) {
    if (Miller_Rabin(n)) {push(n, 1); return;}
    for (k = 97; !(p = Pollard_Rho(n, k)); --k);
    factor(p); n /= p;
    }
        else {
    i = (c[n] ? c[n] : n);
    for (j = 0; !(n % i); n /= i, ++j);
    push(i, j);
    }
}