数学知识

常用代码模板4——数学知识

数论

试除法判定质数

#include <iostream>

using namespace std;

int n, a;

bool is_prime(int n) {
    if (n <= 1) return false;
    for (int i = 2; i <= n / i; ++ i) // 此处不要写成i*i<=n或者i<=sqrt(n)
        if (n % i == 0) return false; // 前者会溢出后者较慢
    return true;
}

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d", &a);
        puts(is_prime(a) ? "Yes" : "No");
    }
}

分解质因数

#include <iostream>

using namespace std;

int n, a;

void divide(int n) {
    for (int i = 2; i <= n / i; ++ i) {
        if (n % i == 0) {
            int s = 0;
            while (n % i == 0) n /= i, s ++;
            printf("%d %d\n", i, s);
        }
    } // n至多有一个大于根号n的质因子
    if (n > 1) printf("%d 1\n", n);
    puts("");
}

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d", &a);
        divide(a);
    }
}

筛质数

#include <iostream>

using namespace std;

const int N = 1e6 + 10;
int n;
int prime[N], cnt;
bool st[N];

void get_primes(int n) {
    for (int i = 2; i <= n; ++ i) {
        if (!st[i]) prime[cnt ++] = i;
        for (int j = 0; prime[j] <= n / i; ++ j) {
            st[prime[j] * i] = true;
            if (i % prime[j] == 0) break;
        }
    }
} // 一个数只用其最小质因数筛掉

int main() {
    scanf("%d", &n);
    get_primes(n);
    printf("%d\n", cnt);
}

试除法求约数

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

using namespace std;

int n, a;

vector<int> get_divisors(int n) {
    vector<int> res;
    for (int i = 1; i <= n / i; ++ i) {
        if (n % i == 0) {
            res.push_back(i);
            if (i != n / i) res.push_back(n / i);
        }
    }
    sort(res.begin(), res.end());
    return res;
}

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d", &a);
        auto res = get_divisors(a);
        for (auto i : res) printf("%d ", i);
        puts("");
    }
}

约数个数

#include <iostream>
#include <unordered_map>

using namespace std;

typedef long long LL;

const int MOD = 1e9 + 7;
int n, a;
unordered_map<int, int> primes;
LL res = 1;

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d", &a);
        for (int i = 2; i <= a / i; ++ i)
            while (a % i == 0) a /= i, primes[i] ++;
        if (a > 1) primes[a] ++;
    }
    for (auto i : primes) res = res * (i.second + 1) % MOD; // 约数个数其实就是质因数的组合数
    printf("%lld\n", res);
}

约数之和

#include <iostream>
#include <unordered_map>

using namespace std;

typedef long long LL;

const int MOD = 1e9 + 7;
int n, a;
unordered_map<int, int> primes;
LL res = 1;

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d", &a);
        for (int i = 2; i <= a / i; ++ i)
            while (a % i == 0) a /= i, primes[i] ++;
        if (a > 1) primes[a] ++;
    }
    for (auto i : primes) {
        int fi = i.first, se = i.second;
        LL t = 1;
        while (se --) t = (t * fi + 1) % MOD;
        res = res * t % MOD;
    }
    printf("%lld\n", res);
}

最大公约数

#include <iostream>

using namespace std;

int n, a, b;

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

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d%d", &a, &b);
        printf("%d\n", gcd(a, b));
    }
}

快速幂

int q_pow(int a, int b, int c) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = (LL)ans * a % c;
        a = (LL)a * a % c;
        b >>= 1;
    }
    return ans;
}

快速幂求逆元

#include <iostream>

using namespace std;

typedef long long LL;

int n, a, p;

int q_pow(int a, int b, int k) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = (LL)ans * a % k;
        a = (LL)a * a % k;
        b >>= 1;
    }
    return ans;
}

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

线性逆元

#include <iostream>
#include <cstdio>

using namespace std;

typedef long long LL;

const int N = 20000528 + 10;
int n, p;
int inv[N];

int main() {
    scanf("%d%d", &n, &p);
    inv[1] = 1;
    for (int i = 2; i <= n; ++ i)
        inv[i] = (LL)(p - p / i) * inv[p % i] % p;
    for (int i = 1; i <= n; ++ i)
        printf("%d\n", inv[i]);
    return 0;
}

拓展欧几里得

int exgcb(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcb(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

线性同余方程

#include <iostream>

using namespace std;

int n, a, b, m, x, y;

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

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%d%d%d", &a, &b, &m);
        int d = exgcd(a, m, x, y);
        if (b % d == 0) printf("%d\n", x * 1ll * b / d % m);
        else printf("impossible\n");
    }
}

中国剩余定理

1

#include <iostream>

using namespace std;

typedef long long LL;

int n;
bool has_answer = true;

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

int main() {
    cin >> n;
    LL a1, m1;
    cin >> a1 >> m1;

    for (int i = 1; i < n; ++ i) {
        LL a2, m2;
        cin >> a2 >> m2;

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

        if ((m2 - m1) % d) {
            has_answer = false;
            break;
        }

        k1 *= (m2 - m1) / d;
        LL t = a2 / d;
        k1 = (k1 % t + t) % t;

        m1 = a1 * k1 + m1;
        a1 = abs(a1 / d * a2);
    }

    cout << (has_answer ? (m1 % a1 + a1) % a1 : -1) << endl;

    return 0;
}

2

#include <iostream>
#include <cmath>

using namespace std;

typedef long long LL;

int n, p, e, i, d;
LL a[4]; int w[] = {0, 23, 28, 33};

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;
}

LL ChRe() {
    int n = 1;
    for (int i = 1; i <= 3; ++ i) a[i] %= w[i], n *= w[i];
    LL ans = 0;
    for (int i = 1; i <= 3; ++ i) {
        int m = n / w[i];
        LL x, y;
        exgcd(m, w[i], x, y);
        ans = (ans + m * x * a[i]) % n;
    }
    ans = (ans % n + n) % n;
    ans -= d;
    while (ans <= 0) ans += n;
    return ans;
}

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

    while (n --) {
        LL Case = 0;
        while (scanf("%lld%lld%lld%lld", &a[1], &a[2], &a[3], &d) != EOF) {
            if (a[1] == -1 && a[2] == -1 && a[3] == -1 && d == -1) break;
            LL ans = ChRe();
            printf("Case %lld: the next triple peak occurs in %lld days.\n", ++ Case, ans);
        }
    }
}

组合计数

公式：(a 选 b)
$$ C ^ b _ a = { {a!} \over {(a - b)! b!} } $$

递推式：
$$ C ^ b _ a = C ^ {b - 1} _ a + C ^ {b - 1}_{a - 1} \\\ C ^ 0 _ a = 1 \\\ C ^ b _ a = 0 (a < b) \\\ \sum ^ {n} _ {i = 0} {C ^ n _ i} = 2 ^ n $$
性质：
$$ C ^ n _ m = C ^ {m - n} _ m $$
数据范围较小时：a,b <= 2e3

#include <iostream>

using namespace std;

const int N = 2e3 + 10, MOD = 1e9 + 7;
int c[N][N];
int n, a, b;

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

int main() {
    init();
    scanf("%d", &n);
    while (n --) {
        scanf("%d%d", &a, &b);
        printf("%d\n", c[a][b]);
    }
}

a,b较大时：a,b <= 1e7

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 1e6 + 10, MOD = 1e9 + 7;
int fact[N], infact[N];
int n, a, b;

int q_pow(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = (LL)res * a % MOD;
        a = (LL)a * a % MOD;
        b >>= 1;
    }
    return res;
}

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

int main() {
    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] * q_pow(i, MOD - 2) % MOD;
    }

    scanf("%d", &n);
    while (n --) {
        scanf("%d%d", &a, &b);
        printf("%d\n", C(a, b));
    }
}

a,b 巨大时：a,b<=1e18 p <=1e5

#include <iostream> // 卢卡斯定理

using namespace std;

typedef long long LL;
int n, p;
LL a, b;

int q_pow(int a, int b, int p) {
    int res = 1;
    while (b) {
        if (b & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        b >>= 1;
    }
    return res;
}

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

    int res = 1;
    for (int i = 1, j = a; i <= b; ++ i, -- j) {
        res = (LL)res * j % p;
        res = (LL)res * q_pow(i, p - 2, p) % p;
    }
    return res;
}

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

int main() {
    scanf("%d", &n);
    while (n --) {
        scanf("%lld%lld%d", &a, &b, &p);
        printf("%d\n", lucas(a, b, p));
    }
}

朴素做法：高精

#include <vector>
#include <iostream>

using namespace std;

const int N = 5010;
int primes[N], cnt;
bool st[N];
int sum[N];
int a, b;
vector<int> res = {1};

void get_primes(int n) {
    for (int i = 2; i <= n; ++ i) {
        if (!st[i]) primes[cnt ++] = i;
        for (int j = 0; primes[j] <= n / i; ++ j) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break; 
        }
    }
}

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

vector<int> mul(vector<int> &a, int b) {
    vector<int> c;
    int t = 0;
    for (int i = 0; i < a.size() || t; ++ i) {
        if (i < a.size()) t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    while (c.size() > 1 && c.back() == 0) c.pop_back();
    return c;
}

int main() {
    scanf("%d%d", &a, &b);
    get_primes(a);
    for (int i = 0; i < cnt; ++ i) {
        int p = primes[i];
        sum[i] = get(a, p) - get(a - b, p) - get(b, p);
    }

    for (int i = 0; i < cnt; ++ i)
        for (int j = 0; j < sum[i]; ++ j)
            res = mul(res, primes[i]);

    for (int i = res.size() - 1; i >= 0; -- i) printf("%d", res[i]);
    puts("");
}

卡特兰数
$$ ktl _ n = \sum ^ {n - 1} _ {i = 0} {ktl _ i * ktl _ {n - 1 - i}} \\\ ktl _ 1 = ktl _ 0 = 1 \\\ ktl _ n = {C ^ n _ {2n} \over {n + 1}} $$

#include <iostream>

using namespace std;

typedef long long LL;

const int MOD = 1e9 + 7;

int q_pow(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = (LL)res * a % MOD;
        a = (LL)a * a % MOD;
        b >>= 1;
    }
    return res;
}

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

    int a = 2 * n, b = n;
    int res = 1;

    for (int i = a; i > a - b; -- i) res = (LL)res * i % MOD;
    for (int i = 1; i <= b; ++ i) res = (LL)res * q_pow(i, MOD - 2) % MOD;
    res = (LL)res * q_pow(n + 1, MOD - 2) % MOD;

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

高斯消元

#include <iostream>
#include <cmath>

using namespace std;

const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];

bool lf_is_same(double a, double b) {
    return fabs(a - b) < eps;
}

int gauss() {
    int c, r; // c--列 r--行
    for (c = 0, r = 0; c < n; ++ c) {
        int t = r; // t--绝对值最大的行
        for (int i = r; i < n; ++ i)
            if (fabs(a[i][c]) > fabs(a[t][c])) t = i;
        // 如果是0就continue
        if (lf_is_same(a[t][c], 0)) continue;
        // 将绝对值最大行放在最上方
        for (int i = c; i < n + 1; ++ i) swap(a[t][i], a[r][i]);
        for (int i = n; i >= c; -- i) a[r][i] /= a[r][c]; // 将改行第一个系数变为1
        // 将下面所有行的第一个系数变为0
        for (int i = r + 1; i < n; ++ i)
            if (!lf_is_same(a[i][c], 0))
                for (int j = n; j >= c; -- j) a[i][j] -= a[r][j] * a[i][c];
        // 行数++
        r ++;
    }
    // 如果行数<n
    if (r < n) {
        for (int i = r; i < n; ++ i)
            if (!lf_is_same(0, a[i][n])) return 2; // 如果出现0=非0则无解
        return 1; // 如果出现0=0
    }

    for (int i = n - 1; i >= 0; -- i)
        for (int j = i + 1; j < n; ++ j) a[i][n] -= a[j][n] * a[i][j]; // 将a[i][n]转换成方程的解
    return 0;
}

int main() {
    scanf("%d", &n);
    for (int i = 0; i < n; ++ i)
        for (int j = 0; j <= n; ++ j)
            scanf("%lf", &a[i][j]);

    int t = gauss();
    if (!t) for (int i = 0; i < n; ++ i) printf("%.2lf\n", a[i][n]);
    else if (t == 1) puts("Infinite group solutions");
    else puts("No solution");
}

容斥原理

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 20;
int n, m;
int p[N];

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 0; i < m; ++ i) scanf("%d", &p[i]);

    int res = 0;
    for (int i = 1; i < 1 << m; ++ i) { // 不能从i=0开始因为一定要选一个集合
        int t = 1, s = 0; // t表示所选质数乘积 s表示所选集合个数
        for (int j = 0; j < m; ++ j)
            if (i >> j & 1) {
                s ++;
                if ((LL)t * p[j] > n) {
                    t = -1; break;
                }
                t *= p[j];
            }
        if (t == -1) continue;
        // 奇数个+ 偶数个-
        if (s & 1) res += n / t;
        else res -= n / t;
    }

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