高斯消元法、组合数学、卡特兰数

高斯消元法：

循环下列过程：

1. 找到当前列绝对值最大的一行
2. 将最大行换到最上面去
3. 将最大行第一个数变为 1
4. 将下面所有行对应的数变为 0

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

using namespace std;

const int N = 110;
const double eps = 1e-8;
int n;
double a[N][N];
void out(){
    for(int i = 0; i < n; i ++){
        for(int j = 0; j < n; j ++){
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}
int gauss(){  // 高斯消元，答案存于a[i][n]中，0 <= i < n
    int c, r;
    for (c = 0, r = 0; c < n; c ++ ){
        int t = r;
        for (int i = r; i < n; i ++ )  // 找绝对值最大的行
            if (fabs(a[i][c]) > fabs(a[t][c]))
                t = i;
        if (fabs(a[t][c]) < eps) continue; // 当前列全为0 
        for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);  // 将绝对值最大的行换到最顶端
        for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];  // 将当前行的首位变成1
        for (int i = r + 1; i < n; i ++ )  // 用当前行将下面所有行的当前列消成0
            if (fabs(a[i][c]) > eps)
                for (int j = n; j >= c; j -- ) // 行首是倍数，所以倒着来
                    a[i][j] -= a[r][j] * a[i][c];
        // out();
        r ++ ;
    }

    if (r < n){
        for (int i = r; i < n; i ++ )
            if (fabs(a[i][n]) > eps)
                return 2; // 无解
        return 1; // 有无穷多组解
    }

    for (int i = n - 1; i >= 0; i -- )
        for (int j = i + 1; j < n; j ++ )// 求解第 i 个解需要用到第 j 个解和 a[i][j]
            a[i][n] -= a[i][j] * a[j][n];
    return 0; // 有唯一解
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i ++ )
        for (int j = 0; j < n + 1; j ++ )
            cin >> a[i][j];
    int t = gauss();
    if(t == 0){
        for (int i = 0; i < n; i ++ ){
            if (fabs(a[i][n]) < eps) a[i][n] = 0;  // 去掉输出 -0.00 的情况
            printf("%.2f\n", a[i][n]);
        }
    }
    else if (t == 1) puts("Infinite group solutions");
    else puts("No solution");
    return 0;
}

组合计数：

$C_a^b=C_{a-1}^b+C_{a-1}^{b-1}$

// 预处理出 C(a, b) O(n^2)
#include <iostream>
using namespace std;
const int N = 2010, mod = 1e9 + 7;
int c[N][N];
void init(){
    for(int i = 0; i < N; i ++)
        for(int j = 0; j <= i; j ++)
            if(j == 0) c[i][j] = 1;
			else c[i][j] = (c[i-1][j] + c[i-1][j-1]) % mod;
}
int main(){
    init();
    int n; cin >> n;
    while(n --){
        int a, b; cin >> a >> b;
        cout << c[a][b] << endl;
    }
}

预处理阶乘：$fact(n)=n!\%mod$ ，$infact(n)=(n!)^{-1}\%mod$ ，

$$C_a^b\%m=fact(a)infact(b-a)infact(b)\%m$$

逆元的递推关系：$(n!)^{-1}=(n(n-1)!)^{-1}=n^{-1}((n-1)!)^{-1}$

// 另一种预处理，更大范围 O(n log n)
#include <iostream>
using namespace std;
using ll = long long;
const int N = 100010, mod = 1e9 + 7;
int fact[N], infact[N];
int qmi(int a, int k, int p){
    int res = 1;
    while(k){
        if(k & 1) res = (ll) res * a % p;
        k >>= 1;
        a = (ll)a * a % p;
    }
    return res;
}
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] * qmi(i, mod - 2, mod) % mod;
    }
    int n; cin >> n;
    while(n --){
        int a, b; cin >> a >> b;
        cout << (ll)fact[a] * infact[b] % mod * infact[a - b] % mod << endl;
        // cout << (ll)(fact[a] * infact[b]) % mod * infact[a - b] % mod << endl;
        // 错，不能加括号，会先溢出而不是先类型转换
    }
}

Lucas 定理：

$$C_a^b\equiv C_{a\%p}^{b\%p}\cdot C_{a/p}^{b/p}\mod p$$

// 更大大的范围 
#include <iostream>
#include <algorithm>
using namespace std;
using ll = long long;
int p;
int qmi(int a, int k){
    int res = 1;
    while(k){
        if(k & 1) res = (ll) res * a % p;
        k >>= 1;
        a = (ll)a * a % p;
    }
    return res;
}
int C(int a, int b){
    int res = 1;
    for(int i = 1, j = a; i <= b; i ++, j --){
        res = (ll)res * j % p;
        res = (ll)res * qmi(i, p - 2) % p;
    }
    return res;
}
int lucas(ll a, ll b){
    if(a < p && b < p) return C(a, b);
    return (ll) C(a % p, b % p) * lucas(a / p, b / p) % p;
}
int main(){
    int n; cin >> n;
    while(n --){
        ll a, b; cin >> a >> b >> p;
        cout << lucas(a, b) << endl;
    }
}

$n!$ 中质因子 $p$ 的个数：

$$[\frac{n}{p}]+[\frac{n}{p^2}]+[\frac{n}{p^3}]+...+$$

// 不带模，算精确结果，需要高精度
#pragma GCC optimize(2) // O2 优化
#include <iostream>
#include <vector>
using namespace std;
const int N = 5010;
int primes[N], cnt;
int sum[N]; // sum[i] 表示质因子primes[i]的个数
int st[N];
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] = 1;
            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; //上一位进位
    } // 去掉前导0, 1234 x 0 = 0000
    while(C.size()>1 && C.back()==0) C.pop_back();
    return C;
}
int main(){
    int a, b; cin >> a >> b;
    get_primes(a);
    for(int i = 0; i < cnt; i ++){
        int p = primes[i];
        sum[i] = get(a, p) - get(b, p) - get(a - b, p);
    }
    vector<int> res;
    res.push_back(1);
    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 --) cout << res[i];
}

卡特兰数：

$$\frac{C_{2n}^n}{n+1}$$

#include <iostream>
using namespace std;
using ll = long long;
const int mod = 1e9 + 7;
int qmi(int a, int k, int p){
    int res = 1;
    while(k){
        if(k & 1) res = (ll) res * a % p;
        k >>= 1;
        a = (ll)a * a % p;
    }
    return res;
}
int main(){
    int n; cin >> 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 * qmi(i, mod - 2, mod) % mod;
    res = (ll)res * qmi(n + 1, mod - 2, mod) % mod;
    cout << res;
}

Q.E.D.