[AcWing]，数学知识(三)，推导高斯消元，组合数

https://www.acwing.com/activity/content/punch_the_clock/11/

高斯消元

高斯消元解线性方程

$$\left\{\begin{array}{c}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ \dots \\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n=b_m\end{array}\right\}$$

初等行（列）变换

1. 把某一行乘一个非$0$的数 (方程的两边同时乘上一个非$0$数不改变方程的解)
2. 交换某两行 (交换两个方程的位置)
3. 把某行的若干倍加到另一行上去 （把一个方程的若干倍加到另一个方程上去）

算法步骤

依次枚举每一列$c$

1. 找到当前列中，绝对值最大的一行元素，记录改行为$t$
2. 使用规则2，将$t$行交换到当前的第$r$行
3. 使用规则1，将第$r$行第$c$列的数置为$1$
4. 使用规则3，将底下的其他行中，第$c$列的数置为$0$

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;

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

/*
    return 0;   存在解
    return 1;   有无穷解
    return 2;   没有解
*/
int gauss() {
    int r,c;
    for(r = 0, c = 0; c < n; c++) {
        int t = r;
        // 寻找当前列中，绝对值最大的一行
        for(int i = r + 1; i < n; i++)
            if(fabs(f[i][c]) > fabs(f[t][c]))   t = i;
        
        if(fabs(f[t][c]) < eps) continue;
        // 交换t行 到第r行
        for(int i = c; i <= n; i++) swap(f[t][i],f[r][i]);
        // 更新f[r][c] 为1
        for(int i = n; i >= c; i--) f[r][i] /= f[r][c];
        // 消除当前列 下面的数为 0
        for(int i = r + 1; i < n; i++) {
            if(fabs(f[i][c]) < eps) continue;
            for(int j = n; j >= c; j--)
                f[i][j] -= f[r][j] * f[i][c];
        }
        
        r++;
    }
    
    if(r < n) {
        // 无解的话，后面最后一列的结果不为0
        for(int i = r; i < n; i++) 
            if(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[j][n] * f[i][j];
    
    return 0;
}

int main() {
    cin >> n;
    for(int i = 0; i < n; i++)
        for(int j = 0; j <= n; j++)
            cin >> f[i][j];
    
    int ret = gauss();
    
    if(ret == 0) {
        for(int i = 0; i < n; i++) printf("%.2lf\n",f[i][n]);
    } else if(ret == 1) {
        puts("Infinite group solutions");
    } else {
        puts("No solution");
    }
    return 0;
}

组合数

组合数 I

$$c[i][j]=c[i-1][j]+c[i-1][j-1]$$
就相当于，要求是在$i$个物品中选$j$个物品，那么他的选法有两种：

1. 选择第$i$个物品，那么剩余的选法就是$c[i-1][j-1]$
2. 不选第$i$个物品，那么剩余的选法就是$c[i-1][j]$

#include <iostream>
using namespace std;
const int N = 2010,mod = 1e9 + 7;
int c[N][N];

int main() {
    int n;
    cin >> n;
    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;
    
    while(n --) {
        int a,b;
        cin >> a >> b;
        cout << c[a][b] << endl;
    }
    return 0;
}

求组合数 II

$${\mathrm{C}}_a^b=\frac{a!}{(a-b)!\times b!}$$
从数据范围来看，我们无法直接开辟$100000\times 100000$大小的二维空间，所以需要对原本的公式进行压缩，统计每一个数的阶乘

但是，$\frac{a(\mathrm{m}\mathrm{o}\mathrm{d}n)}{b(\mathrm{m}\mathrm{o}\mathrm{d}n)}\ne \frac{a}{b}(\mathrm{m}\mathrm{o}\mathrm{d}n)$，所以说，需要对分母中的每个数求一个逆元，又因为$mod=10^9+7$是一个质数，所以可以使用快速幂来求逆元，将除法变为乘法

#include <iostream>
using namespace std;

typedef long long LL;
const int N = 1e5 + 10, 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() {
    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] * qmi(i,mod-2,mod) % mod;
    }
    while( n -- ) {
        int a,b;
        cin >> a >> b;
        cout << (LL)fact[a] * infact[a-b] % mod * infact[b] % mod << endl;
    }
    return 0;
}

组合数 III

$Lucas$定理

$${\mathrm{C}}_a^b\equiv {\mathrm{C}}_{a\%p}^{b\%p}\times {\mathrm{C}}_{\frac{a}{p}}^{\frac{b}{p}}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

#include <iostream>
using namespace std;

typedef long long LL;
int n,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() {
    cin >> n;
    while(n --) {
        LL a,b;
        cin >> a >> b >> p;
        cout << lucas(a,b) << endl;
    }
    return 0;
}

求组合数 IV

$${\mathrm{C}}_a^b=\frac{a!}{b!\times (a-b)!}$$
$$\Rightarrow =\frac{a\times (a-1)\times \cdots \times (a-b+1)}{b\times (b-1)\times \cdots \times 1}$$
$$\Rightarrow =p_1^{{\alpha }_1}\times p_2^{{\alpha }_2}\times \cdots \times p_k^{{\alpha }_k}$$
其中，$p_1{p}_2\cdots p_k$为质数，${\alpha }_1{\alpha }_2\cdots {\alpha }_k$为对应的质数的出现次数

#include <iostream>
#include <vector>
using namespace std;

const int N = 5010;
int sum[N];
int primes[N],cnt;
bool st[N];
int a,b;

void get_prime(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;
        }
    }
}

// n! 中，p出现的次数
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> res;
    int t = 0, i = 0;
    while(i < a.size() || t) {
        if(i < a.size()) t += a[i++] * b;
        res.push_back(t%10);
        t/=10;
    }
    return res;
}

int main() {
    cin >> a >> b;
    get_prime(a);   // 筛选出所有比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];
    return 0;
}