高斯消元法，组合数

1 高斯消元法解线性方程组

在O($n^3$)时间内解决n元线性方程组
解有三种情况：无解，无穷解，唯一解
初等行列变换: 转换成最简阶梯形矩阵

1. 把某一行乘一个非零的数
2. 交换某两行
3. 把某行的若干被加到另一行上

• 完美阶梯形 ——唯一解
• 0 = 非零 ——无解
• 0 = 0 的方程——无穷解

算法步骤：

1. 枚举每一列c，找到绝对值最大的一行
2. 将该行换到最上面去
3. 将这一行的一个数变成1
4. 把下面所有行的当前列消成零

模板

// a[N][N]是增广矩阵
int gauss()
{
    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;

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

        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 ++ )
            a[i][n] -= a[i][j] * a[j][n];

    return 0; // 有唯一解
}

例题

AcWing 883. 高斯消元解线性方程组

/*
## 高斯消元解线性方程组
学过线性代数就会知道运用高斯消元法就是将行列式矩阵转换成行阶梯型矩阵
算法步骤：
1. 首先从列开始枚举（线性代数中也是从前往后化简）
2. 找到当前列上绝对值最大的数（负数可以将这一行乘-1转化成正的所以要是绝对值）
3. 将这一行移动打最上面，并将第一个数化为1
4. 将其他行乘一个数使得这一列为0
`要注意的是每次循环时初始的列和行不同，相当于呈|——右下移动`
*/
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;

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

double a[N][N];
int n;

// 返回1有唯一解，2有无穷多解，0无解
int gauss()
{
    int r, c;//由于每次循环都不会用到上一层换到上面的行了和列所以行和列的值不能改变
    for (r = 0, c = 0; c < n; c ++)//在做行阶梯转变时，只需要将系数化为阶梯形即可，b[i] 不进行，这就是为什么有n + 1列确只循环到n列
    {
        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;

        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，要从后往前边，不然第一个数变了后面就变不了了

        //开始将下面的这列都化为0，后面的数也要变,同理要从后向前变
        for (int i = r + 1; i < n; i ++)
            if (fabs(a[i][c]) > eps)//当前行已经为零就不能变了，不然结果错误
                for (int j = n; j >= c; j --)
                    a[i][j] -= a[r][j] * a[i][c];//当前行的第一个数是a[i][c]
        r ++;
    }

    if (r < n)//方程个数小于未知量的个数，是无解或者无穷解
    {
        for (int i = r; i < n; i ++)
            if (fabs(a[i][n]) > eps)
                return 0;
        return 2;
    }
    //根据阶梯矩阵的性质：第i个位置量的答案就是bi
    for (int i = n - 1; i >= 0; i --)//行从下到上
        for (int j = i + 1; j < n; j ++)//下面行的位置量的答案已经知道利用消解
            a[i][n] -= a[i][j] * a[j][n];// 这个未知量的答案*参数

    return 1;
}

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 >> a[i][j];//如果用scanf读入是有用double类型 scanf("%lf", &a[i][j]);

    int flag = gauss();

    if (flag == 0) printf("No solution");
    else if (flag == 2) printf("Infinite group solutions");
    else{
        for (int i = 0; i < n; i ++)
            printf("%.2lf\n", a[i][n]);
    }
    return 0;
}

2 组合数

如何选相应的算法

• 有10万组询问时 1 <= b <= a <= 2000这个时候用了递推
• 有1万组询问时 1 <= b <= a <= $10^5$时就要预处理了
• 有20组询问时 1<= b <= a <= $10^{18}$,1 <= p <= $10^5$

递推法求组合数

$C\left(\genfrac{}{}{0ex}{}{b}{a}\right)$ = $C\left(\genfrac{}{}{0ex}{}{b}{a-1}\right)$ + $C\left(\genfrac{}{}{0ex}{}{b-1}{a-1}\right)$

模板

时间复杂度：O(n^2)

// c[a][b] 表示从a个苹果中选b个的方案数
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;

例题

AcWing 885. 求组合数 I

在组合的性质中有一条：
$$C\left(\genfrac{}{}{0ex}{}{b}{a}\right)=C\left(\genfrac{}{}{0ex}{}{b-1}{a-1}\right)+C\left(\genfrac{}{}{0ex}{}{b}{a-1}\right)$$
利用这条性质递推处理整个数组

#include <iostream>
#include <cstring>
#include <algorithm>
#define int long long

using namespace std;

const int mod = 1e9 + 7;
const int N = 2010;

int n;
int c[N][N];

signed main()
{
    for (int i = 0; i <= N - 10; i ++)
        for (int j = 0; j <= i; j ++)
        {
            if (!j) c[i][j] = 1;
            else c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
        }
    
    scanf("%d", &n);
    while (n -- )
    {
        int a, b;
        scanf("%d%d", &a, &b);
        
        printf("%lld\n", c[a][b]);
    }
    return 0;
}

通过预处理逆元的方式求组合数

由于$$\frac{a}{b}modm\ne \frac{amodm}{bmodm}$$
所以要使用逆元预处理出b的逆元在进行乘法运算防止溢出
$$a\ast b^{-1}modm=(amodm\ast b^{-1}modm)modm$$

模板

时间复杂度：O(n logn)

首先预处理出所有阶乘取模的余数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;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}

// 预处理阶乘的余数和阶乘逆元的余数
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;
}

例题

AcWing 886. 求组合数 II

利用的是组合公式
$$\left(\genfrac{}{}{0ex}{}{b}{a}\right)=\frac{a!}{(a-b)!\ast b!}$$
但是将%mod加在每个阶乘上就是不正确的这种拆%的方式不适合除法，所以我们要转换成乘法
a / b = a * x
这里x 就是b的逆元
在费马小定理中mod是质数的话b的逆元就是$b^{mod-2}$

// 不能通过#define int long long 直接简单省事
// 中间过程要加上(LL) 来处理不然最后会W/A
#include <iostream>
#include <cstring>
#include <algorithm>

using namespace std;
const int mod = 1e9 + 7;
const int N = 100010;
typedef long long LL;

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

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

LL C(int a, int b)
{
    int f_a = 1, f_b = 1, f_ab = 1;
    //在每个处理步骤上都要加(LL)否则会W/A
    for (int i = 1; i <= a; i ++) f_a = (LL)f_a * i % mod;
    for (int i = 1; i <= b; i ++) f_b = (LL)f_b * i % mod;
    for (int i = 1; i <= a - b; i ++) f_ab = (LL)f_ab * i % mod;
    //这种方法由于没有进行预处理，在处理多个数据是会大量重复计算导致TLE
    
    return (LL) f_a * qmi(f_b, mod - 2) % mod * qmi(f_ab, mod - 2) % mod;
}

signed main()
{
    scanf("%d", &n);
    
    fact[0] = infact[0] = 1;
    
    for (int i = 1; i < N; i ++) 
    {
        fact[i] = (LL)fact[i - 1] * i % mod;
        infact[i] = qmi(fact[i], mod - 2);//事实证明infact先阶乘后取逆元仍正确
    }
    while (n -- )
    {
        int a, b;
        scanf("%d%d", &a, &b);
        
        printf("%lld\n", (LL)fact[a] * infact[b] % mod * infact[a - b] % mod);
        // printf("%lld\n", C(a, b));
    }

    return 0;
}

Lucas定理

$$C\left(\genfrac{}{}{0ex}{}{b}{a}\right)\equiv C\left(\genfrac{}{}{0ex}{}{bmodp}{amodp}\right)\ast C\left(\genfrac{}{}{0ex}{}{b/p}{a/p}\right)(modp)$$

模板

时间复杂度:O($p\ast logn\ast logp$)

int C(int a, int b, int p)  // 通过定理求组合数C(a, b)
{
    if (a < b) return 0;

    LL x = 1, y = 1;  // x是分子，y是分母
    for (int i = a, j = 1; j <= b; i --, j ++ )
    {
        x = (LL)x * i % p;
        y = (LL) y * j % p;
    }

    return x * (LL)qmi(y, p - 2, p) % p;
}

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

例题

AcWing 887. 求组合数 III

特点a, b的范围非常的大必须用long long 来存用lucas定理
利用lucas定理
$$C\left(\genfrac{}{}{0ex}{}{b}{a}\right)=C\left(\genfrac{}{}{0ex}{}{b/p}{a/p}\right)\ast C\left(\genfrac{}{}{0ex}{}{b\%p}{a\%p}\right)\%p$$
其中求普通的组合数用定义法：
$$C\left(\genfrac{}{}{0ex}{}{b}{a}\right)=\frac{a\ast a-1\ast a-2\ast \dots \dots \ast a-b+1}{b\ast b-1\ast b-2\ast \dots \dots \ast 1}$$

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

using namespace std;
typedef long long LL;

int n;

int qmi(int a, int b, int p)
{
    LL 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 (a < b) return 0;
    
    LL x = 1, y = 1;
    for (int i = a, j = 1; j <= b; j ++, i --)//分子是到a - b + 1共乘b次
    {
        x = (LL) x * i % p;
        y = (LL) y * j % p;
    }
    return (LL)x * qmi(y, p - 2, p) % p;
}

int lucas(LL a, LL b, int p)
{
    if (a < p && b < p) return C(a, b, p); // 当a, b都在p的范围内时就可以返回了，不然没有递归返回条件MLE

    return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;// %p肯定在p的范围内所以直接用C减少递归，如果/的话就不一定
}

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

高精度表示组合数的值

AcWing 888. 求组合数 IV

一个组合数可能很大，爆long long 但是我们不用%成小数，而是直接求出这个数

首先我们知道任意一个数都可以表示成若干质数相乘的形式
所以对于$$C\left(\genfrac{}{}{0ex}{}{b}{a}\right)=\frac{a!}{b!\ast (a-b)!}$$
将这个分数阶乘分解成若干质因数相乘的形式

所以这里用到了一个阶乘分解的知识点 阶乘分解
假设阶乘分解其中一个质因数p, 在这三个阶乘中的个数是s_a, s_b, s_ab
则p实际的次方是a_a - s_b - s_ab
最后利用高精度将他们乘起来即可

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

using namespace std;

const int N = 5010;

int a, b;
int primes[N], cnt = 0;
bool st[N];
int sum_p[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] = true;
            if (i % primes[j] == 0) break;
        } 
    }
}

int get(int num, int p)
{
    int s = 0, t = num;
    while (t) s += t / p, t /= p;
    return s;
}

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.back() == 0 && C.size() > 1) C.pop_back();
    return C;
}

int main()
{
    scanf("%d%d", &a, &b);
    
    get_primes(a);
    
    for (int i = 0; i < cnt; i ++)
        sum_p[i] = get(a, primes[i]) - get(b, primes[i]) - get(a - b, primes[i]);
    
    vector<int> res;
    res.push_back(1);//res就是高精度乘法中的大数，所以要给它一个1
    
    for (int i = 0; i < cnt; i ++)
        for (int j = 0; j < sum_p[i]; j ++)
            res = mul(res, primes[i]);
    
    for (int i = res.size() - 1; i >= 0; i --)
        printf("%d", res[i]);
    
    return 0;

}

参考文献

数学知识（三）