0%

算法基础（12）

发表于 更新于 分类于 本文字数： 8.1k
高斯消元，求组合数，卡特兰数

高斯消元

高斯消元解线性方程组：线性代数里面很基础的一个知识点，算法步骤：

枚举每一列c：

  1. 找到当前列绝对值最大的一行
  2. 用初等行变换(2) 把这一行换到最上面（未确定阶梯型的行，并不是第一行）
  3. 用初等行变换(1) 将该行的第一个数变成 1（其余所有的数字依次跟着变化）
  4. 用初等行变换(3) 将下面所有行的当前列的值变成 0

高斯消元解线性方程组（Acwing 883）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

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

int n;
double 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 + 1; 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[j][n] * a[i][j];
    
    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++)  printf("%.2lf\n", a[i][n]);
    }
    else if(t == 1)  puts("Infinite group solutions");
    else puts("No solution");
    
    return 0;
}

高斯消元解异或线性方程组（Acwing 884）

思路和用高斯消元解普通的线性方程组相同，只需在相应的步骤换成异或运算即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include<iostream>
#include<algorithm>
using namespace std;

const int N = 110;

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

int gauss()
{
    int c, r;
    for(c = 0, r = 0; c < n; c ++)   //找到当前列中绝对值最大的一行，本题中最大值即是1
    {
        int t = r;
        for(int i = r; i < n; i ++)
            if(a[i][c])  t = i;
        
        if(!a[t][c])  continue;
        
        for(int i = c; i <= n; i ++)  swap(a[r][i], a[t][i]);  //把这一行换到最上面
        for(int i = r + 1; i < n; i ++)
            if(a[i][c])
                for(int j = n; j >= c; j --)
                    a[i][j] ^= a[r][j];
                    
        r ++;
    }
    
    if(r < n)
    {
        for(int i = r; i < n; i ++)
            if(a[i][n])     //0 等于 非0
                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;
}

int main()
{
    cin >> n;
    for(int i = 0; i < n; i ++)
        for(int j = 0; j < n + 1; j ++)
            cin >> a[i][j];
            
    int res = gauss();
    
    if(res == 0)
        for(int i = 0; i < n; i ++)  cout << a[i][n] << endl;
    else if(res == 1)  puts("Multiple sets of solutions");
    else  puts("No solution");
    
    return 0;
}

求组合数

求组合数I（Acwing 885）：

给定n组询问，每组询问给定两个整数$a,b$，请你输出$C^b_a\, mod \, (10^9+7)$的值。

数据范围：$1 \le n \le 10000, 1 \le b \le a \le 2000$

关于组合数的相关知识可以参考：排列组合；组合数学

组合数的公式：

本题中$a$和$b$的数据范围是$1 \le b \le a \le 2000$，需要计算的$C_a^b$总数不超过400万，因此我们可以先处理出所有$C^b_a$的值，它有下面的递推公式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include<iostream>
using namespace std;
using namespace std;

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

int c[N][N];   //c[i][j] 表示 C_i^j，即从i个不同的东西从抽出j个的组合数

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();
    
    int n;
    scanf("%d", &n);
    
    while(n --)
    {
        int a, b;
        scanf("%d%d", &a, &b);
        printf("%d\n", c[a][b]);
    }
    
    return 0;
}

求组合数II（Acwing 886）

给定n组询问，每组询问给定两个整数$a,b$，请你输出$C^b_a\, mod \, (10^9+7)$的值。

数据范围：$1 \le n \le 10000, 1 \le b \le a \le 10^5$

第二题和第一题的区别是$a$和$b$的数据范围是$1 \le b \le a \le 10^5$，没办法将$C_a^b$全预处理完。那我们可以预处理出所有的$a!$，记为fact[i]，又$(a/b) \,mod \,p \ne (a \, mod \, p)/(b \, mod\, p)$，需要将除法转化为乘法，即转化为计算逆元$a/b \equiv a x(mod \, x)$，记为infact[i]，那么$C_a^b$就等于`fact[a] infanct[b-a] * infact[b]`。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include<iostream>
using namespace std;

typedef long long LL;
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;
        a = (LL) a * a % p;
        k >>= 1;
    }
    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;
    scanf("%d", &n);
    
    while(n--)
    {
        int a, b;
        scanf("%d%d", &a, &b);
        printf("%d\n", (LL) fact[a] * infact[b] % mod * infact[a - b] % mod);  //注意这里乘两个数之后就要取一次模，防止溢出LL
    }
    
    return 0;
}

求组合数III（Acwing 887）

给定n组询问，每组询问给定两个整数$a,b,p$，其中$p$是质数，请你输出$C^b_a\, mod \, p$的值。

数据范围：$1 \le n \le 20, 1 \le b \le a \le 10^{18}, 1 \le p \le 10^5$。

第三题的查询数$n$很小，但$a, b$的值都爆大有$10^{18}$，这时可以使用卢卡斯定理：

这时的时间复杂度为$O(\log_p N \cdot p \cdot \log p)=O(p \log N \log p)$。

关于卢卡斯定理的证明可以参考：卢卡斯定理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<iostream>
using namespace std;

typedef long long LL;

int p;

int qmi(int a, int k)   //快速幂，利用快速幂求逆元
{
    int res = 1;
    while(k)
    {
        if( k & 1)  res = (LL) res * a % p;    //计算a! / (a - b)! (mod p) 
        a = (LL) a * a % p;    ////计算 1 / b! (mod p)
        k >>= 1;
    }
    return res;
}

int C(int a, int b)   //计算C_a^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;   //注意这里要用LL 存a, b，数据范围是1到10^18
        cin >> a >> b >> p;
        cout << lucas(a, b) << endl;
    }
    
    return 0;
}

输入$a,b$，求$C^b_a$的值。注意结果可能很大，需要使用高精度计算。

数据范围：$1 \le b \le a \le 5000$

第四题从定义出发，不要求结果取模，而是用高精度表示，我们需要实现高精度乘法和高精度除法，但是直接计算效率比较低，一般是要先将$a, b$分解质因数，这样只需要高精度乘法就可以了。

如何计算$a!$中$p_i$的次数，可以用下面的公式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include<iostream>
#include<vector>
using namespace std;

const int N = 5010;

int primes[N], cnt;   //cnt 存质因数的个数
int sum[N];   //存a分解质因数后，a!中p_i的指数
bool 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] = true;
            if(i % primes[j] == 0) break;
        }
    }
}


int get(int a, int p)   //求a分解质因数后，a!中p的指数
{
    int res = 0;
    while (a)
    {
        res += a / p;
        a /= p;
    }
    return res;
}

vector<int> mul(vector<int> a, int b)    //高精度乘法，一个很大的数，乘一个较小的数
{
    vector<int> c;
    int t = 0;
    for(int i = 0; i < a.size(); i ++)
    {
        t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    
    while (t)
    {
        c.push_back(t % 10);
        t /= 10;
    }
    return c;     //这里不需要处理前导0，因为本题中不会乘0
}


int main()
{
    int a, b;
    cin >> a >> b;
    
    get_primes(a);   //求a的质因数
    
    for(int i = 0; i < cnt; i ++)
    {
        int p = primes[i];
        sum[i] = get(a, p) - get(a - b, p) - get(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 --)  printf("%d", res[i]);
    puts("");
    
    return 0;
}

卡特兰数

满足条件的01序列（Acwing 889）

给定$n$个0和$n$个1，它们将按照某种顺序排成长度为$2n$的序列，求它们能排列成的所有序列中，能够满足任意前缀序列中0的个数都不少于1的个数的序列有多少个。输出的答案对$10^9+7$取模。

数据范围：$1 \le n \le 10^5$。

如$n=3$时，序列可以是：000111，001101， 001011， 010011， 010101。

这个问题可以抽象成另一个问题，即从原点走路径的问题 ，从(0,0)走到(6,6)，把每种序列转化为一种路径，0表示右走一格，1表示向上走一格。这道题目要求任意前缀序列中0的个数都不少于1的个数，对应到路径问题上，即是路径上任意一个位置都要满足$x \ge y$（在红色边下面），也就是任意一条路径不能经过红边。

从(0,0)走到(6,6)的路径一共有$C_{12}^6$种走法（12步取6步向上），要减去所有经过红边的路径，对于每条经过红边的路径，取其与红边相交的第一个点，后面部分的路径对红色边做轴对称，其终点(6,6)一定对称到点(5,7)，因此所有经过红边的路径都可以转化到一条从(0,0)到点(5,7)的路径，那么合法的路径总数为$C_{12}^6-C_{12}^5$，即$C_{2n}^n-C_{2n}^{n-1}$。

这就是著名的卡特兰数，$H_n=C_{2n}^n-C_{2n}^{n-1}$，可参考：卡特兰数。很多问题的方案数都是卡特兰数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include<iostream>
using namespace std;

typedef long long LL;
const int mod = 1e9 +7;

 //用快速幂求逆元，这里的mod是质数；若mod不是质数，只能用扩展欧几里得算法求逆元
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;
}
int main()
{
    int n;
    cin >> n;
    
    int a = 2 * n, b = n;
    int res = 1;
    
    //用下面两个循环计算C_2n^n
    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;  //乘 1/ (n + 1)
     
    cout << res << endl;
    
    return 0;
}