一点多项式相关
都学不会 qwq
随时更新，主要给自己看。
一. 拉格朗日插值
给出 \(n\) 个点值，求出多项式各项系数。
设第 \(i\) 个点的 \(x\) 值为 \(x_i\)，\(y\) 值为 \(y_i\)。
\[f(x)=\sum\limits_{i=1}^ny_i\cdot\prod_{i\ne j}\dfrac{x-x_j}{x_i-x_j}
\]
解释就不放了，变换一下就能看明白。
时间复杂度 \(O(n^2)\)
代码：
P4781
点击查看代码
#include<bits/stdc++.h>
#define ld long double
#define ll long long
using namespace std;
const int p=998244353;
int n,k,x[2005],y[2005],ans;
int ksm(int x,int mx) {
int anss=1;
while(mx) {
if(mx&1) anss=1ll*anss*x%p;
x=1ll*x*x%p;
mx>>=1;
}
return anss;
}
int main() {
ios_base::sync_with_stdio(false);
cin.tie(NULL); cout.tie(NULL);
cin>>n>>k;
for(int i=1;i<=n;++i) cin>>x[i]>>y[i];
for(int i=1;i<=n;++i) {
int pre=1,suf=1;
for(int j=1;j<=n;++j) {
if(i==j) continue;
pre=1ll*pre*(x[i]-x[j]+p)%p;
suf=1ll*suf*(k-x[j]+p)%p;
}
ans=(ans+1ll*y[i]*suf%p*ksm(pre,p-2))%p;
}
cout<<ans<<"\n";
return 0;
}
二. 快速傅里叶变换
1. DFT
较为快速的求解加法卷积的方法，一般用于加速乘法。
本质是以 \(O(n\log(n))\) 的复杂度将系数表达式映射到点值表达式，并在点值乘法之后将点值表达式映射到系数表达式。
1.1 引理
首先，先看一个欧拉公式
\[e^{ix}=cosx+isinx
\]
若取 \(x=2\pi\)
\[w=e^{\tfrac{2\pi i}{n}}
\]
有几个引理：
- \(w_n^k=w_{2n}^{2k}\)
- \(w_n^k=w_n^{k+tn}(t\in \mathbb{Z})\)
- \(w_n^k=-w_n^{k+\frac{n}{2}}(2\mid n)\)
- \(\sum_{i=0}^{n-1}(w_n^k)^i=0(n>0\land k\nmid n)\)
\(w\) 的几个引理都是比较容易证明的，这里不再多说（如果需要可以看 花花的博客）。
1.2 分治
对于原多项式，\(a\) 为各项系数，有
\[f(w_n^{k})=\sum_{i=0}^{n-1}a_i\cdot w_n^{kj}
\]
构建一个范德蒙德矩阵
\[\begin{bmatrix}(w_n^0)^0&\cdots&(w_n^0)^{n-1}\\\vdots&\ddots&\vdots\\(w_n^{n-1})^0&\cdots&(w_n^{n-1})^{n-1}\end{bmatrix}\times\begin{bmatrix}a_0\\\vdots\\a_{n-1}\end{bmatrix}=\begin{bmatrix}f(w_n^0)\\\vdots\\f(w_n^{n-1})\end{bmatrix}
\]
算法复杂度为 \(O(n\log{n})\)，容易想到需要用分治算法。
考虑将原式按多项式系数奇偶分裂为两个多项式
\[f(x)=f_0(x^2)+xf_1(x^2)
\]
\[f_0(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}
\]
\[f_1(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1}
\]
这里需要注意，为了方便分治，\(n\) 应当补齐到 \(2\) 的幂次。
对于 \(0\le k<\frac{L}{2}\)，有:
\[f(w_L^{k})=f_0(w_L^{2k})+w_L^kf_1(w_L^{2k})
\]
\[f(w_L^{k})=f_0(w_k^{2k+L})+w_k^{k+\frac{L}{2}}f_1(w_k^{2k+L})
\]
由上面的引理可以得到：
\[f(w_L^{k})=f_0(w_{\frac{L}{2}}^k)+w_L^kf_1(w_{\frac{L}{2}}^k)
\]
\[f(w_k^{k+\frac{L}{2}})=f_0(w_{\frac{L}{2}}^k)+w_L^kf_1(w_{\frac{L}{2}}^k)
\]
用这种方法，我们只需要求出 \(L/2\) 次多项式在 \(L/2\) 次单位根处的取值，怎么分治就非常明了了。
1.3 蝴蝶变换
递归实现无疑很慢，那么试图用递推实现，这里就需要引入蝴蝶变换。
设 \(p\) 为系数原位置，很容易理解，第 \(i\) 次递归时的方向决定了下标二进制下 \(w-i\) 位的取值，而它又被 \(\left\lfloor\dfrac{p}{2^{i-1}}\right\rfloor\) 影响，可以得出其位置 \(r\) 为 \(p\) 二进制翻转后的结果。
也可以递推求 \(r\)，\(r_i\) 的第 \(0\) 到 \(w-2\) 位可以由 \(r_{i>>1}\) 得到，因为它们的这几位相同，第 \(w-1\) 位只需要判断 \(i\) 的奇偶性即可。
具体流程：先进行蝴蝶变换，然后枚举区块大小，枚举区块起点，枚举遍历整个区块。
2. IDFT
接下来是逆变换。
显然，只要乘上一个逆矩阵就可以了。这个我不太会严谨证明 qaq，之后会了填坑。注意到两次运算其实十分相似，我们可以用一个函数解决。
\[\dfrac{1}{n}\begin{bmatrix}(w_n^{-0})^0&\cdots&(w_n^{-0})^{n-1}\\\vdots&\ddots&\vdots\\(w_n^{-(n-1)})^0&\cdots&(w_n^{-(n-1)})^{n-1}\end{bmatrix}
\]
PS：由于浮点数精度误差，最后结果四舍五入。
代码：
P3803
点击查看代码
#include<bits/stdc++.h>
#define ld long double
#define ll long long
using namespace std;
const ld pi=acos(-1);
int n,m,L,r[1<<21];
struct comp
{
ld a,b;
comp operator +(comp &x){ return {a+x.a,b+x.b}; }
comp operator -(comp &x){ return {a-x.a,b-x.b}; }
comp operator *(comp &x){ return {a*x.a-b*x.b,a*x.b+b*x.a}; }
}f[1<<21],g[1<<21];
void FFT(comp *f,int type)
{
for(int i=1;i<L;++i)
{
r[i]=(r[i>>1]>>1)+(i&1?L>>1:0);
if(i<r[i]) swap(f[i],f[r[i]]);
}
for(int k=1;k<L;k<<=1)
{
comp wk={cos(pi/k),sin(pi/k)};
wk.b*=type;
static comp w[1<<21];
w[0]={1,0};
for(int i=1;i<k;++i) w[i]=w[i-1]*wk;
for(int i=0;i<L;i+=k<<1)
for(int j=0;j<k;++j)
{
comp x=f[i|j],y=w[j]*f[i|j|k];
f[i|j]=x+y,f[i|j|k]=x-y;
}
}
if(type==-1) for(int i=0;i<L;++i) f[i].a/=L;
}
int main()
{
ios_base::sync_with_stdio(false);
cin.tie(NULL);
cout.tie(NULL);
cin>>n>>m;
for(int i=0;i<=n;++i) cin>>f[i].a;
for(int i=0;i<=m;++i) cin>>g[i].a;
L=1;
while(L<=n+m) L<<=1;
FFT(f,1),FFT(g,1);
for(int i=0;i<L;++i) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=0;i<=n+m;++i) cout<<(int)(f[i].a+0.5)<<(i<n+m?' ':'\n');
return 0;
}
3. 一点常数优化：
3.1 三次变两次优化：
复数平方可得：
\[(a+bi)^2=a^2-b^2+2abi
\]
所以将另一个多项式存进虚部，直接平方，然后取出虚部的一半。
代码：
点击查看代码
#include<bits/stdc++.h>
#define ld long double
#define ll long long
using namespace std;
const ld pi=acos(-1);
int n,m,L,r[1<<21];
struct comp
{
ld a,b;
comp operator +(comp &x){ return {a+x.a,b+x.b}; }
comp operator -(comp &x){ return {a-x.a,b-x.b}; }
comp operator *(comp &x){ return {a*x.a-b*x.b,a*x.b+b*x.a}; }
comp operator /(ld x){ return {a/x,b/x}; }
}f[1<<21];
void FFT(comp *f,int type)
{
for(int i=1;i<L;++i)
{
r[i]=(r[i>>1]>>1)+(i&1?L>>1:0);
if(i<r[i]) swap(f[i],f[r[i]]);
}
for(int k=1;k<L;k<<=1)
{
comp wk={cos(pi/k),sin(pi/k)};
wk.b*=type;
static comp w[1<<21];
w[0]={1,0};
for(int i=1;i<k;++i) w[i]=w[i-1]*wk;
for(int i=0;i<L;i+=k<<1)
for(int j=0;j<k;++j)
{
comp x=f[i|j],y=w[j]*f[i|j|k];
f[i|j]=x+y,f[i|j|k]=x-y;
}
}
if(type==-1) for(int i=0;i<L;++i) f[i]=f[i]/L;
}
int main()
{
ios_base::sync_with_stdio(false);
cin.tie(NULL);
cout.tie(NULL);
cin>>n>>m;
for(int i=0;i<=n;++i) cin>>f[i].a;
for(int i=0;i<=m;++i) cin>>f[i].b;
L=1;
while(L<=n+m) L<<=1;
FFT(f,1);
for(int i=0;i<L;++i) f[i]=f[i]*f[i];
FFT(f,-1);
for(int i=0;i<=n+m;++i) cout<<(int)(f[i].b/2+0.5)<<(i<n+m?' ':'\n');
return 0;
}
