数学讲课

数学

数论

初

埃氏筛，线性筛，数论函数，卷积

反演

例题

UVA11426

P5518 莫比乌斯反演基础练习题

筛法--杜教筛

试求：(displaystyle S(n)=sum_{i=1}^nf(i))

因为

(displaystyle{egin{aligned}&sum_{i=1}^n(f*g)(i)\=&sum_{i=1}^nsum_{d|i}f(d)g(dfrac id)\=&sum_{d=1}^ng(d)sum_{i=1}^{lfloorfrac nd floor}f(i)\=&sum_{d=1}^ng(d)S(lfloordfrac nd floor)end{aligned}})

故将最后拆开再移项，

(displaystyle S(n)=g(1)S(n)=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(d)S(lfloordfrac nd floor))

若可以快速得知 g 的每一项与 f*g 的前缀和，则可通过递归快速求解

注意long long和预处理

(O(n^{frac 34}))

筛法--Min_25

试求同上

要求：

1. 可快速求(f(prime^k))
2. (exist f'(x),forall pin primef(p)=f'(p))，且该函数完全积性，可快速求前缀和

令 (p_k​) 为第k大质数，tot 为n内质数个数，m(x) 为x最小质因子，m(1)=0

预处理：质数

令 (G(n,k)=displaystylesum_{i=1}^nf'(i)[iin prime|m(i)>p_k]​)

则有转移：(G(n,k)=egin{cases}G(n,k-1)&p_k^2>n\G(n,k-1)-f'(p_k)(G(lfloordfrac n{p_k} floor,k-1)-displaystylesum_{i=1}^{k-1}f'(p_i))&p_k^2le nend{cases})

若令L为满足 (p_m^2>n​) 的最小整数，则 (G(n,L)=displaystylesum_{i=1}^{tot}f(p_i)​)

求解

令 (S(n,k)=displaystylesum_{i=2}^nf(i)[m(i)ge p_k]​)

有转移：(displaystyle S(n,k)=S(n,k+1)+f(p_k)+sum_{i=1}^{p_k^{i+1}le n}(f(p_k^{i+1})+f(p_k^i)S(lfloordfrac n{p_k^i} floor,k+1))​)

等价于：(displaystyle S(n,k)=sum_{i=k}^{tot}f(p_i)+sum_{i=k}^{tot}sum_{j=1}^{p_i^jle n}(f(p_i^{j+1})+f(p_i^j)S(lfloordfrac n{p_i^j} floor,i+1)))

答案即 S(n,1)+1

另外两点：

1. 由于G与S第一维太大，而形式均为 (lfloordfrac nx floor)，故用 x 或 (lfloordfrac nx floor) 将下标离散化
2. 注意求S时的返回条件，能剪一点枝

筛法--orz

Legendre , Meissel-Lehmer , Deléglise-Rivat

如果有人能教教我？

阶与原根

阶

由欧拉定理知，对 (ain mathbb{Z})，(pinmathbb{N}^{*})，若 (gcd(a,p)=1)，则 (a^{varphi(p)}equiv 1pmod p)

故存在满足同余式 (a^n equiv 1 pmod p) 的最小正整数 n，称作 a 模 p 的阶，记作 (delta_p(a))

性质

1. 设 (pinmathbb{N}^{*})，(a,binmathbb{Z})，gcd(a,p)=gcd(b,p)=1，则 (delta_p(ab)=delta_p(a)delta_p(b)) 的充要条件是 (gcd(delta_p(a),delta_p(b))=1)

证明：

• 必要性

由 (a^{delta_p(a)}equiv b^{delta_p(b)}equiv1pmod p) 可知 ((ab)^{operatorname{lcm}(delta_p(a),delta_p(b))}equiv1pmod p)

故 (delta_p(ab)|operatorname{lcm}(delta_p(a),delta_p(b)))

而 (delta_p(a)delta_p(b)=delta_p(ab))

即 (gcd(delta_p(a),delta_p(b))=dfrac{delta_p(a)delta_p(b)}{operatorname{lcm}(delta_p(a),delta_p(b))}=1)。

• 充分性

由 ((ab)^{delta_p(ab)}equiv1pmod p)

可知 (1equiv(ab)^{delta_p(ab)delta_p(b)}equiv a^{delta_p(ab)delta_p(b)}pmod p)

故 (delta_p(a)|delta_p(ab)delta_p(b))

而 (gcd(delta_p(a),delta_p(b))=1)，故 (delta_p(a)|delta_p(ab))

同理 (delta_p(b)|delta_p(ab))

所以 (delta_p(a)delta_p(b)|delta_p(ab))

另一方面，有

((ab)^{delta_p(a)delta_p(b)}equiv(a^{delta_p(a)})^{delta_p(b)} imes(b^{delta_p(b)})^{delta_p(a)}equiv 1 pmod p)

故 (delta_p(ab)|delta_p(a)delta_p(b))

综上，(delta_p(ab)=delta_p(a)delta_p(b))

2. 设 (kinmathbb N,pinmathbb N^*,ainmathbb Z)，gcd(a,p)=1，则 (delta_p(a^k)=dfrac{delta_p(a)}{gcd(delta_p(a),k)})

证明：

因

为 (a^{kdelta_p(a^k)}=(a^k)^{delta_p(a^k)}equiv1pmod p)

故 (delta_p(a)|kdelta_p(a^k)Rightarrowdfrac{delta_p(a)}{gcd(delta_p(a),k)}|delta_p(a^k))

同时 ((a^k)^{frac{delta_p(a)}{gcd(delta_p(a),k)}}=(a^{delta_p(a)})^{frac{k}{gcd(delta_p(a),k)}}equiv 1 pmod p)

故 (delta_p(a^k)|dfrac{delta_p(a)}{gcd(delta_p(a),k)})

综上得证

原根

原根：设 (p in mathbb{N}^*)，(ain mathbb{Z})，若 gcd(a,p)=1，且 (delta_p(a)=varphi(p))，则称 a 为模 p 的原根

原根判定定理

若对于 (varphi(p)) 任何大于1且不为自身的因数 x，都有 (g^{varphi(p)/x} otequiv 1pmod p)，则 g 是模 p 的原根

证明：假设存在一个 (t<varphi(p)) 使得 (g^tequiv 1pmod{p}) 且 (forall iinleft[1,m ight]:a^{frac{varphi(p)}{d_{i}}} otequiv 1pmod{p})

由裴蜀定理得，一定存在一组 k,x 满足 (kt=xvarphi(p)+gcd(t,varphi(p)))

故有(g^{gcd(t,varphi(p))}equiv g^{xvarphi(p)+gcd(t,varphi(p))}equiv g^{kt}equiv1pmod p)，得(gcd(t,varphi(p))|varphi(p))

而 (gcd(t,varphi(p))leqslant t<varphi(p))

故 (gcd(t,varphi(p))) 必整除 ({frac{varphi(p)}{d_{1}}},{frac{varphi(p)}{d_{2}}},ldots,{frac{varphi(p)}{d_{m}}}) 中至少一个，设 (gcd(t,varphi(p))mid {frac{varphi(p)}{d_{i}}})，则 (g^{frac{varphi(p)}{d_{i}}}equiv g^{gcd(t,varphi(p))}equiv 1pmod{p})

故假设不成立，原命题成立

原根个数

若一个数 p 有原根，则它原根的个数为 (varphi(varphi(p)))

证明：若 p 有原根 g，则 (delta_p(g^k)=dfrac{delta_p(g)}{gcd(delta_p(g),k)}=dfrac{varphi(p)}{gcd(varphi(p),k)})

所以若 (gcd(k,varphi(p))=1)，则有：(delta_p(g^k)=varphi(p))，即 (g^k) 也是模 p 的原根。

而满足 (gcdig(varphi(p),kig)=1) 且 (1le klevarphi(p)) 的 k 有 (varphi(varphi(p))) 个。所以原根就有 (varphi(varphi(p))) 个

原根存在定理

一个数 m 存在原根当且仅当 (m=2,4,p^{alpha},2p^{alpha})，其中 p 为奇素数，(alphain mathbb N^*)

证明：

1. m=2,4：显然正确
2. (m=p^alpha)

先证明一个引理：存在模 p 的原根 g，使得 (g^{p-1} otequiv1pmod {p^2})

任取 p 的一个原根 x，若 x 符合条件则 g=x

否则 g=x+p

显然 x+p 也是 p 的原根，同时符合条件：

(egin{aligned}(x+p)^{p-1}&equiv C_{p-1}^0x^{p-1}+C_{p-1}^1px^{p-2}\&equiv x^{p-1}+p(p-1)x^{p-2}\&equiv 1-px^{p-2}\& otequiv 1 pmod {p^2}end{aligned})

接下来证明对任意 (alphainmathbb{N}^{*})，g 是模 (p^{alpha}) 的原根

首先证明，对任意 (etainmathbb{N}^{*})，都 (exist p mid k_eta ,g^{varphi(p^eta)}=1+p^{eta} imes k_{eta})

(eta=1) 时，由 g​ 的选取可知结论成立

设上式对 (eta) 时成立，则 (egin{aligned}g^{varphi(p^{eta+1})}&=(g^{varphi(p^{eta})})^{p}\&=(1+p^{eta} imes k_{eta})^p\&= 1+p^{eta+1} imes(k_{eta}+p imes q) end{aligned})

其中 q 为余项（至少有因数 (p^{2eta})）带来的贡献，故命题对 (eta+1) 成立

由数学归纳法，对任意 (etainmathbf{N}^{*}) 命题成立

记 (delta=delta_{p^alpha}(g))，则由欧拉定理，可知 (delta|varphi(p^alpha)= p^{alpha-1}(p-1))

而 g 为模 p 的原根，所以可设 (delta=p^{eta-1}(p-1) (1leq etaleq alpha))

由之前的结论：(g^{varphi(p^{eta})} otequiv 1pmod {p^{eta+1}}Rightarrow g^{delta} otequiv 1pmod {p^{eta+1}})

结合 (g^{delta}equiv 1pmod {p^alpha}) 知 (eta+1>alpha)。

综上，(eta=alpha)，即：(delta_{p^alpha}(g)=p^{alpha-1}(p-1)=varphi(p^alpha))

从而，g 是模 (p^{alpha}) 的原根

3. (m=2p^alpha)

证明：设 g 是模 (p^{alpha}) 的原根，则 (g+p^{alpha}) 也是模 (p^{alpha}) 的原根

在 g 和 (g+p^{alpha}) 中有一个是奇数，设这个奇数是 G，则 (gcd(G,2p^{alpha})=1)

由 (G^{delta_{2p^{alpha}}(G)}equiv 1pmod {2p^{alpha}}) 得 (G^{delta_{2p^alpha}(G)}equiv1pmod{p^{alpha}})

利用 G 为模 (p^{alpha}) 的原根可知 (varphi(p^{alpha})middelta_{2p^{alpha}}(G))

而由欧拉定理，(delta_{2p^{alpha}}(G)midvarphi(2p^{alpha}))

结合 (varphi(p^{alpha})=varphi(2p^{alpha})) 可知 G 为模 (2p^{alpha}) 的原根

4. (m e 2,4,p^alpha,2p^alpha)

对于 (m=2^{alpha})，(alphainmathbb{N}^{*},alphageq 3)，则对任意奇数 a=2k+1 均有

(egin{aligned}a^{2^{alpha-2}}&=(2k+1)^{2^{alpha-2}}\&equiv 1+C_{2^{alpha-2}}^1(2k)+C_{2^{alpha-2}}^{2}(2k)^{2}\&equiv1+2^{alpha-1}k+2^{alpha-1}(2^{alpha-2}-1)k^2\&equiv 1+2^{alpha-1}(k+(2^{alpha-2}-1)k)\&equiv 1 pmod mend{aligned})

其中最后一步用到 (k) 与 ((2^{alpha-2}-1)k) 同奇偶，故其和为偶数

若 (m) 不是 (2) 的幂，且 (m) 为符合题目条件的数，则可设 (m=rt)，这里 (2<r<t) 且 (gcd(r,t)=1)

此时，若 gcd(a,m)=1，由欧拉定理可知 (a^{varphi(r)}equiv1pmod r,a^{varphi(t)}equiv1pmod t)

而 n>2 时，(varphi(n)) 为偶数，故 (a^{frac12varphi(r)varphi(t)}equiv1pmod{rt})

进而 (delta_m(a)leqdfrac12varphi(r)varphi(t)=dfrac12varphi(rt)=dfrac12varphi(m)<varphi(m))

由原根定义可得模 (m) 的原根不存在

二次剩余

一个数 a，如果不是 p 的倍数且模 p 同余于某个数的平方，则称 a 为模 p 的二次剩余。而一个不是 p 的倍数的数 b，不同余于任何数的平方，则称 b 为模 p 的二次非剩余。

对二次剩余求解，也就是对常数 a 解下面的这个方程：(x^2 equiv a pmod p)

通俗一些，可以认为是求模意义下的开方

这里只讨论 p 为奇素数 的求解方法，将会使用 Cipolla 算法。

解的数量

对于 (x^2 equiv n pmod p)，能满足“n 是模 p 的二次剩余的”n 一共有 (frac{p-1}2) 个（0 不包括在内），二次非剩余有 (frac{p-1}{2}) 个。

证明

x=0 对应了 n=0 的特殊情况，因此我们只用考虑 (x in [1,frac{p-1}{2}]) 的情况

显然有 ((p-x)^2 equiv x^2 pmod p)，那么当 (x in [1,frac{p-1}{2}]) 我们可以取到所有解

接下来我们只需要证明当 (xin[1,frac{p-1}{2}]) 时 (x^2 mod p) 两两不同

运用反证法，假设存在不同的两个整数 (x,y in [1,frac{p-1}{2}]) 且 (x^2 equiv y^2 pmod p)

则有 ((x+y)(x-y) equiv 0 pmod p)

显然 (-p<x+y<p,-p<x-y<p,x+y eq 0,x-y eq 0)，故假设不成立，原命题成立

勒让德记号

(left(frac{n}{p} ight)=egin{cases}1,\,&p mid n ext{且}n ext{是}p ext{的二次剩余}\-1,\,&p mid n ext{且}n ext{不是}p ext{的二次剩余}\0,\,&pmid nend{cases})

通过勒让德符号可以判断一个数 n 是否为二次剩余，具体判断 n 是否为 p 的二次剩余，需要通过欧拉判别准则来实现

欧拉判别准则

(left(frac{n}{p} ight)equiv n^{frac{p-1}{2}}pmod p)

若 n 是二次剩余，当且仅当 (n^{frac{p-1}{2}}equiv 1pmod p​)

若 n 是二次非剩余，当且仅当 (n^{frac{p-1}{2}}equiv -1pmod p​)

证明

由于 p 为奇素数，那么 p-1 是一个偶数，根据费马小定理 (n^{p - 1} equiv 1 pmod{p}​)，那么就有

((n^{frac{p-1}{2}}-1)cdot(n^{frac{p-1}{2}}+1)equiv 0 pmod p)

其中 p 是一个素数，所以 (n^{frac{p-1}{2}}-1​) 和 (n^{frac{p-1}{2}}+1​) 中必有一个是 p 的倍数，

因此 (n^{frac{p-1}{2}}mod p) 必为 1 或 -1

p 是一个奇素数，所以关于 p 的原根存在

设 a 是 p 的一个原根，则存在 (1 leqslant j leqslant p-1) 使得 (n=a^j)

于是就有(a^{jfrac{p-1}{2}}equiv1pmod p)

a 是 p 的一个 原根，因此 a 模 p 的指数是 p-1，于是 p-1 整除 (frac{j(p-1)}{2})，这说明 j 是一个偶数

令 (i=dfrac{j}{2})，就有 ((a^i)^2=a^{2i}=n)，即 n 是模 p 的二次剩余

Cipolla 算法

找到一个数 (a) 满足 (a^2-n) 是 非二次剩余，至于为什么要找满足非二次剩余的数，在下文会给出解释。
这里通过生成随机数再检验的方法来实现，由于非二次剩余的数量为 (frac{p-1}{2})，接近 (frac{p}{2})，所以期望约 2 次就可以找到这个数。

建立一个＂复数域＂，并不是实际意义上的复数域，而是根据复数域的概念建立的一个类似的域。
在复数中 (i^2=-1)，这里定义 (i^2=a^2-n)，于是就可以将所有的数表达为 (A+Bi) 的形式，这里的 (A) 和 (B) 都是模意义下的数，类似复数中的实部和虚部。

在有了 (i) 和 (a) 后可以直接得到答案，(x^2equiv npmod p) 的解为 ((a+i)^{frac{p+1}{2}})。

证明

• 定理 1：((a+b)^pequiv a^p+b^ppmod p)

[egin{aligned} (a+b)^p &equiv sum_{i=0}^{p}mathrm C_p^i a^{p-i}b^i \ &equiv sum_{i=0}^{p}frac{p!}{(p-i)!i!}a^{p-i}b^i \ &equiv a^p+b^ppmod p end{aligned} ]

可以发现只有当 (i=0) 和 (i=p) 时由于没有因子 (p) 不会因为模 (p) 被消去，其他的项都因为有 (p) 因子被消去了。

• 定理 2：(i^pequiv -ipmod p)

[egin{aligned} i^p &equiv i^{p-1} cdot i \ &equiv (i^2)^{frac{p-1}{2}}cdot i \ &equiv (a^2-n)^{frac{p-1}{2}}cdot i \ &equiv -i pmod p end{aligned} ]

• 定理 3：(a^pequiv a pmod p) 这是 费马小定理 的另一种表达形式

有了这三条定理之后可以开始推导

[egin{aligned} x &equiv (a+i)^{frac{p+1}{2}} \ &equiv ((a+i)^{p+1})^{frac{1}{2}} \ &equiv ((a+i)^pcdot (a+i))^{frac{1}{2}} \ &equiv ((a^p+i^p)cdot(a+i))^{frac{1}{2}} \ &equiv ((a-i)cdot(a+i))^{frac{1}{2}} \ &equiv (a^2-i^2)^{frac{1}{2}} \ &equiv (a^2-(a^2-n))^{frac{1}{2}} \ &equiv n^{frac{1}{2}}pmod p end{aligned} ]

( herefore xequiv (a+i)^{frac{p+1}{2}} equiv n^{frac{1}{2}}pmod p)

生成函数

普通生成函数

序列 a 的普通生成函数（ordinary generating function，OGF）定义为形式幂级数：(F(x)=sum_{n}a_n x^n)

换句话说，如果序列 a 有通项公式，那么它的普通生成函数的系数就是通项公式。

基本运算

考虑两个序列 a,b 的普通生成函数，分别为 F(x),G(x)，则有 (F(x)pm G(x)=sum_n (a_npm b_n)x^n)

因此 (F(x)pm G(x)) 是序列 (langle a_npm b_n angle) 的普通生成函数。

考虑乘法运算，也就是卷积：(displaystyle F(x)G(x)=sum_n x^n sum_{i=0}^na_ib_{n-i}​)

因此 F(x)、G(x) 是序列 (langledisplaystylesum_{i=0}^n a_ib_{n-i} angle) 的普通生成函数。

封闭形式

在运用生成函数的过程中，我们不会一直使用形式幂级数的形式，而会适时地转化为封闭形式以更好地化简。

例如 (langle 1,1,1,cdots angle) 的普通生成函数 (F(x)=sum_{nge 0}x^n)，发现 F(x)x+1=F(x)

那么解这个方程得到 (F(x)=dfrac{1}{1-x})

这就是 (sum_{nge 0}x^n) 的封闭形式

考虑等比数列 (langle 1,p,p^2,p^3,p^4,cdots angle) 的生成函数 (F(x)=sum_{nge 0}p^nx^n)，有

(egin{aligned}F(x)px+1 &=F(x)\F(x) &=frac{1}{1-px}end{aligned})

等比数列的封闭形式与展开形式是常用的变换手段

牛顿二项式定理

我们重新定义组合数的运算：(inom{r}{k}=frac{r^{underline{k}}}{k!}quad(rinmathbf{C},kinmathbf{N}))

注意 r 的范围是复数域
在这种情况下，对于 (alphainmathbf{C})，有 ((1+x)^{alpha}=sum_{nge 0}inom{alpha}{n}x^n)

二项式定理其实是牛顿二项式定理的特殊情况

指数型生成函数

序列 a 的指数生成函数（exponential generating function，EGF）定义为形式幂级数：(hat{F}(x)=sum_{n}a_n frac{x^n}{n!})

基本运算

指数生成函数的加减法与普通生成函数是相同的，也就是对应项系数相加。

对于两个序列 a,b，设它们的指数生成函数分别为 (hat{F}(x),hat{G}(x))，则

(egin{aligned}hat{F}(x)hat{G}(x)&=sum_{ige 0}a_ifrac{x^i}{i!}sum_{jge 0}b_jfrac{x^j}{j!}\&=sum_{nge 0}x^{n}sum_{i=0}^na_ib_{n-i}frac{1}{i!(n-i)!}\&=sum_{nge 0}frac{x^{n} }{n!}sum_{i=0}^ninom{n}{i}a_ib_{n-i}end{aligned})

因此 (hat{F}(x)hat{G}(x)) 是序列 (langledisplaystylesum_{i=0}^n inom{n}{i}a_ib_{n-i} angle) 的指数生成函数

封闭形式

同样考虑指数生成函数的封闭形式

序列 (langle 1,1,1,cdots angle) 的指数生成函数是 (sum_{nge 0}frac{x^n}{n!}=e^x)。因为你将 (e^x) 在原点处泰勒展开就得到了它的无穷级数形式

类似地，等比数列 (langle 1,p,p^2,cdots angle​) 的指数生成函数是 (sum_{nge 0}frac{p^nx^n}{n!}=e^{px}​)

指数生成函数与普通生成函数

如何理解指数生成函数？我们定义序列 (a​) 的指数生成函数是 (F(x)=sum_{nge 0}a_nfrac{x^n}{n!}​)，但 (F(x)​) 实际上也是序列 (langle frac{a_n}{n!} angle​) 的普通生成函数

也就是说，不同的生成函数只是对问题理解方式的转变

排列与圆排列

长度为 n 的排列数的指数生成函数是 (hat{P}(x)=sum_{nge 0}frac{n!x^n}{n!}=sum_{nge 0}x^n=frac{1}{1-x})

圆排列的定义是把 (1,2,cdots,n) 排成一个环的方案数，也就是说旋转后的方案的等价的（但翻转是不等价的）

n 个数的圆排列数显然是 (n-1)!，因此 n 个数的圆排列数的指数生成函数是

(hat{Q}(x)=sum_{nge 1}frac{(n-1)!x^n}{n!}=sum_{nge 1}frac{x^n}{n}=-ln(1-x)=lnleft( frac{1}{1-x} ight))，也就是说 (exp hat{Q}(x)=hat{P}(x))

但这只是数学层面的推导。我们该怎样直观理解：圆排列数的 EGF 的 exp 是排列数的 EGF？

一个排列，是由若干个置换环构成的。例如 p=[4,3,2,5,1]​ 有两个置换环：

（也就是说我们从 (p_i) 向 i​ 连有向边）

而不同的置换环，会导出不同的排列

如我将第二个置换环改成

那么它对应的排列就是 [5,3,2,1,4]​

即长度为 n​ 的排列的方案数是

1. 把 (1,2,cdots,n) 分成若干个集合
2. 每个集合形成一个置换环

的方案数

而一个集合的数形成置换环的方案数显然为该集合大小的圆排列方案数。因此长度为 n 的排列的方案数就是：把 (1,2cdots,n​) 分成若干个集合，每个集合的圆排列方案数之积

这就是多项式 exp 的直观理解。

推广之

• 如果 n 个点带标号生成树的 EGF 是 (hat{F}(x))，那么 n 个点带标号生成森林的 EGF 就是 (exp hat{F}(x))

  直观理解为，将 n 个点分成若干个集合，每个集合构成一个生成树的方案数之积

• 如果 n 个点带标号无向连通图的 EGF 是 (hat{F}(x))，那么 n 个点带标号无向图的 EGF 就是 (exp hat{F}(x))

  后者计算易得 (exp hat{F}(x)=sum_{nge 0}2^{inom{n}{2}}frac{x^n}{n!})，因此要计算前者再多项式 ln 即可

记 (mathcal{P}) 表示素数集合。

狄利克雷生成函数

对于无穷序列 (f_1, f_2, ldots)，定义其狄利克雷生成函数（Dirichlet series generating function，DGF）为：( ilde{F}(x) = sum_{ige 1}dfrac{f_i}{i^x})

如果序列 f 满足积性（积性函数）：(forall iperp j, ; f_{ij} = f_i f_j​)，那么其 DGF 可以由质数幂处的取值表示：( ilde{F}(x) = prod_{pin mathcal{P}} left(1 + dfrac{f_p}{p^x} + dfrac{f_{p^2}}{p^{2x}} + dfrac{f_{p^3}}{p^{3x}} + cdots ight)​)

对于两个序列 f,g​，其 DGF 之积对应的是两者的狄利克雷卷积序列的 DGF：

( ilde{F}(x) ilde{G}(x) = sum_{i} sum_{j}dfrac{f_i g_j}{(ij)^x} = sum_{i} dfrac{1}{i^x}sum_{d | i} f_d g_{frac{i}{d}}​)

常见积性函数的 DGF

DGF 最适合用于研究与积性函数的狄利克雷卷积相关的问题，因此首先我们要了解常见积性函数的 DGF

黎曼函数

序列 (langle1, 1, 1, ldots angle) 的 DGF 是 (sum_{ige 1}frac{1}{i^x} = zeta(x))，(zeta) 是黎曼函数。

由于其满足积性，因此我们可以得到 (langle1, 1, 1, ldots angle) 的 DGF 的另一种形式：

(zeta(x) = prod_{pinmathcal{P}} left(1 + dfrac{1}{p^x} + dfrac{1}{p^{2x}} + ldots ight) = prod_{pin mathcal{P}} dfrac{1}{1-p^{-x}})

莫比乌斯函数

对于莫比乌斯函数 (mu)，它的 DGF 定义为 ( ilde{M} (x) = prod_{pin mathcal{P}}left(1 - dfrac{1}{p^x} ight) = dfrac1{prod_{pin mathcal{P}}sum_{ege0}p^{-es}}= dfrac{1}{zeta(x)})

欧拉函数

对于欧拉函数 (varphi)，它的 DGF 定义为 ( ilde{Phi}(x) = prod_{pinmathcal{P}} left(1 + dfrac{p-1}{p^x} + dfrac{p(p-1)}{p^{2x}} + dfrac{p^2(p-1)}{p^{3x}} + ldots ight) = prod_{pin mathcal{P}}dfrac{1-p^{-x}}{1-p^{1-x}}= dfrac{zeta(x-1)}{zeta(x)})

幂函数

对于函数 (I_k (n) = n^k)，它的 DGF 定义为 ( ilde{I_k} (x) = prod_{pinmathcal{P}} left(1 + dfrac{p^k}{p^x} + dfrac{p^{2k}}{p^{2x}} + ldots ight) = prod_{pin mathcal{P}} dfrac{1}{1-p^{k-x}} = zeta(x-k))

根据这些定义，容易推导出 (varphi ast 1 = I)，(ast) 表示狄利克雷卷积，因为 ( ilde{Phi}(x)zeta(x) = zeta(x-1))。

其他函数

对于约数幂函数 (sigma_k(n) = sum_{d|n}d^k​)，它的 DGF 可以表示为狄利克雷卷积的形式：( ilde S(x) = zeta(x-k)zeta(x)​)。

对于 (u(n) = |mu(n)|)（无平方因子数），它的 DGF 为 ( ilde{U}(x) = prod_{pin mathcal{P}} (1+p^{-x}) = frac{zeta(x)}{zeta(2x)})。

相关应用

DGF 的应用主要体现在构造积性序列的狄利克雷卷积序列。研究方向通常是质数处的取值。

例如在杜教筛的过程中，要计算积性序列（积性函数在正整数处的取值构成的序列）(f) 的前缀和，我们需要找到一个积性序列 (g) 使得 (fast g) 和 (g) 都可以快速求前缀和。那么我们可以利用 DGF 推导这一过程。

以洛谷 3768 简单的数学题为例，我们要对 (f_i = i^2varphi(i)) 构造一个满足上述条件的积性序列 (g)。由于 (f) 是积性的，考虑其 DGF

[ ilde{F}(x) = prod_{p in mathcal{P}} left(1 + sum_{kge 1} frac{p^{3k-1}(p-1)}{p^{kx}} ight) = prod_{pin mathcal{P}} frac{1-p^{2-x}}{1-p^{3-x}} = frac{zeta(x-3)}{zeta(x-2)} ]

因此 ( ilde{F}(x)zeta(x-2) = zeta(x-3)​)。而 (zeta(x-2)​) 对应的积性函数为 (I_2​)，所以令 (g = I_2​) 即可。这样有 (fast g = I_3​)，两者都是可以快速计算前缀和的。

组合数学

群论

群

给定一个集合(G={a,b,c,dots})和集合(G)上的二元运算，并满足:

1. 封闭性：(forall a,bin G, exists cin G, a*b=c)

2. 结合律：(forall a,b,cin G, (a*b)*c=a*(b*c))

3. 单位元：(exists ein G,forall ain G,a*e=e*a=a)

4. 逆元：(forall ain G,exists bin G, a*b=b*a=e,)记(b=a^{-1})

则称集合G在运算(*)之下是一个群，简称G是群。

置换，置换群，循环

(Z_k) ((K)不动置换类)

设(G)是(1dots n)的置换群,若(K)是(1dots n)中某个元素,(G)中使(K)保持不变的置换的全体,记以(Z_k),叫做(G)中使(K)保持不动的置换类

(E_k)(等价类)

设(G)是(1dots n)的置换群,(K)是(1dots n)中某个元素,(K)在(G)作用下的轨迹，记作(E_k),即(K)在(G)的作用下所能变化成的所有元素的集合

轨道-稳定核定理:

(|E_k|*|Z_k|=|G|)

(Burnside)引理

互异的组合状态的个数(displaystyle L=dfrac{1}{|G|}sum_{i=1}^{|G|}D(a_i))((D(a_i)) 表示在置换(a_i)下不变的元素的个数)

(Pacute{o}lya)定理

设G是p个对象的一个置换群，用m种颜色涂染p个对象，gi的循环节数为c(gi)，则不同染色方案为(displaystyle L=dfrac{1}{|G|}sum_{i=1}^{|G|}m^{c_i})

多项式

FFT

优化--系数 ( ightarrow) 点值

前置知识

单位根：

在复平面上画一个以(0,0)为圆心的单位圆，将该圆n等分

规定(1,0)为第零个等分点，逆时针标号，就得到了第0到n-1个等分点，也就是第0到n-1个n次单位根，记做(omega_n^0,omega_n^1,cdots,omega_n^{n-1})

几个性质与证明

1. (omega_n^k=e^frac{2pi ik}{n})

   根据欧拉公式 (e^{i heta}=cos heta+i*sin heta) 可得

2. (omega_{dn}^{dk}=omega_n^k)

   (omega_{dn}^{dk}=e^frac{2pi idk}{dn}=e^frac{2pi ik}{n}=omega_n^k)

3. 设 (omega_n^k=a+bi) ，则 (omega_n^{-k}=a-bi)

   将 (omega_n^k=e^frac{2pi ik}{n}) 带入即可证明

4. (omega_n^{k+frac n2}=-omega_n^k)

   (omega_n^{k+frac n2}=omega_n^k*omega_n^frac n2=-omega_n^k)

DFT

DFT主要的思想是在求一部分值的时候求出来另外一部分，分治进行

比如我们现在要把多项式 (A(x)=displaystylesum_{i=0}^na_ix^i) 在
(x_1,x_2,cdots,x_n) 处的取值求出

将 A 的系数按照下标的奇偶性分类，即

(displaystyle A_0(x)=a_0+a_2x^1+a_4x^2+cdots, A_1(x)=a_1+a_3x^1+a_5x^2+cdots)

所以 (A(x)=A_0(x^2)+x*A_1(x^2))

分别将 (omega_n^k) 与 (omega_n^{k+frac n2}) 带入得：
(egin{cases}A(omega_n^k)&=&A_0(omega_n^{2k})+omega_n^kA_1(omega_n^k)\A(omega_n^{ k+frac n2})&=&A_0(omega_n^{2k})-omega_n^{2k}A_1(omega_n^{2k})end{cases})

发现只有 (A_1) 的常数项不同，所以我们可以仅递归计算每一个 k 的
(A_0(omega_n^k),A_1(omega_n^k)) ，并用这两项计算出 (A(omega_n^k)) 与 (A(omega_n^{k+frac n2}))

这样每次把问题规模减半，时间复杂度是 (T(n)=2T(dfrac n2)+O(n)=O(nlog n))

优化--点值 ( ightarrow) 系数

前置知识

1. 前面所有

2. 矩阵

IDFT

用 DFT 将系数表示为点值后，还需用同样优秀的复杂度变回去，即 IDFT 了

先证明一个单位根的性质： (displaystyledfrac1nsum_{i=0}^{n- 1}omega_n^{v*i}=[v\%n==0])

由于 (omega_n^{v*i}=omega_n^{v*(i-1)}*omega_n^v)

所以 ({omega_n^{v*i}}) 是个等比数列

1. (v\%n=0)

   (displaystylesum_{i=0}^{n-1}omega_n^{v*i}=sum_{i=0}^{n-1} 1^i=n)

2. (v\%n ot=0)

   根据等比数列求和公式,(displaystylesum_{i=0}^{n-1}omega_n^{v*i}=dfrac{1-omega_n^{n*v}}{1-omega_n^v}=dfrac{1-1^v}{1-omega_n^v}=0)

假设我们要求 c=a*b​：

(displaystyle{egin{aligned}c_i=&sum_{j=0}^{i-1}a_j imes b_{i-1-j}\c_i=&sum_{p=0}sum_{q=0}a_p imes b_q[(p+q-i)\%n=0]\nc_i=&sum_{p=0}sum_{q=0}a_p imes b_qsum_{j=0}omega_n^{(p+q-i)j}\nc_i=&sum_{j=0}omega_n^{(-i)j}(sum_{p=0}omega_n^{pj}a_p)(sum_{q=0}omega_n^{qj}b_q)end{aligned}})

设 (f_a(j)=sum_{i=0}omega_n^{ji}a_i , f_a^{-1}(j)=sum_{i=0}omega_n^{(-j)i}a_i)

则

(displaystyle{egin{aligned}nc_i=&sum_{j=0}omega_n^{(-i)j}f_a(j)f_b(j)\nc_i=&sum_{j=0}omega_n^{(-i)j}f_c(j)\nc_i=&f_{f_c}^{-1}(i)end{aligned}})

而我们发现， (f_a(j)) 就是多项式 a 在 (omega_n^j) 处的点值，也就是说
(f_a) 就是 a 在 DFT 后的结果

所以 (f_a^{-1}) 就是我们的 IDFT

所以我们发现 IDFT 的过程其实和 DFT 的过程极为相似，但有两点不同

1. 带入的值从 (omega_n^j) 变为 (omega_n^{-j})

2. 在最后要将所有的 (f_{f_c}^{-1}(i)) 除以 n 才是真正的 (c_i)

小优化

发现在递归时，会浪费很多时间

而我们递归的原因是为了将系数奇偶分类，偶在前，奇在后

所以如果我们提前模拟将系数全部分好类，再用 for 循环代替递归，就可以节省很多时间

有一个性质：最后的序列其实就是原序列二进制反转了一下：

考虑我们在第 i 次分组时，参照的是二进制第 i 位的奇偶，确定的是位置的二进制第 (len-i+1) 位(也就是 0 向左，1 向右)，刚好奇偶性相同

所以就有了上面的性质

现在我们可以 O(n) 的完成奇偶分类，通过非递归快速实现FFT

const double PI=acos(-1.0);
struct Complex{//复数结构体
	double x,y;
	Complex(double _x=0.0,double _y=0.0){x=_x;y=_y;}
	Complex operator - (const Complex &b)const{return Complex(x-b.x,y-b.y);}
	Complex operator + (const Complex &b)const{return Complex(x+b.x,y+b.y);}
	Complex operator * (const Complex &b)const{return Complex(x*b.x-y*b.y,x*b.y+y*b.x);}
}x1[N],x2[N];
inline void change(Complex y[],int len){//二进制反转
	int i,j,k;
	for(i=1,j=len/2;i<len-1;i++){
		if(i<j) swap(y[i],y[j]);
		k=len/2;
		while(j>=k) j=j-k;k=k/2;//模拟减法的进位
		j+=k;//进位
	}
}
inline void FFT(Complex y[],int len,int on){//on 在 1 时是 DFT，-1 时是IDFT
	change(y,len);
	for(int h=1;h<=len;h<<=1){//模拟递归
		Complex wn(cos(2*PI/h),sin(on*2*PI/h));//注意在递归时每次自变量会平方
		for(int j=0;j<len;j+=h){
			Complex w(1,0);
			for(int k=j;k<j+h/2;k++){
				Complex u=y[k],t=w*y[k+h/2];
				y[k]=u+t;y[k+h/2]=u-t;
				w=w*wn;
			}
		}
	}
	if(on==-1) for(int i=0;i<len;i++) y[i].x/=len;
}
int str1[N],str2[N],sum[N];
signed main(){
	int len1=read()+1,len2=read()+1,len=1;
	for(int i=0;i<len1;i++) str1[i]=read();
	for(int i=0;i<len2;i++) str2[i]=read();
	while(len<=(len1+len2)) len<<=1;//注意在递归进行时，每次lim会减半，中点两侧的值都需要计算。所以我把lim变为2的整次幂，方便递归与计算
	for(int i=0;i<len1;i++) x1[i]=Complex(str1[i],0);
	for(int i=0;i<len2;i++) x2[i]=Complex(str2[i],0);
	FFT(x1,len,1);FFT(x2,len,1);
	for(int i=0;i<len;i++) x1[i]=x1[i]*x2[i];
	FFT(x1,len,-1);
	for(int i=0;i<len;i++) sum[i]=(int)(x1[i].x+0.5);
	len=len1+len2-1;
	for(int i=0;i<len;i++) cout<<sum[i]<<" ";
	cout<<"
";
}

NTT

在某些时候，我们需要求模 p 意义下的卷积

前置知识

1. 阶

   如果 gcd(a,p)=1，那么对于方程 (a^requiv1pmod p) 来说，根据欧拉定理 (a^{varphi(p)}equiv1pmod p)，一定存在解 (r|varphi(p))

   最小的 r 称为 a 关于 p 的阶，记作 (delta_p(a))

   求阶只能从小到大枚举

2. 原根

   在模 p 意义下的 0~p-1 次幂各不相同，取遍 [0,p-1]，也就是说 (delta_p(g)=varphi(p)) (p为质数)

   关于原根的几个性质：

   1. 注意到所有的素数都有原根

   2. 当正整数 p 有原根时，原根的个数为 (varphi(varphi(p)))

   3. 令 p 的原根是 g，那么 (g,g^2,cdots,g^{varphi(p)}) 构成了模 p 简化剩余系

实现

先求出 p 的原根 g，方法如下

从小到大枚举 (din [2,p-1])，令 (a_i) 为 (varphi(p)) 的素因子，若 (d^{frac{varphi(p)}{a_i}}equiv1pmod p) 都不成立，那么 d 就是模 p 的一个原根

这种方法的证明：

假设对于 (d^{frac{varphi(p)}{a_i}}equiv1pmod p) 均不成立，但存在 (k<varphi(p))，使得 (dkequiv1pmod p)

令 (q=gcd(k,varphi(p)))，那么 q<p 并且 x 一定整除 (dfrac{varphi(p)}{a_i}) 中的一个，不妨设它整除 (dfrac{varphi(p)}{a_i})

考虑等式 (kn+varphi(p)m=q)，容易发现一定存在正整数解

由于 (d^qequiv d^{kn+φ(p)m}≡1pmod p)，而 (q|dfrac{varphi(p)}{a_i})，与假设矛盾

得证

在求出原根后，可以发现，(g^frac{p-1}n) 与 (w_n) 的性质类似

所以我们可以用 (g^frac{p-1}n) 来代替 (w_n)

和FFT几乎相同

PS:这种方法只在 (p=a*2^k+1) 的情况下才成立

关于任意模数 NTT：不会写

三模数 NTT

首先挑选 3 个 (a*2^k+1) 形式的模数，进行3次 NTT，求出3组答案

然后将每一个数用中国剩余定理合并即可

inline void NTT(int y[],int len,int on){
	Change(y,len);
	for(int h=1,mul,u,t;h<=len;h<<=1){
		if(on==1) mul=ksm(G,(mod-1)/h);
		else mul=ksm(Ginv,(mod-1)/h);
		for(int j=0,w;j<len;j+=h){
			w=1;
			for(int k=j;k<j+(h>>1);k++){
				u=y[k];t=(w*y[k+(h>>1)])%mod;
				y[k]=(u+t)%mod;y[k+(h>>1)]=(u-t+mod)%mod;
				w=(w*mul)%mod;
			}
		}
	}
	if(on==-1){
		int mul=ksm(len,mod-2);
		for(int i=0;i<len;i++) y[i]=(y[i]*mul)%mod;
	}
}

多项式求逆

主要思想是递归处理，每次减半

首先我们在 (mod x^{n}) 处的逆元可以费马小定理直接求

现在假设我们现在求出来了 (F*Gequiv1pmod {x^{lceilfrac n2 ceil}}) , 要
求 (F*G'equiv 1pmod {x^n})

那么

(egin{aligned}G'&equiv G&&pmod {x^{lceilfrac n2 ceil}}\(G'-G)^2&equiv0&&pmod {x^n}\F*(G'-G)^2&equiv0&&pmod {x^n}\F*G'^2-2F*G*G'+F*G^2&equiv0&&pmod {x^n}\G'-2G+F*G^2&equiv0&&pmod {x^n}\G'&equiv 2G-F*G^2&&pmod {x^n}end{aligned})

那么现在我们就可以 (O(nlog n)) 的时间递归求逆了

void INV(int x[],int inv[],int len){
	if(len==1){inv[0]=ksm(x[0],mod-2);return ;}
	INV(x,inv,(len+1)>>1);
	int lim=1;
	while(lim<(len<<1)) lim<<=1;
	for(int i=0;i<len;i++) c[i]=x[i];
	for(int i=len;i<lim;i++) c[i]=0;//一定要注意各种清空
	NTT(c,lim,1);NTT(inv,lim,1);
	for(int i=0;i<lim;i++) inv[i]=(inv[i]*((2-c[i]*inv[i]%mod)%mod+mod)%mod)%mod;
	NTT(inv,lim,-1);
	for(int i=len;i<lim;i++) inv[i]=0;//一定要注意各种清空
}

多项式求ln

前置知识：求导(高中数学)

假设我们现在有多项式 (G)，需要求 (LNequiv ln(G)pmod{x^n})

由于多项式的求导、求积可以以 (O(n)) 的时间复杂度方便的做到，所以考虑对等式
两侧同时求导

那么
(egin{aligned}LN&equiv ln(G)&&pmod{x^n}\LN'&equiv dfrac{G'}{G}&&pmod{x^n}\end{aligned})

对于(dfrac1{G})，我们通过多项式求逆 (O(nlog n)) 解决

对于(G')，我们可以 (O(n)) 解决

相乘，再求积即可

void LN(int x[],int ln[],int len){
	if(len==1){ln[0]=0;return ;}
	INV(x,ln,len);
	int lim=1;
	while(lim<(len<<1)) lim<<=1;
	for(int i=0;i<len-1;i++) d[i]=(x[i+1]*(i+1))%mod;
	for(int i=len-1;i<lim;i++) d[i]=0;
	NTT(ln,lim,1);NTT(d,lim,1);
	for(int i=0;i<lim;i++) ln[i]=(ln[i]*d[i])%mod;
	NTT(ln,lim,-1);
	for(int i=lim-1;i>=0;i--) ln[i+1]=ln[i]*ksm(i+1,mod-2)%mod;
	ln[0]=0;
}

多项式求exp

前置知识

1. 泰勒展开

   (displaystyle f(x)=sum_{i=0}^{infty}dfrac{f^{(i)}(x_0)}{i!}(x-x_0)^i) ，其中 (x ightarrow x_0)
   主要用于将一个函数表示为多项式形式
   由于 (x ightarrow x_0) ，所以 (dfrac{f^{(i)}(x_0)}{i!}(x-x_0)^i ightarrow 0)
   因此有时可以将 f(x) 近似为 (displaystylesum_{i=0}^ndfrac{f^{(i)}(x_0)}{i!}(x-x_0)^i)

2. 牛顿迭代

   给定一个多项式函数 F(X)，求一个多项式 X 使得 (F(X)equiv0pmod{x^n})
   假设已经迭代得到 (F( X_0)equiv0pmod{x^{lceilfrac n2 ceil}})，考虑将 (F(X)) 在 (X_0) 处展开，那么 (egin{aligned}sum_{i=0}^ndfrac{F^{(i)}(X_0)}{i!}(X-X_0)^i&equiv0&&pmod{x^n}end{aligned})
   发现当 i>1 时，有 ((X-X_0)^iequiv0pmod{x^n})
   所以有
   (egin{aligned}F(X_0)+F'(X_0)(X-X_0)&equiv0&&pmod{x^n}\X&equiv X_0-dfrac{F(X_0)}{F'(X_0)}&&pmod{x^n}end{aligned})

假设我们现在有多项式 G，需要求 (Xequiv e^Gpmod{x^n})

变形可得

(egin{aligned}ln(X)-G&equiv0&&pmod{x^n}end{aligned})

设 (F(X)=ln(X)-G)，代入牛顿迭代，得

(egin{aligned}X&equiv X_0-dfrac{ln(X_0)-G}{dfrac{X'}{X}}&&pmod{x^n}\X&equiv X_0-dfrac{ln(X_0)*X-GX}{X'}&&pmod{x^n}end{aligned})

void EXP(int x[],int exp[],int len){
	if(len==1){
		exp[0]=1;return ;
	}
	EXP(x,exp,(len+1)>>1);
	LN(exp,e,len);
	int lim=1;
	while(lim<(len<<1)) lim<<=1;
	for(int i=0;i<len;i++) e[i]=((x[i]-e[i]+(i==0))%mod+mod)%mod;
	for(int i=len;i<lim;i++) e[i]=0;
	NTT(exp,lim,1);NTT(e,lim,1);
	for(int i=0;i<lim;i++) exp[i]=(exp[i]*e[i])%mod;
	NTT(exp,lim,-1);
	for(int i=len;i<lim;i++) exp[i]=0;
	for(int i=0;i<lim;i++) e[i]=0;
}

分治NTT

假设有一个形式为 (f[i]=displaystylesum_{j=1}^if[i-j]g[j])

它显然不是卷积，但是差别仅在于等号的右侧也出现了 f

处理方法类似 CDQ 分治

假设已经求出了 f(0) 到 (f(frac n2))，要计算它们对未知部分 (f(frac n2+1)) 到 f(n) 的贡献

发现一个已知的 f(x)，对任意一个 f(i)(i>x) 的贡献都只有 f(x)g(i-x)

于是贡献函数 F 就长这样了:

(F(x)=displaystylesum_{i=1}^xf(i)g(x-i))

其中 f 只填充 [l,mid]​ 的部分，因为其他地方的 f=0

g 只填充 [1,r-l]​ 的部分，因为只算对于 [mid,r] 的贡献

void CDQ(int l,int r){
	if(l==r) return ;
	int mid=(l+r)>>1;
	CDQ(l,mid);
	for(int i=0;i<r-l+1;i++) A[i]=B[i]=0;
	for(int i=l;i<=mid;i++) A[i-l]=f[i];
	for(int i=1;i<r-l+1;i++) B[i]=g[i];
	NTT(A,r-l+1,1);NTT(B,r-l+1,1);
	for(int i=0;i<r-l+1;i++) A[i]=(A[i]*B[i])%mod;
	NTT(A,r-l+1,-1);
	for(int i=mid+1;i<=r;i++) f[i]=(f[i]+A[i-l])%mod;
	CDQ(mid+1,r);
}
signed main(){
	while(lim<n) lim<<=1;//先扩展为2^x的形式，方便处理
	CDQ(0,lim-1);
}

Lagrange插值

CRT

若(m_1,m_2,cdots,m_n)是两两互质的正整数，(M=prod_{i=1}^n{m_i})，(M_i=M/m_i) , (t_i)是线性同余方程(M_it_iequiv 1 (mod m_i))的一个解

对于任意的n个整数(a_1,a_2,cdots,a_n)，则同余方程组： (egin{cases}x≡a_1(mod m_1)\x≡a_2(mod m_2)\ cdots cdots\x≡a_n(mod m_n)\end{cases}) 有整数解(x=displaystylesum_{i=1}^n M_it_i a_i).并且在(mod M) 意义下有唯一解。

证明

因为(M_i=M/m_i) 是除(m_i)之外所有模数的倍数，所以(forall k ot=i , M_it_ia_iequiv 0 (mod m_k))

又因为(M_it_ia_iequiv a_i (mod m_i))
所以 (x=displaystylesum_{i=1}^n M_it_i a_i)

结论

CRT给出了模数两两互质的线性同余方程组的一个特解。方程组的通解可以表示为 (x+kM (kinmathbb Z))。有些题目要求我们求出最小的非负整数解，只需把 (x) 对 (M) 取模，并让x落在 ([1,M)) 内即可。

正文

对于一个多项式 (f(x))，其图像在坐标系内经过 (n) 个点 ((x_i,y_i))

我们考虑对于此多项式的“CRT”：

构造 (n) 个多项式 (g_i(n)).对于第 (i) 个多项式，对于(forall k ot= i ，g_i(x_k)=0)，而 (g_i(x_i)=1)，即 (g_i(x_i)*y_i=y_i)

则可得 (g_i(x)=dfrac{(x-x_1)(x-x_2)cdots(x-x_{i-1})(x-x_{i+1})cdots(x-x_n)}{(x_i-x_1)(x_i-x_2)cdots(x_i-x_{i-1})(x_i-x_{i+1})cdots(x_i-x_n)})

所以 (f(x)=displaystylesum_{i=1}^n g_i(x)*y_i)

多项式快速插值