多项式全家桶

$\text{warning: }$ 目前为止本博客仍是半半半半半成品。

以下是本文参考以及推荐阅读的相关博客$\mathrm{\%}\mathrm{\%}\mathrm{\%}$
$p.s.$ 博客都写得很好！没懂是蒟蒻个人问题

点击查看相关博客$Link$

OI-wiki 快速傅里叶变换：难得感觉有点抽象。

OI - wiki 拉格朗日插值：配合其他博客食用更佳。

傅里叶变换(FFT)学习笔记：从前置芝士讲起的详细学习笔记。

小学生都能看懂的FFT！！！：比较通俗并且不枯燥的讲解。

快速傅里叶变换(FFT)详解：证明很详细。

多项式初步与傅里叶变换（Fourier Transform） : 清晰但是初看有一丢没看懂，可以康康迭代部分。

多项式 I：拉格朗日插值与快速傅里叶变换：感觉已经算是全家桶了（）通过这篇学会了拉插！

从拉插到快速插值求值：超级详细的插值笔记。

---

概念

• $\mathrm{DFT}:$ 系数表达转换为点值表达。

• $\mathrm{IDFT}:$ 点值表达转换为系数表达。（ $\mathrm{DFT}$的逆运算 ）

---

拉格朗日插值

考虑构造 $f_i(x)$ ， 使得 $f_i(x_i)=1$ 且 $f_i(x_j)=0(j\ne i)$。
首先考虑满足后者，则有: $f_i(x)=\prod _{i\ne j}(x-x_j)$。
再考虑前者，则在两边同乘 $\frac{y_i}{f_i(x_i)}$ ，可以得到：$f_i(x)=y_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$。
多项式 $f(x)$ 即为：

$$\begin{array}{rl}f(x)& =\sum _{i=1}^n{y}_i{f}_i(x)\\ & =\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\end{array}$$

注意最后一起算分母的逆元，而不是分别求解。
模版题：luogu P4781

点击查看代码

#include<bits/stdc++.h>
#define cs const
#define il inline
#define LL long long
#define pc(i) putchar(i)
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k) for(int i=(j);(i)>=(k);--(i))
using namespace std;
mt19937 srd(time(0));
cs int inf=0x3f3f3f3f,Mod=998244353,N=2003;
int FL,CH;
template<typename T> bool in(T&a)
{
    for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
        if(CH=='-')FL=-1;
    for(a=0;isdigit(CH);CH=getchar())
        a=a*10+CH-'0';
    return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
void wt(int x){if(x<0)x=-x,pc('-');if(x>9)wt(x/10);pc(x%10|48);}
int x[N],y[N],n,k; LL ans;
il int qwq(cs int x,cs int y){return (x-y+Mod)%Mod;}
il int qpow(int b,int p)
{
	int r=1;
	while(p>0)
	{
		if(p&1) r=1ll*r*b%Mod;
		p>>=1,b=1ll*b*b%Mod;
	}
	return r;
}
il int inv(cs int x){return qpow(x,Mod-2);} 
signed main()
{
	in(n,k); fo(i,1,n) in(x[i],y[i]);
	fo(i,1,n) 
	{
		int s1=y[i]%Mod,s2=1; 
		fo(j,1,n) if(i!=j)
			s1=1ll*s1*qwq(k,x[j])%Mod,s2=1ll*s2*qwq(x[i],x[j])%Mod;
		ans=(ans+1ll*s1*inv(s2)%Mod)%Mod;
	}
	wt(ans);
    return 0;
}

横坐标是连续整数的拉格朗日插值

省流：可以 $O(n)$。
横坐标为 $1,\dots ,n$的插值公式：

$$\begin{array}{rl}f(x)& =\sum _{i=1}^n{y}_i\prod _{i\ne j}\frac{x-x_j}{x_i-x_j}\\ & =\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-j}{i-j}\\ & =\sum _{i=1}^n{y}_i\frac{\frac{\prod _{j=1}^n(x-j)}{x-i}}{(-1{)}^{n-i}\cdot (i-1)!\cdot (n-i)!}\\ & =\sum _{i=1}^n{y}_i\cdot \frac{\prod _{j=1}^n(x-j)}{(x-i)\cdot (-1{)}^{n-i}\cdot (i-1)!\cdot (n-i)!}\end{array}$$

例题：CF622F The Sum of the k-th Powers

题目大意：求 $\sum _{i=1}^n{i}^k,n\le {10}^9,k\le {10}^6$

答案是一个 $k+1$ 次的多项式，找 $k+2$ 个点带进去拉插，证明在这里找。
预处理阶乘和前后缀积。

点击查看代码

#include<bits/stdc++.h>
#define cs const
#define il inline
#define LL long long
#define pc(i) putchar(i)
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k) for(int i=(j);(i)>=(k);--(i))
using namespace std;
mt19937 srd(time(0));
cs int inf=0x3f3f3f3f,Mod=1e9+7,N=1e6+7;
int FL,CH;
template<typename T> bool in(T&a)
{
    for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
        if(CH=='-')FL=-1;
    for(a=0;isdigit(CH);CH=getchar())
        a=a*10+CH-'0';
    return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
void wt(int x){if(x<0)x=-x,pc('-');if(x>9)wt(x/10);pc(x%10|48);}
int n,k,ans,fac[N],L[N],R[N],y,num,den;
il int qpow(int b,int p)
{
	int re=1;
	while(p>0)
	{
		if(p&1) re=1ll*re*b%Mod;
		p>>=1,b=1ll*b*b%Mod;
	}
	return re;
}
il int inv(cs int x){ return qpow(x,Mod-2); }
signed main()
{
	in(n,k),L[0]=R[k+3]=fac[0]=1;
	fo(i,1,k+2) L[i]=1ll*L[i-1]*(n-i)%Mod,fac[i]=1ll*fac[i-1]*i%Mod;
	of(i,k+2,1) R[i]=1ll*R[i+1]*(n-i)%Mod;
	fo(i,1,k+2)
	{
		y=(y+qpow(i,k))%Mod,num=1ll*L[i-1]*R[i+1]%Mod,
		den=1ll*fac[i-1]*((k-i)&1?-1:1)*fac[k+2-i]%Mod;
		ans=(ans+1ll*y*num%Mod*inv(den)%Mod)%Mod; 
	}
	wt((ans+Mod)%Mod); 
    return 0;
}

---

$\mathrm{FFT}$

复数乘法

$$\begin{array}{rl}& (a+bi)\ast (c+di)\\ =& ac+adi+bci+bd{i}^2\\ =& ac+adi+bci-bd\\ =& (ac-bd)+(bc+ad)i\end{array}$$

单位根

• ${\omega }_n^k=\cos k\times \frac{2\pi }{n}+i\sin k\times \frac{2\pi }{n}$
• ${\omega }_{2n}^{2k}={\omega }_n^k$
• ${\omega }_n^0={\omega }_n^n=1$

快速傅里叶变换

$A(x)=a_0+a_1\times x+a_2\times x^2+a_3\times x^3+a_4\times x^4+\cdots +a_{n-2}\times x^{n-2}+a_{n-1}\times x^{n-1}$

按照下标的奇偶性分类：

$A(x)=(a_0+a_2\times x^2+a_4\times x^4+\cdots +a_{n-2}\times x^{n-2})+(a_1\times x+a_3\times x^3+a_5\times x^5\dots a_{n-1}\times x^{n-1})$

设：

$A_1(x)=a_0+a_2\times x+a_4\times x^2+\cdots +a_{n-2}\times x^{\frac{n}{2}-1}$

$A_2(x)=a_1+a_3\times x+a_5\times x^2+\cdots +a_{n-1}\times x^{\frac{n}{2}-1}$

则有： $A(x)=A_1(x^2)+x{A}_2(x^2)$

将 ${\omega }_n^k(k<\frac{n}{2})$ 代入：$A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})=A_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k)$

将 ${\omega }_n^{k+\frac{n}{2}}$ 代入： $A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{A}_2({\omega }_n^{2k+n})=A_1({\omega }_n^{2k})-{\omega }_n^k{A}_2({\omega }_n^{2k})$

快速傅里叶逆变换

点值表示法转化为系数表示法。

$(y_0,y_1,\dots ,y_{n-1})$ 为 $(a_0,a_1,\dots ,a_{n-1})$ 的点值表示法。

设有一向量 $(c_0,c_1,\dots ,c_{n-1})$ 满足 $c_k=\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i$ ，即多项式 $B(x)=y_0+y_1\times x+y_2\times x^2+\cdots +y_{n-1}\times x^{n-1}$ 在 ${\omega }_n^0,{\omega }_n^{-1},{\omega }_n^{-2},\dots ,{\omega }_n^{-(n-1)}$ 处的点值表示。

$$\begin{array}{rl}{c}_k& =\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j)({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i)({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i)\\ & =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i\\ & =\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)\end{array}$$

设 $S(x)=\sum _{i=0}^{n-1}{x}^i$，将 ${\omega }_n^k$ 代入：$S({\omega }_n^k)=1+({\omega }_n^k)+({\omega }_n^k{)}^2+\cdots +({\omega }_n^k{)}^{n-1}$

当 $k\ne 0$ 时，同乘 ${\omega }_n^k$ ：${\omega }_n^{k}S({\omega }_n^k)={\omega }_n^k+({\omega }_n^k)+({\omega }_n^k{)}^2+({\omega }_n^k{)}^3+\cdots +({\omega }_n^k{)}^n$

上下两式相减可得：

$$\begin{array}{rl}& {\omega }_n^{k}S({\omega }_n^k)-S({\omega }_n^k)=({\omega }_n^k{)}^n-1\\ & S({\omega }_n^k)=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}\\ & S({\omega }_n^k)=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}\\ & S({\omega }_n^k)=\frac{1-1}{{\omega }_n^k-1}\end{array}$$

因此，$k\ne 0$时，$S({\omega }_n^k)=0$，而 $k=0$ 时 $S({\omega }_n^0)=n$ 。

对于

$$c_k=\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)$$

当 $j\ne k$ 时，值为 $0$ ，而 $j=k$ 时值为 $n$。

因此，$c_k=n{a}_k$ ，$a_k=\frac{c_k}{n}$。

迭代实现

$$\begin{array}{}\underset{―}{a_0{a}_1{a}_2{a}_3}|\underset{―}{a_4{a}_5{a}_6{a}_7}\\ \underset{―}{a_0{a}_2}\left|\underset{―}{a_4{a}_6}\right|\underset{―}{a_1{a}_3}\mid \underset{―}{a_5}{a}_7\\ \underset{―}{a_0}\left|\underset{―}{a_4}\right|\underset{―}{a_2}\left|\underset{―}{a_6}\right|\underset{―}{a_1}\left|\underset{―}{a_5}\right|\underset{―}{a_3}|\underset{―}{a_7}\end{array}$$

$$\begin{array}{}\underset{―}{000,001,010,011}|\underset{―}{100,101,110,111}\\ |\underset{―}{000,010}|\underset{―}{100,110}|\underset{―}{001,011}|\underset{―}{101,111}\\ \underset{―}{000}|\underset{―}{100}|\underset{―}{010}|\underset{―}{110}|\underset{―}{001}|\underset{―}{101}|\underset{―}{011}|\underset{―}{111}\end{array}$$

在原序列中 $i$ 与 $\frac{i}{2}$ 的关系是 ： $i$ 可以看做是 $\frac{i}{2}$ 的二进制上的每一位左移一位得来，那么在反转后的数组中就需要右移一位，同时特殊处理一下奇数。
记为 $rev[i]=(rev[i>>1]>>1)|((rev[i]\mathrm{\&}1)<<(bit-1)),bit=\lceil {\log }_{2}n\rceil$。

$\mathrm{FFT}$常数优化论文
P3803 【模板】多项式乘法（FFT）

点击查看代码

#include<bits/stdc++.h>
#define cs const
#define il inline
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k),for(int i=(j);(i)>=(k);--(i))
#define LL long long
using namespace std;
cs int N=3e6+9;
cs double pi=acos(-1);
int FL,CH;
template<typename T>bool in(T&a)
{
    for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
        if(CH=='-')FL=-1;
    for(a=0;isdigit(CH);CH=getchar())
        a=a*10+CH-'0';
    return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
struct Complex
{
    double x,y;
    Complex operator+(const Complex &a) const
    { return {x+a.x,y+a.y}; }
    Complex operator-(const Complex &a) const
    { return {x-a.x,y-a.y}; }
    Complex operator*(const Complex &a) const
    { return {x*a.x-y*a.y,x*a.y+y*a.x}; }
}A[N],B[N];
int rev[N],bit,tot,n,m;
il void FFT(Complex a[],int inv)
{
    fo(i,0,tot-1) if(i<rev[i]) swap(a[i],a[rev[i]]); //求出要迭代的序列
    for(int mid=1;mid<tot;mid<<=1)//待合并区间的长度的一半
    {
        Complex w1=(Complex){cos(pi/mid),inv*sin(pi/mid)};//单位根
        for(int i=0;i<tot;i+=mid*2)
        {
            Complex wk=(Complex){1,0};
            for(int j=0;j<mid;++j,wk=wk*w1)//枚举左半部分
            {
                Complex x=a[i+j],y=wk*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;
            }
        }
    }
}
signed main()
{
    in(n,m);
    fo(i,0,n) scanf("%lf",&A[i].x);
    fo(i,0,m) scanf("%lf",&B[i].x);
    while((1<<bit)<n+m+1) ++bit; //bit=ceil(log_2^n)
    tot=1<<bit;
    fo(i,0,tot-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    FFT(A,1),FFT(B,1);
    fo(i,0,tot-1) A[i]=A[i]*B[i];
    FFT(A,-1);
    fo(i,0,n+m) printf("%d ",(int)(A[i].x/tot+0.5));
    return 0;
}

P1919 【模板】A*B Problem 升级版（FFT 快速傅里叶变换）

点击查看代码

#include<bits/stdc++.h>
#define cs const
#define il inline
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k) for(int i=(j);(i)>=(k);--(i))
#define LL long long
using namespace std;
cs int N=3e6+9;
cs double pi=acos(-1);
int FL,CH;
template<typename T>bool in(T&a)
{
    for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
        if(CH=='-')FL=-1;
    for(a=0;isdigit(CH);CH=getchar())
        a=a*10+CH-'0';
    return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
struct Complex
{
    double x,y;
    Complex operator+(cs Complex &a) cs
    { return {x+a.x,y+a.y}; }
    Complex operator-(cs Complex &a) cs
    { return {x-a.x,y-a.y}; }
    Complex operator*(cs Complex &a) cs
    { return {a.x*x-y*a.y,x*a.y+y*a.x}; }
}A[N],B[N];
char ca[N],cb[N];
int rev[N],bit,tot,n,m,ans[N],id;
il void FFT(Complex a[],cs int inv)
{
    fo(i,0,tot-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int mid=1;mid<tot;mid<<=1)
    {
        Complex w1=(Complex){cos(pi/mid),inv*sin(pi/mid)};
        for(int i=0;i<tot;i+=mid*2)
        {
            Complex wk=(Complex){1,0};
            for(int j=0;j<mid;++j,wk=wk*w1)
            {
                Complex x=a[i+j],y=wk*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;
            }
        }
    }
}
signed main()
{
    scanf("%s%s",ca,cb);
    n=strlen(ca)-1,m=strlen(cb)-1;
    fo(i,0,n) A[i].x=(ca[n-i]^48);
    fo(i,0,m) B[i].x=(cb[m-i]^48);
    while((1<<bit)<n+m-1) ++bit; 
    tot=1<<bit;
    fo(i,0,tot-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    FFT(A,1),FFT(B,1);
    fo(i,0,tot-1) A[i]=A[i]*B[i];
    FFT(A,-1); 
    for(int i=0,x=0;i<tot||x;++i)
        x+=A[i].x/tot+0.5,ans[++id]=x%10,x/=10;
    while(id>2&&(!ans[id-1])) --id;
    of(i,id-1,1) printf("%d",ans[i]);
    return 0;
}