【算法笔记】ACM数论基础模板

本文中默认：using LL = long long;

几个定理

唯一分解定理

对于一个大于1的正整数$n$，$n$可以分解质因数为：
$$n={\Pi }_{i=1}^k{p}_i^{a_i}=p_1^{a_1}\times p_2^{a_2}\times ...\times p_k^{a_k}$$

$n$的质因子$p$的个数：
$$\sum _{i=1}^n\frac{N}{p^i}$$

约数定理：
$$d(n)=\prod _{i=1}^k(a_i+1)=(a_1+1)(a_2+1)...(a_k+1)$$

鸽巢原理（抽屉原理）

桌上有十个苹果，要把这十个苹果放到九个抽屉里，无论怎样放，我们会发现至少会有一个抽屉里面放不少于两个苹果。这一现象就是我们所说的“抽屉原理”。 抽屉原理的一般含义为：“如果每个抽屉代表一个集合，每一个苹果就可以代表一个元素，假如有n+1个元素放到n个集合中去，其中必定有一个集合里至少有两个元素。” 抽屉原理有时也被称为鸽巢原理。

麦乐鸡定理

给定两个互质的数 $n,m$ ，定义 $x=a\ast n+b\ast m（a\ge 0,b\ge 0）$，当 $x>n\ast m-n-m$ 时，该式子恒成立。

哥德巴赫猜想

任何一个大于 $5$ 的整数都可写成三个质数之和；任何一个大于 $2$ 的偶数都可写成两个素数之和。

容斥原理

对于$n$个集合$A_1,A_2,A_3,...A_n:$$$(A_1\cup A_2\cup A_3\cup ...\cup A_n)=\sum _{0<i\le n}{A}_i-\sum _{0<i<j\le n}(A_i\cap A_j)+\sum _{0<i<j<k\le n}(A_i\cap A_j\cap A_j)+...+(-1{)}^{n+1}(A_1\cap A_2\cap A_3\cap ...\cap A_n)$$（加减加减加减…）

例题

给定一个整数 $n$ 和 $m$ 个不同的质数 $p_1,p_2,...,p_m$，请你求出 1 ∼ $n$ 中能被 $p_1,p_2,...,p_m$ 中的至少一个数整除的整数有多少个。

二进制枚举解

LL n, m;
LL p[N];

int main(){
    ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
    cin >> n >> m;
    for (int i = 0; i < m; i ++ )
        cin >> p[i];
    LL ans = 0;
    for (int i = 1; i < (1 << m); i ++ ){
        LL t = 1, cnt = 0;
        for (int j = 0; j < m; j ++ ){
            if (i >> j & 1){
                cnt ++ ;
                t *= p[j];
                if (t > n){
                    t = -1;
                    break;
                }
            }
        }
        if (t != -1){
            if (cnt & 1) ans += n / t;
            else ans -= n / t;
        }
    }
    cout << ans << endl;
    return 0;
}

dfs解

LL n, res;
int m;
LL p[N];

void dfs(int x, LL s, LL odd) {
    if (x == m) {
        if (s != 1) ans += odd * (n / s);
        return;
    }
    dfs(x + 1, s, odd);
    if (s <= n / p[x]) dfs(x + 1, s * p[x], -odd);
}

int main() {
    ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
    cin >> n >> m;
    for (int i = 0; i < m; i++) {
        cin >> p[i];
    }
    res = 0;
    dfs(0, 1, -1);
    cout << res << endl;
    return 0;
}

裴蜀定理

如果 $a$ 和 $b$ 是两个整数，则存在整数 $x,y$ 使得：
$$a\cdot x+b\cdot y=\gcd (a,b)$$

例题

给定一个序列 $a$，找到一个序列 $x$，使得 ${\sum }_{i=1}^n{a}_i{x}_i$ 最小。

代码

LL n, a, ans;
LL gcd(LL a, LL b){
    return b ? gcd(b, a % b) : a;
}
int main(){
    cin >> n;
    for (int i = 0; i < n; i ++ ){
        cin >> a;
        if (a < 0) a = -a;
        ans = gcd(ans, a);
    }
    cout << ans << endl;
    return 0;
}

最大公约数、最小公倍数

最大公约数

欧几里得算法（辗转相除法）

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

或者直接

__gcd(a, b);

最小公倍数

$a和b的最小公倍数=\frac{a\times b}{a和b的最大公约数}$

防止越界，先除再乘。

int lcm(int a, int b){
	return a / gcd(a,b) * b ;
}

质数

试除法判断质数

见"C语言入门讲义"

bool is_prime(int x){
    if(x < 2) return false;
    for(int i = 2; i <= x / i; i++){
    	if(x % i == 0) return false;
    }
    return true;
}

分解质因数

分解出$n$的$k$个质因数：$n=i_1^{s_1}\times i_2^{s_2}\times ...\times i_k^{s_k}$

void divide(int x){
    for(int i = 2; i <= x / i; i++){
    	if(x % i == 0){
            int s = 0;
            while(x % i == 0){
				x /= i;
				s++;
			}
            cout << i << ' ' << s << endl; // i:底数 s:指数
        }
    }
    if(x != 1) cout << x << ' ' << 1 << endl;
    cout << endl;
}

筛质数

将$1$~$n$中所有的质数筛到$primes$数组中

朴素筛法（埃氏筛法）

vector<int> primes;
bool not_prime[N]; // not_prime[i] = true 表示 i 不是质数

void get_primes(int n){
    for(int i = 2; i <= n; i++){
        if(not_prime[i]) continue;
        primes.push_back(i);
        for(int j = i + i; j <= n; j += i) not_prime[j] = true;
    }
}

线性筛法（欧拉筛法）

vector<int> primes;
bool not_prime[N]; // not_prime[i] = true 表示 i 不是质数

void get_primes(int n){
	for(int i = 2; i <= n; i++){
		if(!not_prime[i]) primes.push_back(i);
		for(int p : primes){
			if(i * p > n) break;
			not_prime[i * p] = true;
			if(i % p == 0) break;
		}
	}	
}

约数

试除法求约数

vector<int> get_divisors(int x){
    vector<int> res;
    for(int i = 1; i <= x / i; i ++)
        if(x % i == 0){
            res.push_back(i);
            if(i != x / i) res.push_back(x / i);
        }
    sort(res.begin(), res.end());
    return res; // 返回一个vector,里面是n的所有约数（因子）
}

求约数个数

一个数求约数个数

int get_divisors(int x){
    unordered_map<int, int> primes; // 存x的质因数和其对应的次数
    for(int i = 2; i <= x / i;i++){
        while(x % i == 0){
            x /= i;
            primes[i]++;
        }
    }
    if(x > 1) primes[x]++; // 如果最后剩下一个大于 1 的数，说明它也是一个质因数
    LL res = 1;
    for(auto p : primes) res = res * (p.second + 1); // 约数定理
    return res; 
}

求1~n所有数的约数个数

$O(nlogn)$筛法

int cnt[N]; // cnt[i] : i 的约数个数

for(int i = 1; i <= n; i++){
    for(int j = i; j <= n; j += i){
    	cnt[j]++;
    }
}

$O(n)$筛法

bool not_prime[N];
int d[N], num[N]; // d[i] 表示 i 的约数个数，num[i] 表示 i 的最后一个质因子的幂（从 i 分解质因数得到的最后一个质因子 j，i = a * j^num）
vector<int> primes; // 筛质数

void get_divisors(int n){
    d[1] = 1;
    num[1] = 0;
    for(int i = 2; i <= n; i++){
        if(!not_prime[i]){
            primes.push_back(i);
            num[i] = 1;
            d[i] = 2;
        }
        for(int p : primes){
            if(i * p > n) break;
            not_prime[i * p] = true;
            if(i % p == 0){
                num[i * p] = num[i] + 1;
                d[i * p] = d[i] / (num[i] + 1) * (num[i * p] + 1);
                break;
            } 
            else{
                num[i * p] = 1;
                d[i * p] = d[i] * 2;
            }
        }
    }
}

约数之和

一个数求约数之和

LL get_divisor_sum(int x){
    unordered_map<int, int> primes; // 存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 p : primes){
        LL t = 1, sum = 1;
        for(int i = 0; i < p.second; i++){
            t *= p.first;
            sum  += t;
        }
        res *= sum;
    }
    return res; // res表示x的所有约数的和
}

求1~n所有数的“约数之和”

$O(nlogn)$筛法

LL sum[N]; // sum[i]表示i的约数和

for(int i = 1; i <= n; i++){
    for(int j = i; j <= n; j += i){
    	sum[j] += i;
    }
}

$O(n)$筛法

bool not_prime[N];
int num[N]; // num[i] 表示 i 的最后一个质因子的幂次（i = a * p^num）
vector<int> primes;
LL sd[N], sp[N]; // sd[i] 表示 i 的所有约数之和, sp[i] 表示 i 最后一个质因子的贡献项（1 + p + p^2 + ... + p^k）

void get_divisors(int n){
    sd[1] = 1;
    num[1] = 0;
    sp[1] = 1;
    for(int i = 2; i <= n; i++){
        if(!not_prime[i]){
            primes.push_back(i);
            num[i] = 1;
            sp[i] = 1 + i;
            sd[i] = 1 + i;
        }
        for(int p : primes){
            if(i * p > n) break;
            not_prime[i * p] = true;
            if(i % p == 0){
                num[i * p] = num[i] + 1;
                sp[i * p] = sp[i] * p + 1;
                sd[i * p] = sd[i] / sp[i] * sp[i * p];
                break;
            } 
            else{
                num[i * p] = 1;
                sp[i * p] = 1 + p;
                sd[i * p] = sd[i] * sp[i * p];
            }
        }
    }
}

欧拉函数

$1\sim N$ 中与 $N$ 互质的数的个数被称为欧拉函数，记为 $\phi (N)$。

若在算术基本定理中，$N=p_1^{a_1}{p}_2^{a_2}\dots p_m^{a_m}$，则：

$$\phi (N)=N\times \frac{p_1-1}{p_1}\times \frac{p_2-1}{p_2}\times \cdots \times \frac{p_m-1}{p_m}$$

定义求一个数的欧拉函数

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; // res即为phi(x)的值
}

筛法求1~n每个数的欧拉函数

vector<int> primes;
bool not_prime[N];
int euler[N]; // eular[i]表示phi(i)的值

void get_eulers(int n){
    euler[1] = 1;
    for(int i = 2; i <= n; i++){
        if(!not_prime[i]){
            primes.push_back(i);
            euler[i] = i - 1;
        }
        for(int p : primes){
            LL t = 1LL * i * p;
            if(t > n) break;
            not_prime[t] = true;
            if(i % p == 0){
                euler[t] = euler[i] * p;
                break;
            } 
            else{
                euler[t] = euler[i] * (p - 1);
            }
        }
    }
}

快速幂

快速幂求a的b次方

不带取模

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

带取模

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

拓展欧几里得算法

原理

给定 $n$ 对正整数 $a_i,b_i$，对于每对数，求出一组 $x_i,y_i$，使其满足 $a_i\times x_i+b_i\times y_i=\gcd (a_i,b_i)$。

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

int main(){
    int n;
    cin >> n;
    while(n --){
        int a, b;
        cin >> a >> b;
        int x, y;
        exgcd(a, b, x, y);
        cout << x << ' ' << y << endl; 
    }
    return 0;
}

拓展欧几里得算法求解线性同余方程

给定 $n$ 组数据 $a_i,b_i,m_i$，对于每组数求出一个 $x_i$，使其满足 $a_i\times x_i\equiv b_i(\mathrm{m}\mathrm{o}\mathrm{d}{m}_i)$，如果无解则输出 impossible。

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

int main(){
    int n;
    cin >> n;
    while(n--){
        int a, b, m;
        cin >> a >> b >> m;
        int x, y;
        int d = exgcd(a, m, x, y);
        if(b % d) cout << "impossible" << endl;
        else cout << (LL)b / d * x % m << endl;
    }
    return 0;
}

逆元

$b$的乘法逆元，记作：$b^{-1}$。
$$\frac{a}{b}(\mathrm{m}\mathrm{o}\mathrm{d}p)=a\times b^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
（将题中的运算除法变成运算乘法）

快速幂求逆元

适用条件：模数$p$一定要是质数，并且$n$和$p$互质，即$gcd(n,p)=1$

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

LL inv(int n, int p){ // n 在模 p下的逆元
	return qmi(n, p - 2, p);
}

拓展欧几里得算法求逆元

适用条件：模数$p$可以是任何数，只要满足$n$和$p$互质，即$gcd(n,p)=1$即可。

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

LL inv(LL n, int p){ // n 在模 p下的逆元
    LL x, y;
    LL d = exgcd(n, p, x, y);
    if(d == 1) return (x % p + p) % p;
    else return -1;
}

预处理1~n的逆元

LL inv[N]; // inv[i] : i 在模 p 下的逆元

inv[1] = 1;
for(int i = 2; i <= n; i++){
    inv[i] = (p - (p / i) % p) * inv[p % i] % p;
}

中国剩余定理

中国剩余定理求解线性同余方程组（所有a彼此互质）

给定 $2n$ 个整数 $a_1,a_2,\dots ,a_n$ 和 $m_1,m_2,\dots ,m_n$，求一个最小的非负整数 $x$，满足
$\forall i\in [1,n],x\equiv m_i(\mathrm{m}\mathrm{o}\mathrm{d}{a}_i)$。

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

LL inv(LL n, LL p){
    LL x, y;
    LL d = exgcd(n, p, x, y);
    if(d != 1) return -1;
    return (x % p + p) % p;
}

LL CRT(LL a[], LL m[], int n){
    LL mm = 1;
    for(int i = 0; i < n; i++){
        mm *= m[i];
    }
    LL res = 0;
    for(int i = 0; i < n; i++){
        LL t = mm / m[i];
        LL inv_t = inv(t, m[i]);
        res = (res + t * inv_t % mm * (a[i] % mm)) % mm;
    }
    if(res < 0) res += mm;
    return res;
}

中国剩余定理求解线性同余方程组（所有a不一定互质）

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

LL exCRT(LL a[], LL m[], int n){
    LL M = m[0], A = a[0];
    for(int i = 1; i < n; i++){
        LL x, y;
        LL d = exgcd(M, m[i], x, y);
        if((a[i] - A) % d != 0) return -1;
        LL t = m[i] / d;
        x = x * ((a[i] - A) / d);
        x = (x % t + t) % t;
        A += x * M;
        M = M / d * m[i];
        A %= M;
    }
    A = (A % M + M) % M;
    return A;
}

组合数

一组测试数据

提供一组测试数据：$C_{132}^{66}=$ 377’389’666’165’540’953’244’592’352’291’892’721’700，模数为 $998244353$ 时为 $241^{\prime }200^{\prime }029$；$10^9+7$ 时为 $598375978$。

几个性质

组合数的两种表示

$$\left(\genfrac{}{}{0ex}{}{n}{k}\right)\iff C_n^k$$

递推式

$$C_n^m=C_n^{m-1}+C_{n-1}^{m-1}$$

阶乘式

$$C_n^m=\frac{n!}{m!\times (n-m)!}$$

卢卡斯定理

对于素数 $p$，有
$${\mathrm{C}}_n^k\equiv {\mathrm{C}}_{\lfloor n/p\rfloor }^{\lfloor k/p\rfloor }\times {\mathrm{C}}_{n\mathrm{m}\mathrm{o}\mathrm{d}p}^{k\mathrm{m}\mathrm{o}\mathrm{d}p}(\mathrm{m}\mathrm{o}\mathrm{d}p).$$

其中，当 $n<k$ 时，组合数 ${\mathrm{C}}_n^k$ 规定为 $0$。

递推求组合数

适用条件：$C_a^b$中的$a$和$b$都比较小，二维数组$C[a][b]$要能存的下

例如：$1\le n\le 10^5,1\le b\le a\le 2000$

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

int C[N][N]; // C[a][b] 即为C_a^b

void init()  {
    for(int i = 0; i < N; i++){
        for(int j = 0; j <= i; j++){
            if(j == 0) C[i][j] = 1;
            else C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
        }
    }
}

预处理阶乘求组合数

用快速幂预处理阶乘和逆元求组合数。

适用条件：$a$和$b$的范围适中，二维数组存不下但可以开一维。

比如：$1\le n\le 10^5,1\le b\le a\le 10^5$

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

LL inv(LL n, int p){
	return qmi(n, p - 2, p);
}

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

LL Cnk(int n, int k){ // Cnk即为C_n^k
	if(k < 0 || k > n) return 0;
    return (LL)fact[n] * infact[k] % mod * infact[n - k] % mod;
}

void init(){
	fact[0] = 1;
	for(int i = 1; i < N; i++){
        fact[i] = fact[i - 1] * i % mod;
    }
    infact[N - 1] = inv(fact[N - 1]);
    for(int i = N - 1; i > 0; i--){
        infact[i - 1] = infact[i] * i % mod;
    }
}

卢卡斯定理求组合数

适用情况：$a$和$b$特别大，用数组存不下，无法预处理，且查询较少，每次直接用定理求组合数不会超时。（模数必须为质数）

比如：$1\le n\le 20，1\le b\le a\le 10^{18}$

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

int Cnk(int a, int b, int p){ // 较小数（0~p-1）的组合数
    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) % p;
    }
    return res;
}

// c[a][b] = c[a % p][b % p] * c[a / p][b / p]
int lucas(LL a, LL b, int p){ // lucas(a, b, p) 即为C_a^b % p
    if(a < p && b < p) return C(a, b, p);
    return (LL)Cnk(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}

卢卡斯定理 + 预处理阶乘求组合数

适用情况：$a$和$b$特别大，用数组存不下，无法预处理，且查询多，每次都求一遍组合数会超时。（模数必须为质数）

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

LL fact[N], infact[N];

void init(){
	fact[0] = 1;
	for(int i = 1; i < p; i++){
        fact[i] = fact[i - 1] * i % mod;
    }
    infact[p - 1] = inv(fact[p - 1]);
    for(int i = p - 1; i > 0; i--){
        infact[i - 1] = infact[i] * i % mod;
    }
}

LL Cnk(int n, int k){ // 较小数（0~p-1）的组合数
	if(k < 0 || k > n) return 0;
    return (LL)fact[n] * infact[k] % mod * infact[n - k] % mod;
}

// c[a][b] = c[a % p][b % p] * c[a / p][b / p]
int lucas(LL a, LL b, int p){ // lucas(a, b, p) 即为C_a^b % p
    if(a < p && b < p) return C(a, b, p);
    return (LL)Cnk(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}