模板——数学类

发表于 分类于 ， 阅读次数： 作者 Pomelorin

模板——数学类

快速幂，质数线性筛，质因数分解，欧几里得（gcd），扩展欧几里得（exgcd），求欧拉函数 \(\varphi(x)\)，欧拉定理，扩展欧拉定理，中国剩余定理（CRT），在任意模意义下的逆元，同余最短路，组合数，Lucas 定理，卡特兰数，递推求逆元，拉格朗日插值，康托展开，高斯消元，矩阵快速幂。

快速幂

1
2
3
4
5
6
7
8
9
10
ll qpow(ll a, ll b, ll p) {
    ll ret = 1;
    while (b) {
        if (b % 2)
            ret = ret * a % p;
        a = a * a % p;
        b /= 2;
    }
    return ret;
}

质数线性筛

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
vector<ll> primes;
ll notprime[M];
void eular() {
    notprime[0] = notprime[1] = 1;
    for (ll i = 2; i <= M - 10; i++) {
        if (notprime[i] == 0) {
            primes.push_back(i);
        }
        for (auto j: primes) {
            if (i * j > M - 10)
                break;
            notprime[i * j] = 1;//i*j 的最小质因数是 j
            if (i % j == 0) {
                break;
            }
        }
    }
}

质因数分解

此方法 \(O(\sqrt{n})\)。

值域不大时可用最小质因数来质因数分解，\(O(\log n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
vector<ll> primes;
vector<ll> fact(ll x) {
    vector<ll> ret;
    for (auto p: primes) {
        if (p * p > x)
            break;
        if (x % p == 0) {
            ret.push_back(p);
            while (x % p == 0) {
                x /= p;
            }
        }
    }
    if (x != 1)
        ret.push_back(x);
    return ret;
}

欧几里得（gcd）

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

扩展欧几里得（exgcd）

求出 \(ax + by = \gcd(a, b)\) 的一组特解 \(x_0,y_0\)。

通解：\(x=x_0+kb/\gcd,y=y_0-ka/\gcd\)。

1
2
3
4
5
6
7
8
9
10
11
ll exgcd(ll a, ll b, ll &x, ll &y) {
    ll ret = 0;
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    ret = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return ret;
}

求欧拉函数 \(\varphi(x)\)

直接求：质因数分解。设 \(n = \prod_{i=1}^{s}p_i^{k_i}\)，其中 \(p_i\) 是质数，有 \(\varphi(n) = n \times \prod_{i = 1}^s{\dfrac{p_i - 1}{p_i}}\)。

线性筛：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
void init() {
    notprime[0] = notprime[1] = 1;
    for (ll i = 2; i <= M - 10; i++) {
        if (notprime[i] == 0) {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (auto j: primes) {
            if (i * j > M - 10)
                break;
            notprime[i * j] = 1;
            if (i % j != 0) {
                phi[i * j] = phi[i] * phi[j];//积性函数
            }
            if (i % j == 0) {
                phi[i * j] = phi[i] * j;//根据“直接求”
                break;
            }
        }
    }
}

欧拉定理

与欧拉函数紧密相关的一个定理就是欧拉定理。其描述如下：

若 \(\gcd(a, m) = 1\)，则 \(a^{\varphi(m)} \equiv 1 \pmod{m}\)。

扩展欧拉定理

当然也有扩展欧拉定理，用于处理一般的 \(a\) 和 \(m\) 的情形。

\[ a^b\equiv \begin{cases} a^{b\bmod\varphi(m)},\,&\gcd(a,\,m)=1\\ a^b,&\gcd(a,\,m)\ne1,\,b<\varphi(m)\\ a^{b\bmod\varphi(m)+\varphi(m)},&\gcd(a,\,m)\ne1,\,b\ge\varphi(m) \end{cases} \pmod m \]

中国剩余定理（CRT）

包含在任意模意义下的逆元求法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
typedef long long ll;
const ll N = 10 + 10;
ll n, a[N], m[N], M = 1, Mi, ti, ans = 0;

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

ll inv(ll x, ll mod) {//在任意模意义下的逆元
    ll y;
    exgcd(x, mod, x, y);
    return (x % mod + mod) % mod;
}

int main() {
    scanf("%lld", &n);
    for (ll i = 1; i <= n; i++) {
        scanf("%lld%lld", m + i, a + i);
        M *= m[i];
    }
    for (ll i = 1; i <= n; i++) {
        Mi = M / m[i];
        ti = inv(Mi, m[i]);
        ans = (ans + a[i] * Mi * ti % M) % M;
    }
    printf("%lld\n", ans);
    return 0;
}

同余最短路

当出现形如「给定 \(n\) 个整数，求这 \(n\) 个整数能拼凑出多少的其他整数（\(n\) 个整数可以重复取）」，以及「给定 \(n\) 个整数，求这 \(n\) 个整数不能拼凑出的最小（最大）的整数」，或者「至少要拼几次才能拼出模 \(K\) 余 \(p\) 的数」的问题时可以使用同余最短路的方法。

同余最短路利用同余来构造一些状态，可以达到优化空间复杂度的目的。

类比差分约束方法，利用同余构造的这些状态可以看作单源最短路中的点。同余最短路的状态转移通常是这样的 \(f(i+y) = f(i) + y\)，类似单源最短路中 \(f(v) = f(u) +edge(u,v)\)。

给定 \(x，y，z，h\)，对于 \(k \in [1,h]\)，有多少个 \(k\) 能够满足 \(ax+by+cz=k\)。（\(0\leq a,b,c\)，\(1\le x,y,z\le 10^5\)，\(h\le 2^{63}-1\)）

不妨假设 \(x < y < z\)。

令 \(d_i\) 为只通过 \(x\) 和 \(y\)，需满足 \(p\bmod x = i\) 能够达到的最低楼层 \(p\)，即 \(x\) 和 \(y\) 操作后能得到的模 \(x\) 下与 \(i\) 同余的最小数，用来计算该同余类满足条件的数个数。

可以得到两个状态：

  • \(i \xrightarrow{y} (i+y) \bmod x\)

  • \(i \xrightarrow{z} (i+z) \bmod x\)

注意通常选取一组 \(a_i\) 中最小的那个数对它取模，也就是此处的 \(x\)，这样可以尽量减小空间复杂度（剩余系最小）。

那么实际上相当于执行了最短路中的建边操作：

add(i, (i+y) % x, y)

add(i, (i+z) % x, z)

接下来只需要求出 \(d_0, d_1, d_2, \dots, d_{x-1}\)，只需要跑一次最短路就可求出相应的 \(d_i\)。

Lucas 定理

含有组合数模板。

对于素数 \(p\)，有 \[ \binom{n}{k}\equiv \binom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}\binom{n\bmod p}{k\bmod p}\pmod p. \]

其中，当 \(n<k\) 时，二项式系数 \(\dbinom{n}{k}\) 规定为 \(0\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
typedef long long ll;
const ll N = 200000 + 10;
ll n, m, p, mul[N], inv[N];

void init() {
    mul[0] = mul[1] = 1;
    for (ll i = 2; i <= n + m; i++)
        mul[i] = mul[i - 1] * i % p;
    inv[1] = 1;
    for (ll i = 2; i <= p - 1; i++)
        inv[i] = ((-p / i * inv[p % i]) % p + p) % p;
    return;
}

ll C(ll n, ll m) {
    if (n < m)
        return 0;
    return mul[n] * inv[mul[m] * mul[n - m] % p] % p;
}

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

int main() {
    ll T;
    scanf("%lld", &T);
    while (T--) {
        scanf("%lld%lld%lld", &n, &m, &p);
        init();
        printf("%lld\n", Lucas(n + m, n, p));
    }
    return 0;
}

卡特兰数

含有递推求逆元。

Catalan 数满足如下递推关系：

\[ C_n = \begin{cases} 1, & n = 0, \\ \sum_{i=0}^{n-1} C_{i}C_{n-1-i}, & n > 0. \end{cases} \] Catalan 数 \(C_n\) 的递推关系有着天然的递归结构：规模为 \(n\) 的计数问题 \(C_n\)，可以通过枚举分界点，分拆为两个规模分别为 \(i\) 和 \((n-1-i)\) 的子问题。这一递推关系使得 Catalan 数广泛出现于各类具有类似递归结构的问题中。

Catalan 数有如下常见的表达式：

\[ C_n = \frac{1}{n+1}\binom{2n}{n} = \dfrac{(2n)!}{n!(n+1)!},~ n\ge 0. \]

\[ C_n = \binom{2n}{n} - \binom{2n}{n+1},~n \ge 0. \]

\[ C_n = \frac{(4n-2)}{n+1}C_{n-1},~ n > 0,~ C_0 = 1. \]

递推求卡特兰数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int p = 998244353; //模
int inv[510]; //逆元
int f[510]; //卡特兰数
void init() {
    inv[1] = 1;
    for (int i = 2; i <= 501; i++) {
        inv[i] = (p - p / i) * inv[p % i] % p;
    }
    f[1] = 1;
    for (int i = 2; i <= 500; i++) {
        f[i] = f[i - 1] * (4 * i - 2) % p * inv[i + 1] % p;
    }
    return;
}

拉格朗日插值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
scanf("%lld%lld", &n, &k);
for (ll i = 1; i <= n; i++)
    scanf("%lld%lld", x + i, y + i);

for (ll i = 1; i <= n; i++) {
    ll tmpa = y[i] % mod, tmpb = 1;
    for (ll j = 1; j <= n; j++) {
        if (i != j) {
            tmpa = tmpa * (k - x[j]) % mod;
            tmpb = tmpb * (x[i] - x[j]) % mod;
        }
    }
    ans += tmpa * inv(tmpb) % mod;
    ans = (ans + mod) % mod;
}

printf("%lld", ans % mod);

康托展开

对于一个排列 \(\{p_n\}\)，它的排名 \[ ans=\sum_{i=1}^na_i(n-i)! \] 其中 \(a_i\) 表示在当前未出现的元素中的排名。\((n-i)!\) 表示后面要填的数字进行全排列的结果。

要求 \(a_i\) 的值实际上就是求逆序对个数，可使用树状数组解决。

可重集全排列：\(n!/\prod c_i\)，\(c_i\) 为重数。

高斯消元

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
typedef double db;
const int N = 100 + 10;
const db eps = 0.000001;
int n;
char ERR[] = "No Solution";
db p[N][N];

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n + 1; j++)
            scanf("%lf", &p[i][j]);
    for (int i = 1; i <= n; i++) {
        int maxn = i;
        for (int j = i; j <= n; j++) {
            if (fabs(p[j][i]) < fabs(p[maxn][i]))
                continue;
            maxn = j;
        }
        for (int j = 1; j <= n + 1; j++)
            swap(p[i][j], p[maxn][j]);
        if (fabs(p[i][i]) <= eps) {
            printf("%s\n", ERR);
            return;
        }
        for (int j = i + 1; j <= n + 1; j++)
            p[i][j] /= p[i][i];
        for (int j = 1; j <= n; j++) {
            if (i == j)
                continue;
            for (int k = i + 1; k <= n + 1; k++)
                p[j][k] -= p[j][i] * p[i][k];
        }
    }
    for (int i = 1; i <= n; i++)
        printf("%.2lf\n", p[i][n + 1]);
    return 0;
}

矩阵快速幂

用于矩阵加速。例如：\(a_i=xa_{i−1}+ya_{i−2}, S_i=\sum a_i^2\)。

\[ \begin{pmatrix} S_i \\ x^2{a_i}^2+y^2{a_{i-1}}^2+2xya_ia_{i-1} \\ {a_i}^2 \\ x{a_i}^2+ya_ia_{i-1} \end{pmatrix}= \begin{pmatrix} 1&1&0&0 \\ 0&x^2&y^2&2xy \\ 0&1&0&0 \\ 0&x&0&y \end{pmatrix}\times \begin{pmatrix} S_{i-1} \\ {a_i}^2 \\ {a_{i-1}}^2 \\ {a_ia_{i-1}} \end{pmatrix} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
typedef long long ll;
const ll N = 5, mod = 1000000007;

struct Matrix //方阵
{
    ll M[N][N], n;
    Matrix(ll type, ll s) //若 type 不为 0，构造一个单位矩阵
    {
        memset(M, 0, sizeof(M));
        n = s;
        if (type)
            for (ll i = 1; i <= N - 1; i++)
                M[i][i] = 1;
    }
    Matrix operator*(const Matrix &N) //矩阵乘法
    {
        Matrix ret(0, n);
        for (ll i = 1; i <= n; i++)
            for (ll k = 1; k <= n; k++)
                for (ll j = 1; j <= n; j++) //j,k 互换可能会更快
                    ret.M[i][j] = (ret.M[i][j] + M[i][k] * N.M[k][j] % mod) % mod;
        return ret;
    }
    Matrix operator^(ll k) //矩阵快速幂
    {
        Matrix ret(1, n), tmp = *this;
        while (k) {
            if (k % 2)
                ret = ret * tmp;
            tmp = tmp * tmp;
            k /= 2;
        }
        return ret;
    }
    void operator^=(ll k) { *this = *this ^ k; }
};