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;
}