多項式全家桶——Part.3 多項式求逆、除法、開根號
多項式全家桶正式進入正片。
有一說一，感覺之前寫的其實都是寫不重要或特別基礎的東東。
現在才是一些比較有難度的東東。
多項式求逆
介紹
多項式求逆即為求一個多項式的乘法逆元。
它長這樣：
\[A(x)*B(x)\equiv 1(mod\ x^n) \]
然後多項式\(B(x)\)就叫做多項式\(A(x)\)的逆元。
這玩意和普通的乘法逆元有一點不同的就是，它的模是一個\(x^n\)，其意義為把結果中次數大於等於n的項都給忽略掉。
過程
引理一： 考慮一個乘法逆元滿足：\(A(x)*B(x)\equiv 1(mod\ x^n)\)，那麼我們可以得到另一個柿子：\(A(x)*B(x)\equiv 1(mod\ x^m)\)
證明？把卷積形式寫出來就知道了。
考慮設\(B(x)\)表示\(A(x)\)在模\(x^n\)意義下的逆元。
設\(G(x)\)表示\(A(x)\)在模\(x^{\lceil \frac n 2 \rceil}\)意義下的逆元。
由引理得：
\[B(x)\equiv G(x)(mod\ x^{\lceil \frac n 2 \rceil})\\B(x)-G(x)\equiv 0(mod\ x^{\lceil \frac n 2 \rceil})\]
開方（這一步我也不能嚴謹證明，但是打表發現這是對的）
\[B(x)^2+G(x)^2-2*B(x)*G(x)\equiv 0(mod\ x^n) \]
乘個\(A(x)\)進去，然後由於是模\(x^n\)意義下，所以可以消去一些\(B(x)\):
\[B(x)+A(x)*G(x)^2-2*G(x)\equiv 0(mod\ x^n) \\B(x)\equiv G(x)(2-A(x)*G(x))(mod\ x^n)\]
這時，就可以用\(G(x)\)來求出\(B(x)\)。具體過程就遞迴求解即可（當然不覺得麻煩的話打分治FFT也是可以），中間需要用到NTT。
學習資料：
https://blog.csdn.net/a_forever_dream/article/details/102483602
https://www.cnblogs.com/LanrTabe/p/11359952.html
https://www.cnblogs.com/guoshaoyang/p/11296025.html
程式碼
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const long long mo=998244353;
const int maxn=300010;
long long a[maxn],b[maxn],rev,A[maxn],G[maxn];
int jl,n,rec[maxn],w[maxn],up,dep;
long long qsm(long long a,long long b)
{
long long t=1;
long long y=a;
while (b>0)
{
if ((b&1)==1) t=t*y%mo;
y=y*y%mo;
b/=2;
}
return t;
}
void FFT(long long a[],int n,int inv)
{
if (n==1) return;
int mid=n/2;
for (int i=0;i<n;i++)
{
rec[i]=(rec[i>>1]>>1)|((i&1)<<(dep-1));
if (i<rec[i]) swap(a[i],a[rec[i]]);
}
for (int len=2;len<=n;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<n;k+=len)
{
int p,q;
q=w[j*(n/len)]*a[k+mid]%mo,p=a[k];
a[k+mid]=(p-q+mo)%mo;
a[k]=(p+q+mo)%mo;
}
}
}
}
void solve(int n)
{
if (n==1)
{
b[0]=qsm(a[0],mo-2);
}
else
{
double op=n;
int oq=ceil(op/2);
solve(oq);
up=1;dep=0;
while (up<=n+jl) up=up*2,dep++;
w[0]=1;
rev=qsm(3,(mo-1)/up);
for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
for (int i=0;i<up;i++)
{
G[i]=0;
A[i]=a[i];
if (i<oq) G[i]=b[i];
}
FFT(A,up,1);
FFT(G,up,1);
for (int i=0;i<up;i++)
{
A[i]=G[i]*(2-A[i]*G[i]%mo+mo)%mo;
}
for (int i=0;i<=up/2;i++)
{
swap(w[i],w[up-i]);
}
FFT(A,up,-1);
for (int i=0;i<n;i++) b[i]=A[i]*qsm(up,mo-2)%mo;
int kk=0;
}
}
int main()
{
scanf("%d",&n);
jl=n;
for (int i=0;i<n;i++)
{
scanf("%lld",&a[i]);
}
solve(n);
for (int i=0;i<n;i++)
{
printf("%lld ",b[i]);
}
printf("\n");
}
多項式除法
介紹
這玩意是求一個長成這樣的柿子的：
\[A(x)=B(x)*D(x)+R(x) \]
其中\(A(x)\)和\(B(x)\)已經給出，且\(A(x)\)次數為n，\(B(x)\)次數為m，求\(D(x)\)和\(R(x)\)
也就是除法運算中求商和餘數，話說如果上過初中就應該明白長除法。
過程
引理二：
考慮一個反轉操作為將一個多項式的係數調換，用數學形式表達即為：
\[Pr(x)=x^n*P(\frac 1 x)=\sum_{i=0}^n a_{n-i}*x^i \]
證明顯然。
那麼現在就推柿子：
\[A(x)=B(x)D(x)+R(x)\\A(\frac 1 x)=B(\frac 1 x)D(\frac 1 x)+R(\frac 1 x)\\x^nA(\frac 1 x)=x^nB(\frac 1 x)D(\frac 1 x)+x^nR(x)\\x^nA(\frac 1 x)=x^mB(\frac 1 x)x^{n-m}D(\frac 1 x)+x^{n-m+1}x^{m-1}R(x) \]
然後再引入引理：
\[Ar(x)=Br(x)Dr(x)+x^{n-m+1}Rr(x) \]
由於我們的\(Dr(x)\)的次數是小於等於\(n-m\)的，所以如果對原式取模只要大於之即可，而\(A(x)\)、\(B(x)\)都已知，所以不必理會。
\[Ar(x)\equiv Br(x)Dr(x)(mod\ x^{n-m+1}) \]
然後我們就可以直接求\(Br(x)\)的逆元，再乘\(Ar(x)\)即可得到\(Dr(x)\)，再把係數反過來即可求出\(D(x)\)。
當然，求出\(D(x)\)之後，再求\(R(x)\)就很簡單了。
學習資料
https://www.cnblogs.com/knife-rose/p/12056617.html
https://www.cnblogs.com/owenyu/p/6724611.html
程式碼
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const long long mo=998244353;
const int maxn=300010;
long long a[maxn],b[maxn],ar[maxn],br[maxn],d[maxn],D[maxn],r[maxn],rev,A[maxn],G[maxn];
int jl,n,m,rec[maxn],w[maxn],up,dep;
long long qsm(long long a,long long b)
{
long long t=1;
long long y=a;
while (b>0)
{
if ((b&1)==1) t=t*y%mo;
y=y*y%mo;
b/=2;
}
return t;
}
void FFT(long long a[],int n,int inv)
{
if (n==1) return;
int mid=n/2;
for (int i=0;i<n;i++)
{
rec[i]=(rec[i>>1]>>1)|((i&1)<<(dep-1));
if (i<rec[i]) swap(a[i],a[rec[i]]);
}
for (int len=2;len<=n;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<n;k+=len)
{
int p,q;
q=w[j*(n/len)]*a[k+mid]%mo,p=a[k];
a[k+mid]=(p-q+mo)%mo;
a[k]=(p+q+mo)%mo;
}
}
}
}
void solve(int n)
{
if (n==1)
{
d[0]=qsm(br[0],mo-2);
}
else
{
double op=n;
int oq=ceil(op/2);
solve(oq);
up=1;dep=0;
while (up<=n+jl) up=up*2,dep++;
w[0]=1;
rev=qsm(3,(mo-1)/up);
for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
for (int i=0;i<up;i++)
{
G[i]=0;
A[i]=br[i];
if (i<oq) G[i]=d[i];
}
FFT(A,up,1);
FFT(G,up,1);
for (int i=0;i<up;i++)
{
A[i]=G[i]*(2-A[i]*G[i]%mo+mo)%mo;
}
for (int i=0;i<=up/2;i++)
{
swap(w[i],w[up-i]);
}
FFT(A,up,-1);
for (int i=0;i<n;i++) d[i]=A[i]*qsm(up,mo-2)%mo;
int kk=0;
}
}
int main()
{
scanf("%d%d",&n,&m);
for (int i=0;i<=n;i++)
{
scanf("%d",&a[i]);
ar[n-i]=a[i];
}
for (int i=0;i<=m;i++)
{
scanf("%d",&b[i]);
br[m-i]=b[i];
}
jl=m;
solve(n-m+1);
up=1;dep=0;
while (up<=2*n-m) up=up*2,dep++;
w[0]=1;
rev=qsm(3,(mo-1)/up);
for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
FFT(ar,up,1);
FFT(d,up,1);
for (int i=0;i<up;i++)
{
d[i]=d[i]*ar[i]%mo;
}
for (int i=0;i<=up/2;i++)
{
swap(w[i],w[up-i]);
}
FFT(d,up,-1);
for (int i=0;i<=n-m;i++) d[i]=d[i]*qsm(up,mo-2)%mo;
int j=0;
for (int i=n-m;i>=0;i--)
{
printf("%d ",d[i]);
D[j++]=d[i];
}
printf("\n");
up=1;dep=0;
while (up<=n) up=up*2,dep++;
w[0]=1;
rev=qsm(3,(mo-1)/up);
for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
FFT(b,up,1);
FFT(D,up,1);
for (int i=0;i<up;i++)
{
b[i]=b[i]*D[i]%mo;
}
for (int i=0;i<=up/2;i++)
{
swap(w[i],w[up-i]);
}
FFT(b,up,-1);
for (int i=0;i<n;i++) r[i]=(a[i]-b[i]*qsm(up,mo-2)%mo+mo)%mo;
for (int i=0;i<m;i++)
{
printf("%d ",r[i]);
}
printf("\n");
}
多項式開根
介紹
這玩意是求一個這樣的東東：
\[B(x)^2\equiv A(x)(mod\ x^n) \]
其中\(A(x)\)次數為\(n-1\)，要求一個\(B(x)\)。注意！由於是模意義下的多項式，所以我們的\(B(x)\)的次數不一定為\(\frac{n-1}2\)。（廢話）
同樣考慮用什麼遞迴分治的思想來搞。
前置芝士
這裡就需要一點額外的芝士了。
就這題而言，有極大的概率會用到二次剩餘。
那麼來康康怎麼搞。
首先，二次剩餘的定義為：
\[\exists 整數x，使得x^2\equiv n(mod\ p)則稱n為p的二次剩餘 \]
再定義一個東東：
\[(\frac n p)=\left\{\begin{matrix}1(n為p二次剩餘) & & \\ -1(n不為p的二次剩餘) & & \\0(p|n) \end{matrix}\right.\]
這個東東叫做勒讓德符號（聽說還有個叫做高斯記號的東東）。當p為奇質數的時候（\(p=2\)的時候每個整數就都是其二次剩餘了），還有個性質：
\[(\frac n p)=n^{\frac {p-1} 2} \]
那麼我們再引入兩個東東：（現在所有的p都為奇質數）
定理一： 對於\(x^2\equiv n(mod\ p)\)中有\(\frac {p-1}2\)個n是p的二次剩餘。
證明？我不費！
定理二： 首先我們根據定理一，可以考慮random一個數\(a\in[0,p-1]\)。設\(w=a^2-n\)。若找到一個\(w\)滿足：\((\frac w p)=-1\)或者寫成這樣：\(w^{\frac{p-1} 2}\equiv -1\)，則\((a+\sqrt w)^{\frac{p+1} 2}\)為一組x值。
證明？這個我費！
\[(a+\sqrt w)^p=\sum_{i=0}^pa^i\sqrt w^{p-i}*C_p^i \]
在模意義下這玩意只有當\(i=0\)或\(i=p\)時有值。
那麼得到：
\[(a+\sqrt w)^p\equiv a^p+\sqrt w^p \]
根據費馬小定理：\(a^{p-1}\equiv 1\)
根據條件：\(w^{\frac{p+1} 2}\equiv 1\)
所以可以得到：
\[a^p+\sqrt w^p\equiv a-\sqrt w \]
那麼把原式乘上這玩意即可得到：
\[(a+\sqrt w)^{p+1}\equiv(a+\sqrt w)(a-\sqrt w)\equiv a^2-w\equiv n \]
證畢。
那麼得到定理後，我們可以去計算x值了。
一點注意的地方就是在跑快速冪時，我們需要分實數部分和整數部分來計算，因為\(\sqrt w\)我們並不能很好表示，所以要帶進去一起算。
程式碼：待填。QWQ
過程
n=1
既然我們知道了二次剩餘，那麼我們可以快速得到\(n=1\)時的情況。
\[B0^2\equiv A0(mod\ x) \]
其實有時候並不需要什麼二次剩餘，因為可能良心出題人給的\(A0\)是個比較好看的數，比如說1。
n>1
現在才是正軌。
考慮原式：
\[B(x)^2\equiv A(x)(mod\ x^n) \]
用同樣的套路設\(G(x)\)表示\(A(x)\)在模\(x^{\lceil \frac n 2 \rceil}\)意義下的開方。
因此得到：
\[G(x)^2\equiv A(x)(mod\ x^{\lceil \frac n 2 \rceil})\\G(x)^2-A(x)\equiv 0(mod\ x^{\lceil \frac n 2 \rceil})\\(G(x)^2-A(x))^2\equiv 0(mod\ x^n)\\G(x)^4-2G(x)^2A(x)+A(x)^2\equiv 0(mod\ x^n)\\G(x)^4+2G(x)^2A(x)+A(x)^2\equiv 4G(x)^2A(x)(mod\ x^n)\\(G(x)^2+A(x))^2\equiv4G(x)^2A(x)(mod\ x^n)\\ \frac{(G(x)^2+A(x))^2}{4G(x)^2}\equiv A(x)(mod\ x^n)\]
再看到原式，可以得到：
\[B(x)=\frac{G(x)^2+A(x)}{2G(x)} \]
愉快遞迴求解吧~
學習資料：
百度百科
https://blog.csdn.net/stevensonson/article/details/85845334
https://www.cnblogs.com/zhangleo/p/11010302.html
https://www.cnblogs.com/xiefengze1/p/9373272.html
程式碼
話說遞迴版的常數是真滴大，洛谷反正我是TLE了QWQ。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const long long mo=998244353;
const long long inv2=499122177;
const int maxn=300010;
int a[maxn],b[maxn],d[maxn],rev,A[maxn],G[maxn],D[maxn];
int jl,n,rec[maxn],w[maxn],up,dep;
int jia(int a,int b)
{
return (a+b)%mo;
}
int cheng(int a,int b)
{
return 1ll*a*b%mo;
}
inline long long qsm(long long a,long long b)
{
long long t=1;
long long y=a;
while (b>0)
{
if ((b&1)==1) t=t*y%mo;
y=y*y%mo;
b/=2;
}
return t;
}
inline void FFT(int a[],int n,int inv)
{
if (n==1) return;
int mid=n/2;
for (int i=0;i<n;i++)
{
rec[i]=(rec[i>>1]>>1)|((i&1)<<(dep-1));
if (i<rec[i]) swap(a[i],a[rec[i]]);
}
for (int len=2;len<=n;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<n;k+=len)
{
int p,q;
q=cheng(w[j*(n/len)],a[k+mid]),p=a[k];
a[k+mid]=jia(p+mo,-q)%mo;
a[k]=jia(p,q)%mo;
}
}
}
}
inline void solve(int n)
{
if (n==1)
{
d[0]=qsm(b[0],mo-2);
}
else
{
double op=n;
int oq=ceil(op/2);
solve(oq);
up=1;dep=0;
while (up<=n+jl) up=up*2,dep++;
w[0]=1;
rev=qsm(3,(mo-1)/up);
for (register int i=1;i<=up;i++) w[i]=cheng(w[i-1],rev)%mo;
for (register int i=0;i<up;i++)
{
G[i]=0;A[i]=0;
if (i<n) A[i]=b[i];
if (i<oq) G[i]=d[i];
}
FFT(A,up,1);
FFT(G,up,1);
for (register int i=0;i<up;i++)
{
A[i]=cheng(G[i],(2-cheng(A[i],G[i])+mo)%mo);
}
for (register int i=0;i<=up/2;i++)
{
swap(w[i],w[up-i]);
}
FFT(A,up,-1);
for (register int i=0;i<n;i++) d[i]=cheng(A[i],qsm(up,mo-2));
for (register int i=n;i<=up;i++) d[i]=0;
int kk=0;
}
}
inline void solve1(int n)
{
if (n==1)
{
b[0]=1;
}
else
{
double op=n;
int oq=ceil(op/2);
solve1(oq);
solve(n);
up=1;dep=0;
while (up<=n+jl) up=up*2,dep++;
w[0]=1;
rev=qsm(3,(mo-1)/up);
for (register int i=1;i<=up;i++) w[i]=cheng(w[i-1],rev);
for (register int i=0;i<up;i++)
{
A[i]=0;
if (i<n) A[i]=a[i];
}
FFT(A,up,1);FFT(b,up,1);FFT(d,up,1);
for (register int i=0;i<up;i++) b[i]=cheng(jia(b[i],cheng(A[i],d[i])),inv2);
for (register int i=0;i<=up/2;i++) swap(w[i],w[up-i]);
FFT(b,up,-1);
for (register int i=0;i<n;i++) b[i]=cheng(b[i],qsm(up,mo-2));
for (register int i=n;i<=up;i++) b[i]=0;
int kk=0;
}
}
int main()
{
scanf("%d",&n);
jl=n;
for (int i=0;i<n;i++)
{
scanf("%d",&a[i]);
}
solve1(n);
for (int i=0;i<n;i++)
{
printf("%d ",b[i]);
}
printf("\n");
}