【整理】算法竞赛数学相关学习笔记

部分内容来自牛客相关课程课件

这个笔记想写得稍微详细一点，部分内容来自我自己之前的博客: OI 里的数学内容整理（提高组）。

其中一些公式是 AI 从 PPT 上识别的，格式混乱，还请谅解。

整数分解与筛法

欧几里得算法

用于求解两个数的最大公因数。

$gcd (a, b) = gcd (b, a mod b)$ 。

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

扩展欧几里得算法与裴蜀定理

用于求解 $a x + b y = gcd (a, b), a, b \in N$ 的一组整数解 $(x, y)$ 。

 int exgcd(int a,int b,int &x,int &y){//注意：结果x,y可能是负数,但是输入的a,b需要均为正数！！！
    if(b==0){x=1;y=0;return a;}
    int GCD=exgcd(b,a%b,x,y);
    int tmp;tmp=x;
    x=y;y=tmp-a/b*y;
    return GCD;
}

也可以用于求逆元：

求 $a (\mod b)$ 意义下的逆元等价于求一个 $x$ 满足 $a x \equiv 1 (\mod b)$ ，由于逆元存在， $a$ 与 $b$ 互质，故 $gcd (a, b) = 1$ ，所以即求 $a x + b y = 1 = gcd (a, b)$ 的解 $x$ ，显然可以用扩欧解决。

例题：青蛙的约会

裴蜀定理（前置定理）：对于 $a x + b y = gcd (a, b), a, b \in N$ ，一定存在一组整数 $x, y$ 使得该式成立。

裴蜀定理扩展：对于 $a_{1} x_{1} + a_{2} x_{2} \dots + a_{n} x_{n} = gcd (a_{1}, a_{2} \dots a_{n}), a_{i} \in N$ ，一定存在一组整数 $x_{1}, x_{2} \dots x_{n}$ 使得该式成立。

算术基本定理（整数的分解）

任何一个大于 $1$ 的正整数 $N$ 都能唯一分解为有限个质数的乘积，即： $N = p_{1}^{c_{1}} p_{2}^{c_{2}} \dots p_{m}^{c_{m}}$ 。

埃氏筛

 bool not_prime[N];
int tot;
int pri[N];
void Sieve_of_Eratosthenes(int nn){
    for(int i=2;i<=nn;i++){
        if(!not_prime[i]){
            pri[++tot]=i;
            for(int j=2*i;j<=n;j+=i){
                not_prime[j]=1;
            }
        }
    }
}

时间复杂度： $O (n \log (\log n))$ 。

例题：Prime Distance

欧拉筛

 bool not_prime[N];
int tot;
int pri[N],phi[N];
void Euler_Sieve(int nn){
    for(int i=2;i<=nn;i++){
        if(!not_prime[i]) pri[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&i*pri[j]<=nn;j++){
            not_prime[pri[j]*i]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
}

时间复杂度： $O (n)$ 。其中的 $p h i$ 数组表示欧拉函数，将在后面的内容中提及。

例题：Sum of Consecutive Prime Numbers

质因数分解

暴力 $O (\sqrt{n})$ ，预处理所有质数可以做到单次询问 $O (\frac{\sqrt{n}}{\log n})$ 。

例题：Prime Land

Stein 算法（大整数的 GCD）

如果 $a = 0$ 或 $b = 0$ ，则 $(0, b) = b, (a, 0) = a$
如果 $a, b$ 均为偶数，则 $(a, b) = (\frac{a}{2}, \frac{b}{2})$
如果 $a$ 是偶数， $b$ 是奇数，则 $(a, b) = (\frac{a}{2}, b)$
如果 $a$ 是奇数， $b$ 是偶数，则 $(a, b) = (a, \frac{b}{2})$
如果 $a, b$ 均为奇数，不妨设 $a > b$ 则 $(a, b) = (b, a - b)$

时间复杂度为 $O ((\log (max a, b))^{2})$ ，这个算法的应用场合是高精度，因为高精度取模速度很慢，而减法和除以 $2$ 则较快。

【补充】数论分块

对于一些整除类问题，如果出现了 $⌊ \frac{n}{i} ⌋$ 型，可以考虑通过数论分块将计算范围缩减到 $\sqrt{n}$ 级别。

 void division_block(ll n){
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        //do something
        //在 l~r 这个范围中，所有的 n/i 得到的数值都是一样的
    }
}

同余与模

逆元

若 $a b \equiv 1 (\mod m)$ ，则称 $a, b$ 在 $\mod m$ 意义下互为逆元。

可以通过扩展欧几里得算法求出。若 $m$ 是质数，也可以通过费马小定理求出（稍后介绍）。

但如果要求出 $n$ 个数的逆元，有以下线性求法：

 ll inv[N];
void init_inv(int nn,int p){
    inv[1]=1;
    for(int i=2;i<=nn;i++){
        inv[i]=(ll)(p-p/i)*inv[p%i]%mod;
    }
}

该算法利用了 $⌊ \frac{p}{i} ⌋ \times i + p % i = p$ 即 $⌊ \frac{p}{i} ⌋ \times i + p % i \equiv 0 (\mod p)$ 这一公式推导而得。

推导过程：

$⌊ \frac{p}{i} ⌋ \times i + p % i \equiv 0 (\mod p)$
$⌊ \frac{p}{i} ⌋ \times i \times (p % i)^{- 1} + 1 \equiv 0 (\mod p)$
$⌊ \frac{p}{i} ⌋ \times (p % i)^{- 1} + i^{- 1} \equiv 0 (\mod p)$
$i^{- 1} \equiv - ⌊ \frac{p}{i} ⌋ \times (p % i)^{- 1} (\mod p)$
$i^{- 1} \equiv (p - ⌊ \frac{p}{i} ⌋) \times (p % i)^{- 1} (\mod p)$

模板题：乘法逆元

例题：Alternating Sum

欧拉函数

记 $φ (n)$ 为 $\leq n$ 的正整数中与 $n$ 互质的数的个数（因此 $φ (1) = 1$ ）。

通式： $φ (n) = n \prod_{i = 1}^{n} (1 - \frac{1}{p_{i}})$ ，其中 $p_{1}, p_{2} \dots p_{n}$ 为 $n$ 的所有不重复的质因数。

记 $n = p_{1}^{α_{1}} p_{2}^{α_{2}} \dots p_{k}^{α_{k}}$ ，则也有 $φ (n) = (p_{1} - 1) p_{1}^{α_{1} - 1} \times (p_{2} - 1) p_{2}^{α_{2} - 1} \times \dots \times (p_{k} - 1) p_{k}^{α_{k} - 1}$ 。

容易发现，当 $p$ 为一个质数时，有 $φ (p) = p - 1$ 。

可以通过欧拉筛线性求出所有欧拉函数（见上文），但如果只需要求一项欧拉函数，也可以暴力求。

 ll phi(int x){
    int ret=1,ori=x;
    for(int i=2;i*i<=ori;i++){
        int cnt=0;
        while(x%i==0){
            x/=i;cnt++;
            if(cnt==1) ret*=(i-1);
            else ret*=i;
        }
    }
    if(x>1) ret*=(x-1);
    return ret;
}

欧拉定理和费马小定理

欧拉定理：对于任意互素的 $a$ 和 $p$ ，有 $a^{φ (p)} \equiv 1 (\mod p)$ 。

欧拉定理推论：若 $a$ 和 $p$ 互素，则有 $a^{b} \equiv a^{b \mod φ (p)} (\mod p)$ 。

费马小定理：若 $p$ 为素数且 $a$ 和 $p$ 互素，则有 $a^{p - 1} \equiv 1 (\mod p)$ 。

费马小定理求逆元：若 $p$ 为质数，由 $a^{p - 1} \equiv 1 (\mod p)$ ，故 $a \times a^{p - 2} \equiv 1 (\mod p)$ ，故 $a^{p - 2}$ 即为 $a (\mod p)$ 意义下的逆元，显然可以快速幂解决。

同余方程组与扩展中国剩余定理

求解 $x$ ，满足方程组 $x \equiv a_{i} (\mod m_{i})$ ，这就是同余方程组。

利用 $x = y_{1} m_{1} + a_{1} = y_{2} m_{2} + a_{2}$ ，结合扩展欧几里得算法可以完成同余方程组的合并，建议读者自己推导一遍式子。

 int msc(int a,int b,int mod){//慢速乘，用于防止乘法溢出
    if(a<0&&b<0) a=-a,b=-b;
    else if(b<0) swap(a,b);
    int ret=0;
    while(b){
        if(b&1) ret=(ret+a)%mod;
        a=(a+a)%mod;b>>=1;
    }
    return ret;
}
int gcd(int a,int b){
    return b==0?a:gcd(b,a%b);
}
int lcm(int a,int b){
    return a/gcd(a,b)*b;
}
int exgcd(int a,int b,int &x,int &y){//注意：结果x,y可能是负数,但是输入的a,b需要均为正数！！！
    if(b==0){x=1;y=0;return a;}
    int GCD=exgcd(b,a%b,x,y);
    int tmp;tmp=x;
    x=y;y=tmp-a/b*y;
    return GCD;
}
void merge_equation(int &A,int &B,int a,int b){
    int y1,y2;
    int d=exgcd(B,b,y1,y2);
    if((a-A)%d!=0){A=-1;return;}
    y1=msc(y1,(a-A)/d,b/d);
    int m=lcm(B,b);
    A=((B*y1+A)%m+m)%m;
    B=m;
}
pi excrt(int nn,int a[],int b[]){//x=a(mod b)
    int A,B;
    for(int i=1;i<=nn;i++){
        if(i==1) A=a[i],B=b[i];
        else merge_equation(A,B,a[i],b[i]);
        if(A==-1) return mk(-1,-1);
    }
    return mk(A,B);
}

模板题：扩展中国剩余定理

例题：Biorhythms

中国剩余定理

求解 $x$ ，满足方程组 $x \equiv a_{i} (\mod m_{i})$ ，其中 $m_{i}$ 两两互质。

令 $M = m_{1} m_{2} \dots m_{n}, M_{i} = M / m_{i}$ 。
显然 $(M_{i}, m_{i}) = 1$ ，所以 $M_{i}$ 关于模 $m_{i}$ 的逆元存在。把这个逆元设为 $t_{i}$ 。
于是有： $M_{i} t_{i} \equiv 1 (\mod m_{i}), M_{i} t_{i} \equiv 0 (\mod m_{j}) (j \neq i)$ 。
进一步： $a_{i} M_{i} t_{i} \equiv a_{i} (\mod m_{i}), a_{i} M_{i} t_{i} \equiv 0 (\mod m_{j}) (j \neq i)$ 。

因此解为 $x \equiv a_{1} M_{1} t_{1} + a_{2} M_{2} t_{2} + \dots + a_{n} M_{n} t_{n} (\mod M)$ ，且在 $mod M$ 意义下是唯一解。

扩展欧拉定理

$a^{b} \equiv {\begin{cases} a^{b} & b < φ (m) \\ a^{b \mod φ (m) + φ (m)} & b \geq φ (m) \end{cases} (\mod m)$

例题：Power Tower

在题目中，你可以通过修改快速幂的代码来帮助使用扩展欧拉定理，具体如下：

 ll ksm(ll a,ll b,int p){
    ll ret=1;
    while(b){
        if(b&1){ret=ret*a;if(ret>=p) ret=ret%p+p;}
        a=a*a;if(a>=p) a=a%p+p;
        b>>=1;
    }
    return ret;
}

该代码计算的是 ${\begin{cases} a^{b} \mod p & a^{b} < p \\ a^{b} \mod p + p & a^{b} \geq p \end{cases}$ 的值。

威尔逊定理

$p$ 是质数等价于 $(p - 1)! \equiv - 1 (\mod p)$ 。即两者互为充要条件。

排列组合与容斥原理

约定符号

下降幂： $n^{\underset{―}{r}} = n \cdot (n - 1) \cdot (n - 2) \cdot \dots \cdot (n - r + 1)$
上升幂： $n^{\bar{r}} = n \cdot (n + 1) \cdot (n + 2) \cdot \dots \cdot (n + r - 1)$
阶乘： $n! = n \cdot (n - 1) \cdot \dots \cdot 1, 0! = 1$

注： $n^{\underset{―}{r}} = (n - r + 1)^{\bar{r}}$

加法和乘法原理

加法原理：做某件事情有几种选择，每种选择的方案数之和就是做这件事情的方案数。
乘法原理：做某件事情分为几步，每步的方案数是独立的，则它们的积就是做这件事情的方案数。

排列与组合

排列数： $P_{n}^{m} / A_{n}^{m} = \frac{n!}{m!}$

组合数： $C_{n}^{m} / (\binom{n}{m}) = \frac{n!}{(n - r)! r!}$

容斥原理

基础形式

设 $S$ 是一个有限集， $A_{1}, A_{2}, \dots, A_{n}$ 是 $S$ 的 $n$ 个子集，则
$| S - ⋃_{i = 1}^{n} A_{i} | = \sum_{i = 0}^{n} (- 1)^{i} \sum_{1 \leq j_{1} < j_{2} < \dots < j_{i} \leq n} | ⋂_{k = 1}^{i} A_{j_{k}} |$

符号形式

设 $S$ 是一个有限集， $a_{1}, a_{2}, \dots, a_{n}$ 是 $n$ 种性质。

记 $N (a_{i})$ 为 $S$ 中有 $a_{i}$ 性质的元素的数量。特殊的，记 $N (1) = | S |$ 。
记 $N (1 - a_{i})$ 为 $S$ 中没有 $a_{i}$ 性质的元素的数量。
$N (a_{i_{1}} a_{i_{2}} \dots a_{i_{k}})$ 为 S 中同时有 $a_{i_{1}}, a_{i_{2}}, \dots, a_{i_{k}}$ 性质的元素的数量。
记 $N (a \pm b) = N (a) \pm N (b)$ 。

则容斥原理可以写成
$N ((1 - a_{1}) (1 - a_{2}) \dots (1 - a_{n})) = \sum_{i = 0}^{n} (- 1)^{i} \sum_{1 \leq j_{1} < j_{2} < \dots < j_{i} \leq n} N (a_{j_{1}} a_{j_{2}} \dots a_{j_{i}})$
如果需要有其中一些性质，而不需要其他的性质，可以写作
$N (a_{1} \dots a_{x} (1 - a_{x + 1}) \dots (1 - a_{x + n})) = \sum_{i = 0}^{n} (- 1)^{i} \sum_{x + 1 \leq j_{1} < j_{2} < \dots < j_{i} \leq x + n} N (a_{1} \dots a_{x} a_{j_{1}} \dots . a_{j_{i}})$
例题：大水题

最后的晚餐 (dinner)

第二类斯特林数

有 $n$ 个互不相同的球，放到 $k$ 个互不区分的盒子里，每个盒子里至少要有一个球，求方案数。

记号 $S_{2} (n, k)$ , ${\binom{n}{m}}$
递推式 ${\binom{n}{k}} = {\binom{n - 1}{k - 1}} + k {\binom{n - 1}{k}}$
容斥 ${\binom{n}{k}} = \frac{1}{k!} \sum_{i = 0}^{k} (- 1)^{i} (\binom{k}{i}) (k - i)^{n}$
重要的公式 $n^{m} = \sum_{k = 0}^{m} {\binom{m}{k}} n^{\underset{―}{k}}$

积性函数

定义与求法

积性函数如果函数 $f : N \to R$ , 满足对于任意一对互质的正整数 $p, q$ , 都有 $f (p q) = f (p) f (q)$ , 则称 $f$ 为积性函数。

使用质因数分解求解单个 $f (n)$

 int calc_f(int x,int y){//计算 f(x^y)
    //return y+1;
}
ll get_f(ll nn) {
    ll ans=1;
    for (int i=2; i*i<=nn;i++){
        int cnt=0;
        while (nn%i==0) cnt++,nn/=i;
        if(cnt!=0) ans*=calc_f(i,cnt);
    }
    if (nn>1) ans*=f(nn,1);
    return ans;
}

使用欧拉筛求解多个 $f (n)$

 int f[N];
bool not_prime[N];
int tot;
int pri[N],cnt[N];//cnt 表示最小的质因子的次数
int calc_f(int x,int y){//计算 f(x^y)
    //return y+1;
}
void Euler_Sieve(int nn){
    f[1]=1;//根据情况定初始值
    for(int i=2;i<=nn;i++){
        if(!not_prime[i]) pri[++tot]=i,f[i]=calc_f(i,1),cnt[i]=1;
        for(int j=1;j<=tot&&i*pri[j]<=nn;j++){
            not_prime[pri[j]*i]=1;
            if(i%pri[j]==0){
                cnt[i*pri[j]]=cnt[i]+1;
                f[i*pri[j]]=f[i]/calc_f(pri[j],cnt[i])*calc_f(pri[j],cnt[i]+1);
                break;
            }
            cnt[i*pri[j]]=1;
            f[i*pri[j]]=f[i]*calc_f(pri[j],1);
        }
    }
}

模板题：【线性筛求积性函数】

例题：华华给月月出题

莫比乌斯反演

前置定义（莫比乌斯函数）： $μ (n)$

$μ (n) = {\begin{cases} (- 1)^{k}, & if n 无平方因子且有 k 个质因子 \\ 0, & if n 有平方因子 \end{cases}$

举例： $μ (1) = 1, μ (2) = - 1, μ (4) = μ (8) = 0$ .

莫比乌斯反演形式一

设 $f : N \to R, g : N \to R$ 是两个函数，则
$f (n) = \sum_{d ∣ n} g (d) \Leftrightarrow g (n) = \sum_{d ∣ n} μ (\frac{n}{d}) f (d)$

莫比乌斯反演形式二

设 $f : N \to R, g : N \to R$ 是两个函数，且存在正整数 $N$ , 对于所有 $n > N$ , $f (n) = g (n) = 0$ , 则
$f (n) = \sum_{\begin{matrix} n | m \\ m \leq N \end{matrix}} g (m) \Leftrightarrow g (n) = \sum_{\begin{matrix} n | m \\ m \leq N \end{matrix}} μ (\frac{m}{n}) f (m)$
例题：序列

Sum of gcd of Tuples (Hard)

常见的积性函数

单位函数： $ϵ (n) = [n = 1]$
常数函数： $I (n) = 1$
恒等函数： $I d (n) = n$
幂函数： $I d_{k} (n) = n^{k}$
约数函数： $σ (n) = \sum_{d | n} d$
约数个数函数： $σ_{0} (n) = \sum_{d ∣ n} d^{0}$
欧拉函数： $ϕ (n) = \sum_{a = 1}^{n} [(a, n) = 1]$ 或 $φ (n) = \sum_{a = 1}^{n} [(a, n) = 1]$
莫比乌斯函数： $μ (n) = {\begin{cases} 1 & : n = 1 \(- 1)^{k} & : n = p_{1} p_{2} \dots p_{k} \0 & : o t h e r w i s e \end{cases}$

其中 $I (n), I d (n), I d_{k} (n)$ 为完全积性函数，即不要求 $a, b$ 互质，仍有 $f (a) f (b) = f (a b)$ 。

狄利克雷卷积

设 $f : N \to R, g : N \to R$ 是两个函数，则它们的狄利克雷卷积为
$(f * g) (n) = \sum_{d | n} f (d) g (\frac{n}{d})$

定理

如果 $f (n), g (n)$ 是积性函数，则 $h (n) = (f * g) (n)$ 也是积性函数。
$f=g1\Leftrightarrow g=f\mu$

常用公式

$ε = μ * 1$
$i d = φ * 1$
$φ = μ * i d$

杜教筛

对于 $M (n) = \sum {i= 1}^n\mu ( i) $,$ ( n\leq 10^{11}) $，有$ M (n)=1-\sum{i=2}^nM (\lfloor\frac ni\rfloor) $。$ \lfloor\frac ni\rfloor $的取值最多有$ \sqrt n$ 个，可以用整除分块得出。

例题：LCMs

矩阵与高斯消元

这部分的内容基本是线性代数的内容，故简略叙述。

矩阵的定义、乘法与快速幂

矩阵的定义：

 struct Mat{
    ll a[MATN][MATN];
    Mat(){
        for(int i=0;i<MATN;i++){
            for(int j=0;j<MATN;j++){
                a[i][j]=0;
            }
        }
        for(int i=0;i<MATN;i++) a[i][i]=1;
    }
    Mat(ll a1,ll a2,ll a3,ll a4){
        a[0][0]=a1;a[0][1]=a2;
        a[1][0]=a3;a[1][1]=a4;
    }
}a[N];

矩阵的乘法：

 Mat operator *(Mat x,Mat y){
    Mat z;
    for(int i=0;i<MATN;i++){
        for(int j=0;j<MATN;j++){
            z.a[i][j]=inf;
        }
    }
    for(int i=0;i<MATN;i++){
        for(int j=0;j<MATN;j++){
            for(int k=0;k<MATN;k++){
                z.a[i][j]=min(z.a[i][j],x.a[i][k]+y.a[k][j]);
            }
        }
    }
    return z;
};

矩阵快速幂：将普通快速幂中的运算对象改为矩阵即可，不过要注意左乘还是右乘，这在矩阵中是不同的。

矩阵快速幂能够有效优化先行递推，可以用于动态 DP 或者是一些数据结构与矩阵的结合问题。

例题：Magic Gems

Sasha and Array

线性方程组

高斯消元

高斯消元是求解线性方程组的基本方法，也可以用来求解异或方程组（甚至更简单）。

 for(int i=1;i<=n;i++){
    int maxx=i;
    for(int j=i+1;j<=n;j++){
        if(fabs(a[j][i])>fabs(a[maxx][i])) maxx=j;
    }
    if(i!=maxx) swap(a[i],a[maxx]);
    if(fabs(a[i][i])<eps){
        printf("No Solution\n");
        return 0;
    }
    for(int j=i+1;j<=n+1;j++) a[i][j]/=a[i][i];
    for(int j=1;j<=n;j++){
        if(j==i) continue;
        for(int k=i+1;k<=n+1;k++){
            a[j][k]-=a[j][i]*a[i][k];
        }
    }
}

例题：EXTENDED LIGHTS OUT

行列式与逆矩阵

行列式：我们可以用高斯消元的方法将矩阵化为上三角矩阵，从而快速求出行列式。

逆矩阵：构造出增广矩阵 $W = [A | E]$ ，然后使用高斯消元的方法（行变换）将 $A$ 化成对角矩阵，再化成单位矩阵，则此时增广矩阵的右半部分即为 $A$ 的逆矩阵。

线性基

向量空间的基：向量空间中最大的线性无关组称为向量空间的一组基。

线性基：一般指 $Z_{2}^{k}$ 的一个子空间的基。应用中，常常是把若干个数看成它们的二进制分解，二进制的第 $k$ 位即为其对应向量的第 $k$ 维的值。

通过线性基，可以求一个集合 $S$ 中取一个子集异或得到所有数的数量，或是求一个集合 $S$ 中取一个子集异或可以得到的最大 / 小值。

例题：xor 序列

多项式与生成函数

生成函数

多项式 $A (x) = \sum_{i = 0}^{n} a_{i} x^{i}$
形式幂级数 $A (x) = \sum_{i \geq 0} a_{i} x^{i}$
常生成函数：一个数列 $a_{n}$ 对应的常生成函数为 $A (x) = \sum_{n \geq 0} a_{n} x^{n}$
形式幂级数 $A (x)$ 的逆元： $A (x) B (x) = 1$
逆元存在的条件： $[x^{0}] A (x) \neq 0$
暴力计算的方法：递推
Catalan 数：
- 生成函数： $C (x) = \frac{1 - \sqrt{1 - 4 x}}{2 x}$
- 通项公式： $c_{n} = \frac{1}{n + 1} (\binom{2 n}{n})$
一个数列 $a_{n}$ 对应的指数生成函数为 $A (x) = \sum_{n \geq 0} a_{n} \frac{x^{n}}{n!}$

FFT（快速傅里叶变换）

FFT 通常包含 DFT 和 IDFT。

DFT（离散傅里叶变换）：将多项式 $A (x) = a_{0} + a_{1} x + \dots + a_{n - 1} x^{n - 1}$ 转化为其点值形式 $(ω_{n}^{k}, A (ω_{n}^{k})), (k = 0, 1, \dots, n - 1) 。$

IDFT 是 DFT 的逆过程，同样可以使用 FFT 解决。

对于两个多项式相乘的问题，通过分别 DFT -> 将结果依次相乘 -> IDFT 的步骤可以在 $O (n \log n)$ 的时间复杂度内解决。

不过由于 FFT 使用的是虚数运算，故精度较低。

NTT（快速数论变换）

NTT 和 FFT 一样，用于解决两个多项式相乘的问题，不过 NTT 是在模数意义下进行的，其要求模数是满足 $p = r 2^{l} + 1$ 的质数 $p$ ，常见的模数有 $65537, 998244353, 1004535809$ 等。

拉格朗日插值定理

$n$ 个点值 $(x_{i}, y_{i}), (1 \leq i \leq n)$ , 满足 $x_{i} \neq x_{j}, (i \neq j)$ , 它们唯一确定一个 $n - 1$ 次多项式 $f (x)$ ,
$f (x) = \sum_{i = 1}^{n} y_{i} \prod_{j \neq i} \frac{x - x_{j}}{x_{i} - x_{j}}$

暴力实现：先计算 $g (x) = \prod_{i = 1}^{n} (x - x_{i})$ , 再求 $n$ 次 $g (x) / (x - x_{i})$ , 按照权值和将它们相加。 $O (n^{2})$
特殊情况 1: 只需要求一个 $f (x_{0})$ 。 $O (n^{2})$
特殊情况 2: 只需要求一个 $f (x_{0})$ , 且 $x_{i} = i$ 。 $O (n)$

阶与原根

阶：若正整数 $m, a$ , 满足 $(a, m) = 1$ , 则使得 $a^{n} \equiv 1$ (mod $m)$ 的最小的正整数 $n$ 称为 $а$ 模 $m$ 的阶，记作 $δ_{m} (a)$ 。
原根：若 $δ_{m} (a) = φ (m)$ , 则称 $a$ 为 $m$ 的一个原根。
假设 $(a, m) = 1$ , $δ = δ_{m} (a)$ , 则
- $a^{0}, a^{1}, \dots, a^{δ - 1}$ 在模 $m$ 意义下两两不同
- $a^{γ} \equiv a^{γ^{'}} (\mod m) \Leftrightarrow γ \equiv γ^{'} (\mod δ)$
- $δ ∣ φ (m)$
原根的存在定理：只有模 $2, 4, p^{a}, 2 p^{a}$ 存在原根
原根的判定定理：设 $m > 1$ , $g$ 为正整数且 $(g, m) = 1$ 。则 $g$ 是 $m$ 的原根当且仅当对于任意 $φ (m)$ 的质因子 $q_{i}, g^{\frac{φ (m)}{q_{i}}} ≢ 1 (\mod m)$
指标：对于质数 $p$ , 假设 $g$ 是 $p$ 的一个原根，则 $g^{0}, g^{1}, \dots, g^{p - 2}$ 在模 $p$ 意义下是 $1, 2, \dots, p - 1$ 的
一个排列。假设对于 $1 \leq x \leq p - 1$ 有 $g^{c} \equiv x (\mod p)$ , 则称 $x$ 的指标为 $c$ , 记作
$ind (x) = c_{0}$

代码

FFT 和 NTT 代码来自于 A999999。

FFT

 #include <algorithm>
#include <iostream>
#include <vector>
#include <cmath>
const double PI = std::acos(-1);
class Complex{
public:
    Complex(double _real = 0, double _virtual = 0);
    double getReal()const;
    double getVirtual()const;
    Complex operator +(const Complex& b)const;
    Complex operator -(const Complex& b)const;
    Complex operator *(const Complex& b)const;
private:
    double x, y;
};
Complex::Complex(double _real, double _virtual): x(_real), y(_virtual){}
double Complex::getReal()const{ return this->x;}
double Complex::getVirtual()const{ return this->y;}
Complex Complex::operator +(const Complex& b)const{
    return Complex(this->x + b.x, this->y + b.y);
}
Complex Complex::operator -(const Complex& b)const{
    return Complex(this->x - b.x, this->y - b.y);
}
Complex Complex::operator *(const Complex& b)const{
    return Complex(this->x * b.x - this->y * b.y, this->x * b.y + this->y * b.x);
}
/* recursive FFT
void FFT(std::vector<Complex>& A, int flag = 1){
    int n = A.size();
    if(n == 1) return;
    std::vector<Complex> A1(n >> 1), A2(n >> 1);
    for(int i = 0; i < (n >> 1); ++i)
        A1[i] = A[i << 1],
        A2[i] = A[i << 1 | 1];
    FFT(A1, flag); FFT(A2, flag);
    Complex w1(std::cos(2.0 * PI / n), std::sin(2.0 * PI / n) * flag), wk(1, 0);
    for(int k = 0; k < (n >> 1); wk = wk * w1, ++k)
        A[k] = A1[k] + A2[k] * wk,
        A[k + (n >> 1)] = A1[k] - A2[k] * wk;
}
*/
// make reversed binary representation array
std::vector<int> makerev(const int& len){
    std::vector<int> ans;
    ans.emplace_back(0);
    ans.emplace_back(len >> 1);
    int l = 0;
    while((1 << l) < len) l++;
    for(int i = 2; i < len; ++i)
        ans.emplace_back(ans[i >> 1] >> 1 | (i & 1) << (l - 1));
    /*
        ans[i >> 1] is the reversed representation of i >> 1
        (i >> 1) << 1 = i, so in reversed representation we need ans[i >> 1] >> 1
        if i & 1 == 1, then the MSB of reversed representation should be 1
        that is (i & 1) << (l - 1)
    */
    return ans;
}
// iterative FFT
void FFT(std::vector<Complex>& A, const int& flag = 1){
    static std::vector<int> rev;
    int n = A.size();
    if(int(rev.size()) != n)
        rev.clear(),
        rev = makerev(n);
    for(int i = 0; i < n; ++i)
    if(rev[i] > i)
        std::swap(A[i], A[rev[i]]);
    for(int len = 2, m = 1; len <= n; m = len, len <<= 1){
        Complex w1(std::cos(2.0 * PI / len), std::sin(2.0 * PI / len) * flag), wk;
        for(int l = 0, r = len; r <= n; l += len, r += len){
            wk = Complex(1, 0);
            for(int k = l; k < l + m; wk = wk * w1, ++k){
                Complex x = A[k] + A[k + m] * wk,
                        y = A[k] - A[k + m] * wk;
                A[k] = x;
                A[k + m] = y;
            }
        }
    }
}
signed main(int argc, char ** argv){
    std::cin.tie(nullptr)->sync_with_stdio(false);
    int n, m, len = 1;
    std::cin >> n >> m;
    while(len <= (n + m)) len <<= 1;// to make the length a power of 2
    std::vector<Complex> A(len), B(len);
    for(int i = 0, x; i <= n; ++i){
        std::cin >> x;
        A[i] = Complex(x, 0);
    }
    for(int i = 0, x; i <= m; ++i){
        std::cin >> x;
        B[i] = Complex(x, 0);
    }
    FFT(A); FFT(B);
    for(int i = 0; i < len; ++i)
        A[i] = A[i] * B[i];
    FFT(A, -1);
    for(int i = 0; i <= (n + m); ++i)
        std::cout << int(A[i].getReal() / len + 0.5) << ' ';
    std::cout << std::endl;
    return 0;
}

NTT

 #include <algorithm>
#include <iostream>
#include <vector>
typedef long long i64;
/// quick power
/// @brief 
/// @tparam MOD 
/// @param x 
/// @param y 
/// @return 
template<i64 MOD>
i64 qpow(i64 x, i64 y){
    i64 ans = 1;
    while(y){
        if(y & 1) ans = ans * x % MOD;
        x = x * x % MOD;
        y >>= 1;
    }
    return ans;
}
// make reversed binary representation array
/// @brief  
/// @param rev 
/// @param len 
void makerev(std::vector<int>& rev, const int& len){
    rev.resize(len);
    rev[0] = 0;
    rev[1] = len >> 1;
    int l = 0;
    while((1 << l) < len) l++;
    for(int i = 2; i < len; ++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
    /*
        rev[i >> 1] is the reversed representation of i >> 1
        (i >> 1) << 1 = i, so in reversed representation we need rev[i >> 1] >> 1
        if i & 1 == 1, then the MSB of reversed representation should be 1
        that is (i & 1) << (l - 1)
    */
}
// iterative NTT
/// @brief 
/// @tparam MOD 
/// @param A 
/// @param flag 
template<i64 MOD>
void NTT(std::vector<i64>& A, const int& flag = 1){
    static std::vector<int> rev;
    static const i64 prt = 3, invprt = qpow<MOD>(prt, MOD - 2);
    int n = A.size();
    if(int(rev.size()) != n)
        makerev(rev, n);
    for(int i = 0; i < n; ++i)
    if(rev[i] > i)
        std::swap(A[i], A[rev[i]]);
    for(int len = 2, m = 1; len <= n; m = len, len <<= 1){
        i64 w1 = qpow<MOD>(flag == 1? prt: invprt, (MOD - 1)/ len), wk;
        for(int l = 0, r = len, k; r <= n; l += len, r += len){
            for(k = l, wk = 1; k < l + m; wk = wk * w1 % MOD, ++k){
                i64 x = A[k], y = A[k + m] * wk % MOD;
                A[k] = (x + y) % MOD;
                A[k + m] = (x - y + MOD) % MOD;
            }
        }
    }
}
const i64 mod = 998244353;
signed main(int argc, char ** argv){
    std::cin.tie(nullptr)->sync_with_stdio(false);
    int n, m, len = 1;
    std::cin >> n >> m;
    while(len <= (n + m))
        len <<= 1;// to make the length a power of 2
    std::vector<i64> A(len), B(len);
    for(int i = 0; i <= n; ++i)
        std::cin >> A[i];
    for(int i = 0; i <= m; ++i)
        std::cin >> B[i];
    NTT<mod>(A); NTT<mod>(B);
    for(int i = 0; i < len; ++i)
        A[i] = A[i] * B[i] % mod;
    NTT<mod>(A, -1);
    i64 invlen = qpow<mod>(len, mod - 2);
    for(int i = 0; i <= (n + m); ++i)
        std::cout << A[i] * invlen % mod << ' ';
    std::cout << std::endl;
    return 0;
}