[数论第三节]高斯消元法/求组合数/卡特兰数

手机

顺络电子：董事长部分股权办理股票质押业务

深圳7月二手住宅成交2259套，中介称近期咨询客户开始增加

家电

[数论第三节]高斯消元法/求组合数/卡特兰数

2023-08-10 22:59:05 来源：博客园

(资料图)

高斯消元

求解含有n个未知数，n个方程的多元线性方程组 O(n^3)
初等行变换：
- 某行乘以一个非零数
- 交换两行
- 某行加上另一行的若干倍
利用初等行变换将方程组化为上三角矩阵
解的情况：
- 完美阶梯型：唯一解
- 非完美阶梯型：
  - 0 == 非0：无解
  - 0 == 0：无穷解
步骤:
- 枚举每一列
  - 找到这一列系数的绝对值最大的一行
  - 将这一行与第一行交换
  - 将改行的第一个数变成一(方程两边同乘某数)
  - 把下面所有行的当前列的系数消成0(某行加上第一行的若干倍)

代码：

const int N = 110;const double esp = 1e-6; //x < esp,则x = 0,否则 x != 0,由于浮点数精度问题,0可能是0.000001double a[N][N];int n;int gauss(){int r, c;for(r = 1, c = 1; c <= n; ++ c){//遍历每一列int t = r;for(int i = r; i <= n; ++ i) //找当前列的系数最大的行if(fabs(a[t][c]) < fabs(a[i][c]))t = i; //记录最大行if(fabs(a[t][c]) < esp) continue; //最大行系数为0说明该列系数均为0for(int i = c; i <= n + 1; ++ i) swap(a[r][i], a[t][i]);//否则交换第一行与系数最大行for(int i = n + 1; i >= c; -- i) a[r][i] /= a[r][c]; //将第一行当前列系数变为1for(int i = r + 1; i <= n; ++ i) //将下面所有行的当前列的系数变为0if(fabs(a[i][c]) > esp) //当前列系数非0for(int j = n + 1; 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 + 1]) > esp) return 0; //无解return 1; //无数解}for(int i = n; i >= 1; -- i) //从后往前求解for(int j = i + 1; j <= n; ++ j)a[i][n + 1] -= a[i][j] * a[j][n + 1];return 2; //唯一解}int main(){cin >> n;for(int i = 1; i <= n; ++ i)for(int j = 1; j <= n + 1; ++ j)cin >> a[i][j];int t = gauss();if(t == 1) cout << "Infinite group solutions";else if(t == 0) cout << "No solution";else if(t == 2){for(int i = 1; i <= n; ++ i) printf("%.2f\n", a[i][n + 1]);}return 0;}

求组合数

\(C_a^b = \frac{a(a-1)(a-2)···(a-b+1)}{b!}=\frac{a!}{b!(a-b)!}\)
一般直接用公式求组合数容易超时
故可以通过预处理的方式，用递推式求出所有可能的组合数
递推式：\(C_a^b = C_{a-1}^{b-1} + C_{a-1}^{b}\) \(O(N^2)\)
当a，b <= 2e3，询问n = 1e4次时

代码：

//询问10000次//a,b <= 2000const int N = 2e3 + 10;const int mod = 1e9 + 7;int c[N][N];int n;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 - 1] + c[i - 1][j]) % mod;}

当a，b <= 1e5，询问n = 1e4次，模数p = 1e9 + 7时

开二维数组会爆内存，于是利用公式，预处理出阶乘以及阶乘的逆\(O(Nlog(N))\)

代码：

typedef long long LL;const int N = 1e5 + 10;const int mod = 1e9 + 7;int fact[N], infact[N];int n;//快速幂int qmi(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;}//预处理各阶乘以及阶乘的逆void init(){    fact[0] = infact[0] = 1;    for(int i = 1; i < N; ++ i){        fact[i] = (LL)fact[i - 1] * i % mod; //递推式n! = (n - 1)! * n;        infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod) % mod; //递推式(n!)^-1 = ((n - 1)!)^-1 * n^-1，取模意义下的逆就是逆元    }}int main(){    cin >> n;    init();    while(n -- ){        int a, b;        cin >> a >> b;        cout << (LL)fact[a] * infact[b] % mod * infact[a - b] % mod << endl;//两个int最大值相乘不会爆longlong，但是三个相乘会爆    }    return 0;}

当a，b <= 1e18，询问n = 2e5次，模数p不定时

两数相乘会爆longlong，利用卢卡斯定理预处理组合数 \(O(plog(N)log(p))\)
定理：\(C_a^b\equiv C_{a \mod p}^{b \mod p}·C_{a/p}^{b/p} \pmod p\)

代码：

typedef long long LL;int p;int n;//快速幂求逆元int qmi(int a, int b){    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){ //求组合数C_a^b    if(a < b) return 0; //可以不要    int x = 1, y = 1; //分子分母    for(int i = 0; i < b; ++ i){        x = (LL)x * (a - i) % p; //求分子        y = (LL)y * (i + 1) % p; //求分母    }    return (LL)x * qmi(y, p - 2) % p;}//卢卡斯定理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;}

当a，b <= 5e5，询问一次，不取模时

此时组合数较大，需要用高精度存储结果
阶乘中素数的个数：
- \(v_p(n!) = \sum_{i=1}^{\infty} \left\lfloor \frac{n}{p^i} \right\rfloor\)
- 证明：将 \(n!\) 记为 \(1\times 2\times \cdots \times p\times \cdots \times 2p\times \cdots \times \lfloor n/p\rfloor p\times \cdots \times n\) 那么其中 \(p\) 的倍数有 \(p\times 2p\times \cdots \times \lfloor n/p\rfloor p=p^{\lfloor n/p\rfloor }\lfloor n/p\rfloor !\) 然后在 \(\lfloor n/p\rfloor !\) 中继续寻找 \(p\) 的倍数即可，这是一个递归的过程

代码：

const int N = 5e3 + 10;int primes[N], cnt;int st[N];int sum[N];//线性筛，1-n中质数就是n!中的质因子void get_primes(int a){    for(int i = 2; i <= a; ++ i){        if(!st[i]) primes[ ++ cnt] = i;        for(int j = 1; primes[j] <= a / i; ++ j){            st[primes[j] * i] = 1;            if(i % primes[j] == 0) break;        }    }}//求a!中质数p的个数int get(int a, int p){    int res = 0;    while(a){        res += a / p;        a /= p;    }    return res;}//高精度乘法vector mul(vector &A, int b){    vector 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(){int a, b;    cin >> a >> b;    get_primes(a);    for(int i = 1; i <= cnt; ++ i){        int p = primes[i];        sum[i] = get(a, p) - get(b, p) - get(a - b, p); //a!中质因子个数-b!中质因子个数-(a-b)!中质因子个数    }        vector res;    res.push_back(1);        for(int i = 1; i <= cnt; ++ i) //累乘        for(int j = 1; j <= sum[i]; ++ j)            res = mul(res, primes[i]);        for(int i = res.size() - 1; i >= 0; -- i)//输出        cout << res[i];        return 0;}

卡特兰数

\(C_{2n}^n-C_{2n}^{n-1}=\frac{1}{n+1}·C_{2n}^n\)
用于解决序列排列等问题
AcWing 889. 满足条件的01序列 - AcWing

代码：

#include using namespace std;typedef long long LL;const int N = 1e5 + 10;const int p = 1e9 + 7;int n;int qmi(int a, int k){    int res = 1;    while(k){        if(k & 1) res = (LL)res * a % p;        a = (LL)a * a % p;        k >>= 1;    }    return res;}int main(){    cin >> n;        int res = 1;        int a = 2 * n, b = n;    int x = 1, y = 1;    for(int i = 0; i < b; ++ i){        x = (LL)x * (a - i) % p;        y = (LL)y * (i + 1) % p;    }        res = (LL)x * qmi(y, p - 2) % p;    res = (LL)res * qmi(n + 1, p - 2) % p;    cout << res;        return 0;}

关键词：