多项式各种算法学习笔记
1.FFT(快速傅里叶变换)
1.前置技能
复数:
基本表示法及性质:
\(i\)是虚数单位
1.坐标(代数)形式：
当b为0是z为实数，当a为0时为纯虚数
注:复数包括实数和虚数，虚数下有纯虚数
虚数z对应了复平面上的一点(a,b)
运算法则：
设复数\(z_1,z_2,z_1=a+bi,z_2=c+di\)
加法:\(z_1+z_2=(a+c)+(b+d)i\)
减法:\(z_1-z_2=(a-c)+(b-d)i\)
乘法:\(z_1*z_2=(ac-bd)+(bc+ad)i\)
除法:\(\dfrac{z_1}{z_2}=\dfrac{(ac+bd)+(bc-ad)i}{(c^2+d^2)}\)
2.三角形式
\(\theta\)是复数\(z\)的幅角,\(r\)是该复数的模长
这种形式下的乘法和除法运算更方便，通过和角公式，对于复数\(z_1=r_1(cos\theta_1+isin\theta_1),z_2=r_2(cos\theta_2+isin\theta_2)\)
那么：
几何意义:相当于把该复数 拉长/缩短 到另一个复数模长 倍/分之一 ,然后 顺时针/逆时针 旋转另一个复数的幅角大小的角度
于是有了如下非常有用的公式:
3.指数形式
于是我们知道了:\(e^{i\theta}=cos\theta+isin\theta\)
可以发现:$$e^{i\pi}=-1$$
(优美)
这个算乘除法就更好算了:
有:
上面的那个公式就可写成这样：
4.单位复数根
学FFT最重要的就这个了吧
设有如下方程:
这方程的复数根\(z\)为\(n\)次单位根,通常记为\(w\)
这样的根\(w\)有n个，也就是说\(n\)次单位根有\(n\)个,记为\(w_k (k=0,1,2,\dots n-1)\)
其中:
不难发现其实\(n\)次单位复数根就是把复平面上的单位圆 \(n\) 等分后的那些\(n\)等分点
一些次数单位的单位根举例:
1次单位根: \(1\)
2次单位根: \(1,-1\)
3次单位根: \(1,\frac{-1+\sqrt{3}i}{2},\frac{-1-\sqrt{3}i}{2}\)
\(\dots\)
其实就是个解二元n次方程
可以发现1是任意次的单位复数根，-1是任意偶数次单位复数根
某些引理:
1.消去引理
把复数的指数形式带进去就可以了
说人话就是不同次数的单位根可以互相转化
2.折半引理
假设\(n\)是大于0的偶数:
说人话就是n次单位复数根的前后两半平方后是对应相等的
3.求和引理
对于大于1的整数n和小于等于n的整数k有:
这就是个等比数列求和
2.步入正题
常规的一个最高次数为n-1的多项式的表示形式是系数表示法,如:
一共有n项
按照朴素的多项式乘法(卷积),就是每一项两两相乘，复杂度为\(O(n^2)\)
如果我们把多项式看成一个函数，我们取图像上的n个点来表示这个函数也即该多项式，这样的表示法叫做点值表示法
对于两个n-1次多项式,由于我们最后卷积出来的多项式是2n-2次的,如果我们知道了卷积后的多项式函数图像上的至少2n-1个点，那么我们就可以确定这个多项式了
所以有如下想法：
现在原来的两个多项式上选取好x值相同的点(个数为原来两多项式次数的和加1)，用\(O(n)\)的时间将选的点的y值相乘，得到的值用某种方法(待定系数法)转化为多项式系数的形式
如果选取的x都是n次单位复数根的k次方，那么以上两个过程就是\(DFT\)和\(IDFT\)了
1.\(DFT\)
离散傅里叶变换：
对于\(k\in [0,n-1],\)和n-1次多项式\(A(x),\)定义:
这个叫做离散傅里叶变换，记做\(y=DFT_n(a)\)
朴素求\(DFT\),是\(O(n^2)\)的
2.\(IDFT\)
逆离散傅里叶变换:
就是\(DFT\)的逆运算，用于求出多项式的系数a，记为\(DFT^{-1}\)
本人太弱不会证，丢个式子:
由于写的时候系数表达式和点值表达式是共用的数组，所以写起来和\(DFT\)没什么差别
3.\(FFT\)
快速傅里叶变换:
是一种快速求出\(DFT\)和\(DFT^{-1}\)的算法，利用了单位复数根的优良性质
我们先列出朴素求\(DFT\)的步骤,为了方便这里先不妨假设多项式次数为2的幂:
1.求出n次单位复数根的幂:\(w_n^0,w_n^1,w_n^2.....w_n^{n-1}\)
2.代入多项式\(A(x)\),求得以下式子:
根据折半引理,我们可以发现:
对于一个\(A(w_n^k)\),\((k<\frac{n}{2})\)它的每一个偶数次方项与\(A(w_n^{k+\frac{n}{2}})\)是对应相同的
以n=4的情况来考虑:
就是要求:
把每一个多项式的偶数次数项的系数提出来组成的多项式记为\(A^{[0]}(x)\),奇数项记为\(A^{[1]}(x)\)
那么有:
奇数项就没有这么优美的性质,但是我们可以提出来一个x使得它变成偶数项
所以我们现在考虑把一个多项式奇偶分组:
显然会有如下结论:
就是
折半引理化简一下:
那么我们知道如果把上面的分成两半,根据折半引理:前后两半对应的只有后面要乘上的单位复数根不一样,所以其实我们已经把问题的规模缩小了一半了,因为我们只需求解如下东西:
\(P.S.\)
其实我们还可以发现:
所以说其实只是后面的系数的符号不一样
再来观察一下现在要求的东西,可以发现如果把整个子问题也进行奇偶分组的话,每一组其实要求的是长度为原来一半的\(DFT\),如果采用递归进行处理,那么递归分组一次之后
要求的子问题是:
于是我们只要计算这4个子问题原\(DFT\)的两边就都可以直接算出来了
并且对于该子问题,我们也只需求解\(A^{[0]}(w_2^0)和A^{[1]}(w_2^0)\),另一半也是和上面的一样由这一半推出,因为这是一个长度为原来一半的\(DFT\)
我们先看看递归版应该怎么实现:
1.不停将要求的\(DFT\)进行奇偶分组并递归下去
2.如果只有一个元素,显然这时\(x=w_1^0=1\),直接返回当前的系数就可以了
3.合并答案:
把每层要求的东西写出来:
根据上面的公式就很容易知道每个元素是由下一层的哪一个得来的了,这很像个蝴蝶,于是被成为蝴蝶操作
还没有理解的话可以看看8的情况,建议自己手画一下,标出贡献来源,这样就很容易理解了
从下往上蝴蝶的跨度不断增大,倒数第2层就是相邻的两个元素宽度半径为1,然后往上不断倍增
,每一次只是改变两个位置对应的数,所以其实也可以写成非递归的
但是我们需要最后分组后的最下面一层
其实每个数最后在的位置是他的二进制反序数,所以是一个Rader排序(二进制反转),这个并不是非常重要,记一下就差不多了
贴个非递归的板子:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdlib>
using namespace std;
inline int read()
{
int x=0;char ch=getchar();int t=1;
for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=-1;
for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
return x*t;
}
typedef long long ll;
typedef double db;
const db PI=acos(-1);
struct Complex{//手写复数
db x,y;
Complex(db x1=0.0,db y1=0.0){x=x1,y=y1;}
inline Complex operator +(const Complex &b)const{return Complex(x+b.x,y+b.y);}
inline Complex operator -(const Complex &b)const{return Complex(x-b.x,y-b.y);}
inline Complex operator *(const Complex &b)const{return Complex(x*b.x-y*b.y,x*b.y+y*b.x);}
inline void init(){x=read();}
inline void out(){printf("%lf %lf\n",x,y);}
};
typedef Complex C;
const int N=4e6+10;
int r[N],l;int n,m;
C a[N],b[N];
inline void FFT(C *P,int f)
{
for(register int i=0;i<n;++i) if(i<r[i]) swap(P[i],P[r[i]]);
for(register int i=1;i<n;i<<=1){//枚举每一层蝴蝶操作的跨度
C W(cos(PI/i),f*sin(PI/i));//同层所需要的单位根是相同的(是蝴蝶操作直径次单位根)
for(register int p=i<<1,j=0;j<n;j+=p){
//两倍当前蝴蝶操作半径为一组(这里是多项式系数奇偶分组后的组数)
C w(1,0);
for(register int k=0;k<i;++k,w=w*W){
//每组内有蝴蝶操作直径个子DFT,但是蝴蝶操作每次处理两个,只要枚举到半径的长度
C X=P[j+k],Y=w*P[j+k+i];
P[j+k]=X+Y,P[j+k+i]=X-Y;
}
}
}
}
int main()
{
n=read();m=read();
for(register int i=0;i<=n;++i) a[i].init();
for(register int i=0;i<=m;++i) b[i].init();
m+=n;l=0;
for(n=1;n<=m;n<<=1) ++l;//补成2的幂
for(register int i=1;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);//Rader
FFT(a,1); FFT(b,1);
for(register int i=0;i<n;++i) a[i]=a[i]*b[i];
FFT(a,-1);
for(register int i=0;i<=m;++i) printf("%d ",(int)(a[i].x/n+0.5));//IDFT最后除n
return 0;
}
注意事项:
1.由于FFT有大量实数运算,不仅常数大,精度也会有很大问题,转化为整数时要四舍五入。
2.注意蝴蝶操作的式子不要记错, Y 要乘上单位复数根!
2.NTT(快速数论变换)
1.NTT模数NTT
在模意义下有如下结论:
其中\(g\)是\(p\)的原根,然后就直接把这个代换进FFT模板就可以了,加上模意义下的基本操作
常用NTT模数以及其原根:
\(998244353=119*2^{23}+1\;\ \ ------3\)
\(1004535809=479*2^{21}+1\ ------3\)
注意,虽然
\(1000000007\)的原根是\(5\)
但是我们必须要保证 NTT 的时候数组的长度 \(N\) 始终是 \((P-1)\) 的因数,所以 \(10^9+7\) 并不能用来直接做 \(NTT\) 模数\(!!\)
主体代码:
inline int fpow(int x,int k){int ret=1;for(;k;k>>=1,x=1ll*x*x%mod) if(k&1) ret=1ll*ret*x%mod;return ret;}
inline void Mod(int&x){if(x>=mod) x-=mod;else if(x<0) x+=mod;return;}
const int mod=998244353;
const int phi=998244352;
const int G=3;
inline void NTT(int *A,int n,int f)
{
for(int i=1;i<n;++i) if(rader[i]>i) swap(A[i],A[rader[i]]);
for(int i=1;i<n;i<<=1){
int W=fpow(G,phi/(i<<1));
if(f==-1) W=fpow(W,mod-2);
for(int p=i<<1,j=0;j<n;j+=p){
int w=1;
for(int k=0;k<i;++k,w=1ll*w*W%mod){
int X=A[j|k],Y=1ll*w*A[j|k|i]%mod;
Mod(A[j|k]=X+Y);Mod(A[j|k|i]=X-Y);
}
}
}
if(f==-1) for(int i=0;i<n;++i) A[i]=1ll*A[i]*inv[n]%mod;
return;
}
注意事项
1.如果确定计算结果不会超过你的 NTT 模数,那么可以用 NTT 来代替 FFT 进行运算
2.任意模数NTT
如果真的就是让你用 \(10^9+7\) 等模数来\(NTT\),怎么办
那么就要用 \(MTT\)(任意模数\(NTT\))了
这里用的是拆系数的方法
假设模数为 \(P\),两个多项式分别为 \(A(x),B(x)\)
我们取 \(M=\sqrt P\),并将多项式的系数\(K\)分为 \(K/M,K\mod M\) 两部分,那么可以把每一个多项式拆分为两个系数都不超过\(M\)的多项式,分别\(DFT和IDFT\),最后还原出系数并对原模数取模即可,这样的话最后的系数最大可能达到 \(nP\) 的级别,要确保\(long\ long\)能存下
主体代码:
inline void MTT(int n,int *A,int *B,int *C){
for(int i=0;i<n;++i) {A1[i]=Co(A[i]/P,0);A2[i]=Co(A[i]%P,0);B1[i]=Co(B[i]/P,0);B2[i]=Co(B[i]%P,0);}
for(int i=1;i<n;i<<=1)
for(int k=0;k<i;++k)
w[n/i*k]=Co(cos(1.00*k*PI/i),sin(1.00*k*PI/i));
FFT(n,A1,1);FFT(n,A2,1);FFT(n,B1,1);FFT(n,B2,1);
for(int i=0;i<n;++i){
Co X=A1[i]*B1[i],Y=A2[i]*B1[i]+A1[i]*B2[i],Z=A2[i]*B2[i];
A1[i]=X,A2[i]=Y,B1[i]=Z;
}
FFT(n,A1,-1);FFT(n,A2,-1);FFT(n,B1,-1);
int P2=1ll*P*P%mod;
for(int i=0;i<n;++i){
ll X=(ll)(A1[i].x+0.5)%mod,Y=(ll)(A2[i].x+0.5)%mod,Z=(ll)(B1[i].x+0.5)%mod;
C[i]=0;C[i]=X*P2%mod;upd(C[i],Y*P%mod);upd(C[i],Z);
}
return;
}
3.多项式的一些基本操作
1.多项式求逆
求\(G(X)\)使$$F(x)G(x)\equiv 1\ (mod\ x^n)$$
其含义是无论当\(x\)取何值时,该式子总成立,那么如果该多项式有逆,那么他的逆是唯一的
采用倍增法,假设我们已经知道了 \(F(x)H(x)\equiv 1\ (mod\ x^{\lceil \frac{n}{2} \rceil})\)
那么因为显然有:
上下相减并除掉 \(F(x)\)
记\(D(x)=G(X)-H(x)\)那么必须满足\(D(x)\)的\(0\sim {\lceil \frac{n}{2}\rceil-1}\)次方项的系数均为 0
把\(D(x)\)平方,那么显然小于 \(n\) 次方项的系数均为 0, 因为其中必定含有一个次数小于 \(\lceil \frac{n}{2} \rceil-1\)的项,那么显然会有
之后回代并在两边乘上\(F(x)\),变形可得:
所以不断递归就行了,如果模数不是 NTT 模数,那么需要 MTT
边界就是当 n=1 的时候直接求逆元
主体代码:
void Poly_inv(int n){
if(n==1){return (void)(B[0]=fpow(A[0],mod-2));}
Poly_inv((n+1)>>1);
int m=n<<1;
int len,up;for(len=1,up=-1;len<=m;len<<=1,++up);
rader[0]=0;for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
for(int i=0;i<n;++i) C[i]=A[i];for(int i=n;i<len;++i) C[i]=0;
NTT(len,C,1);NTT(len,B,1);
for(int i=0;i<len;++i) B[i]=(2ll*B[i]-1ll*C[i]*B[i]%mod*B[i]%mod+mod)%mod;
NTT(len,B,-1);
for(int i=n;i<len;++i) B[i]=0;
return;
}
2.多项式除法&取模
求解:
给定次数分别为\(n,m\)的多项式\(F(x),G(x)\),求出次数为 \(n-m\) 的多项式\(Q(x)\)和次数小于\(m\)的多项式\(R(x)\),使得:
可以发现\(Q(x)\)即为多项式除法的商,而\(R(x)\)为余数
首先你可能会想,除法不就是乘上多项式的逆吗,所以直接求逆即可。
但是你会发现好像这个东西没办法下手，因为有次数限制和余数,并且这个不是多项式求逆的式子的形式。
所以我们必须要找到一种方法使得能够把式子转化为能和多项式求逆搭上边的。
那么做法就是,定义\(F_R(x)=x^nF(\dfrac{1}{x})\),这个显然就是把 \(F(x)\) 的系数翻转了一下。
那么我们把原式中的\(x\)用 \(\dfrac{1}{x}\) 来替换,并且两边乘上一个\(x^n\)这样可以得到:
于是:
这下终于能用同余式来表示原式子了,并且未知多项式还能够少掉一个:
那么:
这下只要求出多项式\(G_R(x)\)的逆然后乘法就得到了\(Q_R(x)\),翻转回去就可以进一步通过多项式乘法得到\(R(x)\)
主体代码:
int D[N];
void Poly_divide(int *A,int *B,int n,int m){
reverse(A,A+1+n);reverse(B,B+1+m);
Poly_inv(B,C,D,n-m+1);
Set(D,0);for(int i=0;i<=n-m;++i) D[i]=A[i];
int r=(n-m)<<1,len,up;
for(len=1,up=-1;len<=r;len<<=1,++up);
for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
NTT(C,len,1);NTT(D,len,1);
for(int i=0;i<len;++i) C[i]=1ll*C[i]*D[i]%mod;
NTT(C,len,-1);for(int i=n-m+1;i<len;++i) C[i]=0;
reverse(A,A+1+n);reverse(B,B+1+m);reverse(C,C+1+n-m);
for(int i=0;i<=n-m;++i) printf("%d ",C[i]);
puts("");
for(len=1,up=-1;len<=n;len<<=1,++up);
for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
NTT(C,len,1);NTT(B,len,1);
for(int i=0;i<len;++i) C[i]=1ll*C[i]*B[i]%mod;
NTT(C,len,-1);
for(int i=0;i<m;++i) Mod(C[i]=A[i]-C[i]),printf("%d ",C[i]);
puts("");
}
3.多项式求 ln
给定 \(A(x)\),求\(B(x)\)使得
求对数函数直接两边求导就可以了:
所以:
只需要会求导,求不定积分,多项式求逆就行了
主体代码:
inline void Direv(int *A,int n){for(int i=1;i<n;++i) A[i-1]=1ll*A[i]*i%mod;A[n-1]=0;return;}
inline void Inter(int *A,int n){for(int i=n-2;~i;--i) A[i+1]=1ll*A[i]*inv[i+1]%mod;A[0]=0;return;}
inline void Poly_ln(int *A,int *B,int n){
Poly_inv(A,B,n);Set(D,0);//D是辅助数组
for(int i=0;i<n;++i) D[i]=A[i];
Direv(D,n);int len,up;int m=n<<1;
for(len=1,up=-1;len<=m;len<<=1,++up);
for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
NTT(D,len,1);NTT(B,len,1);
for(int i=0;i<len;++i) B[i]=1ll*B[i]*D[i]%mod;
NTT(B,len,-1);for(int i=n;i<len;++i) B[i]=0;
Inter(B,n);
return;
}
4.多项式 exp
已知\(A(x)\),求\(B(x)\)使得:
首先你需要学会牛顿迭代,这是个不断求导并迭代来求函数零点的方法。
就是要求:$$G(F(x))\equiv 0\ (mod\ x^n)$$
假设我们已经知道了:
把\(G(F(x))\),在\(F_0(x)\)处泰勒展开:
我们易知\(F(x)-F_0(x)\)的低于\(x^{\lceil\frac{n}{2}\rceil}\)的项的系数为0,而我们的运算在模\(x^n\)意义下进行,那么大于等于2次的展开项都可以直接视作0,因为\((x^{\lceil\frac{n}{2}\rceil})^2\geq x^n\)
于是:
而\(G(F(x))\equiv 0\ (mod\ x^n)\),故可得:
回到多项式 \(exp\),两边取\(ln\)后可得:
直接把 \(F(x)=B(x),G(F(x))=ln\ B(x)-A(x)\)代入牛顿迭代公式:
可得:
倍增求出这个东西,需要 多项式求\(ln\) 的所有内容(注意求ln的数组要清空)
主体代码:
void Poly_exp(int *A,int *B,int n){
if(n==1) return (void)(B[0]=1);
Poly_exp(A,B,(n+1)>>1);Poly_ln(B,C,n);
int m=n<<1,len,up;
for(len=1,up=-1;len<=m;len<<=1,++up);
for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
for(int i=0;i<n;++i) Mod(C[i]=A[i]-C[i]);for(int i=n;i<len;++i) C[i]=0;
Mod(++C[0]);
NTT(C,len,1);NTT(B,len,1);
for(int i=0;i<len;++i) B[i]=1ll*B[i]*C[i]%mod;
NTT(B,len,-1);
for(int i=n;i<len;++i) B[i]=0;
return;
}
\(P.S.\): 一个常错点,当运算在 \(mod\; x^n\) 下进行时,一定不能合并点值相乘的步骤,必须一次次来,并把多余的项赋成 0
5.多项式多点求值
求解多项式 \(F(x)\) 在 \(x_1,x_2,x_3\dots x_m\) 处的点值
预备知识:
余式定理:
一个多项式函数 \(F(x)\) 在 \(x_0\) 处的点值是 \(F(x)\mod (x-x_0)\)
根据多项式取模的原理，因为除式是一个一次时，那么结果只有常数项。
理解可以从我们奥数中学习过的因式定理出发(尽管因式定理是余式定理的一种特殊情况)。
因式定理指出，一个多项式\(F(x)\) 含有因式 \((x-x_0)\) 当且仅当 \(F(x_0)=0\)，这与我们因式分解法解方程的根是一样的原理。
那么我们可以知道 \(F(x)\mod (x-x_0)=0\)，令\(G(x)=F(x)+y_0\)，由于只是增加了一个常数项，我们可以得到:
\(G(x_0)=F(x_0)+y_0=y_0\)
\(G(x)\mod (x-x_0)=(F(x)+y_0)\mod (x-x_0)=F(x)\mod (x-x_0)+y_0=y_0\)
故余式定理成立。
接下来就是多点求值的方法了。
我们知道 \(F(x) \mod (x-x_1)=F(x) \mod (x-x_1)(x-x_2) \mod (x-x_1)\) ，因为之后取模的多项式是上次取模的多项式的一个因式，显然结果不受上次取模的影响。
这样我们可以发现多项式的次数可以降低而不影响最终的点值。
那么我们定义 \(P_{l,r}(x)=\prod_{i=l}^r (x-x_i)\)
然后考虑分治求解点值。按照线段树那样dfs。
对于当前多项式 \(F_{l,r}(x)\)。
往左走就令 \(F_{l,mid}(x)=F_{l,r}(x)\mod P_{l,mid}(x)\)
往右走就是 \(F_{mid+1,r}(x)=F_{l,r}(x)\mod P_{mid+1,r}(x)\)
这样每次往下走一步多项式的长度都减半,次数之和就是 \(O(nlogn)\) ，总体复杂度就是\(O(nlog^2n)\)
当然中间的常数巨大...
这样我们先用 分治+多项式乘法 来求出每一个线段树上节点的 \(P(x)\)， 然后递归分治进取求解点值，到了最底下的点值就是其 \(F(x)\) 的常数项了。
为了卡常可以当当前区间长度较小的话改用暴力直接计算点值。
主体code:
#define ls (u<<1)
#define rs (u<<1|1)
int POOL[MAXM],cnt=0;
int *P[MAXN],*F[MAXN],X[N];int ans[N];
inline int* Neospace(int len){int*ret=&POOL[cnt];cnt+=len;return ret;}
inline int Calc(int*F,int n,const int x){int X=1,ret=0;for(int i=0;i<n;++i) Inc(ret,(ll)F[i]*X%mod),X=(ll)X*x%mod;return ret;}
inline void Divide(int u,int l,int r){P[u]=NULL;
if(l==r) {P[u]=Neospace(2);P[u][0]=mod-X[l];P[u][1]=1;return;}
static int L[MAXN],R[MAXN];int mid=(l+r)>>1;int LS=ls,RS=rs;
Divide(LS,l,mid);Divide(RS,mid+1,r);
int n=r-l+2,nl=mid-l+2,nr=r-mid+1;
for(int i=0;i<nl;++i) L[i]=P[LS][i];
for(int i=0;i<nr;++i) R[i]=P[RS][i];
Mul(L,R,L,nl,nr);P[u]=Neospace(n);
for(int i=0;i<n;++i) P[u][i]=L[i];
return;
}
inline void Poly_Evaluation(int u,int l,int r){//注意最开始的多项式也要取一次模
if(r-l+1<=200) {for(int i=l;i<=r;++i) ans[i]=Calc(F[u],r-l+1,X[i]);return;}// 长度小的暴力计算
int mid=(l+r)>>1,LS=ls,RS=rs;
int n=r-l+1,nl=mid-l+1,nr=r-mid;
static int R[MAXN];
F[LS]=Neospace(nl),F[RS]=Neospace(nr);
Poly_Mod(F[u],P[LS],__,R,n-1,nl);
for(int i=0;i<nl;++i) F[LS][i]=R[i];
Poly_Mod(F[u],P[RS],__,R,n-1,nr);
for(int i=0;i<nr;++i) F[RS][i]=R[i];
Poly_Evaluation(LS,l,mid);
Poly_Evaluation(RS,mid+1,r);
return;
}
void Solve(){
int n,m;init(n),init(m);++n,++m;
static int A[MAXN];
for(int i=0;i<n;++i) init(A[i]);
Input_Array(X,1,m);Divide(1,1,m-1);
F[1]=Neospace(m-1);
Poly_Mod(A,P[1],__,F[1],n-1,m-1);
Poly_Evaluation(1,1,m-1);
for(int i=1;i<m;++i) printf("%d\n",ans[i]);
}
6.多项式快速插值
问题:
给定 \(n\) 个点 \((x_i,y_i)\) ，求出这些点确定的一个 \(n-1\) 次多项式。
前置知识:
拉格朗日插值与多项式多点求值。
考虑拉格朗日插值:
复杂度是 \(O(n^2)\)
考虑快速计算后面这一部分。
我们定义 \(\pi(x)=\prod_{i=1}^n(x-x_i)\)，把式子变形一下就是:
显然第一个东西下面的分式是一个 \(\frac{0}{0}\) 的形式，因为 \(x\rightarrow x_i时 ,(x-x_i)\rightarrow 0 , \pi(x)\rightarrow 0\)
根据洛必达法则求函数极限: $$\lim_{x\rightarrow x_i} \frac{\pi(x)}{(x-x_i)}=\lim_{x\rightarrow x_i} \frac{\pi'(x)}{(x-x_i)'}$$
因为 \((x-x_i)'=1,\)所以结果就是 \(\pi'(x_i)\)
那么原式就变成了:
前面的那部分就可以直接多点求值先全部算出来了。
然后就是处理后面的部分了，这个东西也可以类似分治的求。
设 \(F_{l,r}(x)\) 表示 \([l,r]\) 内的点插值出来的多项式，那么这个东西可以由左右两边合并得到，具体来说是:
推导略。
主体code:
inline void Poly_Interpolate(int u,int l,int r){// 多项式快速插值
if(l==r) {G[u]=Neospace(1);G[u][0]=val[l];return;}
int mid=(l+r)>>1,LS=ls,RS=rs,n=r-l+1,nl=mid-l+1,nr=r-mid;
Poly_Interpolate(LS,l,mid);
Poly_Interpolate(RS,mid+1,r);
static int L[MAXN],R[MAXN];
G[u]=Neospace(n);
Mul(G[LS],P[RS],L,nl,nr+1);
Mul(G[RS],P[LS],R,nr,nl+1);
for(int i=0;i<n;++i) G[u][i]=Sum(L[i],R[i]);
return;
}
inline void Solve_Interpolation(int*_X,int*Y,int*ans,int n){
cnt=0;for(int i=1;i<=n;++i) X[i]=_X[i];
Divide(1,1,n);static int Val[N];val=Val;int hcnt=cnt;
F[1]=Neospace(n);for(int i=0;i<=n;++i) F[1][i]=P[1][i];
Direv(F[1],n+1);Poly_Evaluate(1,1,n);
Clear(POOL,hcnt,cnt);cnt=hcnt;
for(int i=1;i<=n;++i) Val[i]=(ll)Y[i]*fpow(Val[i],mod-2)%mod;
Poly_Interpolate(1,1,n);for(int i=0;i<n;++i) ans[i]=G[1][i];
return;
}
7.多项式三角函数
求 \(F(x)\) 使得 \(F(x)=sin(A(x))\) 或 \(F(x)=cos(A(x))\)
(并不知道这玩意有什么鸟用)
学过复数就知道: \(e^{i\theta}=cos\theta+isin\theta\)
然后就代入就行了:
三角函数诱导公式...
所以:
\(i\) 是个虚数，我们把它放在模意义下就行了，因为 \(i^2=-1\) 所以 \(i\) 是 -1 的二次剩余，求出来是 \(86583718\)
然后 \(exp\) 就没了。
代码咕咕咕。
8.多项式反三角函数
(这个就更加不知道有什么鸟用了)
求 \(F(x)\) 使得 \(F(x)=arcsin(A(x))\) 或 \(F(x)=arctan(A(x))\)
求导和积分随便化一下就弄出来了，你需要一个反三角函数求导公式。
直接扔出最后的式子了:
板子一顿上就行了。
代码依然咕咕咕。
4.FWT(快速沃尔什变换)
1.概述
我们知道假设有两个多项式 \(A,B\),他们的多项式卷积可以写成以下形式(中间打个乘号的话表示的是按位相乘,没有则为卷积)
但假设现在有一个问题是给你两堆集合,求出在两堆中各选出一个集合使得他们的并集是某个集合的方案数
这个问题其实就是要求出以下形式的卷积:
\(x^i\)的系数即为集合为\(i\)时的答案
所以\(FWT\)就是把普通多项式乘法卷积中的\(j+k=i\)变成了\(j\oplus k=i\),其中 \(\oplus\)表示一个位运算的算法,也就是:
现在就要找到一种方法能够快速求出以上卷积。
仿照\(FFT\)的方法,我们希望找到一种变换\(FWT\),能够使得:
以上操作类似于\(FFT\)点值相乘,复杂度为O(n),这样的化我们通过找到一种\(FWT\)的逆变换\(IFWT\)就可以还原出结果多项式了
复杂度就在于求出 \(FWT\)
一下开始讲三种位运算的卷积,注意之后的下标可以将其理解为一个集合
2.或卷积
由于\(j|k=i\) ,显然的 \(i\subseteq j\)且\(k\subseteq i\)
通过分析(xia cai),我们感觉可以构造这样的\(FWT\)形式(\(FWT(A)_i\)表示多项式\(A\)经过变换后的第\(i\)项):
这样对于 \(C_i=\sum_{j|k=i}A_j*B_k\),可得:
我们可以发现对于任意一对 \(p,q\),当且仅当\(j\)枚举到 \(p|q\)时会被计算一次贡献,而我们显然有\(p\subseteq j\),\(q\subseteq j\),因此式子可以改写为:
显然可以用分配率提出来\(A_p\)
万事大吉,这样就找到一种优秀的\(FWT\)形式了
剩下的就是看怎么快速求出这个\(FWT\)和\(IFWT\)了
如果你去看很多其他人的博客,他们会通过分治思想,多项式拼接什么的把非巨佬以外的人绕的云里雾里,然后再甩给你一个代码。
但是你没发现或运算下的\(FWT\)就是个子集的前缀和吗,所以直接用高维前缀和求就行了。
对,就是这段代码:
for(int i=0;i<m;++i)
for(int j=0;j<n;++j){if(j>>i&1) Inc(A[j],A[j^bin[i]]);}
\(FWT\)的精髓在于构造优秀的变换形式,求解变换就是高维前缀和。
至于为什么可以写成和\(FFT\)代码高度相似,是因为那个是直接枚举已经确定下来的\(1\)左右两边的0/1情况,理论上能少一半的常数,但是会多一重循环(说不定常数更大了),也就是长成这样:
for(int i=1;i<n;i<<=1)
for(int p=i<<1,j=0;j<n;j+=p)
for(int k=0;k<i;++k) f? Inc(A[j|k|i],A[j|k]):Dec(A[j|k|i],A[j|k]);//这是模意义下加减法
后面那个减法就是 \(IFWT\)了,就是把高位前缀和的操作反着来一边就行了,把子集的和减去就得到原本自己的值
这就是或卷积了。
3.与卷积
与和或运算非常相似,所以找\(FWT\)的思想一样,直接把子集前缀和换成超集前缀和就行了:
后面的推导是一样的,求解\(FWT\)也一样高位前缀和解决,代码如下:
for(int i=1;i<n;i<<=1)
for(int p=i<<1,j=0;j<n;j+=p)
for(int k=0;k<i;++k) f? Inc(A[j|k],A[j|k|i]):Dec(A[j|k],A[j|k|i]);
仅仅只是把操作的对应数组下标前后交换了一下,其他一模一样
4.异或卷积
用\(\oplus\)表示异或
异或的话要自己瞎找就有点难了只能直接放结论了:
不过这个代码写起来就真的可以和\(FFT\)一模一样了,代码:
for(int i=1;i<n;i<<=1)
for(int p=i<<1,j=0;j<n;j+=p)
for(int k=0;k<i;++k){
int X=A[j|k],Y=A[j|k|i];
Mod(A[j|k]=X+Y);Mod(A[j|k|i]=X-Y);
if(!f) A[j|k]=1ll*A[j|k]*inv%mod,A[j|k|i]=1ll*A[j|k|i]*inv%mod;//inv是2的逆元,不是模意义下就直接除以2
}
5.小结
\(FWT\)的主要难点(其实就由这两部分构成),一是构造出合适的\(FWT\)变换形式,而是能够快速求出\(FWT\)及\(IFWT\)
我的总体的\(FWT\)代码:
const int mod=998244353;
const int inv=499122177;
inline void Inc(int &x,int y){x+=y;if(x>=mod) x-=mod;}
inline void Dec(int &x,int y){x-=y;if(x < 0) x+=mod;}
inline void Mod(int &x){if(x>=mod) x-=mod;else if(x<0) x+=mod;return;}
inline void FWT(int *A,int n,int f,int opt){
if(opt==1) {
for(int i=1;i<n;i<<=1) for(int p=i<<1,j=0;j<n;j+=p) for(int k=0;k<i;++k) f? Inc(A[j|k|i],A[j|k]):Dec(A[j|k|i],A[j|k]);
}
else if(opt==2){
for(int i=1;i<n;i<<=1) for(int p=i<<1,j=0;j<n;j+=p) for(int k=0;k<i;++k) f? Inc(A[j|k],A[j|k|i]):Dec(A[j|k],A[j|k|i]);
}
else if(opt==3){
for(int i=1;i<n;i<<=1)
for(int p=i<<1,j=0;j<n;j+=p)
for(int k=0;k<i;++k){
int X=A[j|k],Y=A[j|k|i];
Mod(A[j|k]=X+Y);Mod(A[j|k|i]=X-Y);
if(!f) A[j|k]=1ll*A[j|k]*inv%mod,A[j|k|i]=1ll*A[j|k|i]*inv%mod;
}
}
return;
}
