好的繼續一點一點啃吧
Todo List
- “鴿了”的實現
- 多項式復合
- 多項式復合逆
- 常系數非齊次線性遞推
一.數學前置
1.1 二項式反演
等價於
Proof
證明：參考
顯然只要證明前推后即可。
考慮二式的右邊。
得證。
1.2 拉格朗日反演
則
Proof
證明：
設
所以
而當\(i\neq n\)時，顯然對左式沒有貢獻，可得
注意到
所以
即
證畢。
1.3 單位根反演
設 \(\omega\) 表示 \(n\) 次本原單位根。
我們有
Proof
證明：
當 \(n|k\) 時，有 \(\omega^{ik}=1\)，顯然成立！
當 \(n\nmid k\) 時，有
得證！
1.4 斯特林數
先給出斯特林數的定義：
第一類斯特林數： \(\begin{bmatrix}n\\k\end{bmatrix}\) 表示將 \(n\) 個兩兩不同的元素，划分為 \(k\) 個非空圓排列的方案數。
第二類斯特林數： \(\begin{Bmatrix}n\\k\end{Bmatrix}\) 表示將 \(n\) 個兩兩不同的元素，划分為 \(k\) 個非空子集的方案數。
有如下遞推式：
有用的式子：
斯特林反演：
等價於
通項公式：
這個可以用二項式反演或者直接組合意義容斥來證
二.進階多項式算法
2.1 Bluestein's Algorithm
orz myy！
注：原論文似乎有點問題，如果發現這里也有問題歡迎指出。
\(\text{Bluestein's Algorithm}\) 支持任意長度 \(\text{CZT}\)，即給出多項式 \(B(x)\)，計算
考慮如何將其化為普通卷積的形式。
下面我們用一種膜法，注意到這么一個恆等式：
正確性顯然。
通常稱為 \(\text{Bluestein}\) 的第一個拆乘積的方法。
代入上面的式子：
發現是個卷積，於是直接 NTT 就好了……
嗎？
我們發現指數不是整數，這意味着只有在 \(c\) 有二次剩余的時候才可以，但是通常不存在。
有沒有不依賴二次剩余的方法呢？
我們發現了另一個式子：
通常稱為 \(\text{Bluestein}\) 的第二個拆乘積的方法。
於是
這就是一個卷積的形式了，直接 NTT 即可。
復雜度 \(O(n\log n)\)。
Code
vec CZT(vec f,int len,int c,int m){
vec A(len),B(len+m),g(m);
for(int i=0;i<len;i++)A[i]=1LL*qpow(c,P-1-F(i))*f[i]%P;
for(int i=0;i<len+m;i++)B[i]=qpow(c,F(i));
reverse(A.begin(),A.end());
A=mul(A,B,len*2+m);
for(int i=0;i<m;i++)g[i]=1LL*qpow(c,P-1-F(i))*A[i+len-1]%P;
return g;
}
2.2 多項式多點求值
多項式多點求值，即求多項式 \(A(x)\) 在點 \(x_1,x_2,\cdots,x_n\) 上分別取到的值。
Sol 1
常用做法。
構造函數
設 \(A(x)=P_0(x)D(x)+R(x)\)，且 \(\deg R<\deg P_0=\lfloor\frac{n}{2}\rfloor\)。
那么對於 \(\forall x\in \{x_1,x_2,\cdots,x_{\lfloor\frac{n}{2}\rfloor}\}\)，都有 \(A(x)=R(x)\)，於是我們就將一個 \(\deg =n\) ，有 \(n\) 個點的問題，轉化為一個 \(\deg <\lfloor\frac{n}{2}\rfloor\) ，有 \(\lfloor\frac{n}{2}\rfloor\) 個點的問題。
另一半點值與上面同理。
因此分治 NTT 即可，注意計算 \(P_0(x)\) 和 \(P_1(x)\) 同樣可以用一次分治 NTT 預處理。
設計算的復雜度為 \(T(x)\)，則 \(T(x)=2T(\frac{x}{2})+O(n\log n)\)，此處 \(O(n\log n)\) 的部分是算 \(D(x),R(x)\) 的復雜度，容易發現預計算復雜度也是 \(T(x)\)。
由主定理可得總復雜度為 \(O(n\log^2 n)\)。
實際應用時可以通過用暴力秦九韶進行剪枝來優化常數。
Code
inline vec mod(vec f,vec g,int n,int m){
f.resize(n),g.resize(m);
vec ff(n),gg(n),iv(n),t(n);
for(int i=0,t=n-1;i<n;i++,t--)ff[i]=f[t];
for(int i=0,t=m-1;i<m;i++,t--)gg[i]=g[t];
int len=n-m+1;
for(int i=len;i<n;i++)ff[i]=gg[i]=0;
getinv(gg,iv,len);
ff=mul(ff,iv,len<<1);
vec q(n);
for(int i=0,t=len-1;i<len;i++)q[i]=ff[t--];
for(int i=len;i<n;i++)q[i]=0;
t=mul(q,g,n<<1);
vec r(m-1);
for(int i=0;i<m-1;i++)r[i]=(f[i]-t[i]+P)%P;
return r;
}
void build(int k,int l,int r,const vec&a){
if(l==r){
Len[k]=1;
p[k].resize(2);
p[k][0]=P-a[l];
p[k][1]=1;
return;
}
int mid=l+r>>1;
build(k<<1,l,mid,a);
build(k<<1|1,mid+1,r,a);
Len[k]=Len[k<<1]+Len[k<<1|1];
p[k]=mul(p[k<<1],p[k<<1|1],Len[k]+1);
p[k].resize(Len[k]+1);
}
void solve(int k,int l,int r,const vec&a,const vec&A,vec&ans){
if(Len[k]<=500){
int m=Len[k]-1;
for(int i=l;i<=r;i++)
for(int j=m;j>=0;j--)
ans[i]=(1LL*ans[i]*a[i]+A[j])%P;
return;
}
if(l==r){ans[l]=A[0];return;}
int mid=l+r>>1;
vec R;
R=mod(A,p[k<<1],Len[k],Len[k<<1]+1);
solve(k<<1,l,mid,a,R,ans);
R=mod(A,p[k<<1|1],Len[k],Len[k<<1|1]+1);
solve(k<<1|1,mid+1,r,a,R,ans);
}
vec evaluation(vec f,vec a,int n,int m){
build(1,1,m,a);
if(n>m){
f=mod(f,p[1],n,Len[1]+1);
}
vec ans(m+1);
solve(1,1,m,a,f,ans);
return ans;
}
Sol 2
轉置做法。
先簡單解釋一下 轉置原理 的核心思想：
對矩陣 \(\mathbf M\) 和向量 \(\mathbf v\)，為了快速計算 \(\mathbf {Mv}\)，我們可以先考察怎么計算 \(\mathbf M^\mathsf T\mathbf v\)，也就是矩陣的轉置。設 \(\mathbf M\) 可以分解為 \(\mathbf M=\mathbf E_1 \mathbf E_2\cdots \mathbf E_k\)，那么 \(\mathbf M^\mathsf T=\mathbf E_1^\mathsf T \mathbf E_2^\mathsf T\cdots \mathbf E_k^\mathsf T\)。在這些初等矩陣中，它們的功能可以分為兩類：
- \(\mathbf {Ev}\) 的效果是讓向量 \(\mathbf v\) 的第 \(i\) 位乘上 \(c\)，那么其轉置效果不變。
- \(\mathbf{Ev}\) 的效果是讓向量 \(\mathbf v\) 的第 \(i\) 位乘上 \(c\) 然后加到第 \(j\) 位，那么它的轉置的效果就是讓向量的第 \(j\) 位乘上 \(c\) 然后加到第 \(i\) 位。
接下來我們考慮多點求值的過程，其本質是 VanderMonde 矩陣
對應的線性變換 \(\mathbf{V}(x_0,x_1,\dots,x_{n-1})\mathbf v\)，其中 \(\mathbf v\) 是系數向量 \(\begin{bmatrix}v_0\\v_1\\\vdots\\v_{n-1}\end{bmatrix}\)。
根據轉置原理，我們考慮其轉置 \(\mathbf{V}^\mathsf T(x_0,x_1,\dots,x_{n-1})\mathbf v\)：
考慮這個轉置所求的是什么。
第 \(k\) 項是
而為了求這個，我們可以這么做：
- 分治，維護左右兩部分的 \((u_1,L)\) 和 \((u_2,R)\)，也就是答案的分子和分母的部分，然后合並為 \((u_1R+u_2L,LR=P)\)。
- 答案就是 \(\frac{u}{P}=uP^{-1}\)，這里 \(P\) 指的是分母，即 \(\prod_i(1-x_ix)\)。
我們只要算出其轉置，就可以得到新的多點求值做法。在此之前，我們還要考慮一下多項式乘法 \(\mathbf{MUL}(A)\mathbf v\) 的轉置 \(\mathbf {MUL}(A)^\mathsf T\mathbf v\)。
對多項式 \(A(x)=\sum_{i}a_ix^i\)，在多項式乘法中會把每個 \(v_j\) 乘上 \(a_i\) 然后貢獻給 \(i+j\)，因此根據轉置原理，轉置的多項式乘法就是把每個 \(v_{i+j}\) 乘上 \(a_i\) 然后貢獻給 \(i\)。可以看到，這就是一個標准的減法卷積。此外，如果原來的變換沒有溢出部分，那么轉置后的溢出部分我們也要扔掉。
於是，我們就可以得到新的多點求值做法：
- 整體分治（其實就是線段樹的 build），算出每個線段樹節點對應的區間里 \((1-x_ix)\) 的乘積。
- 對要求的 \(m-1\) 次多項式，我們對它計算 \(\mathbf{MUL}(P^{-1}\bmod x^m)^\mathsf T\) 並保留前 \(n\) 位。
- 從線段樹自頂向下遞歸，令下傳的兩個向量分別為 \(\mathbf l=\mathbf {MUL}(R)^\mathsf T\mathbf v\) 和 \(\mathbf r=\mathbf{MUL}(L)^\mathsf T\mathbf v\)。
- 最后得到的葉節點的向量長度均為 \(1\)，對應於該處的點值。
新做法不需要多項式取模，而且常數也有顯著提升。
Code
vec mulT(vec f,vec g){
int n=f.size(),m=g.size()-1;
init(n);
f.resize(lmt),g.resize(lmt);
reverse(g.begin()+1,g.end());
vec ret(lmt);
NTT(f,lmt,1);NTT(g,lmt,1);
for(int i=0;i<lmt;i++)ret[i]=1LL*f[i]*g[i]%P;
NTT(ret,lmt,-1);
ret.resize(n-m+1);
return ret;
}
pair<vec,vec> mul2(vec a,vec b1,vec b2){
int n=a.size(),m1=b1.size()-1,m2=b2.size()-1;
init(n);
a.resize(lmt),b1.resize(lmt),b2.resize(lmt);
reverse(b1.begin()+1,b1.end());
reverse(b2.begin()+1,b2.end());
vec ret1(lmt),ret2(lmt);
NTT(a,lmt,1);NTT(b1,lmt,1);NTT(b2,lmt,1);
for(int i=0;i<lmt;i++)ret1[i]=1LL*a[i]*b1[i]%P,ret2[i]=1ll*a[i]*b2[i]%P;
NTT(ret1,lmt,-1);NTT(ret2,lmt,-1);
ret1.resize(n-m1+1);
ret2.resize(n-m2+1);
return make_pair(ret1,ret2);
}
void build(int n){
lc.assign(n*2,0);
rc.assign(n*2,0);
deg.assign(n*2,0);
p.assign(n*2,{0,0});
q.assign(n*2,{0,0});
pos.resize(n);
int cnt=0;
function<int(int,int)> dfs=[&](int l,int r){
if(l==r){
pos[l]=cnt;
deg[cnt]=1;
return cnt++;
}
int ret=cnt++,mid=l+r>>1;
lc[ret]=dfs(l,mid);
rc[ret]=dfs(mid+1,r);
deg[ret]=deg[lc[ret]]+deg[rc[ret]];
return ret;
};
dfs(0,n-1);
}
vec eval(vec&f,const vec&x){
int n=f.size();
build(n);
rep(i,0,n-1)q[pos[i]]={1,P-x[i]},deg[pos[i]]=1;
Rep(i,n*2-2,0)if(lc[i])q[i]=mul(q[lc[i]],q[rc[i]],deg[i]+1);
f.resize(n<<1);
vec g;
int len=q[0].size()+1;
getinv(q[0],g,len);
g.resize(len);
p[0]=mulT(f,g);
rep(i,0,n*2-1)if(lc[i])tie(p[lc[i]],p[rc[i]])=mul2(p[i],q[rc[i]],q[lc[i]]);
vec ans(n);
rep(i,0,n-1)ans[i]=p[pos[i]][0];
return ans;
}
vec evaluation(vec f,vec x){
int n=f.size(),m=x.size();
f.resize(max(n,m));
x.resize(max(n,m));
vec ret=eval(f,x);
ret.resize(m);
return ret;
}
2.3 多項式快速插值
多項式快速插值，即給定 \(n\) 個點 \((x_i,y_i)\)，你需要求出一個多項式 \(f(x)\)，使得 \(\deg f<n\) 且 \(\forall i,f(x_i)=y_i\)。
首先根據拉格朗日插值公式可得
這時如果你直接模擬那么復雜度是 \(O(n^2)\) 的，我們考慮如何優化。
先考慮系數部分的
不難發現如果設
那么
（最后一步是導數的定義）
那么使用多點求值就可以做到 \(O(n\log^2n)\)。
下設 \(\frac{y_i}{P'(x_i)}=v_i\)。
考慮和多點求值類似的做法，我們還是構造函數
設
那么容易發現
分治 NTT 即可。
設這部分復雜度為 \(T(x)\)，則 \(T(x)=2T(\frac{x}{2})+O(n\log n)\)，此處 \(O(n\log n)\) 的部分是計算 \(f(x)\) 的復雜度。由於 \(P(x)\) 在多點求值時已經算過，所以不需再算。
由主定理可得總復雜度為 \(O(n\log^2 n)\)。
講一些卡常的技巧：注意到在遞推的時候大部分都是先 \(n\) 大小 IDFT 再 \(2n\) 大小 DFT，可以考慮直接維護 DFT 結果，可以顯著優化常數。
Code
鴿了
2.4 普通多項式轉下降冪多項式
給定普通多項式
求下降冪多項式
使得 \(F(x)=G(x)\)。
首先有
所以
所以
多點求值預處理 \(F(t)\) ，然后卷積即可。
復雜度 \(O(n\log ^2n)\)。
Code
鴿了
2.5 下降冪多項式乘法
給定 \(n\) 次下降冪多項式 \(A(x)\) 和 \(m\) 次下降冪多項式 \(B(x)\)，求一個 \(n+m\) 次的下降冪多項式 \(F(x)\)，使 \(F(x)=A(x)B(x)\)。
對於下降冪多項式 \(F(x)\)，考慮它的點值的 EGF
當 \(F(x)=x^\underline c\) 時
設
則
其中
也就是說，\(F^\#\) 就是把 \(F\) 當成普通多項式，乘上 \(e^x\) 的結果。
同理可得其逆運算就是 \(e^{-x}\)。
點值 EGF 點乘即可，\(e^x\) 和 \(e^{-x}\) 可以手算。
復雜度 \(O(n\log n)\)。
Code
void FDT(vec&A,int len,int tp){
init(len<<1);
vec ee(lmt);
if(A.size()<lmt)A.resize(lmt);
if(tp==-1)for(int i=0;i<lmt;i++)A[i]=1LL*A[i]*invfac[i]%P;
for(int i=0;i<len;i++){
if(tp==-1&&i&1)ee[i]=P-invfac[i];
else ee[i]=invfac[i];
}
for(int i=len;i<lmt;i++)ee[i]=A[i]=0;
NTT(A,lmt,1);NTT(ee,lmt,1);
for(int i=0;i<lmt;i++)A[i]=1LL*A[i]*ee[i]%P;
NTT(A,lmt,-1);
if(tp==1)for(int i=0;i<lmt;i++)A[i]=1LL*A[i]*fac[i]%P;
for(int i=0;i<lmt;i++)ee[i]=0;
}
vec mulDown(vec f,vec g,int len){
FDT(f,len,1);FDT(g,len,1);
for(int i=0;i<len;i++)f[i]=1LL*f[i]*g[i]%P;
FDT(f,len,-1);
return f;
}
2.6 下降冪多項式轉普通多項式
給定下降冪多項式
求普通多項式
使得 \(F(x)=G(x)\)。
根據上題結論，下降冪多項式的點值是與 \(e^x\) 的卷積。
那我們只要求出 \(F(x)\) 在 \(0,1,\dots,n-1\) 的點值，然后用多點求值就可以得到原多項式。
復雜度 \(O(n\log^2 n)\)。
2.7 多項式復合
2.8 多項式復合逆
2.9 快速沃爾什變換
給定序列 \(A\) 和 \(B\)，求
其中 \(\oplus\) 是位運算
嘗試采取和 DFT 相類似的手法，考慮另一個數列變換 \(\operatorname{FWT}(A)\)，稱之為沃爾什變換。
為了能做到快速進行乘法，我們需要使得
其中 \(\cdot\) 表示點乘。
下面分與、或、異或三種進行討論。
2.9.1 或
顯然可得
構造
所以
下面考慮如何快速計算和逆運算。
和 FFT 一樣，將 \(a\) 拆成下標最高位為 \(0\) 與下標最高位為 \(1\) 兩個序列，分別記為 \(a_0\) 和 \(a_1\)。
所以可得
這里 \(\operatorname{merge}\) 表示拼接，\(+\) 表示對應位置相加。
分治就好了。
逆運算：根據上面那個式子我們同樣可以得到 \(a=\operatorname{merge}(a_0,a_1-a_0)\)，就可以了。
2.9.2 與
完全同理地，我們定義
推理過程略。
2.9.3 異或
這里我們定義運算 \(\otimes\)，使得 \(x\otimes y=\operatorname{popcount}(x\& y)\bmod 2\)。
因此我們有 \((i\otimes j)\oplus(i\otimes k)=i\otimes(j\oplus k)\)。
那么定義
所以
所以
以上復雜度均為 \(O(n\log n)\)。
2.10 常系數齊次線性遞推
常系數齊次線性遞推，是指給定 \(f[1,2,\dots,k]\) 和 \(a[0,1,\dots,k-1]\)，且對 \(\forall m\ge k\)，都有
求 \(a_n\) 的值。
先考慮一般是如何處理這個問題的。
考慮矩陣
那么
所以
使用矩陣快速冪就可以做到 \(O(k^3\log n)\) 的復雜度。
考慮如何優化這個過程，由於我們只要求 \(a_n\)，所以我們只關心右式的第一項即可。
下記 \(m=n-k+1\)，且
現在假設我們不知道從哪找來一個數組 \(c\)，使得
則
下考慮如何構造 \(c\)。
考慮構造一個
其中 \(Q,G,R\) 是一些多項式，且使得 \(\deg R<\deg G\)。
令 \(\deg G=k\) 且 \(G(A)=0\)，有：
因此我們如果構造出了 \(G\)，就直接求出 \(x^n\bmod G\) 就可以了，使用快速冪即可，復雜度為 \(O(k\log k\log n)\) 的。
因此只要取
就行。
結論：如果遞推系數為 \(\{a\}\)，那么 \(g_{k-i}=-a_i,g_k=1\)。
Proof
先定義特征值：如果 $$ (\lambda I-A)v=0 $$ 那么 $\lambda$ 是 $A$ 的特征值，$v$ 是 $A$ 的特征向量。再定義特征多項式：
是特征多項式。
根據 Cayley-Hamilton 定理可得，\(G(A)=0\)，因此 \(G\) 就是我們要求的。
得證。
復雜度 \(O(k\log k\log n)\)。