数论基础知识和模板

质数筛

用于快速处理 1～n 中所有素数的算法

因为依次遍历判断每一个数是否质数太慢，所以把一些明显不能质数的筛出来

普通筛法，对于每个整数，删除掉其倍数。

bool vis[N];//0表示是质数
int pri[N],o; //质数表
void get(int n) {
    vis[1] = true; //1不是素数
    for(int i = 2; i <= n; ++i) {//从2开始
        if(!vis[i])
            pri[++o] = i; //加入质数表
        for(int j = 2 * i; j <= n; j += i)
            vis[j] = 1;//删除到n的整数倍
    }
}

线性筛和埃氏筛-CSDN博客

埃氏筛

在普通筛法的过程中，发现筛 4 的倍数毫无意义，因为 4 的倍数肯定已被 2 的倍数筛过了，所以优化方法是只用质数去筛。

 对于每一个整数，只有他是质数时才去筛

if(!vis[i]){
            primes[ ++ tot] = i; 
            for(int j = 2 * i; j <= n; j += i)
                vis[j] = 1;//进入if里了
        }

线性筛（欧拉

让每个数都被其最小的质因数筛掉，故一个数只被筛一次，复杂度为O(n)

实现就是，对于每个数，我们都将他对整个现在已有的质数表筛一遍，比如3是质数，那么把已有的2相乘，筛掉6，但是并不往下进行，因为比如12，最小质数因数是3，在遍历到4时自然而然就筛掉了，这样每个数就只会被他最小质因数筛掉一次，达到了线性

for(int i = 2; i <= n; ++i) {
    if(!vis[i]) //是素数
        pri[ ++ o] = i;
   
    for(int j=1;j<o;j++)
        if(i*pri[j]<=n)vis[i*pri[j]]=1;
        else break;//
    }
}

快速幂：

1. 基本概念

一个在O(logn)的时间内计算an的魔术技巧小技巧，而暴力的计算需要O(n)的时间，基本思想就是将指数按照二进制表示来分割。例如

ll ksm(ll a,ll b,ll p){
    ll ans=1;
    while(b){
        if(b&1) ans=ans*a%p;

//b & 1：用按位与运算判断 b 的二进制最低位是否为 1（比如 b=5 是 101，b&1 结果为 1；b=4 是 100，b&1 结果为 0 ）。
如果最低位是 1，说明当前这一位对应的幂次需要乘到结果里。
ans = ans * a % p：把当前的 a（对应某一 “二进制位权” 的幂）乘到结果 ans 中，再对 p 取模，保证数值不会溢出，且符合题目对结果取模的要求。

        a=a*a%p;//每次循环，底数 a 自乘一次（即 a = a² % p ），对应指数二进制位的位权翻倍。
比如初始 a=3（对应 3 1），第一次自乘后 a=3²=9（对应 3 2 ），第二次自乘后 a=9²=81（对应 3 4），第三次自乘后 a=81²=6561（对应 3 8 ）…… 以此类推，刚好匹配指数二进制各位的权值

        b>>=1;//后移
    }
    return ans;
}

• a：底数，即要做幂运算的基数。
• b：指数，即幂次。
• p：模数，计算结果要对 p 取模，避免数值过大溢出或满足题目对结果范围的要求。

2. 矩阵快速幂

矩阵快速幂 = 矩阵乘法 + 快速幂。

单位矩阵是对角线全为 1，其他位置均为 0 的矩阵，记为 I，性质：A×I=A。

include <iostream>
#include <cstring>
using namespace std;

typedef long long LL;         // 使用LL作为long long的别名，简化代码
const int MOD = 1e9 + 7;      // 定义模数，防止整数溢出
const int N = 105;            // 定义矩阵的最大尺寸

int n;                        // 矩阵的实际大小
LL k;                         // 矩阵需要乘方的次数

// 矩阵结构体，封装矩阵乘法和单位矩阵设置操作
struct Matrix {
    LL a[N][N];               // 存储矩阵元素
    
    // 默认构造函数，初始化矩阵为全0
    Matrix() {
        memset(a, 0, sizeof a);
    }
    
    // 重载乘法运算符，实现矩阵乘法
    Matrix operator*(const Matrix& b) const {
        Matrix c;
        // 三层循环实现矩阵乘法，注意取模防止溢出
        for (int i = 1; i <= n; i++)
            for (int k = 1; k <= n; k++)
                for (int j = 1; j <= n; j++)
                    c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j]) % MOD;
        return c;
    }
    
    // 设置当前矩阵为单位矩阵（对角线为1，其余为0）
    void setIdentity() {
        for (int i = 1; i <= n; i++) a[i][i] = 1;
    }
};

int main() {
    cin >> n >> k;  // 输入矩阵大小n和幂次k
    
    // 读取原始矩阵
    Matrix base;
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++)
            cin >> base.a[i][j];
    
    // 初始化结果矩阵为单位矩阵，相当于数值计算中的1
    Matrix res;
    res.setIdentity();
    
    // 快速幂算法：计算矩阵的k次幂
    // 原理与数值快速幂相同，但使用矩阵乘法代替数值乘法
    while (k) {
        if (k & 1) res = res * base;  // 如果当前位为1，则乘上当前的base
        base = base * base;            // base平方
        k >>= 1;                       // k右移一位，相当于除以2
    }
    
    // 输出结果矩阵
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++)
            cout << res.a[i][j] << ' ';
        cout << endl;
    }
    
    return 0;
}

  扩展欧几里得算法

 一、核心目标：求解 ax + by = gcd(a, b) 的整数解

扩展欧几里得算法本质是 在求最大公约数（gcd）的同时，找到满足 ax + by = gcd(a, b) 的整数 x 和 y 。
比如：a=5, b=3，gcd 是 1，同时能找到 x=2, y=-3（因为 5×2 + 3×(-3) = 10-9=1 = gcd(5,3) ）。

 二、递归推导：从 gcd(b, a%b) 反推 gcd(a, b) 的解

欧几里得算法的核心是 gcd(a, b) = gcd(b, a%b)（比如 gcd(5,3)=gcd(3,2)=gcd(2,1)=gcd(1,0)=1 ）。
扩展欧几里得在此基础上，递归地从 gcd(b, a%b) 的解，反推 gcd(a, b) 的解，

ll exgcd(ll a, ll b, ll &x, ll &y){ // 传引用修改 x、y
    if(b == 0){ 
        // 递归出口：b=0时，gcd是a，解为x=1, y=0
        x = 1; 
        y = 0; 
        return a; 
    }
    // 递归求解 gcd(b, a%b)，同时得到子问题的解 (x', y')
    ll d = exgcd(b, a % b, y, x); // 注意：这里交换了 y、x 的位置！
    // 子问题的解是 (x', y')，但因为交换了参数，实际拿到的是 y=x', x=y' 
    // 根据推导公式，当前层的 y = x' - (a/b)*y' 
    y -= (a / b) * x; 
    return d; // 返回 gcd(a,b)
}

以 a=5, b=3 为例，逐步看递归和求解过程：

1. 第一层调用：exgcd(5, 3, x, y)

   • b≠0，递归调用 exgcd(3, 5%3=2, y, x)

2. 第二层调用：exgcd(3, 2, y, x)

   • b≠0，递归调用 exgcd(2, 3%2=1, y, x)

3. 第三层调用：exgcd(2, 1, y, x)

   • b≠0，递归调用 exgcd(1, 2%1=0, y, x)

4. 第四层调用：exgcd(1, 0, y, x)

   • b=0，触发递归出口：x=1, y=0，返回 d=1（即 gcd=1 ）

5. 返回第三层：

   • 此时 d=1，且因为参数交换，当前层的 y 是子问题的 x=1，x 是子问题的 y=0
   • 执行 y -= (2/1)*x → y = 0 - 2*1 = -2
   • 返回 d=1，当前层的解是 x=0, y=-2 不对？别急，继续看… 其实这里的 x,y 是相对于当前层 a=2, b=1 的， 哦，因为参数交换后，逻辑需要再仔细看：

   实际第三层的调用是 exgcd(2, 1, y, x)，所以子问题（第四层）的 x=1, y=0 会被赋值给 当 前层的 y 和 x ，即：

   • 当前层（第三层）的 y = 子问题的 x=1
   • 当前层（第三层）的 x = 子问题的 y=0
     然后执行 y -= (2/1)*x → y = 1 - 2*0 = 1
     所以第三层的解是 x=0, y=1，对应等式 2*0 + 1*1 = 1 = gcd(2,1)，正确！

6. 返回第二层：

   • 第二层调用是 exgcd(3, 2, y, x)，子问题（第三层）返回 d=1，且当前层的 y = 第三层的 x=0，x = 第三层的 y=1
   • 执行 y -= (3/2)*x → 3//2=1，所以 y = 0 - 1*1 = -1
   • 此时第二层的解是 x=1, y=-1，对应等式 3*1 + 2*(-1) = 3-2=1 = gcd(3,2)，正确！

7. 返回第一层：

   • 第一层调用是 exgcd(5, 3, x, y)，子问题（第二层）返回 d=1，且当前层的 y = 第二层的 x=1，x = 第二层的 y=-1
   • 执行 y -= (5/3)*x → 5//3=1，所以 y = 1 - 1*(-1) = 2
   • 此时第一层的解是 x=-1, y=2，对应等式 5*(-1) + 3*2 = -5+6=1 = gcd(5,3)，正确！

最终，exgcd(5,3,x,y) 执行后，x=-1, y=2，满足 5*(-1) + 3*2 = 1 。

乘法逆元

1. 基本概念

更准确的来说是模意义下的乘法逆元。单位元：在一个集合中，对于某种运算∗，如果对于任何的集合元素 a，和元素 e 运算，得到还是集合元素 a 本身，则称 e 为这个运算下的单位元。例如加法运算的单位元是 0，乘法的单位元是 1。逆元：在一个集合中，对于某种运算∗，如果任意两个元素的运算结果等于单位元，则称这两个元素互为逆元。加法中 a 的逆元是 - a，乘法中 a 的逆元是a1​即a−1。所以，数学上的乘法逆元就是直观上的倒数，即ax=1，x 即为 a 的逆元。对于模意义下的乘法逆元，即满足ax≡1(modb)时，则称 x 为 a 关于模 b 的逆元。

很容易看出来，这是扩展欧几里得的一种运用 

扩展欧几里得法

ax≡1(modb)⇒ax+by=1

 但是利用以上的欧几里得，我们只能求出最近的一个乘法逆元，比如3m10，得到-3；

此时，应该进行魔术技巧

x = (x % n + n) % n

操作，把不论正负，求取最小正数逆元

费马小定理法

定义：若 p 为质数，gcd(a,p)=1，则ap−1≡1 (mod p)。

证明：

• 因为ax≡1(modb)
• 所以ax≡ab−1(modb)
• 故x≡ab−2(modb)

然后就可以用快速幂求解，时间复杂度为O(logb)。 

ll inv(ll a,ll p){
    return ksm(a,p-2,p);//ksm上面呢
}

线性求逆元

求出 1-n 中每个数的逆元。

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

线性求逆元的递推公式

对于任意一个数 i，它的逆元 inv[i] 可以通过如下递推公式计算得出

inv[i]=(p−⌊p/i⌋)×inv[p%i]%p

推导：

        我们令 k=⌊p/i⌋ 以及 r=p%i，那么显然有 p=k×i+r，也就是 k×i+r≡0(modp)。

        对这个同余式进行变形，在等式两边同时乘以 i−1×r−1，就可以得到 k×r−1+i−1≡0(modp)。

        进一步移项，就能得到 i−1≡−k×r−1(modp)。

        为了保证逆元是正数，我们把公式调整为 i−1≡(p−k)×r−1(modp)，这里的 r 其实就是 p%i。

递推的起始条件
当 i 等于 1 时，它的逆元毫无疑问是 1，即 inv[1]=1，这是整个递推过程的起点。

组合数学 

基本求法

ll res[N][N];//记忆化，避免重复计算
ll C(ll n,ll m){
    if(m ==0 || n == m) return 1;
    if(res[n][m])
        return res[n][m];//调用记忆
    return res[n][m] = C(n - 1, m) + C(n - 1, m - 1);
}

也可以将记忆化演变成动态规划

#include<bits/stdc++.h>
#define ll long long
using namespace std;
int main(){
    ll n, m, p;
    cin>>n>>m>>p;
    ll total = n + m;
    vector<vector<ll>> C(total + 1, vector<ll>(total + 1, 0));        
    // 预处理组合数
    for(int i = 0; i <= total; i++){
        C[i][0] = 1 % p;
        C[i][i] = 1 % p;
        for(int j = 1; j < i; j++){
            C[i][j] = (C[i-1][j] + C[i-1][j-1]) % p;
        }
                
    cout<<C[total][n]<<endl;
    return 0;
}

ps:面对过大数据空间会炸 

Lucas 定理 

ll Lucas(ll n, ll m){
    if(m == 0) return 1;
    return Lucas(n / p, m / p) * C(n % p, m % p) % p;
}

 逆元法求组合（利用费马

由于公式

我们想，能不能直接利用阶乘自己取模后相除

但是当需要计算 C(n,m)modp 时，不能直接对分子和分母分别取模后相除，因为模运算的除法不满足分配律

逆元的作用是将模意义下的除法转换为乘法

（ps，-1不是次方，是逆元的意思）

ll comb(ll n, ll m, ll p) {//这里n上，m下
    if (m < 0 || m > n) return 0;
    
    // 预处理阶乘数组
    vector<ll> fact(n + 1, 1);
    for (ll i = 1; i <= n; i++) {
        fact[i] = fact[i - 1] * i % p;
    }
    
    // 计算逆元：fact[m]^(-1) 和 fact[n-m]^(-1)
    ll inv_m = mod_pow(fact[m], p - 2, p);
    ll inv_nm = mod_pow(fact[n - m], p - 2, p);//费马
    
    // 计算组合数
    return fact[n] * inv_m % p * inv_nm % p;
}

 卢卡斯与逆元结合

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

typedef long long ll;

// 快速幂计算 a^b % p
ll quick_pow(ll a, ll b, ll p) {
    ll res = 1;
    while (b > 0) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

// 计算组合数 C(a, b) % p，其中 a < p，b < p
ll comb(ll a, ll b, ll p) {
    if (b < 0 || b > a) return 0;
    if (b == 0 || b == a) return 1;
    // 预处理阶乘和逆元
    vector<ll> fact(p), inv_fact(p);
    fact[0] = 1;
    for (int i = 1; i < p; ++i) {
        fact[i] = fact[i - 1] * i % p;
    }
    inv_fact[p - 1] = quick_pow(fact[p - 1], p - 2, p);
    for (int i = p - 2; i >= 0; --i) {
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % p;
    }
    return fact[a] * inv_fact[b] % p * inv_fact[a - b] % p;
}

// 卢卡斯定理计算 C(n, m) % p
ll lucas(ll n, ll m, ll p) {
    if (m == 0) return 1;
    return comb(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    
    int T;
    cin >> T;
    while (T--) {
        ll n, m, p;
        cin >> n >> m >> p;
        // 计算 C(n+m, n) % p
        cout << lucas(n + m, n, p) << '\n';
    }
    
    return 0;
}

卡特兰数 

 catalan数-CSDN博客

排列组合与容斥原理

1. 排列数公式

• 从 n 个元素中选 m 个排列：Anm​=n×(n−1)×⋯×(n−m+1)=n!/(n−m)!​。

2. 容斥原理基础

• 两个集合：∣A∪B∣=∣A∣+∣B∣−∣A∩B∣。
• 三个集合：∣A∪B∪C∣=∣A∣+∣B∣+∣C∣−∣A∩B∣−∣A∩C∣−∣B∩C∣+∣A∩B∩C∣。
• 应用场景：计算不满足某些条件的元素个数，如求 1~n 中不被 2、3、5 整除的数的个数。

快速傅里叶变换（FFT）

1. 核心作用

• 加速多项式乘法，将暴力O(n2)的复杂度优化到O(nlogn)。
• 原理：利用复数单位根的性质，将多项式从系数表示转为点值表示，相乘后再逆变换回系数表示。

2. 基本步骤

1. 预处理单位根：计算复数域上的 n 次单位根。
2. FFT 正变换（DFT）：将多项式转换为点值形式。
3. 点值相乘：对应点值相乘得到结果多项式的点值表示。
4. 逆变换（IDFT）：转换回系数形式

 框架（搞不懂，看不懂，不明白，所以就贴了个码）

const double PI = acos(-1);
struct Complex {
    double x, y;
    Complex(double x=0, double y=0) : x(x), y(y) {}
};
Complex operator+(Complex a, Complex b) { return Complex(a.x+b.x, a.y+b.y); }
Complex operator-(Complex a, Complex b) { return Complex(a.x-b.x, a.y-b.y); }
Complex operator*(Complex a, Complex b) { return Complex(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }

void fft(Complex *a, int n, int inv) {
    if (n == 1) return;
    Complex a1[n/2], a2[n/2];
    for (int i = 0; i < n; i++) {
        if (i % 2 == 0) a1[i/2] = a[i];
        else a2[i/2] = a[i];
    }
    fft(a1, n/2, inv); fft(a2, n/2, inv);
    Complex w(1, 0), wn(cos(2*PI/n), inv*sin(2*PI/n));
    for (int i = 0; i < n/2; i++) {
        a[i] = a1[i] + w * a2[i];
        a[i+n/2] = a1[i] - w * a2[i];
        w = w * wn;
    }
    if (inv == -1) {
        for (int i = 0; i < n; i++) a[i].x /= n;
    }
}

// 计算多项式乘法c = a * b，a和b的次数分别为n和m
void multiply(ll *a, ll *b, ll *c, int n, int m) {
    int len = 1;
    while (len < n + m) len <<= 1;
    Complex *fa = new Complex[len], *fb = new Complex[len];
    for (int i = 0; i < n; i++) fa[i] = Complex(a[i], 0);
    for (int i = n; i < len; i++) fa[i] = Complex(0, 0);
    for (int i = 0; i < m; i++) fb[i] = Complex(b[i], 0);
    for (int i = m; i < len; i++) fb[i] = Complex(0, 0);
    fft(fa, len, 1); fft(fb, len, 1);
    for (int i = 0; i < len; i++) fa[i] = fa[i] * fb[i];
    fft(fa, len, -1);
    for (int i = 0; i < n + m; i++) c[i] = (ll)(fa[i].x + 0.5);
    delete[] fa; delete[] fb;
}