数论-2

中国剩余定理

则求

得到这个数为

c++解法

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
#pragma GCC diagnostic error "-std=c++11"
#include <bits/stdc++.h>
#define debug(x) cout << #x << " = " << x << endl

using namespace std;

typedef pair<int, int> PII;
typedef long long LL;

const int N = 1e6 + 5;
//const double eps = 1e-6;
//int dx[] = {1, 0, -1, 0}, dy[] = {0, -1, 0, 1};
// int dx[8] = {-1, -1, -1, 0, 1, 1, 1, 0};
// int dy[8] = {-1, 0, 1, 1, 1, 0, -1, -1};

LL exgcd(LL a, LL b, LL &x, LL &y)
{
    if (b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    LL  d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

LL inline mod(LL a, LL b)
{
    return ((a % b) + b) % b;
}

int main()
{
    int n;
    cin >> n;
    LL a1, m1;
    cin >> a1 >> m1;
    for (int i = 1; i < n; i ++ )
    {
        LL a2, m2, k1, k2;
        cin >> a2 >> m2;
        LL d = exgcd(a1, -a2, k1, k2);
        if((m2 - m1) % d)
        {
            puts("-1");
            return 0;
        }
        k1 = mod(k1 * (m2 - m1) / d, abs(a2 / d));
        m1 = k1 * a1 + m1;
        a1 = abs(a1 / d * a2);
    }
    cout << m1;
    return 0;
}

python解法

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
def exgcd(a,b,x,y):
    if not b:
        x,y = 1,0
        return a,x,y

    d,x,y = exgcd(b, a % b, y, x)
    temp = y
    y = x - a // b * y
    x = temp
    return d,x,y

if __name__ == "__main__":
    n = int(input())
    a1,m1 = map(int,input().split())

    for i in range(n-1):
        a2,m2 = map(int,input().split())

        k1,k2 = 0,0

        d,k1,k2 = exgcd(a1,a2,k1,k2)

        if (m2-m1) % d:
            print(-1)
            exit()

        k1 *= (m2-m1) // d
        t = a2 // d
        k1 = k1 % t
        m1 = a1 * k1 + m1
        a1 = abs(a1 // d * a2)

    print(m1 % a1)

求组合数

Ⅰ求 (递推)

a, b范围是2000，询问次

通过递推进行处理

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
#include <iostream>
#include <algorithm>

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

Ⅱ求 (预处理逆元)

a, 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
#include <iostream>

using namespace std;

const int N = 1e5 + 10, mod = 1e9 + 7;
typedef long long LL;

int fact[N], infact[N];  //  预处理的阶乘及逆元

int q_pow(int a, int k, int mod)
{
    int res = 1;
    while(k)
    {
        if(k&1)  res = (LL)res * a % mod;
        a = (LL)a * a % mod;
        k >>= 1;
    }
    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] * q_pow(i, mod - 2, mod) % mod;
    }
    while (n -- )
    {
        int a, b;
        cin >> a >> b;
        cout << (LL)fact[a] * infact[b] % mod * infact[a - b] % mod << endl;
    }

    return 0;
}

Ⅲ 求 (Lucas定理)

a, b范围是 ，询问次

利用Lucas定理 Lucas 定理用于求解大组合数取模的问题，其中模数必须为素数

其中的求法推导如下 

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
#pragma GCC diagnostic error "-std=c++11"
#include <bits/stdc++.h>
#define debug(x) cout << #x << " = " << x << endl

using namespace std;

typedef pair<int, int> PII;
typedef long long LL;

const int N = 1e6 + 5, mod = 1e9 + 7;
//const double eps = 1e-6;
//int dx[] = {1, 0, -1, 0}, dy[] = {0, -1, 0, 1};
// int dx[8] = {-1, -1, -1, 0, 1, 1, 1, 0};
// int dy[8] = {-1, 0, 1, 1, 1, 0, -1, -1};

int q_pow(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 C(int a, int b, int p)
{
    if(b > a)  return 0;
    int res = 1;
    for (int i = 1, j = a; i <= b; i ++ ,j -- )
    {
        res = (LL)res * j % p;
        res = (LL)res * q_pow(i, p - 2, p) % p;
    }
    return res;
}

int lucas(LL a, LL b, int p)
{
    if (a < p && b < p)  return C(a, b, p);
    else return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}

int main()
{
    int n;
    cin >> n;
    while(n -- )
    {
        LL a, b; 
        int p;
        cin >> a >> b >> p;
        cout << lucas(a, b, p) << endl;
    }
    return 0;
}

Ⅳ 求(高精度)

a, b的范围是5000，求1次，但是要求计算出来所以要用高精度

c++版本

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
#include <iostream>
#include <cstring>
#include <algorithm>

using namespace std;

const int N = 5010;

int primes[N], cnt;
int sum[N];
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 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(); i ++ )
    {
        t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    while(t)
    {
        c.push_back(t % 10);
        t /= 10;
    }
    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(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]);
    
    
    return 0;
}

python版本

利用定义，但是跑起来比c++略慢

1
2
3
4
5
6
7
8
a, b = map(int, input().split())
up = 1
down = 1
for i in range(a - b + 1, a + 1):
    up *= i
for i in range(1, b + 1):
    down *= i
print(up // down)

卡特兰数

求满足条件的01序列

给定n个1，n个0，求它们能排列成的所有序列中，能够满足任意前缀序列中0的个数都不少于1的个数的序列有多少个

下图中，表示从走到的路径，在绿线及以下表示合法，若触碰红线即不合法

结果： 求组合数和逆元即可，注意除以n+1时也要用逆元

img 
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
#include <iostream>
#include <cstring>
#include <algorithm>

using namespace std;

typedef long long LL;

const int mod = 1e9 + 7;

LL q_pow(LL a, LL k)
{
    LL ret = 1;
    while(k)
    {
        if(k&1) ret = ret * a % mod;
        a = a * a % mod;
        k >>= 1;
    }
    return ret;
}

int main()
{
    int n;
    cin >> n;
    LL ans = 1;
    // z
    for(int i = 2 * n; i > n; i -- ) ans = ans * i % mod;
    for(int i = n; i > 0; i -- )  ans = ans * q_pow(i, mod - 2) % mod;
    
    ans = ans * q_pow(n + 1, mod - 2) % mod;
    
    cout << ans << endl;
    return 0;
}

斯特林数

错排问题

n封信，n个信封，求完全放错的方案数?

记方案数为f[n],可以递推来求

f[1] = 0, f[2] = 1

f[i] = (i-1) * (f[i-1] + f[i-2])

推导逻辑：

假设第一封信a占据了b的位置，那么此时b放在哪个信封分两种情况，b放在a位置，或b不放在a位置；

  • 第一类：第一种情况是放在a位置，此时b放在a位置，剩下n-2封信进行错排，方案数是f(n-2)

  • 第二类：第二种情况是b没有去a的位置，那么b可能出现在除α之外的任何位置，b有n-2个位置可以去，不能去a,b位置，其余所有元素都有n-2个位置可以去（a，b位置不能去），这种情况下相当于除a之外的其它元素的错排问题，即n-1个元素的错排问题，方案数是f(n-1)

  • 加法法则：汇总上述分类计数原理，使用加法法则，计算结果是f(n-1)+f(n-2)

  • 乘法法则：对于每一封信都如此，当有n-1封信已经错排时，已经完全错排了

    即f[i] = (i-1) * (f[i - 1] + f[i - 2])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#pragma GCC diagnostic error "-std=c++11"
#include <algorithm>
#include <iostream>
#include <cstring>

using namespace std;

typedef long long LL;
const int N = 1e6 + 5, mod = 1e9 + 7;

LL f[N];

int main()
{
    int n;
    cin >> n;
    f[1] = 0;
    f[2] = 1;
    for (int i = 3; i <= n; i ++ )
        f[i] = (i - 1) * (f[i - 1] + f[i - 2]);
    cout << f[n] << endl;
    return 0;
}

高斯消元

解线性方程组

高斯消元是一种可以可以将增广矩阵化成阶梯型矩阵的算法

适用于求解包含n个方程，n个未知数的多元线性方程组

算法步骤：

  • 找到当前列绝对值最大的一行
  • 将这一行移到最上面
  • 乘以相应的系数将改行的第一个数变成1
  • 相减相消让下面所有的行的当前列变成0

  • 最后再把阶梯型矩阵从下到上回代到第一层即可得到方程的解

复杂度

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
#pragma GCC diagnostic error "-std=c++11"
#include <bits/stdc++.h>
#define debug(x) cout << #x << " = " << x << endl

using namespace std;

typedef pair<int, int> PII;
typedef long long LL;

const int N = 1e3 + 5, mod = 1e9 + 7;
const double eps = 1e-6;  // 用eps来判断是否为0
//int dx[] = {1, 0, -1, 0}, dy[] = {0, -1, 0, 1};

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

int  gauss()
{
    int c, r;  // c代表列col, r代表行row
    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，那么这一行全是0

        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，方程两边同时除以第一个系数
        // 必须要从后往前，不然第一个数会直接变为1，然后不修改其他系数
        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];
                    //  第四步，把第r+1行~n行的第c列元素都消除为0
                    //  用最上面一行和这行做差除去这一行
       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 ++ )
        if(fabs(a[i][n]) < eps) puts("0.00");
        else printf("%.2lf\n", a[i][n]);
    }
    else if (t == 1)  puts("Infinite group solutions");
    else puts("No solution");
    return 0;
}

高斯消元解异或方程组

异或=>不进位加法

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
#pragma GCC diagnostic error "-std=c++11"
#include <algorithm>
#include <iostream>
#include <cstring>

using namespace std;

typedef pair<int, int> PII;
typedef long long LL;

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

int n;
int 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(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])
                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 t = gauss();

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