数学知识算法模板

Ayanami_2147483647 · 2025-02-22 15:43:26 · 算法·理论

试除法判定质数

#include <iostream>

using namespace std;

const int N = 110;

int n;

bool is_prime(int val) {
    if (val < 2) return false;
    for (int i = 2; i <= val / i; ++i) {
        if (val % i == 0) return false;
    }

    return true;
}

int main() {
    scanf("%d", &n);

    while (n--) {
        int val;
        scanf("%d", &val);
        if (is_prime(val)) puts("Yes");
        else puts("No");
    }

    return 0;
}

线性筛质数

Acwing1292.哥德巴赫猜想

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

const int N = 1e6 + 10;

int primes[N], cnt;
bool visited[N];

void init(int val) {
    for (int i = 2; i <= val; ++i) {
        if (!visited[i]) primes[cnt++] = i;
        for (int j = 0; primes[j] * i <= val; ++j) {
            visited[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        } 
    }
}

int main() {
    init(N - 1);

    int n;
    while (cin >> n, n) {
        for (int i = 1; ; ++i) {
            int val1 = primes[i];
            int val2 = n - val1;

            if (!visited[val2]) {
                printf("%d = %d + %d\n", n, val1, val2);
                break;
            }
        }
    }

    return 0;
}

Acwing1293.夏洛克和他的女朋友

每一条边都是从质数连到非质数, 因此是二分图

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

const int N = 1e5 + 10;

int n;
int primes[N], cnt;
bool vis[N];

void init(int val) {
    for (int i = 2; i <= val; ++i) {
        if (!vis[i]) primes[cnt++] = i;

        for (int j = 0; primes[j] * i <= val; ++j) {
            vis[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int main() {
    scanf("%d", &n);

    init(n + 1);
    // 输出需要的颜色的个数
    if (n <= 2) puts("1");
    else puts("2");

    for (int i = 2; i <= n + 1; ++i) {
        if (!vis[i]) printf("1 ");
        else printf("2 ");
    }

    return 0;
}

Acwing195.质数距离

寻找在区间[l, r]之间, 相邻质数差值最大和最小的一对

范围比较小的时候, 直接筛出质数, 然后在区间中进行扫描即可, 但该题数据范围非常大, 因此不能直接扫描

如果一个数是合数, 必然存在一个小于\sqrt{n}的因子

先找出[1, 50000]所有的质因子, 对于全部区间的数字, 如果该数字是合数, 必然存在一个[1, 50000]的质因子, 因此将区间中的质因子的倍数全部删去

时间复杂度O(n\log{}{\log{}{n}})

用范围很小的质数处理范围很大的合数

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e6 + 10;

int primes[N], cnt;
bool vis[N];

void init(int val) {
    memset(vis, false, sizeof vis);
    cnt = 0;

    for (int i = 2; i <= val; ++i) {
        if (!vis[i]) primes[cnt++] = i;
        for (int j = 0; i * primes[j] <= val; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
} 

int main() {
    int l, r;
    while (scanf("%d%d", &l, &r) == 2) {
        init(50000);
        memset(vis, false, sizeof vis);

        // 找到大于等于l的最小的p的倍数, 将p的倍数全部删去
        for (int i = 0; i < cnt; ++i) {
            LL p = primes[i];
            for (LL j = max(2 * p, (l + p - 1) / p * p); j <= r; j += p) {
                vis[j - l] = true;
            }
        }

        // 找到区间中剩余的质数
        cnt = 0;
        for (int i = 0; i <= r - l; ++i) {
            if (!vis[i] && i + l >= 2) primes[cnt++] = l + i;
        }

        if (cnt < 2) puts("There are no adjacent primes.");
        else {
            int minp = 0, maxp = 0;
            for (int i = 0; i < cnt - 1; ++i) {
                int d = primes[i + 1] - primes[i];
                if (d < primes[minp + 1] - primes[minp]) minp = i;
                if (d > primes[maxp + 1] - primes[maxp]) maxp = i;
            }

            printf("%d,%d are closest, %d,%d are most distant.\n", primes[minp], primes[minp + 1], primes[maxp], primes[maxp + 1]);
        }
    }

    return 0;
}

注意: 忘记开long long也会导致段错误

分解质因数之阶乘分解

Acwing197.阶乘分解

• 筛出[1, 10^6]中所有的质数
• 求出所有p的次数

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

const int N = 1e6 + 10;

int primes[N], cnt;
bool vis[N];

void init(int val) {
    for (int i = 2; i <= val; ++i) {
        if (!vis[i]) primes[cnt++] = i;

        for (int j = 0; i * primes[j] <= val; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int main() {
    int n;
    scanf("%d", &n);

    init(n);

    for (int i = 0; i < cnt; ++i) {
        int res = 0;
        int p = primes[i];
        for (int j = n; j; j /= p) res += j / p;

        printf("%d %d\n", p, res);
    }

    return 0;
}

快速幂

Acwing1289.序列的第k个数

将情况分为等差数列和等比数列

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int MOD = 200907;

LL quick_pow(int val, int k) {
    LL res = 1;

    while (k) {
        if (k & 1) res = (LL) res * val % MOD;
        val = (LL) val * val % MOD;
        k >>= 1;
    }

    return res;
}

int main() {
    int cases;

    scanf("%d", &cases);
    while (cases--) {
        int a, b, c, k;
        scanf("%d%d%d%d", &a, &b, &c, &k);
        if (b - a == c - b) {
            int d = b - a;
            LL res = a + (LL) (k - 1) * d % MOD;
            printf("%lld\n", res % MOD);
        }
        else {
            int d = b / a;
            LL res = (LL) a * quick_pow(d, k - 1) % MOD;
            printf("%lld\n", res % MOD);
        }
    }

    return 0;
}

Acwing1290.越狱

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int MOD = 100003;

int m;
LL n;

LL quick_pow(int val, LL k) {
    LL res = 1;

    while (k) {
        if (k & 1) res = (LL) res * val % MOD;
        val = (LL) val * val % MOD;
        k >>= 1;
    }

    return res;
}

int main() {
    scanf("%d%lld", &m, &n);

    LL res = quick_pow(m, n) - (LL) m * quick_pow(m - 1, n - 1) % MOD + MOD;

    printf("%lld\n", res % MOD);

    return 0;
}

快速幂求逆元

Acwing876.快速幂求逆元

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e5 + 10;

int n;
int a, p;

int quick_pow(int val, int k, int mod) {
    int res = 1 % mod;

    while (k) {
        if (k & 1) res = (LL) res * val % mod;
        val = (LL) val * val % mod;
        k >>= 1;
    }

    return res;
}

int main() {
    scanf("%d", &n);
    while (n--) {
        scanf("%d%d", &a, &p);
        if (a % p == 0) puts("impossible");
        else {
            int res = quick_pow(a, p - 2, p);
            printf("%d\n", res);
        }
    }

    return 0;
}

约数

约数个数: (c_1 + 1)(c_2 + 1)(...)(c_k + 1)

• 最大次数30, 因为2^{31} > 2 \times 10 ^ 9
• 所有质因子次数递减, 一定是降序, 因为如果不是降序交换一下一定会变得更小

因为搜索空间非常小, 因此直接DFS

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int primes[9] = {2, 3, 5, 7, 11, 13, 17, 19, 23};

int n;
int cnt, res;

void dfs(int u, int last, int curr, int sum) {
    if (sum > cnt || sum == cnt && curr < res) {
        cnt = sum;
        res = curr;
    }

    if (u == 9) return;

    for (int i = 1; i <= last; ++i) {
        if ((LL) curr * primes[u] > n) break;
        curr *= primes[u];
        dfs(u + 1, i, curr, sum * (i + 1));
    }
}

int main() {
    scanf("%d", &n);

    dfs(0, 30, 1, 1);

    printf("%d\n", res);
    return 0;
}

Acwing200.Hankson的趣味题

给定4个正整数, 有多少个正整数x, gcd(a, x) = b, 并且[c, x] = d, 数据范围很大

[a, b]$最小公倍数等于$\frac {a \times b}{(a, b)}

直接枚举约数时间复杂度是O(\sqrt n), 用分解质因数方法计算效率会高很多

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 50010;

int primes[N], cnt;
bool vis[N];
struct Factor {
    int p, s;
} factors[10];
int fcnt;
// 存储约数
int divisor[N], dcnt;

void init() {
    for (int i = 2; i <= N - 1; ++i) {
        if (!vis[i]) primes[cnt++] = i;

        for (int j = 0; i * primes[j] <= N - 1; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

void dfs(int u, int curr) {
    if (u == fcnt) {
        divisor[dcnt++] = curr;
        return;
    }

    for (int i = 0; i <= factors[u].s; ++i) {
        dfs(u + 1, curr);
        curr *= factors[u].p;
    }
}

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

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    init();

    int cases;
    cin >> cases;

    while (cases--) {
        int a, b, c, d;
        // 求一个x使得(x, a) = b 并且[x, c] = d
        cin >> a >> b >> c >> d;

        // 将d分解质因数
        fcnt = 0;
        int val = d;
        for (int i = 0; primes[i] <= val / primes[i]; ++i) {
            int p = primes[i];
            if (val % p == 0) {
                int s = 0;
                while (val % p == 0) val /= p, s++;
                factors[fcnt++] = {p, s};
            }
        }

        if (val > 1) factors[fcnt++] = {val, 1};

        dcnt = 0;
        // 搜索d的所有约数
        dfs(0, 1);

        // 枚举d的所有约数, 找到符合条件的一个
        int res = 0;
        for (int i = 0; i < dcnt; ++i) {
            int val = divisor[i];
            if (gcd(a, val) == b && (LL) c * val / gcd(c, val) == d) res++;
        }

        cout << res << endl;
    }

    return 0;
}

欧拉函数

定义: 1到n中与n互质的数的个数

线性筛法求欧拉函数

在筛质数的时候同时计算出欧拉函数

• 1的时候, 欧拉函数的值是1
• 当p_j是i的最小质因子, \varphi (i \times p_j) = p_j \times \varphi(i)
• 当p_j不是i的最小质因子, \varphi(i \times p_j) = (p_j - 1)\cdot \varphi(i)

Acwing201.可见的点

Acwing220.最大公约数

求最大公约数是素数的数的个数 将最大公约数是p的数的个数, 等价于

最大公约数是1的数对的个数, 也就是互质的数的个数, 与上一题类似

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e7 + 10;

int n;
int primes[N], cnt;
bool vis[N];
int phi[N];
LL sum[N];

void init(int val) {
    phi[1] = 0;

    for (int i = 2; i <= val; ++i) {
        if (!vis[i]) {
            primes[cnt++] = i;
            phi[i] = i - 1;
        }

        for (int j = 0; i * primes[j] <= val; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                break;
            }
            phi[i * primes[j]] = phi[i] * (primes[j] - 1);
        }
    }

    for (int i = 1; i <= val; ++i) sum[i] = sum[i - 1] + phi[i];
}

int main() {
    scanf("%d", &n);
    init(n);

    LL res = 0;
    // 枚举所有质数
    for (int i = 0; i < cnt; ++i) {
        int p = primes[i];
        res += 2 * sum[n / p] + 1;
    }

    printf("%lld\n", res);

    return 0;
}

扩展欧几里得算法

裴蜀定理: 给定任意正整数a, b, 一定存在整数x, y使得ax + by = (a, b)

扩展欧几里得算法求同余方程

Acwing203.同余方程

求最小正整数解, 也就是求最小正余数

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;

int exgcd(int a, int b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }

    int gcd = exgcd(b, a % b, y, x);
    y -= (LL) (a / b) * x;

    return gcd;
}

int main() {
    int a, b;
    cin >> a >> b;

    LL x, y;
    exgcd(a, b, x, y);

    cout << (x % b + b) % b << endl;

    return 0;
}

同余综合应用

Acwing222.青蛙的约会

求最小正整数解

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>

using namespace std;

typedef long long LL;

LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }

    LL d = exgcd(b, a % b, y, x);
    y -= (a / b) * x;

    return d;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    LL a, b, m, n, len;
    cin >> a >> b >> m >> n >> len;

    LL x, y;
    LL gcd = exgcd(m - n, len, x, y);

    // 如果gcd无法转化为b - a一定无解
    if ((b - a) % gcd) puts("Impossible");
    else {

        x = x * (b - a) / gcd;

        // 求最小的正整数解, 就是求正余数
        LL mod = abs(len / gcd);
        cout << (x % mod + mod) % mod << endl;
    }
    return 0;
}

Acwing.202最幸运的数字

将问题转化为上述形式, 将9L和8的公因子筛掉

证明如下

在将问题转化为欧拉定理, 但是求的是最小的x, 欧拉定理求的不是最小的解, 需要继续推导

结论: 最小的x一定是\varphi(c)的约数

枚举所有约数, 找到满足要求的\varphi(c)的约数

该问题的时间复杂度瓶颈在求约数: O(\sqrt n)

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;

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

// 计算欧拉函数
LL get_euler(LL val) {
    LL res = val;
    for (LL i = 2; i <= val / i; ++i) {
        // 如果i是val的质因子, 根据定义计算
        if (val % i == 0) {
            while (val % i == 0) val /= i;
            res = res / i * (i - 1);
        }
    }
    if (val > 1) res = res / val * (val - 1);

    return res;
}

LL quick_mul(LL val, LL k, LL mod) {
    LL res = 0;

    while (k) {
        if (k & 1) res = (res + val) % mod;
        val = (val + val) % mod;
        k >>= 1;
    }

    return res;
}

LL quick_pow(LL val, LL k, LL mod) {
    LL res = 1;

    while (k) {
        if (k & 1) res = quick_mul(res, val, mod);
        val = quick_mul(val, val, mod);
        k >>= 1;
    }

    return res;
}

int main() {
    int cases = 1;
    LL len;

    while (scanf("%lld", &len), len) {
        LL c = gcd(len, 8);
        c = 9 * len / c;
        LL phi = get_euler(c);

        LL res = 1e18;
        if (c % 2 == 0 || c % 5 == 0) res = 0;
        else {
            for (LL i = 1; i * i <= phi; ++i) {
                // 枚举phi的约数
                if (phi % i == 0) {
                    LL d1 = i;
                    LL d2 = phi / i;

                    if (quick_pow(10, d1, c) == 1) res = min(res, d1);
                    if (quick_pow(10, d2, c) == 1) res = min(res, d2);
                }
            }
        }

        printf("Case %d: %lld\n", cases++, res);
    }

    return 0;
}

中国剩余定理

Acwing1298.曹冲养猪

用扩展欧几里得欧几里得算法求M_i对m_i的逆元

**如果计算出的x是负数, 对M取模***

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 20;

int n;
int arr1[N], arr2[N];

LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }

    LL d = exgcd(b, a % b, y, x);
    y -= (a / b) * x;

    return d;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    LL mul = 1;
    cin >> n;

    for (int i = 1; i <= n; ++i) {
        cin >> arr1[i] >> arr2[i];
        mul *= arr1[i];
    }

    LL res = 0;
    for (int i = 1; i <= n; ++i) {
        LL Mi = mul / arr1[i];
        LL ti, y;
        // 求Mi对mi的逆元
        exgcd(Mi, arr1[i], ti, y);
        // 计算中国剩余定理
        res = (res + (__int128) arr2[i] * Mi * ti) % mul;
    }

    cout << (res % mul + mul) % mul << endl;

    return 0;
}

中国剩余定理扩展

Acwing204.奇怪的整数表达

推导过程

每次将新的方程合并到现有方程当中

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;

int n;

LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }

    LL d = exgcd(b, a % b, y, x);
    y -= (a / b) * x;
    return d;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> n;
    LL x = 0, a1, m1;
    cin >> a1 >> m1;

    // 合并两个式子
    for (int i = 0; i < n - 1; ++i) {
        LL a2, m2;
        cin >> a2 >> m2;

        LL k1, k2;
        LL gcd = exgcd(a1, a2, k1, k2);

        // 如果m2 - m1不能整除最大公约数 无解
        if (abs(m2 - m1) % gcd) {
            x = -1;
            break;
        }

        k1 *= (m2 - m1) / gcd;
        LL mod = a2 / gcd;
        // 因为数据范围的限制, 每次计算新的k都只能用最小正整数解进行迭代
        k1 = (k1 % mod + mod) % mod;

        m1 = a1 * k1 + m1;
        // 新的a是a1, a2的最小公倍数
        a1 = a1 / gcd * a2;
    }

    if (x == -1) cout << -1 << endl;
    // 最后剩余一个式子就是求m1 % a1的正余数
    else cout << (m1 % a1 + a1) % a1 << endl;

    return 0;
}

矩阵乘法

Acwing1303.斐波那契数列的前n项和

因为矩阵乘法具有结合律但是没有交换律, 因此能用快速幂计算

要求保证第一个矩阵的行数和第二个矩阵的列数相等

第一个矩阵的第i行和第j个矩阵的第j列相乘得到结果矩阵的(i,j)位置的值

• 枚举结果矩阵的每一行每一列
• 枚举B的位置

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 3;

int n, mod;

void mul(int res[][N], int a[][N], int b[][N]) {
    static int tmp[N][N];
    memset(tmp, 0, sizeof tmp);

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            for (int k = 0; k < N; ++k) {
                tmp[i][j] = (tmp[i][j] + (LL) a[i][k] * b[k][j]) % mod;
            }
        }
    }

    memcpy(res, tmp, sizeof tmp);
}

int main() {
    cin >> n >> mod;
    int f[3][3] = {1, 1, 1};

//  系数矩阵
    int m[3][3] = {{0, 1, 0},
                   {1, 1, 1},
                   {0, 0, 1}};
    n--;
//  快速幂求矩阵乘法
    while (n) {
        if (n & 1) mul(f, f, m);
        mul(m, m, m);
        n >>= 1;
    }

    cout << f[0][2] % mod << endl;

    return 0;
}

Acwing1304.佳佳的斐波那契数列

为了写代码的方便, 可以将一维矩阵扩充成系数矩阵大小, 但是只有第一行有值

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 4;

int n, mod;

void mul(LL res[][N], LL a[][N], LL b[][N]) {
    static LL tmp[N][N];
    memset(tmp, 0, sizeof tmp);

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            for (int k = 0; k < N; ++k) {
                tmp[i][j] = (tmp[i][j] + (LL) a[i][k] * b[k][j]) % mod;
            }
        }
    }

    memcpy(res, tmp, sizeof tmp);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> n >> mod;
//  fn, fn+1, sn, pn
    LL f[4][4] = {1, 1, 1, 0};
    LL m[4][4] = {{0, 1, 0, 0},
                  {1, 1, 1, 0},
                  {0, 0, 1, 1},
                  {0, 0, 0, 1}};

    int k = n - 1;
    while (k) {
        if (k & 1) mul(f, f, m);
        mul(m, m, m);
        k >>= 1;
    }

    LL res = n * f[0][2] - f[0][3];

    cout << (res % mod + mod) % mod << endl;

    return 0;
}

矩阵乘法优化状态机DP

Acwing1305.GT考试

#include <iostream>
#include <cstring>

using namespace std;

const int N = 25;

int n, len, mod;
char str[N];
int _next[N], m[N][N];

void mul(int res[][N], int a[][N], int b[][N]) {
    static int tmp[N][N];
    memset(tmp, 0, sizeof tmp);

    for (int i = 0; i < len; ++i) {
        for (int j = 0; j < len; ++j) {
            for (int k = 0; k < len; ++k) {
                tmp[i][j] = (tmp[i][j] + a[i][k] * b[k][j]) % mod;
            }
        }
    }

    memcpy(res, tmp, sizeof tmp);
}

int quick_pow(int k) {
    int f[N][N] = {1};

    while (k) {
        if (k & 1) mul(f, f, m);
        mul(m, m, m);
        k >>= 1;
    }

    int res = 0;
//  枚举KMP状态机的每个位置
    for (int i = 0; i < len; ++i) res = (res + f[0][i]) % mod;
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> n >> len >> mod;
    cin >> str + 1;

    for (int i = 2, j = 0; i <= len; ++i) {
        while (j && str[j + 1] != str[i]) j = _next[j];
        if (str[j + 1] == str[i]) j++;
        _next[i] = j;
    }

    for (int i = 0; i < len; ++i) {
        for (char c = '0'; c <= '9'; ++c) {
            int k = i;
            while (k && str[k + 1] != c) k = _next[k];
            if (str[k + 1] == c) k++;
            if (k < len) m[i][k]++;
        }
    }

    cout << quick_pow(n) << endl;

    return 0;
}

组合计数

• 加法原理
• 乘法原理
• 组合数
• 排列数
• 容斥原理

递推法求组合数

数字范围小, 但是询问次数多

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 2010, MOD = 1e9 + 7;

int a, b, n;
int f[N][N];

int main() {
    cin >> n;

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j <= i; ++j) {
            if (!j) f[i][j] = 1;
            else f[i][j] = (f[i - 1][j] + f[i - 1][j - 1]) % MOD;
        }
    }

    while (n--) {
        cin >> a >> b;
        cout << f[a][b] << endl;
    }

    return 0;
}

时间复杂度: O(a \times b)

注意: 枚举j的时候要小于等于i, 正因为这样i = 0的情况被直接跳过不会发生段错误

预处理法求组合数

数字范围大, 但是询问次数不大

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e5 + 10, MOD = 1e9 + 7;

int fact[N], infact[N];

LL quick_pow(int val, int k) {
    LL res = 1;

    while (k) {
        if (k & 1) res = (LL) res * val % MOD;
        val = (LL) val * val % MOD;
        k >>= 1;
    }

    return res;
}

LL C(int a, int b) {
    LL res = (LL) fact[a] % MOD * infact[a - b] % MOD * infact[b] % MOD;
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int n;
    cin >> n;

    // 预处理阶乘和阶乘的逆元
    fact[0] = infact[0] = 1;
    for (int i = 1; i < N; ++i) {
        fact[i] = (LL) fact[i - 1] * i % MOD;
        infact[i] = (LL) infact[i - 1] * quick_pow(i, MOD - 2) % MOD;
    }

    while (n--) {
        int a, b;
        cin >> a >> b;
        cout << C(a, b) << endl;
    }

    return 0;
}

时间复杂度: O(n\log{}{n})

Lucas定理求组合数

询问个数非常小, 但是数据范围超级大, 但是模数不能大最好是质数

时间复杂度: O(p\log_{n}{p})

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e5 + 10;

int fact[N], infact[N];
int mod;

LL quick_pow(LL val, LL k) {
    LL res = 1 % mod;

    while (k) {
        if (k & 1) res = (LL) res * val % mod;
        val = (LL) val * val % mod;
        k >>= 1;
    }

    return res;
}

LL C(int a, int b) {
    if (b > a) return 0;
    LL res = (LL) fact[a] * infact[a - b] % mod * infact[b] % mod;;
    return res;
}

LL lucas(LL a, LL b) {
    if (b == 0) return 1;
    return C(a % mod, b % mod) * lucas(a / mod, b / mod) % mod;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int n;
    cin >> n;
    while (n--) {
        LL a, b;
        cin >> a >> b >> mod;

        fact[0] = infact[0] = 1;
        for (int i = 1; i < N; ++i) {
            fact[i] = (LL) fact[i - 1] * i % mod;
            infact[i] = (LL) infact[i - 1] * quick_pow(i, mod - 2) % mod;
        }

        cout << lucas(a, b) << endl;
    }

    return 0;
}

分解质因数 + 高精度计算组合数

没有模数, 直接计算组合数, 需要结合高精度乘法

如何求n!的质数的个数?

参考阶乘分解

#include <iostream>
#include <algorithm>
#include <cstring>
#include <vector>

using namespace std;

const int N = 5010;

int a, b;
int primes[N], cnt;
bool vis[N];
int s[N];

void init(int val) {
    for (int i = 2; i <= val; ++i) {
        if (!vis[i]) primes[cnt++] = i;
        for (int j = 0; i * primes[j] <= val; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

// 计算阶乘有多少个p
int get(int val, int p) {
    int res = 0;

    while (val) {
        res += val / p;
        val /= p;
    }

    return res;
}

void mul(vector<int> &res, int val) {
    int c = 0;
    for (int i = 0; i < res.size(); ++i) {
        c += res[i] * val;
        res[i] = c % 10;
        c /= 10;
    }

    while (c) {
        res.push_back(c % 10);
        c /= 10;
    }
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> a >> b;

    init(a);
    for (int i = 0; i < cnt; ++i) {
        int p = primes[i];
        s[i] = get(a, p) - get(a - b, p) - get(b, p);
    }

    vector<int> res;
    res.push_back(1);

    // 枚举所有质数
    for (int i = 0; i < cnt; ++i) {
        for (int j = 1; j <= s[i]; ++j) {
            mul(res, primes[i]);
        }
    }

    while (res.size() > 0 && res.back() == 0) res.pop_back();

    while (res.size()) {
        cout << res.back();
        res.pop_back();
    }
    cout << "\n";

    return 0;
}

递推法求组合计数问题

Acwing1307.牡牛和牝牛

题目大意: 当前位置是1, 当前位置到前一个位置1的距离不能小于等于k

将集合按照上一个1的位置进行分类 最坏情况下时间复杂度是O(n^2), 使用前缀和进行优化到O(n)

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 1e5 + 10, MOD = 5000011;

int f[N], s[N];

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int n, k;
    cin >> n >> k;

    f[0] = s[0] = 1;
    // 递推res
    for (int i = 1; i < n; ++i) {
        // 因为i - k - 1可能是负数, 因此当左边界是负数时, 将左边界设置为0
        f[i] = s[max(i - k - 1, 0)];
        // 计算前缀和
        s[i] = (s[i - 1] + f[i]) % MOD;
    }

    cout << s[n] << endl;

    return 0;
}

排列组合问题之方程的解

Acwing1308.方程的解

关于方程的右侧使用后快速幂解决

方程解的数量用隔板法进行计算, 假如三个未知数 那就在小球中放两个隔板, 将整个集合分为三个部分, 每个部分代表一个未知数C_{n - 1}^{k - 1}

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 150;

int k, x;
int f[1000][100][N];

int quick_pow(int val, int k, int mod) {
    int res = 1;

    while (k) {
        if (k & 1) res = res * val % mod;
        val = val * val % mod;
        k >>= 1;
    }

    return res;
}

void add(int res[], int a[], int b[]) {
    int c = 0;
    for (int i = 0; i < N; ++i) {
        c += a[i] + b[i];
        res[i] = c % 10;
        c /= 10;
    }
}

int main() {
    cin >> k >> x;

    int n = quick_pow(x % 1000, x, 1000);
    // 计算C(n - 1, k - 1)
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j <= i && j < k; ++j) {
            if (!j) f[i][j][0] = 1;
            else add(f[i][j], f[i - 1][j], f[i - 1][j - 1]);
        }
    }

    if (!n) cout << 0 << endl;
    else {
        int i = N - 1;
        // 去除前导零
        while (!f[n - 1][k - 1][i]) i--;
        // 输出结果的每一位
        while (i >= 0) cout << f[n - 1][k - 1][i--];
    }

    return 0;
}

组合数排列数求不规则图形

Acwing1309.车的放置

如果是规则图形

• 首先先从n行中选择k行, 保证不在同一行
• 然后考虑列, 第一行的车有m中选择, 第二行的车有k - 1种选择...

将问题拆分成两个规则图形, 因为不知道每个规则图形中放置几个车, 因此可以枚举车在两个规则图形中的数量

对于下面的矩形, 因为上面已经用了i列, 因此是从a + c - i中有顺序的选出k - i个数字

因为组合数排列数都用到阶乘, 因此可以预处理阶乘

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;
const int N = 2010, MOD = 100003;

int fact[N], infact[N];

LL quick_pow(LL val, LL k) {
    LL res = 1;

    while (k) {
        if (k & 1) res = res * val % MOD;
        val = val * val % MOD;
        k >>= 1;
    }

    return res;
}

LL C(int a, int b) {
    if (b > a) return 0;
    LL res = (LL) fact[a] * infact[a - b] % MOD * infact[b] % MOD;
    return res;
}

LL A(int a, int b) {
    if (b > a) return 0;
    LL res = (LL) fact[a] * infact[a - b] % MOD;
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int a, b, c, d, k;
    cin >> a >> b >> c >> d >> k;

    // 求阶乘和阶乘的逆元
    fact[0] = infact[0] = 1;
    for (int i = 1; i < N; ++i) {
        fact[i] = (LL) fact[i - 1] * i % MOD;
        // 因为MOD是质数, 因此使用费马小定理计算逆元
        infact[i] = (LL) infact[i - 1] * quick_pow(i, MOD - 2) % MOD;
    }

    LL res = 0;
    // 枚举上面矩形放置的车的数量
    for (int i = 0; i <= k; ++i) {
        LL up = (LL) C(b, i) * A(a, i) % MOD;
        LL down = (LL) C(d, k - i) * A(a + c - i, k - i) % MOD;
        res = (res + up * down % MOD) % MOD;
    }

    cout << res % MOD << endl;
    return 0;
}

注意: 在计算组合数和排列数时, 一定要添加边界条件

组合数之网格求三角形数量

Acwing1310.数三角形

因为正向枚举所有点的位置不好计算, 因此可以采用补集的思想, 减去所有不合法情况

• 对于斜率大于0的点, 将直线按照横纵坐标差值进行划分, 也就是在矩形中枚举直角三角形, 直角三角形的斜边就是我们要求的线段
• 对于斜边上的两个端点是确定的, 因此需要知道两个端点中在整点的点的数量, 也就是枚举第三个点

为什么是gcd(i, j)能代表线段上在整点的数量?

• 首先将问题转化为简单形式, 将左下角端点移动到原点, 设右上角端点的坐标是(Δx,Δy), 设最大公约数是d
• 那么线段上的点的坐标就能表示为\left( \frac{k \Delta x}{d}, \frac{k \Delta y}{d} \right), 其中k取值0, 1, ..., d
• 算上端点, 一共d + 1的数

为什么是最大公约数, 而不是任意一个公约数?

因为只有最大公约数求出的点的坐标能够代表最小步长, 如果是任意一个公约数, 可能遗漏某些点

图片引用自小呆呆

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;

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

LL C(int n) {
    return (LL) n * (n - 1) * (n - 2) / 6;
}

int main() {
    int row, col;
    cin >> row >> col;

    row++, col++;
    // 所有情况减去斜率为0和斜率是正无穷的情况
    LL res = C(row * col) - (LL)row * C(col) - (LL) col * C(row);

    // 枚举横坐标差值和纵坐标差值
    for (int i = 1; i <= row; ++i) {
        for (int j = 1; j <= col; ++j) {
            // 后面两项是枚举斜边两个端点的取值情况
            res -= 2ll * (gcd(i, j) - 1) * (row - i) * (col - j);
        }
    }

    cout << res << endl;
    return 0;
}

隔板法求不等式解的个数

Acwing1312.序列统计

将原来的范围进行映射, 然后进行差分拆解, 将问题转化为求不等式的解的个数

再将x_i映射到y_i = x_i + 1, 这样保证每个解都是大于等于1的

因为是不等式因此放置k个隔板, 如果是等式放置k - 1个板子

因为枚举k会超时, 因此考虑将原式进行化简

依据组合数递推式进行化简

因为数据范围很大, 但是模数很小, 因此使用lucas定理进行计算组合数

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e6 + 10, MOD = 1e6 + 3;

int cases;
int n, l, r;
int fact[N], infact[N];

int quick_pow(int val, int k) {
    int res = 1;

    while (k) {
        if (k & 1) res = (LL) res * val % MOD;
        val = (LL) val * val % MOD;
        k >>= 1;
    }

    return res;
}

// 处理所有阶乘和阶乘的逆元
void init() {
    fact[0] = infact[0] = 1;

    for (int i = 1; i < N; ++i) {
        fact[i] = (LL) fact[i - 1] * i % MOD;
        infact[i] = (LL) infact[i - 1] * quick_pow(i, MOD - 2) % MOD;
    }
}

int C(int a, int b) {
    if (a < b) return 0;
    return (LL) fact[a] * infact[a - b] % MOD * infact[b] % MOD;
}

int lucas(int a, int b) {
    if (a < MOD && b < MOD) return C(a, b);
    return (LL) C(a % MOD, b % MOD) % MOD * lucas(a / MOD, b / MOD) % MOD;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    init();

    cin >> cases;
    while (cases--) {
        cin >> n >> l >> r;
        int res = (lucas(r - l + n + 1, r - l + 1) - 1 + MOD) % MOD;
        cout << res << endl;
    }

    return 0;
}

因为数量一定是整数, 因此在输出最终结果的时候一定要计算正余数

卡特兰序列

Acwing889.满足条件的01序列

求满足任意前缀中0的个数不少于1的个数的序列的数量

将问题转化为从原点走路径的问题

任意一条路径不能走到红颜色线上!

再将卡特兰数进行转化, \frac{C_{2n}^{n}}{n + 1}

• 卡特兰数递推式, f(n) = f(1) \times (n - 1) + f(2) \times (n - 2) + ...
• 挖掘卡特兰数性质, 任意前缀里面某种东西的个数不能少于另一种东西的个数

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;
const int MOD = 1e9 + 7;

int quick_pow(int val, int k) {
    int res = 1;
    while (k) {
        if (k & 1) res = (LL) res * val % MOD;
        val = (LL) val * val % MOD;
        k >>= 1;
    }

    return res;
}

int C(int a, int b) {
    if (a < b) return 0;

    int val1 = 1;
    for (int i = a; i >= a - b + 1; --i) {
        val1 = (LL) val1 * i % MOD;
    }

    // 计算b阶乘的逆元
    int val2 = 1;
    for (int i = b; i >= 1; --i) {
        val2 = (LL) val2 * quick_pow(i, MOD - 2) % MOD;
    }

    return (LL) val1 * val2 % MOD;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int n;
    cin >> n;
    int res = (LL) C(2 * n, n) * quick_pow(n + 1, MOD - 2) % MOD;
    cout << (res % MOD + MOD) % MOD << endl;

    return 0;
}

组合数学之卡特兰数

Acwing1315.网格

因为此题没有模数, 因此需要第四种求组合数的方法, 也就是分解质因数 + 阶乘分解 + 高精度

图片引用自文章链接

#include <iostream>
#include <algorithm>
#include <cstring>
#include <vector>

using namespace std;

const int N = 10010;

int row, col;
int primes[N], cnt;
bool vis[N];

void init() {
    for (int i = 2; i <= N - 1; ++i) {
        if (!vis[i]) primes[cnt++] = i;

        for (int j = 0; i * primes[j] <= N - 1; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int get(int val, int p) {
    int res = 0;
    while (val) res += val / p, val /= p;
    return res;
}

//高精度乘法
void mul(vector<int> &res, int val) {
    int c = 0;
    for (int i = 0; i < res.size(); ++i) {
        c += res[i] * val;
        res[i] = c % 10;
        c /= 10;
    }

    while (c) res.push_back(c % 10), c /= 10;
}

vector<int> C(int a, int b) {
    vector<int> res;
    if (a < b) return res;

    res.push_back(1);
//  记录每个质数出现的次数
    int tmp[N] = {0};
    for (int i = 0; i < cnt; ++i) {
        int p = primes[i];
        tmp[i] = get(a, p) - get(a - b, p) - get(b, p);
    }

//  计算高精度乘法
    for (int i = 0; i < cnt; ++i) {
        for (int j = 0; j < tmp[i]; ++j) {
            mul(res, primes[i]);
        }
    }

    return res;
}

//高精度减法
void sub(vector<int> &res, vector<int> &val) {
    int c = 0;
    for (int i = 0; i < res.size(); ++i) {
        int tmp = res[i] - (i < val.size() ? val[i] : 0) - c;
        if (tmp < 0) {
            res[i] = tmp + 10;
            c = 1;
        }
        else res[i] = tmp, c = 0;
    }
    while (res.size() > 1 && res.back() == 0) res.pop_back();
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    init();
    cin >> row >> col;

    vector<int> res = C(row + col, row);
    vector<int> val = C(row + col, col - 1);
    sub(res, val);

    for (int i = res.size() - 1; i >= 0; --i) cout << res[i];
    cout << "\n";

    return 0;
}

Acwing1316.有趣的数列

如果发现数列的第三项是5, 大概率是卡特兰序列

对于当前位置是偶数项, 必须满足前面数字中奇数项的个数不能小于偶数项的个数, 如果奇数项个数小于偶数项个数, 必然至少存在某个奇数项是空的, 一定会从当前偶数位置的后面取值, 因为相邻两项的左侧奇数项要求小于等于偶数项, 但是奇数项是从当前偶数项后面选择, 因此一定大于前一个偶数项, 矛盾, 因此奇数项数量大于等于偶数项的数量

因此就将问题转化为了卡特兰序列问题, 因为模数比较大, 因此不能用lucas定理

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 2e6 + 10;

int n, mod;
int primes[N], cnt;
bool vis[N];

void init() {
    for (int i = 2; i < N; ++i) {
        if (!vis[i]) primes[cnt++] = i;
        for (int j = 0; i * primes[j] < N; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int get(int val, int p) {
    int res = 0;

    while (val) {
        res += val / p;
        val /= p;
    }

    return res;
}

int quick_pow(int val, int k) {
    int res = 1 % mod;

    while (k) {
        if (k & 1) res = (LL) res * val % mod;
        val = (LL) val * val % mod;
        k >>= 1;
    }

    return res % mod;
}

int C(int a, int b) {
    if (a < b) return 0;
    int tmp[N] = {0};

    for (int i = 0; i < cnt; ++i) {
        int p = primes[i];
        tmp[i] = get(a, p) - get(a - b, p) - get(b, p);
    } 

    int res = 1;
    for (int i = 0; i < cnt; ++i) {
        if (!tmp[i]) continue;
        res = (LL) res * quick_pow(primes[i], tmp[i]) % mod;
    }
    return res;
}

int main() {
    init();
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> n >> mod;
    int res = C(2 * n, n) - C(2 * n, n - 1);
    cout << (res + mod) % mod << endl;

    return 0;
}

高斯消元解线性方程组

Acwing883.高斯消元解线性方程组

高斯消元过程

• 枚举每一列
• 为了精度问题, 找到绝对值最大的一行, 也就是找主元
• 交换到最上面
• 将该行的第一个数变为1
• 把下面的行的这一列c全部消成0

时间复杂度: O(n^ 3)

高斯消元解的情况

• 完美阶梯形就是恰好一组解
• 左边没有未知数, 右边非0, 无解
• 消完之后不是完美阶梯形 不够n个方程, 有无穷多组解

#include <iostream>
#include <algorithm>
#include <cmath>

using namespace std;

const int N = 110;
const double EPS = 1e-8;

int n;
double f[N][N];

int gauss() {
    int row, col;
    for (row = col = 0; col < n; ++col) {
        int u = row;

        for (int i = row + 1; i < n; ++i) {
            if (fabs(f[i][col]) > fabs(f[u][col])) u = i;
        }

        if (fabs(f[u][col]) < EPS) continue;
        for (int i = col; i <= n; ++i) swap(f[row][i], f[u][i]);

        // 归一化
        for (int i = n; i >= col; --i) f[row][i] /= f[row][col];

        for (int i = row + 1; i < n; ++i) {
            if (fabs(f[i][col]) < EPS) continue;
            for (int j = n; j >= col; --j) {
                f[i][j] -= f[row][j] * f[i][col];
            }
        }
        row++;
    }

    if (row < n) {
        for (int i = row; i < n; ++i) if (fabs(f[i][n]) > EPS) return 2;
        return 1;
    }

    for (int i = n - 1; i >= 0; --i) {
        for (int j = i + 1; j < n; ++j) {
            f[i][n] -= f[i][j] * f[j][n];
        }
    }

    return 0;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> n;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n + 1; ++j) {
            cin >> f[i][j];
        }
    }

    int res = gauss();

    if (res == 2) cout << "No solution" << endl;
    else if (res == 1) cout << "Infinite group solutions" << endl;
    else {
        for (int i = 0; i < n; ++i) printf("%.2lf\n", f[i][n]);
    }

    return 0;
}

高斯消元解线性异或方程组

• 在解决哟异或线性方程组的时候, 异或是不进位的加法, 也可以理解为加法或者减法
• 乘法是与的关系

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 110;

int n;
int f[N][N];

int gauss() {
    int row, col;
    for (row = col = 0; col < n; ++col) {
        int u = row;

        for (int i = row; i < n; ++i) {
            if (f[i][col]) {
                u = i;
                break;
            }
        }

        if (f[u][col] == 0) continue;
        for (int i = col; i <= n; ++i) swap(f[row][i], f[u][i]);

        // 消元
        for (int i = row + 1; i < n; ++i) {
            if (!f[i][col]) continue;
            for (int j = n; j >= col; --j) {
                f[i][j] ^= f[row][j] & f[i][col];
            }
        }
        row++;
    }

    if (row < n) {
        for (int i = row; i < n; ++i) if (f[i][n]) return 2;
        return 1;
    }

    for (int i = n - 1; i >= 0; --i) {
        for (int j = i + 1; j < n; ++j) {
            f[i][n] ^= f[i][j] & f[j][n];
        }
    }

    return 0;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    cin >> n;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n + 1; ++j) {
            cin >> f[i][j];
        }
    }

    int res = gauss();
    if (res == 2) cout << "No solution" << "\n";
    else if (res == 1) cout << "Multiple sets of solutions" << "\n";
    else {
        for (int i = 0; i < n; ++i) cout << f[i][n] << "\n";
    }

    return 0;
}

高斯消元应用

Acwing207.球形空间产生器

欧式空间 + Gauss消元

将方程转化为线性方程, 使用Gauss消元求解 球面方程的展开

对于点 P_i = (x_{i1}, x_{i2}, \dots, x_{in})，

球面方程为：(x_{i1} - c_1)^2 + (x_{i2} - c_2)^2 + \dots + (x_{in} - c_n)^2 = r^2

展开后得到：(x_{i1}^2 - 2x_{i1}c_1 + c_1^2) + (x_{i2}^2 - 2x_{i2}c_2 + c_2^2) + \dots + (x_{in}^2 - 2x_{in}c_n + c_n^2) = r^2

减去第一个点的方程 对于第一个点 P_1 = (x_{11}, x_{12}, \dots, x_{1n})，

方程为: (x_{11}^2 - 2x_{11}c_1 + c_1^2) + (x_{12}^2 - 2x_{12}c_2 + c_2^2) + \dots + (x_{1n}^2 - 2x_{1n}c_n + c_n^2) = r^2

将其他点 P_i 的方程减去第一个点的方程：\begin{aligned} & \left[(x_{i1}^2 - 2x_{i1}c_1 + c_1^2) + (x_{i2}^2 - 2x_{i2}c_2 + c_2^2) + \dots + (x_{in}^2 - 2x_{in}c_n + c_n^2)\right] \\ & - \left[(x_{11}^2 - 2x_{11}c_1 + c_1^2) + (x_{12}^2 - 2x_{12}c_2 + c_2^2) + \dots + (x_{1n}^2 - 2x_{1n}c_n + c_n^2)\right] = 0 \end{aligned}

消去平方项和r^2

整理后得到：(x_{i1}^2 - x_{11}^2) - 2(x_{i1} - x_{11})c_1 + (x_{i2}^2 - x_{12}^2) - 2(x_{i2} - x_{12})c_2 + \dots + (x_{in}^2 - x_{1n}^2) - 2(x_{in} - x_{1n})c_n = 0

将方程整理为线性形式：$ 2(x{i1} - x{11})c1 + 2(x{i2} - x_{12})c2 + \dots + 2(x{in} - x_{1n})cn = (x{i1}^2 - x{11}^2) + (x{i2}^2 - x{12}^2) + \dots + (x{in}^2 - x_{1n}^2)

• 否则就是无解

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

const int N = 40;

int n;
int f[N][N];

// 下标范围1~n
int gauss() {
    int row, col;
    for (row = 1, col = 1; col <= n; ++col) {
        int u = row;

        for (int i = row + 1; i <= n; ++i) {
            if (f[i][col]) {
                u = i;
                break;
            }
        }

        if (!f[u][col]) continue;
        for (int i = col; i <= n + 1; ++i) swap(f[row][i], f[u][i]);

        for (int i = row + 1; i <= n; ++i) {
            for (int j = n + 1; j >= col; --j) {
                f[i][j] ^= f[row][j] & f[i][col];
            }
        }
        row++;
    }

    int res = 1;

    if (row <= n) {
        for (int i = row; i <= n; ++i) {
            if (f[i][n + 1]) return -1;
            res *= 2;
        }
    }

    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int cases;
    cin >> cases;

    while (cases--) {
        memset(f, 0, sizeof f);

        cin >> n;
        for (int i = 1; i <= n; ++i) cin >> f[i][n + 1];
        for (int i = 1; i <= n; ++i) {
            int val;
            cin >> val;
            f[i][n + 1] ^= val;
            f[i][i] = 1;
        }

        int val1, val2;
        while (cin >> val1 >> val2, val1 || val2) {
            f[val2][val1] = 1;
        }

        int res = gauss();
        if (res == -1) cout << "Oh,it's impossible~!!" << "\n";
        else cout << res << "\n";
    }

    return 0;
}

容斥原理

Acwing890.能被整除的数

\begin{array}{c|cccccccccc} i & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \ \hline \left\lfloor \frac{k}{i} \right\rfloor & 5 & 2 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \ \end{array}$

可见 \left\lfloor \frac{k}{i} \right\rfloor 分成了 3 块，我们只需要计算 n \times k 减去每一块的和即可

枚举块的左边界根据左边界和k计算出右边界

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    LL n, k;
    cin >> n >> k;
    LL res = n * k;

    for (int l = 1, r; l <= n; l = r + 1) {
        if (k / l != 0) r = min(k / (k / l), n);
        else r = n;
        res -= (k / l) * (r - l + 1) * (l + r) / 2;
    }

    cout << res << endl;
    return 0;
}

莫比乌斯函数

假设有问题求与n互质的数的个数

其中S_i代表i的倍数的个数

对于当前问题, 每一个数前面的系数就是莫比乌斯函数

如何求莫比乌斯函数?

• 当前数是质数, 只有一个质因子, 函数值是-1
• 如果p_j是i的最小质因子, i \times p_j当中至少包含两个p_j, 函数值是0
• 如果未break, i \times p_j 中的所有质因子都只出现过一次 函数值\times -1

莫比乌斯函数应用

Acwing215.破译密码

首先将问题进行映射

直接做的时间复杂度是O(n ^ 2), 不同的数字的个数是有限的, 参考上面的整除分块, 通过该算法将复杂度优化成O(n\sqrt n)

段内部整除的值是相等的, 因此求莫比乌斯函数的前缀和就能求出这一段的值

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 50010;

int cases;
int a, b, d;
int primes[N], cnt;
bool vis[N];
int mobius[N], s[N];

// 筛质数过程中计算莫比乌斯函数
void init() {
    // 质因子数量0偶数
    mobius[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (!vis[i]) {
            primes[cnt++] = i;
            mobius[i] = -1;
        }

        for (int j = 0; i * primes[j] < N; ++j) {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                mobius[i * primes[j]] = 0;
                break;
            }
            mobius[i * primes[j]] = mobius[i] * -1;
        }
    }

    // 预处理莫比乌斯函数前缀和
    for (int i = 1; i < N; ++i) s[i] = s[i - 1] + mobius[i];
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    init();

    cin >> cases;
    while (cases--) {
        cin >> a >> b >> d;
        a /= d, b /= d;

        int n = min(a, b);
        LL res = 0;
        // 整除分块
        for (int l = 1, r; l <= n; l = r + 1) {
            r = min({n, a / (a / l), b / (b / l)});
            res += (s[r] - s[l - 1]) * (LL) (a / l) * (b / l);
        }

        cout << res << "\n";
    }

    return 0;
}

概率与数学期望

Acwing217.绿豆蛙的归宿

计算路径长度的期望

#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

const int N = 1e5 + 10, EDGE_NUM = N << 1;

int n, edge_num;
int head[N], edge_end[EDGE_NUM], next_edge[EDGE_NUM], w[EDGE_NUM], edge_index;
int deg[N];
double f[N];

void add(int ver1, int ver2, int val) {
    edge_end[edge_index] = ver2, next_edge[edge_index] = head[ver1], w[edge_index] = val, head[ver1] = edge_index++;
}

// 记忆化搜索
double dfs(int u) {
    if (f[u] >= 0) return f[u];

    f[u] = 0;
    for (int i = head[u]; ~i; i = next_edge[i]) {
        int ver = edge_end[i];
        // 每次出边概率都是相等的
        f[u] += (w[i] + dfs(ver)) / deg[u];
    }

    return f[u];
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    memset(head, -1, sizeof head);

    cin >> n >> edge_num;
    while (edge_num) {
        int ver1, ver2, val;
        cin >> ver1 >> ver2 >> val;
        add(ver1, ver2, val);
        deg[ver1]++;
    }

    memset(f, -1, sizeof f);
    dfs(1);
    printf("%.2lf", f[1]);

    return 0;
}

Acwing218.扑克牌

问题描述: 给定四堆扑克牌数量, 求最少的翻的牌的数学期望是多少(如果是大小王可以选择), 因为大小王放置到四个牌堆的概率是相等的, 因此可以转化为dp问题

• 当异或值不是0一定可以通过某种操作使其变为0

#include <iostream>

using namespace std;

const int N = 1e5 + 10;

int main() {
    int n;
    cin >> n;
    int res = 0;
    for (int i = 0; i < n; ++i) {
        int val;
        cin >> val;
        res ^= val;
    }
    if (res == 0) cout << "No" << "\n";
    else cout << "Yes" << "\n";

    return 0;
}

Acwing892.台阶Nim游戏

任选一个台阶拿若干个石子放到下一个台阶

结论: 将奇数台阶上石子当做普通Nim游戏

• 如果当前异或值不是0, 一定可以将当前状态转换为0
• 如果对手从偶数台阶拿石子, 我们就把他拿的放到下一层台阶, 使得对手操作后奇数台阶石子不变
• 如果对手拿的奇数石子, 拿完后的异或值一定不是0, 也就是经典Nim游戏

因此将奇数台阶看作经典Nim游戏是正确的

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 1e5 + 10;

int main() {
    int n;
    cin >> n;

    int res = 0;
    for (int i = 1; i <= n; ++i) {
        int val;
        cin >> val;
        if (i % 2) res ^= val;
    }

    if (res) cout << "Yes" << "\n";
    else cout << "No" << "\n";

    return 0;
}

博弈论之SG函数