基础数论模板

模板

基础的数学模板应该都有了。

namespace easymath
{
typedef long long ll;
const ll mod = 998244353;
const int maxsz = 1e3 + 9;
// inv
ll inv[maxsz];
// sieve
vector<int> prime;
bool vis[maxsz];
int phi[maxsz], pcnt[maxsz], minp[maxsz];
// pcnt : cnt prime factor
// P.S. sum_of_factors[i] = (1LL << pcnt[i]) 
// minp : min prime factor
ll sf[maxsz]; // sum of all factors in [1, n]
// fact
ll fact[maxsz], ifact[maxsz];
ll lcm (ll a, ll b)
{
    return a / __gcd (a, b) * b;
}
ll ksm (ll b, ll k)
{
    ll ret = 1;
    while (k) {
        if (k & 1) ret = ret * b % mod;
        b = b * b % mod;
        k /= 2;
    }
    return ret;
}
void exgcd (ll a, ll b, ll& d, ll& x, ll& y)
{
    if (b == 0) d = a, x = 1, y = 0;
    else {
        exgcd (b, a % b, d, y, x);
        y -= x * (a / b);
    }
}
// inv
ll pinv (ll x) // mod is prime
{
    return ksm (x, mod - 2);
}
ll rinv (ll a, ll m) // relatively prime
{
    ll x, y, d;
    exgcd (a, m, d, x, y);
    if (d == 1) return (x % m + m) % m;
    else return -1;
}
void init_inv (int n, int m) // n : upper bound, m : mod
{
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) {
        inv[i] = (m - m / i) * inv[m % i] % m;
    }
}
// sieve
void prime_sieve (int n)
{
    for (int i = 2; i <= n; ++i) {
        if (!vis[i]) {
            prime.push_back (i);
            pcnt[i] = 1;
            minp[i] = i;
        }
        for (int j = 0; j < prime.size ( ); ++j) {
            if (i * prime[j] >= n) break;
            vis[i * prime[j]] = 1;
            minp[i * prime[j]] = prime[j];
            if (i % prime[j] == 0) {
                pcnt[i * prime[j]] = pcnt[i];
                break;
            }
            pcnt[i * prime[j]] = pcnt[i] + 1;
        }
        // s[i] = (1LL << pcnt[i]);
        // binomial theorem
    }
}
void factor_sum (int n)
{
    // don't use minp[] after this function
    sf[1] = 1, minp[1] = 0;
    for (int i = 2; i <= n; ++i) {
        if (!vis[i]) {
            prime.push_back (i);
            sf[i] = i + 1;
            minp[i] = i;
        }
        for (int j = 0; j < prime.size ( ); ++j) {
            if (i * prime[j] >= n) break;
            vis[i * prime[j]] = 1;
            minp[i * prime[j]] = prime[j];
            if (i % prime[j] == 0) {
                sf[i * prime[j]] = sf[i] * (minp[i] * prime[j] * prime[j] - 1) / (minp[i] * prime[j] - 1);
                minp[i * prime[j]] = minp[i] * prime[j];
                break;
            }
            sf[i * prime[j]] = sf[i] * sf[prime[j]];
            minp[i * prime[j]] = prime[j];
        }
    }
}
ll get_phi (ll n) // euler's totient O(sqrt n)
{
    ll ret = n;
    for (ll i = 2; i * i <= n; ++i)
        if (n % i == 0) {
            ret -= ret / i;
            while (n % i == 0) n /= i;
        }
    if (n > 1) ret -= ret / n;
    return ret;
}
void init_phi (int n)
{
    memset (vis, 1, sizeof (vis));
    int cnt = 0;
    vis[1] = 0;
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (vis[i]) {
            prime.push_back (i);
            phi[i] = i - 1;
        }
        for (int j = 0; j < prime.size ( ); j++) {
            if (i * prime[j] > n) break;
            vis[i * prime[j]] = 0;
            if (i % prime[j])
                phi[i * prime[j]] = phi[i] * phi[prime[j]];
            else {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
        }
    }
}
// fact
void init_fact (int n)
{
    fact[0] = 1;
    for (int i = 1; i <= n; ++i)
        fact[i] = fact[i - 1] * i % mod;
    ifact[n] = pinv (fact[n]);
    for (int i = n - 1; i >= 0; --i)
        ifact[i] = ifact[i + 1] * (i + 1) % mod;
    // 1 / (n!) * n = 1 / (n - 1)!
}
ll C (int n, int m)
{
    return fact[n] * ifact[n - m] % mod * ifact[m] % mod;
}
ll crt (ll r[ ], ll m[ ], int n)
{
    ll p = 1, ret = 0;
    for (int i = 1; i <= n; ++i) p *= r[i];
    for (int i = 1; i <= n; ++i) {
        ll a = p / r[i];
        ll aa = rinv (a, m[i]);
        ret += r[i] * a * aa;
    }
    return ret;
}
}

GCD

欧几里得算法（辗转相除）：

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

最坏复杂度 \(O(n\log n)\)，在求斐波那契数列相邻项的 gcd 时取到。

算术基本定理：因数唯一分解。

EXGCD：求 \(ax+by=\gcd(a,b)\) 的一组解。

根据裴蜀定理，只要 \(a,b\) 不全为零，该方程就存在整数解。

下面来对式子 \(ax+by=\gcd(a,b)\) 做一点变形：\(ax+by=\gcd(a,b)=\gcd(b,a\mod b)\)，有 \(bx+(a\mod b)y=bx+(a-\lfloor \frac a b \rfloor b)y=\gcd(b,a\mod b)\)，即 \(ay+b(x-\lfloor \frac a b \rfloor y)=\gcd(b,a\mod b)\)，这意味着可以递归求解 \((x,y)\)。这个递归过程和欧几里得算法的递归是一样的，在递归边界 \(b=0\) 处，可以发现 \(\gcd(a,0)=a\)，所以 \(x=1,y=0\)。

因此 exgcd 的代码如下：

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

显然，其复杂度并没有变化。

线性筛

其实 \(O(n\log\log n)\) 的埃氏筛如果做好优化也不算很慢，但还是没什么用，所以这里就不写了。核心思想是对每个质数找出所有倍数。

埃氏筛的缺点是重复标记了一些合数，欧拉筛优化了这一点，因此复杂度是线性的。

void init (int n)
{
    for (int i = 2; i <= n; ++i) {
        if (!vis[i]) prime.push_back (i);
        for (int j = 0; j < prime.size ( ); ++j) {
            if (i * prime[j] >= n) break;
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }
}

这是欧拉筛的核心部分，注意第二个 break 的意义，这保证了每个合数都被其最小的质因数筛到。

最小质因数、质因数个数、因数和、因数个数等都可以用线性筛求出。

欧拉函数

欧拉函数 \(\varphi(n)\) 表示在 \([1,n]\) 中与 \(n\) 互质的数的个数。

显然：

• \(\varphi(n)\ge\varphi(1)=1\)
• \(\varphi(p)=p-1\)，\(p\) 为质数

性质：

• 欧拉函数是积性函数，即 \(\gcd(a,b)=1\Rightarrow\varphi(a\times b)=\varphi(a)\times\varphi(b)\)
• \(\large n=\sum_{d|n}\varphi(d)\)（莫比乌斯反演）
• 若 \(n=p^k\)（\(p\) 为质数），则 \(\varphi(n)=p^k -p^{k-1}\)（定义）
• 设 \(\large n=\Pi_{i=1}^m p_i^{k_i}\)（\(p_i\) 为质数），\(\large\varphi(n)=n\Pi_{i=1}^m (1-\dfrac{1}{p_i})\)（算术基本定理）

\(O(\sqrt n)\) 单点求值：

ll get_phi (ll n)
{
    ll ret = n;
    for (ll i = 2; i * i <= n; ++i)
        if (n % i == 0) {
            ret -= ret / i;
            while (n % i == 0) n /= i;
        }
    if (n > 1) ret -= ret / n;
    return ret;
}

线性筛：

void pre ( )
{
    memset (is_prime, 1, sizeof (is_prime));
    int cnt = 0;
    is_prime[1] = 0;
    phi[1] = 1;
    for (int i = 2; i <= 5000000; i++) {
        if (is_prime[i]) {
            prime[++cnt] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= 5000000; j++) {
            is_prime[i * prime[j]] = 0;
            if (i % prime[j])
                phi[i * prime[j]] = phi[i] * phi[prime[j]];
            else {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
        }
    }
}

欧拉定理：

若 \(\large \gcd(a,m)=1\) 则 \(\large a^{\varphi(m)}\equiv 1\pmod m\)。这就是欧拉降幂的原理。

扩展欧拉定理：

\[\large a^b\equiv \begin{cases}a^{b\mod \varphi(p)}\, ,\,\gcd(a,p)=1\\ a^b\, ,\, \gcd(a,p)\ne 1\, , \, b\lt \varphi(p)\\ a^{b\mod {2\times \varphi(p)}}\, ,\, \gcd(a,p)\ne 1\, , \, b\ge \varphi(p)\end{cases}\pmod p \]

费马小定理

若 \(p\) 为质数，\(\gcd(a,p)=1\)，则 \(a^{p-1}\equiv 1\pmod p\)。

乘法逆元

单点求值：exgcd 或快速幂（费马小定理）

注意：快速幂求逆元的依据是费马小定理，因此需要保证模数为质数，而 exgcd 求逆元只需保证互质。

中国剩余定理

\(x\) 满足对模数 \(m_i\) 的模为 \(r_i\)，求最小的 \(x\)。

1. 令 \(\large p=\Pi_{i=1}^n r_i\)
2. 令 \(\large a_i=\dfrac{p}{r_i}\)
3. 令 \(\large a_i' = \text{inv}(a_i,m_i)\)

最后得到

\[\large x = \sum_{i=1}^n r_i a_i a_i' \]