好的一點一點啃吧
Todo List
- 任意模數多項式乘法逆
- 多項式開根（加強版）
- 多項式快速冪 (加強版)
目錄
零.預處理逆元
原來其實沒寫這個。。。但是多項式反三角卡常必須要用，那就寫一下。
記
\[\operatorname{inv}(i)=i^{-1}\pmod P \]
有：
\[\operatorname{inv}(i)=-\left[\frac{P}{i}\right]\operatorname{inv}(P\ \bmod i) \]
Proof
顯然
\[P-\left[\frac{P}{i}\right]i\equiv -\left[\frac{P}{i}\right]i\pmod P \]
所以
\[\operatorname{inv}(i)=\frac{1}{i}=-\frac{\left[\frac{P}{i}\right]}{P-\left[\frac{P}{i}\right]i}=-\frac{\left[\frac{P}{i}\right]}{P\ \bmod i}=-\left[\frac{P}{i}\right]\operatorname{inv}(P\ \bmod i) \]
所以可以遞推，復雜度\(O(n)\)。
Code
ll prep(ll n){
ll lmt=1;
inv[0]=inv[1]=1;
while(lmt<n)lmt<<=1;
for(ll i=2;i<lmt;i++)inv[i]=(P-P/i)*inv[P%i]%P;
return lmt;
}
一.模數為NTT模數的部分
1. 多項式乘法（FFT）
orz myy！！！
給定序列 \(a_i,b_i\)，下標從 \(0\) 開始，求
\[c_i=\sum_{k=0}^{i}a_kb_{i-k} \]
顯然 \(|c|=|a|+|b|-1\)。
首先，我們可以找一個 \(n=2^k\)，使得 \(n>c\)，至於為什么要求 \(n=2^k\) 后面再說。
於是：
\[c_r=\sum_{p,q}[(p+q)\bmod n=r]a_pb_q \]
記 \(\omega\) 為 \(n\) 次本原單位根，則根據單位根反演：
\[[n\bmod k=0]=\frac{1}{n}\sum_{i=0}^{n-1}\omega^{ik} \]
所以
\[\begin{aligned} &[(p+q)\bmod n=r] \cr =&[(p+q-r)\bmod n=0] \cr =&\frac{1}{n}\sum_{k=0}^{n-1}\omega^{k(p+q-r)}\cr =&\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-kr}\omega^{kp}\omega^{kq} \end{aligned} \]
所以
\[\begin{aligned} c_r&=\sum_{p,q}[(p+q)\bmod n=r]a_pb_q\\ &=\sum_{p,q}\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-kr}\omega^{kp}\omega^{kq}a_pb_q\\ &=\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-rk}\sum_{p,q}\omega^{pk}a_p\omega^{qk}a_q\\ &=\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-rk}\sum_p\omega^{pk}a_p\sum_q\omega^{qk}a_q \end{aligned} \]
於是我們發現了這兩個式子：
\[g_m=\sum_{k=0}^{n-1}\omega^{mk}f_k \]
\[f_m=\sum_{k=0}^{n-1}\omega^{-mk}g_k \]
根據單位根反演，這兩式等價。
根據上面求 \(c\) 的式子，我們發現只要對於每一種關系式，知道右邊即可快速求得左邊，就可以解決這個問題。
另一方面，我們發現第一個式子等價於
\[g_m=F(\omega^m) \]
其中 \(F\) 為數列 \(f\) 的 OGF。
所以這個操作相當於，給出一個多項式的系數表示，對於所有 \(k\) 求出這個多項式在 \(\omega^k\) 的值，也就是在系數表示和點值表示之間轉化。
我們稱這樣的一個過程為離散傅里葉變換(\(\text{DFT}\))，反過來就是 \(\text{IDFT}\)。
下面我們介紹一種可以在 \(O(n\log n)\) 時間內完成這個操作的算法，叫做快速傅里葉變換 \(\text{FFT}\)。
下面我們只考慮 \(\text{DFT}\)，\(\text{IDFT}\) 同理。
設 \(\omega_n\) 表示 \(n\) 次本原單位根，則有：
\[\omega_{2n}^{2m}=\omega_n^m \]
\[\omega_{2n}^{m}=\omega_{2n}^{m+n} \]
設 \(A_0(x)\) 表示 \(A(x)\) 偶次項的和，\(A_1(x)\) 表示 \(A(x)\) 奇次項的和，有：
\[\begin{aligned} A(\omega_n^m)&=A_0((\omega_n^m)^2)+\omega_n^mA_1((\omega_n^m)^2)\\ &=A_0(\omega_{\frac{n}{2}}^m)+\omega_n^mA_1(\omega_{\frac{n}{2}}^m) \end{aligned} \]
\[\begin{aligned} A(\omega_n^{m+\frac{n}{2}}))&=A_0((\omega_n^m)^2)+\omega_n^{m+\frac{n}{2}}A_1((\omega_n^m)^2)\\ &=A_0(\omega_{\frac{n}{2}}^m)-\omega_n^mA_1(\omega_{\frac{n}{2}}^m) \end{aligned} \]
我們將上述兩式稱為蝴蝶操作。
於是可以考慮迭代/遞歸，考慮到遞歸常數太大，使用迭代法進行處理。
Code
int lmt,rev[N];
struct C{
double x,y;
C(double a=0,double b=0){x=a,y=b;}
};
C operator+(C a,C b){return C(a.x+b.x,a.y+b.y);}
C operator-(C a,C b){return C(a.x-b.x,a.y-b.y);}
C operator*(C a,C b){return C(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline void init(int n){
lmt=1;int t=0;
while(lmt<n)lmt<<=1,t++;
for(int i=1;i<lmt;i++)rev[i]=(rev[i>>1]>>1)|(i&1)<<(t-1);
}
inline void FFT(vector<C>&A,int lmt,int tp){
for(int i=0;i<lmt;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int mid=1;mid<lmt;mid<<=1){
C Wn(cos(Pi/mid),tp*sin(Pi/mid));
for(int R=mid<<1,j=0;j<lmt;j+=R){
C w(1,0);
for(int k=0;k<mid;k++,w=w*Wn){
C x=A[j+k],y=w*A[j+k+mid];
A[j+k]=x+y,A[j+k+mid]=x-y;
}
}
}
}
vec mul(vec f,vec g,int len){
init(len);
f.resize(lmt),g.resize(lmt);
vector<C>a(lmt),b(lmt);
for(int i=0;i<lmt;i++)a[i].x=f[i],b[i].x=g[i],a[i].y=b[i].y=0;
FFT(a,lmt,1);FFT(b,lmt,1);
for(int i=0;i<lmt;i++)a[i]=a[i]*b[i];
FFT(a,lmt,-1);
vec ret(lmt);
for(int i=0;i<lmt;i++)ret[i]=lround(a[i].x/lmt);
return ret;
}
時間復雜度 \(O(n\log n)\)。
2. 多項式乘法（NTT）
模板同上
注意到 FFT 要用復數的原因是因為他有一個單位根，而我們注意到由於有了模數，所以可以取它的一個原根 \(g\)，並且用 \(g\) 來代替單位根。
具體來說，我們有
\[e^{\frac{i\pi}{mid}}=g^{\frac{P-1}{2mid}} \]
不過這玩意有個限制，就是要求 \(2mid|P-1\) ，否則上式顯然不成立。因此，我們一般來說要求 \(V_2(P-1)\ge 20\) ，此處 \(V_2(x)\) 表示 \(x\) 的因子中 \(2\) 的次冪。
通常我們把滿足這個要求的模數稱為 NTT 模數，典型例子有 \(998244353=1+2^{23}\times 7\times 17\)（\(g=3\)），\(1004535809=1+2^{21}\times 479\)（\(g=3\)）等。
Code
inline int qpow(int a,int k){
int ret=1;
while(k){
if(k&1)ret=1LL*ret*a%P;
a=1LL*a*a%P;
k>>=1;
}
return ret%P;
}
inline void init(int n){
lmt=1;int t=0;
while(lmt<n)lmt<<=1,t++;
for(int i=1;i<lmt;i++)rev[i]=(rev[i>>1]>>1)|(i&1)<<(t-1);
}
inline void NTT(vec&A,int lmt,int tp){
if(A.size()<lmt)A.resize(lmt);
for(int i=0;i<lmt;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int m=1,t=0;m<lmt;m<<=1,t++)
for(int j=0,Wn=pw[t+1];j<lmt;j+=m<<1)
for(int k=0,w=1,x,y;k<m;k++,w=1LL*w*Wn%P)
x=A[j+k],y=1LL*w*A[j+k+m]%P,A[j+k]=(x+y)%P,A[j+k+m]=(x-y+P)%P;
if(tp==1)return;
reverse(A.begin()+1,A.begin()+lmt);
for(int i=0,inv=qpow(lmt,P-2);i<=lmt;i++)A[i]=1LL*A[i]*inv%P;
}
vec mul(vec f,vec g,int len){
init(len);
f.resize(lmt),g.resize(lmt);
NTT(f,lmt,1);NTT(g,lmt,1);
for(int i=0;i<lmt;i++)f[i]=1LL*f[i]*g[i]%P;
NTT(f,lmt,-1);
return f;
}
時間復雜度還是 \(O(n\log n)\)。
3. 多項式乘法逆
首先設我們已經找到了一個多項式 \(G_1(x)\)，使得
\[G_1(x)F(x)\equiv 1\pmod {x^{\lceil\frac{n}{2}\rceil}} \]
則將其與原式相減即可得到
\[G(x)-G_1(x)\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}} \]
平方：
\[G(x)^2-2G(x)G_1(x)+G_1(x)^2\equiv 0\pmod {x^n} \]
兩邊同乘 \(F(x)\)，根據 \(F(x)G(x)\equiv 1\pmod {x^n}\) 可得：
\[G(x)-2G_1(x)+G_1^2(x)F(x)\equiv 0\pmod {x^n} \]
移項則有
\[G(x)\equiv 2G_1(x)-G_1^2(x)F(x)\pmod {x^n} \]
over。
時間復雜度 \(T(n)=T(\frac{n}{2})+n\log n\)，根據主定理可知 \(T(n)=n\log n\)。
Code
void getinv(vec&f,vec&g,int len){
if(f.size()<(len<<2))f.resize(len<<2);
if(g.size()<(len<<2))g.resize(len<<2);
if(len==1){g[0]=qpow(f[0],P-2);return;}
getinv(f,g,len+1>>1);
init(len<<1);
vec c(lmt);
copy(f.begin(),f.begin()+len,c.begin());
NTT(c,lmt,1),NTT(g,lmt,1);
for(int i=0;i<lmt;i++)g[i]=1LL*(2LL-1LL*g[i]*c[i]%P+P)%P*g[i]%P;
NTT(g,lmt,-1);
for(int i=len;i<lmt;i++)g[i]=0;
}
4. 多項式除法
占位。
4.5 多項式求導，積分
並不想說什么。
Code
vec getdev(vec f,int len){
vec g(len-1);
for(int i=1;i<len;i++)g[i-1]=1LL*i*f[i]%P;
return g;
}
vec getinvdev(vec f,int len){
vec g(len+1);
for(int i=1;i<len;i++)g[i]=1LL*f[i-1]*inv[i]%P;
return g;
}
復雜度 \(O(n)\)。
5. 多項式 ln
\[G(x)\equiv \ln F(x)\pmod{x^n} \]
兩邊同時求導
\[G'(x)\equiv \frac{F'(x)}{F(x)}\pmod{x^{n-1}} \]
積回去
\[G(x)\equiv \int \frac{F'(x)}{F(x)}dx\pmod{x^n} \]
over。
復雜度 \(O(n\log n)\).
Code
vec getln(vec f,int len){
vec a=getdev(f,len),b;
getinv(f,b,len);
return getinvdev(mul(a,b,len<<1),len);
}
5.5 多項式牛頓迭代
若
\[G(F(x))\equiv 0\pmod {x^n} \]
\[G(F_0(x))\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}} \]
有
\[F\equiv F_0-\frac{G(F_0)}{G'(F_0)}\pmod{x^n} \]
Proof
根據泰勒展開：
\[G(F)=G(F_0)+G'(F_0)(F-F_0)+G''(F_0)(F-F_0)^2+\cdots \]
因為
\[F\equiv F_0 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]
所以
\[F^2\equiv F_0^2\pmod {x^n} \]
所以
\[0\equiv G(F)\equiv G(F_0)+G'(F_0)(F-F_0) \pmod {x^n} \]
所以
\[F\equiv F_0-\frac{G(F_0)}{G'(F_0)}\pmod{x^n} \]
6 多項式exp
\[G(x)\equiv e^{F(x)}\pmod {x^n} \]
\[\ln G-F\equiv 0 \pmod {x^n} \]
設
\[y=\ln G-F \]
則
\[dy=d(\ln G-F) \]
\[dy=\frac{dG}{G} \]
\[\frac{dy}{dG}=\frac{1}{G} \]
設
\[\ln G_0-F\equiv 0 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]
則根據牛頓迭代公式：
\[G\equiv G_0-\frac{\ln G_0-F}{\frac{1}{G_0}}\equiv G_0(1-\ln G_0+F)\pmod{x^n} \]
於是遞歸即可。
復雜度 \(T(n)=T(\frac{n}{2})+O(n\log n)\)，根據主定理可知 \(T(n)=O(n\log n)\)。
Code
void getexp(vec f,vec&g,int len){
if(g.size()<(len<<2))g.resize(len<<2);
if(len==1){g[0]=1;return;}
getexp(f,g,len+1>>1);
init(len<<1);
vec d=getln(g,len),e(f.begin(),f.begin()+len);
NTT(g,lmt,1),NTT(d,lmt,1),NTT(e,lmt,1);
for(int i=0;i<lmt;i++)g[i]=1LL*(1-d[i]+e[i]+P)*g[i]%P;
NTT(g,lmt,-1);
for(int i=len;i<lmt;i++)g[i]=0;
}
7 多項式快速冪
\[G\equiv F^k\pmod {x^n} \]
\[\ln G\equiv k\ln F\pmod {x^n} \]
\[G\equiv e^{k\ln F}\pmod {x^n} \]
over。
復雜度 \(O(n\log n)\)。
Code
vec getpow(vec f,int len,int k){
vec g(len),h=getln(f,len);
for(int i=0;i<len;i++)h[i]=1LL*h[i]*k%P;
getexp(h,g,len);
return g;
}
8.多項式開根
\[G\equiv \sqrt F\pmod {x^n} \]
\[\ln G\equiv \frac{F}{2}\pmod{x^n} \]
\[G\equiv e^\frac{F}{2}\pmod{x^n} \]
over。
復雜度 \(O(n\log n)\)。
Code
vec getsqrt(vec f,int len){
vec h=getln(f,len),g;
for(int i=0;i<len;i++)h[i]=1LL*h[i]*inv[2]%P;
getexp(h,g,len);
return g;
}
9.多項式三角函數
我們有
\[\begin{cases} e^{ix}=\cos x+i\sin x\\ e^{-ix}=\cos x-i\sin x \end{cases} \]
解得
\[\begin{cases} \cos x=\frac{e^{ix}+e^{-ix}}{2}\\ \sin x=\frac{e^{ix}-e^{-ix}}{2i} \end{cases} \]
又
\[\left(\frac{-1}{P}\right)=(-1)^\frac{P-1}{2}=1 \]
所以存在 \(i_0\in \mathbb Z_P\)，使得 \(i_0^2\equiv -1\pmod P\)。
比如當 \(P=998244353\) 時，\(i=86583718\)。
於是
\[\begin{cases} \cos F=\frac{e^{i_0F}+e^{-i_0F}}{2}\\ \sin F=\frac{e^{i_0F}-e^{-i_0F}}{2i_0} \end{cases} \]
over。
總復雜度 \(O(n\log n)\)。
Q：如果 \(P\) 不是固定的怎么做？
A：暴力枚舉 \(i\) 即可，當然也可以用 Cipolla 算法。
Code
vec sin(vec f,int len){
vec x(len),X,Y,g(len);
for(int i=0;i<len;i++)x[i]=1LL*img*f[i]%P;
getexp(x,X,len),getinv(X,Y,len);
for(int i=0,inv=qpow(img<<1,P-2);i<len;i++)g[i]=1LL*(X[i]-Y[i]+P)%P*inv%P;
return g;
}
vec cos(vec f,int len){
vec x(len),X,Y,g(len);
for(int i=0;i<len;i++)x[i]=1LL*img*f[i]%P;
getexp(x,X,len),getinv(X,Y,len);
for(int i=0;i<len;i++)g[i]=1LL*(X[i]+Y[i])%P*inv[2]%P;
return g;
}
10.多項式反三角函數
和前面多項式\(\ln\)差不多。。。
首先推一發 \(\arcsin\)：
\[\begin{aligned} y&=\arcsin x\\ \sin y&=x\\ \text d(\sin y)&=\text dx\\ \text dy\times \cos y&=\text dx\\ \frac{dy}{dx}&=\frac{1}{\cos y}\\ \frac{dy}{dx}&=\frac{1}{\sqrt{1-\sin y^2}}\\ \frac{dy}{dx}&=\frac{1}{\sqrt{1-x^2}} \end{aligned} \]
返回原題
\[\begin{aligned} G&\equiv \arcsin F\pmod {x^n}\\ G'&\equiv \frac{F'}{\sqrt{1-F^2}}\pmod {x^{n-1}}\\ G&\equiv \int \frac{F'}{\sqrt{1-F^2}}dx\pmod {x^n} \end{aligned} \]
\(y=\arctan x\) 同理可得：
\[G\equiv \int \frac{F'}{1+F^2}dx\pmod {x^n} \]
over。
總復雜度依然是\(O(n\log n)\)。
Code
vec arcsin(vec f,int len){
vec x=getdev(f,len);
init(len<<1);
vec y(lmt);
NTT(f,lmt,1);
for(int i=0;i<lmt;i++)y[i]=(1+P-1LL*f[i]*f[i]%P)%P;
NTT(y,lmt,-1);
for(int i=len;i<lmt;i++)y[i]=0;
vec z=getsqrt(y,len);
y.clear();
getinv(z,y,len);
NTT(x,lmt,1),NTT(y,lmt,1);
for(int i=0;i<lmt;i++)x[i]=1LL*x[i]*y[i]%P;
NTT(x,lmt,-1);
return getinvdev(x,len);
}
vec arctan(vec f,int len){
vec x=getdev(f,len);
init(len<<1);
vec y(lmt),z;
NTT(f,lmt,1);
for(int i=0;i<lmt;i++)y[i]=(1+1LL*f[i]*f[i]%P)%P;
NTT(y,lmt,-1);
for(int i=len;i<lmt;i++)y[i]=0;
getinv(y,z,len);
NTT(x,lmt,1),NTT(z,lmt,1);
for(int i=0;i<lmt;i++)x[i]=1LL*x[i]*z[i]%P;
NTT(x,lmt,-1);
return getinvdev(x,len);
}
二.模數不為NTT模數的部分
1. 多項式乘法（MTT）
orz myy！！！
1.1 9 次NTT的MTT
考慮到值域最大可以達到 \(10^{23}\)，我們可以取 \(3\) 個NTT模數（如 \(998244353\)，\(1004535809\)，\(167772161\)），用 CRT 合並即可。
共計使用 \(9\) 次NTT。
1.2 7 次FFT的MTT
取 \(M_0=[\sqrt M]\)，設 \(F(x)=M_0A(x)+B(x)\)，\(G(x)=M_0C(x)+D(x)\)，所以 \(F(x)G(x)=M_0^2A(x)C(x)+M_0(A(x)D(x)+B(x)C(x))+B(x)D(x)\)。
注意到中間系數為 \(M_0\) 的兩項我們可以直接合並，因此可以省去一次 IDFT，共計 \(7\) 次 DFT。
1.3 DFT 二合一
這里二合一的前置要求是：兩個數列都是實數數列。
我們考慮兩個實多項式\(A(x),B(x)\)，希望能一次DFT同時求出這兩個的DFT。
令
\[P(x)=A(x)+iB(x) \]
\[Q(x)=A(x)-iB(x) \]
為了方便起見，記 \(F_p[k]=P(\omega^k),F_q[k]=Q(\omega^k)\)，並定義 \(X=\frac{2\pi jk}{2L}\)。
所以
\[\begin{aligned} F_p[k]&=A(\omega_{2L}^k)+iB(\omega_{2L}^k)\\ &=\sum_{j=0}^{2L-1}A_j\omega_{2L}^{jk}+iB_j\omega_{2L}^{jk}\\ &=\sum_{j=0}^{2L-1}(A_j+iB_j)(\cos X+i\sin X) \end{aligned} \]
\[\begin{aligned} F_q[k]&=A(\omega_{2L}^k)-iB(\omega_{2L}^k)\\ &=\sum_{j=0}^{2L-1}A_j\omega_{2L}^{jk}-iB_j\omega_{2L}^{jk}\\ &=\sum_{j=0}^{2L-1}(A_j-iB_j)(\cos X+i\sin X)\\ &=\sum_{j=0}^{2L-1}(A_j\cos X+B_j\sin X)+i(A_j\sin X-B_j\cos X)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j\cos X+B_j\sin X)-i(A_j\sin X-B_j\cos X)\right)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j+iB_j)(\cos(-X)+i\sin(-X))\right)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j+iB_j)\omega_{2L}^{-jk}\right)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j+iB_j)\omega_{2L}^{(2L-k)j}\right)\\ &=\operatorname{conj}(F_p[2L-k]) \end{aligned} \]
所以只要一次 \(\text{DFT}\) 即可同時算出 \(F_p,F_q\)，從而
\[\text{DFT}(A[k])=\frac{F_p[k]+F_q[k]}{2} \]
\[\text{DFT}(B[k])=i\frac{F_p[k]-F_q[k]}{2} \]
於是一次DFT即可。
1.4 IDFT 二合一
這里二合一的要求是：IDFT的結果必須是實數數列。
考慮前面兩個式子，我們知道 \(\text{DFT}(A[k]),\text{DFT}(B[k])\) 之后，就可以得到 \(F_p,F_q\)，就可以求出 \(P(x)\)，於是就可以求出 \(A(x),B(x)\)。
於是一次 IDFT 即可。
1.5 4 次 DFT 的 MTT
考慮上面的7次DFT的方法，總共用了4次DFT與3次IDFT。
於是4次DFT兩兩合並壓縮為2次，3次IDFT兩兩合並壓縮為2次。
Q：為什么不能再壓縮了？
A：再壓縮就不符合壓縮的條件了。
總共使用4次\(\text{DFT}\)。
Code
vec MTT(vec f,vec g,int n,int m,int P){
init(n+m);
const int lim=(1<<15)-1;
vector<C>a(lmt),b(lmt),c(lmt),d(lmt);
vec ans(lmt);
for(int i=0;i<n;i++)a[i]=C(f[i]&lim,f[i]>>15);
for(int i=0;i<m;i++)b[i]=C(g[i]&lim,g[i]>>15);
FFT(a,lmt,1),FFT(b,lmt,1);
for(int i=0;i<lmt;i++){
int t=(lmt-i)&(lmt-1);
c[i]=C((a[i].x+a[t].x)*0.5,(a[i].y-a[t].y)*0.5)*b[i];
d[i]=C((a[i].y+a[t].y)*0.5,(a[t].x-a[i].x)*0.5)*b[i];
}
FFT(c,lmt,-1),FFT(d,lmt,-1);
for(int i=0;i<lmt;i++)c[i]=c[i]/lmt,d[i]=d[i]/lmt;
for(int i=0;i<lmt;i++){
long long p=c[i].x+0.5,o=c[i].y+0.5,x=d[i].x+0.5,u=d[i].y+0.5;
ans[i]=1LL*(p%P+((o+x)%P<<15)+(u%P<<30))%P;
}
return ans;
}