前言

在考试中，经常会出现一类问题，它们不涉及很深的算法，却和数学息息相关。这样的问题通常难度不大，也不需要特别的数学知识，只要掌握简单的数理逻辑即可。——《算法笔记》胡凡

质数(素数)

质数又称素数，是指除了1和本身外，不能被其它的任何数整除。

质数的判定

如何判断给定的正整数n是否是质数？

时间复杂度：O(sqrt(n))O(sqrt(n))

1
2
3
4
5
6
7
8
bool isPrime(int n) {
    if (n < 2) return false;
    
    for (int i = 2; i <= n / i; i++) // i <= n / i
        if (n % i == 0) return false;
    
    return true;
}

分解质因数

n中有且仅有一个大于sqrt(n)的质因子。（如果有两个，两者相乘将大于n，明显不符合）

据此，我们可以先枚举小于sqrt(n)的所有质因子

如果最后n仍大于1，那么这个n就是那个大于sqrt(n)的质因子

最坏时间复杂度：O(sqrt(n))O(sqrt(n))

最好时间复杂度：O(log(n))O(log(n))，（若n为2k2^k，那么在i=2时，所有因子就被分解出来了，所以时间应该是log(2k)=klog(2^k)=k）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
void divide(int n) {
    
    for (int i = 2; i <= n / i; i++) 
        if (n % i == 0) {
            int s = 0;
            while (n % i == 0) {
                s++;
                n /= i;
            }
            printf("%d %d\n", i, s);
        }
    
    if (n > 1) printf("%d %d", n, 1);
}

筛质数

筛选出从1-n之间质数的个数

数据量小于等于10610^6时，以下两个算法效率不相上下

埃氏筛法思想

一个数是p，如果2~p-1中的质数都不是它的质因数的话，那么p是质数(以此对算法进行优化)

朴素筛法

O(nlog(n))O(nlog(n))

时间复杂度计算，当n=2时，执行n/2次；当n=3时，执行n/3次

n2+n3+n4+...+nn\frac{n}{2}+\frac{n}{3}+\frac{n}{4}+...+\frac{n}{n}

=n∗(12+13+14+...+1n)=n*(\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n})

=n∗lim⁡x→+∞12+13+14+...+1n=n*\lim_{x \to +\infty}{\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n}}

=n∗ln(n)<n∗log(n)=nlog(n)=n*ln(n)<n*log(n)=nlog(n)

1
2
3
4
5
6
7
8
9
10
11
12
int primes[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 = i + i; j <= n; j += i) st[j] = true;
    }
}

cout << cnt << endl; // 输出质数个数

优化

O(nlog(log(n)))O(nlog(log(n)))。相当于线性的时间复杂度

1~n中有nln(n)\frac{n}{ln(n)}个质数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int primes[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 = i + i; j <= n; j += i) st[j] = true;
        }
    }
}

cout << cnt << endl; // 输出质数个数

线性筛法(*常用)

n只会被最小质因子筛掉

  1. i % primes[j] == 0

    primes[j]primes[j]一定是ii的最小质因子，primes[j]primes[j]一定是primes[j]∗iprimes[j]*i的最小质因子

  2. i % primes[j] != 0

    primes[j]primes[j]一定小于ii的所有质因子，primes[j]primes[j]也一定是primes[j]∗iprimes[j]*i的最小质因子

1
2
3
4
5
6
7
8
9
10
11
12
int primes[N]; // 用于存储质数
bool st[N]; // false时为质数

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; // primes[j]一定是i的最小质因子
        }
    }
}

约数

12的约数为1,12,3,4,2,6

试除法求约数

O(sqrt(n))O(sqrt(n))

1
2
3
4
5
6
7
8
9
10
11
12
vector<int> get_divisors(int n) {
    vector<int> res;
    
    for (int i = 1; i <= n / i; i++) // 只枚举约数里最小的
        if (n % i == 0) {
            res.push_back(i);
            if (i != n / i) res.push_back(n / i); // 避免出现5*5的情况
        }
        
    sort(res.begin(), res.end());
    return res;
}

约数个数

若要求出N的约数个数，首先将N质因数分解，p代表质数，α\alpha为p的次数

N=p1α1+p2α2+p3α3+...+pnαnN=p_1^{\alpha_1}+p_2^{\alpha_2}+p_3^{\alpha_3}+...+p_n^{\alpha_n}

然后求出每个质因数的约数个数，如

p1α1:p10,p11,...,p1α1−1,p1α1p_1^{\alpha_1}:p_1^0,p_1^1,...,p_1^{\alpha_1-1},p_1^{\alpha_1}

所以p1α1p_1^{\alpha_1}共α1+1\alpha_1+1个约数，同理p2α2p_2^{\alpha_2}共α2+1\alpha_2+1个约数……

根据乘法原则，N的约数个数为(α1+1)(α2+1)∗...∗(αn+1)(\alpha_1+1)(\alpha_2+1)*...*(\alpha_n+1)

约定nn个整数aia_i，请你输出这些数的乘积的约数个数，答案对109+710^9+7取模

数据范围

1≤n≤1001 \le n \le 100

1≤ai≤2∗1091 \le a_i \le 2*10^9

输入格式

1
2
3
4
3
2
6
8

输出格式

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

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

int main() {
    cin >> n;
    
    unordered_map<int, int> primes;
    while (n -- ) {
        int x;
        cin >> x;
        
        // 分解质因数
        for (int i = 2; i <= x/i; i++) 
            while (x % i == 0) {
                x /= i;
                primes[i]++;
            }
        // 记得收尾
        if (x > 1) primes[x]++;
    }
    
    LL res = 1;
    // 代约数个数公式
    for (auto prime: primes) res = res * (prime.second + 1) % mod;
    
    cout << res << endl;
    
    return 0;
}

约数之和

若要求出N的约数之和，首先将N质因数分解，p代表质数，α\alpha为p的次数

N=p1α1+p2α2+p3α3+...+pnαnN=p_1^{\alpha_1}+p_2^{\alpha_2}+p_3^{\alpha_3}+...+p_n^{\alpha_n}

约数之和公式为

(p10+p11+...+p1α1)(p20+p21+...+p2α2)...(pn0+pn1+...+pnαn)(p_1^0+p_1^1+...+p_1^{\alpha_1})(p_2^0+p_2^1+...+p_2^{\alpha_2})...(p_n^0+p_n^1+...+p_n^{\alpha_n})

至于为什么是这样计算呢？可以利用乘法分配律进行验证，两两相乘后可以得到N的所有约数(因为这些质因数本就是由约数拆解而来，这步相当于还原了用质数还原了约数)，进行加和就是答案

约定nn个整数aia_i，请你输出这些数的乘积的约数之和，答案对109+710^9+7取模

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

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

int n;

int main() {
    cin >> n;
    
    unordered_map<int,int> primes;
    // 质因数分解
    while (n -- ) {
        int x;
        cin >> x;
        
        for (int i = 2; i <= x/i; i++)
            while (x % i == 0) {
                x /= i;
                primes[i]++;
            }
        if (x > 1) primes[x]++;
    }
    
    LL res = 1;
    // 约数之和公式
    for(auto prime: primes) {
        int p = prime.first, a = prime.second;
        LL t = 1;
        while(a--) t = (t * p + 1) % mod;
        res = (res * t) % mod;
    }
    
    cout << res << endl;
    
    return 0;
}

最大公约数

欧几里得算法，即辗转相除法，是将两个数辗转相除直到余数为0

O(log(n))O(log(n))

1
2
3
int gcd(int a, int b) {
    return b? gcd(b, a % b):a;
}

欧拉函数

在数论，对正整数n，欧拉函数是小于等于n的正整数中与n互质的数的数目

互质：若N个整数的最大公因数是1，则称这N个整数互质。

若要求出N的欧拉函数，首先将N质因数分解，p代表质数，α\alpha为p的次数

N=p1α1+p2α2+p3α3+...+pnαnN=p_1^{\alpha_1}+p_2^{\alpha_2}+p_3^{\alpha_3}+...+p_n^{\alpha_n}

那么欧拉函数公式如下：

ϕ(n)=N∗(1−1p1)∗(1−1p2)∗...∗(1−1pk)\phi(n)=N*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})

证明如下（容斥原理）：

N同时代表个数，首先我们减去N中所有质因数倍数的个数

N−Np1−Np2−...−NpkN-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}

但这样减会出现一个问题，有的数同时是p1、p2的倍数，就会被减去两次，所以我们需要加回来

N−Np1−Np2−...−Npk+Np1∗p2+Np2∗p3+...N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1*p_2}+\frac{N}{p_2*p_3}+...

一些数同时是p1、p2、p3的倍数，那么减去3次的同时，又加上了3次，相当于没有计算，因此

N−Np1−Np2−...−Npk+Np1∗p2+Np2∗p3+...−Np1∗p2∗p3−Np2∗p3∗p4N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1*p_2}+\frac{N}{p_2*p_3}+...-\frac{N}{p_1*p_2*p_3}-\frac{N}{p_2*p_3*p_4}

根据以上式子，可以看出一定的规律，将它们进行整理即可得到欧拉函数

求一个数的欧拉函数

1
2
3
4
5
6
7
8
9
10
11
12
13
int phi(int x)  // 欧拉函数
{
    int res = x;
    for (int i = 2; i <= x / i; i ++ )
        if (x % i == 0)
        {
            res = res / i * (i - 1);
            while (x % i == 0) x /= i;
        }
    if (x > 1) res = res / x * (x - 1);

    return res;
}

筛选法求欧拉

给定一个数n，求出1-n中每个数欧拉函数的和。这里使用到了前面线性筛质数的方法

  1. 若i为质数，那么i的欧拉函数ϕ(i)\phi(i)为i-1

  2. 若i不为质数，求ϕ(primes[j]∗i)\phi(primes[j]*i)分两种情况

    1. 若i % primes[j] == 0，意味着primes[j]是i的一个质因子

      此时ϕ(i)=i∗(1−1p1)∗(1−1p2)∗...∗(1−1pk)\phi(i)=i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})

      则ϕ(primes[j]∗i)=primes[j]∗i∗(1−1p1)∗(1−1p2)∗...∗(1−1pk)\phi(primes[j]*i)=primes[j]*i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})

      对比两个式子发现，后者只比前者多乘了一个primes[j]primes[j]，因此我们得到公式为

      ϕ(primes[j]∗i)=primes[j]∗ϕ(i)\phi(primes[j]*i)=primes[j]*\phi(i)

    2. 若i % primes[j] != 0，意味着意味着primes[j]不是i的一个质因子

      此时ϕ(i)=i∗(1−1p1)∗(1−1p2)∗...∗(1−1pk)\phi(i)=i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})

      则ϕ(primes[j]∗i)=primes[j]∗i∗(1−1p1)∗(1−1p2)∗...∗(1−1pk)∗(1−1pj)\phi(primes[j]*i)=primes[j]*i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})*(1-\frac{1}{p_j})

      对比两个式子，因为此时的primes[j]不是i的一个质因子，所以在求primes[j]*i时，需要比原先式子多乘了primes[j]和(1-1/pj)，化简得最后结果如下

      ϕ(primes[j]∗i)=p[j]∗ϕ(i)∗(1−1pj)=ϕ(i)∗(pj−1)\phi(primes[j]*i)=p[j]*\phi(i)*(1-\frac{1}{p_j})=\phi(i)*(p_j-1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
int get_eulers(int n) {
    
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            phi[i] = i - 1;
        }
        for (int j = 0; primes[j] <= n / i; j++) {
            st[primes[j]*i] = true;
            if (i % primes[j] == 0) {
                phi[i*primes[j]] = primes[j] * phi[i];
                break;
            }
            phi[i*primes[j]] = (primes[j]-1) * phi[i];
        }
    }
    LL res = 0;
    for (int i = 1; i <= n; i++) res += phi[i];
    
    return res;
}

欧拉数论定理

内容

若aa与nn互质（gcd(a,n)=1），则$a^{\phi(n)} \equiv 1 ; mod ; n $

其中ϕ(n)\phi(n)称为对模nn缩系的元素个数

  • 当n为质数时，$a^{n-1} \equiv 1 ; mod ; n $（费马小定理）

证明

1~n中，与n互质的数为a1,a2,a3,...,aϕ(n)a_1,a_2,a_3,...,a_{\phi(n)}（注意这里的a与前面的a不同含义）

将上述的数同时乘以a，得

aa1,aa2,aa3,...,aaϕ(n)aa_1,aa_2,aa_3,...,aa_{\phi(n)}

因a1,a2,a3与n互质并且a与n也互质，所以两个数的乘积也和n互质。将两组分别内部同时相乘，得

a1∗a2∗a3∗...∗aϕ(n)≡aa1∗aa2∗aa3∗...∗aaϕ(n)a_1*a_2*a_3*...*a_{\phi(n)} \equiv aa_1*aa_2*aa_3*...*aa_{\phi(n)}

a1∗a2∗a3∗...∗aϕ(n)≡aϕ(n)(a1∗a2∗a3∗...∗aϕ(n))a_1*a_2*a_3*...*a_{\phi(n)} \equiv a^{\phi(n)}(a_1*a_2*a_3*...*a_{\phi(n)})

约去相同部分，得

aϕ(n)≡1(mod  n)a^{\phi(n)} \equiv 1 (mod \;n)

因此我们得出结论

$a^{\phi(n)} \equiv 1 \ mod \ n $

应用

这个定理可以用来简化幂的模运算。比如计算72227^{222}的个位数，实际是求72227^{222}被10除的余数。7和10互质，且φ(10)=4。由欧拉定理知74≡1 (mod 10)7^4\equiv 1\ (mod\ 10)。所以7222=(74)55∗72≡49≡9 (mod10)7^{222}=(7^{4})^{55}*7^2\equiv 49\equiv 9\ (mod 10)

快速幂

1≤a,k,p≤1091 \le a,k,p \le 10^9

以log(k)log(k)的时间复杂度，快速求出ak mod pa^k \ mod\ p的结果

  • 要快速求出上述的结果，我们需要先预处理出a20,a21...(mod p)a^{2^0},a^{2^1}...(mod\ p)的值
  • 然后将aka^k进行二进制分解，ak=a20∗a21∗...=a20+21+...(mod p)a^k=a^{2^0}*a^{2^1}*...=a^{2^0+2^1+...}(mod\ p)
  • 最后用原先预处理的结果，相乘取模得出答案

此算法主要的时间花在了预处理O(log(k))O(log(k))和k的填二进制拆解O(log(k))O(log(k))上

虽然上述过程有一些复杂，但在最后的算法实现非常简洁

1
2
3
4
5
6
7
8
9
10
11
12
13
// a^k%p
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; 
}

扩展欧几里得算法

求出一组xi,yix_i,y_i，使其满足ai∗xi+bi∗yi=gcd(ai,bi)a_i*x_i+b_i*y_i=gcd(a_i,b_i)

推导如下：

by+(a mod b)x=dby+(a\ mod\ b)x=d

⇒by+(a−[ab]∗b)x=d\Rightarrow by+(a-[\frac{a}{b}]*b)x=d

⇒ax+b(y−[ab]∗x)=d\Rightarrow ax+b(y-[\frac{a}{b}]*x)=d

1
2
3
4
5
6
7
8
9
10
11
int exgcd(int a, int b, int& x, int& y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    } 
    
    int d = exgcd(b, a%b, y, x);
    y -= a / b * x;
    
    return d;
}

线性同余方程

题目描述

给定nn组数据ai,bi,mia_i,b_i,m_i，对于每组数求出一个xix_i，使其满足ai∗xi≡bi(mod mi)a_i*x_i \equiv b_i(mod\ m_i)

推导

ai∗xi≡bi(mod mi)a_i*x_i \equiv b_i(mod\ m_i)

⇒ai∗xi/mi=yi...bi\Rightarrow a_i*x_i / m_i = y_i...b_i

⇒ai∗xi−mi∗yi=bi\Rightarrow a_i*x_i - m_i*y_i=b_i

令yi′=yiy_i\prime = y_i

⇒a∗x+m∗y′=b\Rightarrow a*x + m*y\prime=b

如果b不是最大公约数d的倍数，那么一定无解

根据扩展欧几里得算法，对于任意的整数a和m，必定存在一组x和y'，使得ax+my'等于gcd(a,m)

所以可以求得a∗x+m∗y′=da*x + m* y\prime =d

⇒a∗x∗bd+m∗y′∗bd=b\Rightarrow a*x*\frac{b}{d} + m*y\prime*\frac{b}{d}=b

所以要求解的xi=x*b/d

中国剩余定理

中国剩余定理(Chinese remainder theorem, CRT)

小学数学时大家都学过这样的问题：有一群人，3个人一组，剩1人，5个人一组，剩3人，10人一组，剩7人。

那么有多少人呢？CRT的出现就是用来解决上述问题

以上的问题可以看成一个线性同余方程组

CRT公式：x=(∑i=1nai∗Mi−1∗Mi) mod Mx=(\sum_{i=1}^n a_i*M_i^{-1}*M_i)\ mod \ M

以上问题可以转换成以下的线性同余方程组表达

{x≡a1(mod m1)x≡a2(mod m2)⋯⋯x≡ak(mod mk)\begin{cases} x\equiv a_1(mod \ m_1)\\ x\equiv a_2(mod \ m_2) \\ \cdots \cdots \\ x\equiv a_k(mod \ m_k)\end{cases}

M=m1∗m2∗⋯∗mkM=m_1*m_2*\cdots*m_k

Mi=M/mi=⋯∗mi−1∗mi+1∗⋯M_i=M/m_i=\cdots*m_{i-1}*m_{i+1}*\cdots

此外，还需要求出MiM_i的逆元Mi−1M_i^{-1}，使其满足条件Mi∗Mi−1≡1(mod mi)M_i*M_i^{-1} \equiv 1(mod \ m_i)

而对于本组中的其它i，满足条件Mj∗Mj−1≡0(mod mi)M_j*M_j^{-1} \equiv 0(mod \ m_i)

根据x=ai∗Mi−1∗Mix=a_i*M_i^{-1}*M_i，求出其中的一个值

对每一线性同余方程执行上述过程，最后将结果对M取模，得出最后答案

举个例子，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

即，一个整数除以3余2，除以5余3，除以7余2，求这个整数。

M=3∗5∗7=105M=3*5*7=105

首先求解x≡2(mod 3)x \equiv 2 (mod \ 3)

M1=M/m1=105/3=35M_1=M/m_1=105/3=35

M1−1∗35≡1(mod 3)M_1^{-1}*35 \equiv 1 (mod \ 3)，此步利用扩展欧几里得算法，求出M1−1=2M_1^{-1} = 2

x=ai∗Mi−1∗Mi=2∗2∗35=2∗70=140x=a_i*M_i^{-1}*M_i=2*2*35=2*70=140

同理，可以求出x≡3(mod 5)x \equiv 3 (mod \ 5)和x≡2(mod 7)x \equiv 2 (mod \ 7)

x=15∗2+21∗3+70∗2=233 mod M=233 mod 105=23x=15*2+21*3+70*2=233 \ mod \ M = 233 \ mod \ 105 = 23

这个余数23就是合乎条件的最小数.

扩展CRT

根据题目所述，这里的a1,a2,⋯ ,ana_1,a_2,\cdots,a_n不满足两两互质的条件，所以需要重新推导CRT

{x≡m1(mod a1)x≡m2(mod a2)\begin{cases} x\equiv m_1(mod \ a_1)\\ x\equiv m_2(mod \ a_2) \\ \end{cases}

化简得，

{x=k1a1+m1x=k2a2+m2\begin{cases} x= k_1a_1+m_1\\ x= k_2a_2+m_2 \\ \end{cases}

所以k1a1+m1=k2a2+m2k_1a_1+m_1=k_2a_2+m_2

⇒k1a1−k2a2=m2−m1\Rightarrow k_1a_1-k_2a_2=m_2-m_1，这个式子有解则满足gcd(a1,a2)∣m2−m1gcd(a_1,a_2) | m_2-m_1

运用扩展欧几里得算法，能够求出k1,k2的通解,非齐次=齐次通解+非齐次特解

{k1=k1+ka2dk2=k2+ka1d\begin{cases} k_1= k_1+k\frac{a_2}{d}\\ k_2= k_2+k\frac{a_1}{d} \\ \end{cases}

x=k1a1+m1x=k_1a_1+m_1

⇒x=(k1+ka2d)a1+m1\Rightarrow x=(k_1+k\frac{a_2}{d})a_1+m_1

⇒x=a1k1+m1+ka1a2d\Rightarrow x=a_1k_1+m_1+k\frac{a_1a_2}{d}

其中a1a2/d就是a1a2的最小公倍数，因此

⇒x=a1k1+m1+k[a1a2]\Rightarrow x=a_1k_1+m_1+k[a_1a_2]

通过对比最前面化简后的式子，我们可以把a1k1+m1a_1k_1+m_1看成m，k[a1a2]k[a_1a_2]看成ka

⇒x=x0+ka\Rightarrow x=x_0+ka

以这样的形式重复更新合并后来的线性同余方程，可以得到最终答案

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>
#include <cmath>

using namespace std;

typedef long long LL;
int n;

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


int main() {
    cin >> n;
    
    LL a1, m1;
    cin >> a1 >> m1;
    
    bool has_answer = true;
    for (int i = 0; i < n-1; i ++ ) {
        LL a2, m2;
        cin >> a2 >> m2;
        
        LL k1, k2;
        LL d = exgcd(a1, a2, k1, k2);
        if ((m2-m1) % d) { // 当m2-m1不是d的倍数时，无解
            has_answer = false;
            break;
        }
        k1 *= (m2 - m1)/d; // 欧几里得求出的解是以最小公约数d为结果，转换为m2-m1
        LL t = a2 / d;
        k1 = (k1 % t + t) % t; // 最小正整数解
        
        m1 = a1*k1+m1; // 为下一轮循环计算m1,a1
        a1 = abs(a1 / d * a2); // 不能a1*a2就近相乘，否则会溢出
    }
    
    if (has_answer) cout << (m1 % a1 + a1) % a1 << endl;
    else puts("-1");
    
    return 0;
}

204.表达整数的奇怪方式

高斯消元

求解线性方程组

题目描述

输入一个包含n个方程和n个未知数的线性方程组。方程组的系数为实数，求解这个方程组

数据范围

1≤n≤1001 \le n \le 100

步骤如下，

  • 枚举每一列
    1. 找到该列中绝对值最大的那个数所在的行
    2. 将该行移动到最上面（次上、次次上）
    3. 将该行的第一个数变为1
    4. 用这个数把下面行该列的数消成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
const int N = 110;
const double eps = 1e-6;

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

int gauss() {
    int r, c;
    for (r = 0, c = 0; c < n; c++) {
        int t = r;
        // 1、找到该列的最大值的行号
        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;
        
        // 2、交换其到第r行
        for (int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
        // 3、将该行归一化
        for (int i = n; i >= c; i--) a[r][i] /= a[r][c];
        // 4、用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[i][c] * a[r][j];
                }
        r++;
        
    }
    if (r < n) {
        for (int i = r; i < n; i++)
            if (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; // 唯一解
}

组合计数

排列组合中的组合

组合I

给定两个整数a,b，请你输出CabC_a^b

根据递推式Cab=Ca−1b+Ca−1b−1C_a^b=C_{a-1}^b+C_{a-1}^{b-1}，时间复杂度n2n^2

image-20220305214658795

数据范围

1≤b≤a≤20001 \le b \le a \le 2000

1≤n≤100001 \le n \le 10000

1
2
3
4
5
6
7
8
9
const int N = 2010, mod = 1e9+10;
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; // 递推
}

组合II

预处理出阶乘，用infact(a!)表示a!的逆元（快速幂）

Cab=a!(a−b)!b!=fact(a!)∗infact(b!)∗infact((a−b)!)C_a^b=\frac{a!}{(a-b)!b!}=fact(a!)*infact(b!)*infact((a-b)!)，时间复杂度nlog(n)nlog(n)

数据范围

1≤b≤a≤1051 \le b \le a \le 10^5

根据费马小定理，an−1≡1 mod na^{n-1} \equiv 1\ mod \ n，其中a与n互质，且n为质数

得an∗an−2≡mod na^n*a^{n-2} \equiv mod \ n

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
#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;
        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);
    }
    
    return 0;
}

组合III

数据范围

1≤b≤a≤10181 \le b \le a \le 10^{18}

1≤n≤201 \le n \le 20

1≤p≤1051 \le p \le 10^5

卢卡斯定理

Cab=Ca mod pb mod p∗Ca/pb/p(mod p)C_a^b=C_{a \ mod \ p}^{b \ mod \ p}*C_{a / p}^{b / p}(mod \ p)

参考资料

AcWing算法基础课，当你还在犹豫的时候，别人已经学完了