【简单数论】（模运算，快速幂，乘法逆元，同余，exgcd，gcd，欧拉函数，质数，欧拉筛，埃式筛，调和级数枚举，约数，组合数）

数论

模运算

• $amodb=a-\lfloor a/b\rfloor \times b$

• $nmodp$得到的结果的正负至于被除数$n$有关

模运算的性质：

$(a+b)modm=((amodm)+(bmodm))modm$

$(a-b)modm=((amodm)-(bmodm))modm$ = $((amodm)-(bmodm)+m)modm$

$(a\times b)modm=((amodm)\times (bmodm))modm$

但是除法例外，除法的取模需要用到逆元。

计算减法的时候，通常需要加上模数，防止出现负数。

快速幂

快速求解$a^{b}modp$，可以采用二进制拼凑的思想

对于$b$，它可以看成一个二进制数，可以写成$2$的幂相加的形式，那么$a^b$可以写成$a^{(2^i+2^j+...+2^k)}=a^{2^i}\times a^{2^j}\times ...\times a^{2^k}(i,j,k\in S,其中集合S为b转换成二进制后位置上是1的位置)$

我们只需计算的同时预处理每一项$a^{2^i}$即可

using ll = long long;
ll qmi(ll a, ll b, ll c){
    ll res = 1;
    while(b){
        if(b & 1 == 1) res = res * a % c;
        a = a * a % c;
        b >>= 1;
    }
    return res;
}

乘法逆元

若$a\times x\equiv 1(modb)$，且$a$与$b$互质，我们定义$x$为$a$的逆元，记为$a^{-1}$,$x$可称为$a$在模$b$意义下的倒数

对于$\frac{a}{b}modp$，我们可以求出$b$在$modp$意义下的逆元，然后乘上$a$，再$modp$即可

注意对于数$a$在模$p$意义下的乘法逆元，我们需要确保$a$与$p$互质，此时才有$a$的乘法逆元$a^{-1}$，此条件等价于$p$是质数

费马小定理

$a^{p-1}\equiv 1(modp)$，其中$p$是素数

费马小定理求逆元
$a\times a^{p-2}\equiv 1(modp)$，根据乘法逆元的定义可知，此时$a$的逆元是$a^{p-2}$

ll inv(ll a, ll p){
    return qmi(a, p - 2, p);
}

扩展欧几里得求逆元

根据逆元的定义$ax\equiv 1(modp)$，我们得到方程：$ax-1=yp$（$ax-1$ 是$p$的倍数），则$ax+py=1$

解方程得$xmodp$得到的就是$a$的乘法逆元，同时逆元存在的条件是$gcd(a,p)=1$

求阶乘的乘法逆元

同余

两个整数$a$和$b$对用一个正整数$m$取模后余数相同，则称$a$和$b$对模$m$同余，记作$a\equiv b(modm)$，这意味着$m|(a-b)$

同余的性质：

• 自反性：$a\equiv a(modm)$
• 对称性：若$a\equiv b(modm)$，则$b\equiv a(modm)$
• 传递性：若$a\equiv b(modm),b\equiv c(modm)$，则$a\equiv c(modm)$
• 同加性：若$a\equiv b(modm)$，则$a\pm c\equiv b\pm c(modm)$
• 同乘性：若$a\equiv b(modm)$，则$a\times c\equiv b\times c(modm)$，若$a\equiv b(modm),c\equiv d(modm)$，则$a\times c\equiv b\times d(modm)$
• 同幂性：若$a\equiv b(modm)$，则$a^c\equiv b^c(modm)$
• 不满足同除性，但是有，若$ca\equiv cb(modm)$，则$a\equiv b(mod\frac{m}{gcd(m,c)})$

扩展欧几里得定理

对于给定的两个整数$a$，$b$，必须存在整数$x$与$y$，使得$a\times x+b\times y=gcd(a,b)$

裴蜀定理

方程$a\times x+b\times y=c$有整数解的充分必要条件是$c$是$gcd(a,b)$的倍数

证明：

• 必要性：令$p=gcd(a,b)$，则有$a=p\times a^{\prime }$，$b=p\times b^{\prime }$；则$c=p\times (a^{\prime }x+b^{\prime }y)$，即$gcd(a,b)$的倍数。

• 充分性：考虑$gcd$的欧几里得算法。$a,b$通过$b,a\%b$的方式一直辗转下去，最后会出现$p,0$的形式。原因是$0\le a\%b\le b-1$。对于递归出口$p,0$，方程$a\times p+0\times y=c$是否存在解呢？显然有解，由于充分性，这里$a$取$\frac{c}{p}$即可。那么辗转相除过程中，下面的成立能否推得上面的成立呢？

  • 对于方程$b{x}_1+a\%b\times y_1=c=ax+by$

  • $左边=b{x}_1+(a-\lfloor \frac{a}{b}\rfloor \times b)\times y_1=a{y}_1+b{x}_1-\lfloor \frac{a}{b}\rfloor \times b{y}_1=a{y}_1+b\times (x_1-\lfloor \frac{a}{b}\rfloor \times y_1)=右边=ax+by$

    所以找到一组特解：$x=y_1$，$y=x_1-\lfloor \frac{a}{b}\rfloor \times y_1$，就可以通过最下面的特解推得所有的特解。

  • 对于特解$(x_0,y_0)$，设下一组解是$(x_0+d_1,y_0+d_2)$

    有$a\times (x_0+d_1)+b\times (y_0+d_2)=c$，又，$a\times x_0+b\times y_0=c$，则$a{d}_1+b{d}_2=0$

    故，$\frac{d_1}{d_2}=-\frac{b}{a}=-\frac{b/gcd(a,b)}{a/gcd(a,b)}$

    则$d_1=\frac{b}{gcd(a,b)}$，$d_2=-\frac{a}{gcd(a,b)}$，故一般解为：

    $x=x_0+k\times \frac{b}{gcd(a,b)}$，$y=y_0-k\times \frac{a}{gcd(a,b)}(k\in Z)$

  • 若要求$x$的最小整数解，则可以通过$(x_{0}mod(\frac{b}{gcd(a,b)})+\frac{b}{gcd(a,b)})mod{x}_0$，原因是$x_0$的周期是$\frac{b}{gcd(a,b)}$，所以可通过取模的方式得到结果

扩展欧几里得算法

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

GCD

最大公约数：欧几里得算法，时间复杂度$O(loga)$

$gcd(a,b)=gcd(b,amodb)$

证明如下：

• 令$p=gcd(a,b)$，则有$a=p\times a^{\prime }$，$b=p\times b^{\prime }$，$gcd(a^{\prime },b^{\prime })=1$

  • $amodb=a-\lfloor a/b\rfloor \times b=p\times a^{\prime }-p\times b^{\prime }\times \lfloor a/b\rfloor$

  • 可以看出，$p|b$， $p|(amodb)$

• 那么$p$是否是最大公约数呢？即我们需要证明$gcd(b^{\prime },a^{\prime }-b^{\prime }\times \lfloor a/b\rfloor )=1$

  • 令$gcd(b^{\prime },a^{\prime }-b^{\prime }\times \lfloor a/b\rfloor )=d$，则有$b^{\prime }=d\times b^{\prime \prime }$，$a^{\prime }-b^{\prime }\times \lfloor a/b\rfloor =d\times c$
  • 转换得到$a^{\prime }=d\times c+b^{\prime }\times \lfloor a/b\rfloor =d\times c+d\times b^{\prime \prime }\times \lfloor a/b\rfloor$
  • 可以发现$d|b^{\prime }$， $d|a^{\prime }$，与上述$gcd(a^{\prime },b^{\prime })=1$矛盾，则命题得证。

欧拉函数

$\mathrm{1...}N$中与$N$互质的数的个数，被称为欧拉函数，记作$\varphi (N)$

在唯一分解定理中，$N=p_1^{a_1}\times p_2^{a_2}\times ...\times p_m^{a_m}$，则$\varphi (N)=N\times \frac{p_1-1}{p_1}\times \frac{p_2-1}{p_2}\times ...\times \frac{p_m-1}{p_m}=N\times (1-\frac{1}{p_1})\times (1-\frac{1}{p_2})\times ...\times (1-\frac{1}{p_m})$

证明：

先对$N$分解质因数，得到$N=p_1^{a_1}\times p_2^{a_2}\times ...\times p_m^{a_m}$

在$[1,N]$中，我们先将所以$p_1,p_2,...,p_m$的倍数减掉，再将所有$p_i\times p_j$的倍数的数加上，再减去$p_i\times p_j\times p_k$的倍数的数，以此类推。最后剩下的就是与$N$互质的数的个数。式子表达为$N-(\frac{N}{p_1}+\frac{N}{p_2}+...+\frac{N}{p_m})+(\frac{N}{p_1{p}_2}+\frac{N}{p_1{p}_3}+...)-(\frac{N}{p_1{p}_2{p}_3}+...)...$

该式就是$N\times (1-\frac{1}{p_1})\times (1-\frac{1}{p_2})\times ...\times (1-\frac{1}{p_m})$的展开式

公式法求欧拉函数

#include <bits/stdc++.h>
using namespace std;

void solve() {
    int n;
    cin >> n;
    auto f = [&](int x) -> int {
        int res = x;
        for (int i = 2; i <= x / i; i ++) {
            if (x % i != 0) continue;
            res = res / i * (i - 1);
            while (x % i == 0) x /= i;
        }  
        if (x > 1) res = res / x * (x - 1); 
        return res;
    };
    cout << f(n) << endl;
}

int main() {
    int t;
    cin >> t;
    while (t --) {
        solve();
    }
    return 0;
}

筛法求欧拉函数

对于一个质数$p$来说，$\varphi (p)=p-1$

对于每一个数$primes[j]\times i$来说

• 当$i\%primes[j]=0$时，$i$中的质因子包括$primes[j]$，则有$\varphi (primes[j]\times i)=\varphi (i)\times primes[j]$
• 当$i\%primes[j]\ne 0$时，$i\times primes[j]$中的质因子既包括$i$中的质因子，又包括$primes[j]$，则有$\varphi (primes[j]\times i)=\varphi (i)\times primes[j]\times (1-\frac{1}{primes[j]})$

#include <bits/stdc++.h>s
using namespace std;

#define int long long

const int N = 1000010;
int cnt, primes[N];
bool st[N];
int phi[N];

int getEuler(int n) {
    phi[1] = 1;
    for (int i = 2; i <= n; i ++) {
        if (st[i] == false) {
            primes[cnt ++] = i;
            phi[i] = i - 1;
        }
        for (int j = 0; j < cnt && primes[j] <= n / i; j ++) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                break;
            }
            phi[i * primes[j]] = phi[i] * primes[j] / primes[j] * (primes[j] - 1);
        }
    }
    int res = 0;
    for (int i = 1; i <= n; i ++) {
        res += phi[i];
    }
    return res;
}

signed main() {
    int n;
    cin >> n;
    cout << getEuler(n) << endl;
    return 0;
}

质数

分解质因数

给定一个数$x$，由唯一分解定理可知，$x=p_1^{a_1}\times p_2^{a_2}\times p_3^{a_3}...(p_1<p_2<p_3)$，我们从小到大枚举每一个数，最先整除$x$的除数$n$一定是素数。同时将这个数中所有$n$的因子除干净，那么下次$x$被整除时除数也是素数，从而分解了质因数。

#include <bits/stdc++.h>
using namespace std;

int main() {
    int n;
    cin >> n;
    for (int i = 0; i < n; i ++) {
        int x;
        cin >> x;
        map<int, int> cnt;
        vector<int> res;
        for (int j = 2; j <= x / j; j ++) {
            if (x % j != 0) continue;
            res.push_back(j);
            while (x % j == 0) {
                cnt[j] ++;
                x /= j;
            }
        }
        if (x > 1) {
            res.push_back(x);
            cnt[x] ++;
        }
        for (auto& u : res) {
            cout << u << " " << cnt[u] << endl;
        }
        cout << endl;
    }
    return 0;
}

普通筛法

给定一个整数$x$，并筛掉所有倍数，如$2x,3x,4x...$

筛完后，$vis[i]=false$的为素数

这种方法的时间复杂度是$O(nlogn)$级别

对于每一个外层循环变量$i$，内层循环枚举所有不超过$n$的$i$的倍数，枚举次数是$\lfloor \frac{n}{i}\rfloor$，又$\lfloor \frac{n}{i}\rfloor <\frac{n}{i}$，有调和级数的渐近公式${\sum }_{i=1}^n\frac{1}{i}\approx ln(n+1)+0.577218$，其中$0.577218$是欧拉常数

所以复杂度是$O(nlogn)$级别

虽然这个筛法的时间复杂度很差，但是调和级数枚举引人深思，这种枚举方式非常巧妙，而且复杂度仅有$O(logn)$级别

埃式筛法

在调和级数枚举的筛法中，我们发现枚举4的倍数无意义，因为从这里筛掉的数一定从2的倍数里筛掉了，所以我们想到：只用质数去筛

其复杂度只有$O(nloglogn)$

#include <bits/stdc++.h>
using namespace std;

const int N = 1e6 + 10;
int cnt, primes[N], n;
bool st[N];

int main() {
    cin >> n;
    for (int i = 2; i <= n; i ++) {
        if (st[i]) continue;
        primes[cnt ++] = i;
        for (int j = 2 * i; j <= n; j += i) {
            st[j] = true;
        }
    }
    cout << cnt << endl;
    return 0;
}

欧拉筛

通过对埃氏筛法的分析，我们发现$6=2\times 3$会被质数$2$和$3$各筛一次。一个数有多少个质因子，它就会被筛多少次。那么我们能否让一个数只被筛掉一次呢？

欧拉筛法可以做到并且每个数会被它的最小质因子筛掉

#include <bits/stdc++.h>
using namespace std;

const int N = 1e6 + 10;
int cnt, primes[N], n;
bool st[N];

int main() {
    cin >> n;
    for (int i = 2; i <= n; i ++) {
        if (!st[i]) primes[cnt ++] = i; //是素数
        for (int j = 0; j < cnt && primes[j] <= n / i; j ++) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
    cout << cnt << endl;
    return 0;
}

• 设$x=i\times primes[j]$
• 若$imodprimes[j]\ne 0$，则$i$中没有$primes[j]$这个质因子。又因为$j$是从小到大遍历的，所以$primes[j]<i$，则$primes[j]$是$x$的最小质因子
• 若$imodprimes[j]=0$ ，说明$primes[j]$是$i$的一个质因子。又因为$j$是从小到大遍历的，所以$primes[j]$是$i$的最小质因子。所以$primes[j]$是$x$的最小质因子
• 综上所述，每个数都是被其最小质因子$primes[j]$筛掉且只会被筛掉一次

调和级数枚举模型

内容参考于：[OI&ACM]调和级数枚举倍数模型 - 知乎
给定一个数$n(n\le 10^6)$，求$[1,n]$中所有整数$i$的正因子个数

思路：枚举$[1,n]$中的所有数$i$，对$i$的所有倍数，$i$ 会增加$1$的贡献。考虑调和级数枚举：

int d[1000001] = {0}; // 储存因数倍数，初始全为 0
for (int i = 1; i <= n; i++) // 枚举 [1, n] 的数
{
    for (int j = i; j <= n; j += i) // 枚举 i 的倍数
    {
        d[j]++;
    }
}

时间复杂度为$O(nlogn)$级别

下面的题目都应用了该模型，读者不妨观看思考

题目1：Problem - D - Codeforces

大意：给定一个长度为$n(n\le 10^5)$的数组，每次操作可以任取一些位置上的数，并将该位置上的数染成黑色。然后将这些被染成黑色的位置索引的所有倍数染成绿色，每次操作的得分是该数组中黑色和绿色位置上的数的最大值。由于选择染成黑色的操作有$2^n-1$种，所以我们要求出所有操作的得分的和

选择每个数所能得到的最大值是容易处理的。

vector<int> maxv(n + 1, 0);
for (int i = 1; i <= n; i ++) {
    for (int j = i; j <= n; j += i) {
        maxv[i] = max(maxv[i], a[j]);
    }
}

我们明确$maxv$数组的含义：选择第$i$上的数，将其涂成黑色，并将其索引的倍数涂成绿色，所能得到的最大值。假设有这样的$maxv$数组，$maxv=[1,6,7,3]$，$1$贡献了$1$次，$3$贡献了$2$次，$6$贡献了$4$次，$7$贡献了$8$次。这是因为：对于第一个位置上的数，我们有如下的选择方式$[1],[1,2],[1,3],[1,4],[1,2,3],[1,2,4],[1,3,4],[1,2,3,4]([1,3,4]表示我们选择第1，3，4个数字，其余同理)$，由于$1$是最小的，在如下选择中只贡献了$1$次，就是选择方式为$[1]$的时候。所以我们可以将$maxv$数组排序并用贡献法求出答案

sort(maxv.begin() + 1, maxv.begin() + n + 1);
LL res = 0, pow2 = 1;	
for (int i = 1; i <= n; i ++) {
    res = (res % mod + maxv[i] * pow2 % mod) % mod;
    pow2 = pow2 * 2 % mod;
}

cout << res << endl;

题目2：Problem - D - Codeforces

大意：给定一个数组$a$，其长度为$n$，其中$1\le n\le 10^6$，$1\le a_i\le n$。若不存在$k，(1\le k\le n)$使得$a_k|a_i$，$a_k|a_j$，那么称$(i,j)$是对儿好整数。求好整数的数量。

• 条件$a_k|a_i$，$a_k|a_i$可以转化为$a_k|gcd(a_i,a_j)$。由于值域仅有$10^6$级别，我们可以预处理出$f(x)(1\le x\le 10^6)$，$f(x)$表示$gcd(i,j)=x$好整数的个数。

  • 令$y|gcd(a_i,a_j)$，此时$a_i$，$a_j$是$y$的倍数，故可以枚举$y$的倍数计算出数量$sum$。那么$y|gcd(a_i,a_j)$中$(i,j)$的对数是$C_{sum}^2=\frac{sum\times (sum-1)}{2}$。其中存在两种情况$y=gcd(a_i,a_j)$或$y<gcd(a_i,a_j)$，（也就是$gcd(a_i,a_j)$是$y$的大于$1$倍的倍数）。所以$f(x)=\frac{sum\times (sum-1)}{2}-f(2x)-f(3x)-...$

  for (int i = n; i >= 1; i --) {
      int s = 0, t = 0;
      for (int j = i; j <= n; j += i) {
          s += cnt[j];
          t += f[j];
      }
      f[i] = s * (s - 1) / 2 - t;
  
  }
  

• 在值域之中，我们将$a_m(1\le m\le n)$及其倍数筛掉，剩下的就是$f(x)$中$x$可以取的值。于是$ans={\sum }_{i=1}^{n}f(i)(其中i未被筛掉)$

#include <bits/stdc++.h>
using namespace std;

using LL = long long;

#define int LL

void solve() {
	int n;
	cin >> n;
	vector<int> f(n + 1, 0);
	vector<int> cnt(n + 1, 0);
	for (int i = 0; i < n; i ++) {
		int x;
		cin >> x;
		cnt[x] ++;
	}
	
	for (int i = n; i >= 1; i --) {
		int s = 0, t = 0;
		for (int j = i; j <= n; j += i) {
			s += cnt[j];
			t += f[j];
		}
		f[i] = s * (s - 1) / 2 - t;
		
	}
	
	vector<int> vis(n + 1, false);
	for (int i = 1; i <= n; i ++) {
		if (cnt[i]) {
			for (int j = i; j <= n; j += i) {
				vis[j] = true;
			}
		}
	}
	
	int res = 0;
	for (int i = 1; i <= n; i ++) {
		if (vis[i] == false) {
			res += f[i];
		}
	}
	cout << res << endl;
}

signed main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	int t;
	cin >> t;
	while (t --) {
		solve();
	}
	
	return 0;
}

约数

约数个数

对于一个正整数$x$，由于唯一分解定理，$x=p_1^{a_1}\times p_2^{a_2}\times p_3^{a_3}\times p_4^{a_4}...$

对于$x$的每一个约数$y$，$y$都能写成$y=p_1^{b_1}\times p_2^{b_2}\times p_3^{b_3}\times p_4^{b_4}...(0\le b_1\le a_1,0\le b_2\le a_2,0\le b_3\le a_3,0\le b_4\le a_4)$

所以约数个数为$(a_1+1)\times (a_2+1)\times (a_3+1)\times (a_4+1)\times ...$

约数之和

约数之和$sum=(p_1^0+p_1^1+...+p_1^{a_1})\times (p_2^0+p_2^1+...+p_2^{a_2})\times (p_3^0+p_3^1+...+p_3^{a_3})\times (p_4^0+p_4^1+...+p_4^{a_4})\times ...$

由分配律将上述括号拆开可得：$x_1+x_2+x_3+...+$，其中每一项都是$x$的约数

#include <bits/stdc++.h>
using namespace std;

using LL = long long;

const int mod = 1e9 + 7;

int n;
map<int, int> cnt;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    
    cin >> n;
    for (int i = 0; i < n; i ++) {
        int x;
        cin >> x;
        for (int j = 2; j <= x / j; j ++) {
            if (x % j != 0) continue;
            int t = 0;
            while (x % j == 0) {
                t ++;
                x /= j;
            }
            cnt[j] += t;
        }
        if (x > 1) cnt[x] ++;
    }   
    LL res = 1;
    for (auto& [u, v] : cnt) {
        LL t = 0;
        for (int j = 0; j <= v; j ++) {
            t = (t * u % mod + 1) % mod;
        }
        res = res * t % mod;
    }
    cout << res << endl;
    return 0;
}

组合数

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

可由此式子递推计算组合数

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

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

可以预处理出来$fact[N]$与$infact[N]$

故$C_n^m=fact[N]\times infact[m]\times infact[n-m]$

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
#define endl "\n"
typedef pair<int, int> pII;
const int p = 1e9 + 7;
const int N = 100010;
LL fact[N], infact[N];

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

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

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

void solve() {
	int n, m; cin >> n >> m;
	cout << (fact[n] % p) * (infact[n - m] % p * infact[m] % p) % p << endl;
}


int main() {
	ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
	
	init(100005);
	int t; cin >> t;
	while (t --) {
		solve();
	}
	
	return 0;
}

卢卡斯定理

$C_n^m=C_{\frac{n}{p}}^{\frac{m}{p}}\times C_{nmodp}^{mmodp}$