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

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

整数分解与筛法

欧几里得算法

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

$\gcd(a,b)=\gcd(b,a\bmod b)$。

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

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

用于求解 $ax+by=\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\pmod{b}$ 意义下的逆元等价于求一个 $x$ 满足 $ax\equiv1\pmod{b}$，由于逆元存在，$a$ 与 $b$ 互质，故 $\gcd(a,b)=1$，所以即求 $ax+by=1=\gcd(a,b)$ 的解 $x$，显然可以用扩欧解决。

例题：青蛙的约会

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

裴蜀定理扩展：对于 $a_1x_1+a_2x_2 \dots+a_nx_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)$。其中的 $phi$ 数组表示欧拉函数，将在后面的内容中提及。

例题：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$ 则较快。

【补充】数论分块

对于一些整除类问题，如果出现了 $\lfloor\frac ni\rfloor$ 型，可以考虑通过数论分块将计算范围缩减到 $\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 得到的数值都是一样的
    }
}

同余与模

逆元

若 $ab\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;
    }
}

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

推导过程：

1. $\lfloor \frac{p}{i}\rfloor\times i +p\%i \equiv 0 (\mod p)$
2. $\lfloor \frac{p}{i}\rfloor\times i\times (p\%i)^{-1} +1 \equiv 0 (\mod p)$
3. $\lfloor \frac{p}{i}\rfloor\times (p\%i)^{-1} +i^{-1} \equiv 0 (\mod p)$
4. $i^{-1} \equiv -\lfloor \frac{p}{i}\rfloor\times (p\%i)^{-1} (\mod p)$
5. $i^{-1} \equiv (p-\lfloor \frac{p}{i}\rfloor)\times (p\%i)^{-1} (\mod p)$

模板题：乘法逆元

例题：Alternating Sum

欧拉函数

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

通式：$\varphi(n)=n \prod\limits_{i=1}^{n}(1-\frac{1}{p_i})$，其中 $p_1,p_2\dots p_n$ 为 $n$​ 的所有不重复的质因数。

记 $n=p_1^{\alpha_1}p_2^{\alpha_2}\dots p_k^{\alpha_k}$，则也有 $\varphi(n)=(p_1-1)p_1^{\alpha_1-1}\times(p_2-1)p_2^{\alpha_2-1}\times \dots \times (p_k-1)p_k^{\alpha_k-1}$。

容易发现，当 $p$ 为一个质数时，有 $\varphi(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^{\varphi(p)}\equiv 1\pmod{p}$​​。

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

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

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

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

求解 $x$，满足方程组 $x ≡ a_i \pmod {m_i}$，这就是同余方程组。

利用 $x=y_1m_1+a_1=y_2m_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 ≡ a_i \pmod {m_i}$，其中 $m_i$ 两两互质。

令 $M = m_1m_2 … m_n, M_i = M/m_i$。
显然 $(M_i, m_i) = 1$，所以 $M_i$ 关于模 $m_i$ 的逆元存在。把这个逆元设为 $t_i$。
于是有：$M_it_i ≡ 1\pmod {m_i}, M_it_i ≡ 0\pmod {m_j}(j \ne i)$。
进一步：$a_iM_it_i ≡ a_i \pmod {m_i},a_iM_it_i ≡ 0 \pmod {m_j}(j \ne i)$。

因此解为 $x ≡ a_1M_1t_1 + a_2M_2t_2 + ⋯ + a_nM_nt_n \pmod M$，且在 $\bmod M$​ 意义下是唯一解。

扩展欧拉定理

$$
a^b\equiv\begin{cases}\quad a^b&b<\varphi(m)\newline a^{b\operatorname{mod}\varphi(m)+\varphi(m)}&b\geq\varphi(m)&\end{cases}(\mathrm{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\newline a^b\mod p + p&a^b\geq p&\end{cases}$ 的值。

威尔逊定理

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

排列组合与容斥原理

约定符号

• 下降幂：$n^{\underline{r}}=n\cdot(n-1)\cdot(n-2)\cdot…\cdot(n-r+1)$
• 上升幂：$n^{\bar{r}}=n\cdot(n+1)\cdot(n+2)\cdot…\cdot(n+r-1)$
• 阶乘：$n!=n\cdot(n-1)\cdot…\cdot1,\quad0!=1$

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

加法和乘法原理

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

排列与组合

排列数：$P^m_n/A^m_n=\frac{n!}{m!}$

组合数：$C^m_n/\dbinom{n}{m}=\frac{n!}{(n-r)!r!}$​

容斥原理

基础形式

设 $S$ 是一个有限集，$A_1, A_2,\dots , A_n$ 是 $S$ 的 $n$ 个子集，则
$$
|S-\bigcup_{i=1}^nA_i|=\sum_{i=0}^n(-1)^i\sum_{1\leq j_1<j_2<…<j_i\leq n}|\bigcap_{k=1}^iA_{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}…a_{i_k})$ 为 S 中同时有 $a_{i_1},a_{i_2},…,a_{i_k}$ 性质的元素的数量。
• 记 $N(a\pm b)=N(a)\pm N(b)$。​

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

最后的晚餐(dinner)

第二类斯特林数

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

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

积性函数

定义与求法

积性函数如果函数 $f:\mathbb{N}\to\mathbb{R}$ ,满足对于任意一对互质的正整数 $p,q$, 都有 $f(pq)=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);
        }
    }
}

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

例题：华华给月月出题

莫比乌斯反演

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

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

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

莫比乌斯反演 形式一

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

莫比乌斯反演 形式二

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

Sum of gcd of Tuples (Hard)

常见的积性函数

• 单位函数：$\epsilon(n)=[n=1]$

• 常数函数：$I(n)=1$

• 恒等函数：$Id(n)=n$

• 幂函数：$Id_k\left(n\right)=n^k$

• 约数函数：$\sigma(n)=\sum_{d|n}d$

• 约数个数函数： $\sigma_0\left(n\right)=\sum_{d\mid n}d^0$

• 欧拉函数：$\phi(n)=\sum_{a=1}^n[(a,n)=1]$ 或 $\varphi(n)=\sum_{a=1}^n[(a,n)=1]$

• 莫比乌斯函数：$\mu(n)=\begin{cases}1&: n=1\(-1)^k&: n=p_1p_2\cdots p_k\0&: otherwise\end{cases}$

其中 $I(n),Id(n),Id_k(n)$ 为完全积性函数，即不要求 $a,b$ 互质，仍有 $f(a)f(b)=f(ab)$。

狄利克雷卷积

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

定理

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

常用公式

• $\varepsilon=\mu*1$

• $id=\varphi*1$

• $\varphi=\mu*id$

杜教筛

对于 $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$ 的逆矩阵。

线性基

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

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

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

例题：xor序列

多项式与生成函数

生成函数

• 多项式 $A(x)=\sum_{i=0}^na_ix^i$
• 形式幂级数 $A(x)=\sum_{i\geq0}a_ix^i$

• 常生成函数：一个数列 ${a_n}$ 对应的常生成函数为 $A(x)=\sum_{n\geq0}a_nx^n$

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

FFT（快速傅里叶变换）

FFT 通常包含 DFT 和 IDFT。

DFT（离散傅里叶变换）：将多项式 $A(x)=a_{0}+a_{1}x+…+a_{n-1}x^{n-1}$ 转化为其点值形式 $(\omega_{n}^{k},A(\omega_{n}^{k})),\mathrm{~}(k=0,1,…,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}^ny_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\equiv1$ (mod $m)$ 的最小的正整数 $n$ 称为 $а$ 模 $m$ 的阶，记作 $\delta_m(a)$。
• 原根：若 $\delta_m(a)=\varphi(m)$ ,则称 $a$ 为 $m$ 的一个原根。

• 假设 $(a,m)=1$ ,$\delta=\delta_{m}(a)$ ,则

  • $a^{0},a^{1},…,a^{\delta-1}$ 在模 $m$ 意义下两两不同
  • $a^{\gamma}\equiv a^{{\gamma^{\prime}}}\mathrm{~(mod~}m)\Leftrightarrow\gamma\equiv\gamma^{\prime}\mathrm{~(mod~}\delta)$
  • $\delta\mid\varphi(m)$
• 原根的存在定理：只有模 $2,4,p^a,2p^a$ 存在原根

• 原根的判定定理：设 $m>1$ , $g$ 为正整数且 $(g,m)=1$。则 $g$ 是 $m$ 的原根当且仅当对于任意 $\varphi(m)$ 的质因子 $q_i,g^{\frac{\varphi(m)}{q_i}}\not\equiv1\mathrm{~(mod~}m)$
• 指标：对于质数 $p$,假设 $g$ 是 $p$ 的一个原根，则 $g^0,g^1,…,g^{p-2}$ 在模 $p$ 意义下是 $1,2,…,p-1$ 的
  一个排列。假设对于 $1\leq x\leq p-1$ 有 $g^c\equiv x(\operatorname{mod}p)$,则称 $x$ 的指标为 $c$,记作
  $\operatorname{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;
}