[營業日誌]萌新的多項式入門
同步更新於caijiのBlog
菜雞不會多項式石錘了
多項式工業入門
一些定義
template<const int mod> struct modint{ int x; modint<mod>(int o=0){x=o;} modint<mod> &operator = (int o){return x=o,*this;} modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;} modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;} modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;} modint<mod> &operator ^=(int b){ modint<mod> a=*this,c=1; for(;b;b>>=1,a*=a)if(b&1)c*=a; return x=c.x,*this; } modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;} modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;} modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;} modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;} modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);} template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;} template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;} template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;} template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;} friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;} friend bool operator ==(modint<mod> a,int b){return a.x==b;} friend bool operator !=(modint<mod> a,int b){return a.x!=b;} bool operator ! () {return !x;} modint<mod> operator - () {return x?mod-x:0;} modint<mod> &operator++(int){return *this+=1;} }; const int N=4e6+5; const int mod=998244353; const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718; struct poly{ vector<modint<mod>>a; modint<mod>&operator[](int i){return a[i];} int size(){return a.size();} void resize(int n){a.resize(n);} }; int rev[N]; inline poly one(){poly a;a.a.push_back(1);return a;}
基礎快速變幻
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;} inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));} inline void ntt(poly&a,int k,int typ){ int n=1<<k; for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]); for(int mid=1;mid<n;mid<<=1){ modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1)); for(int r=mid<<1,j=0;j<n;j+=r){ modint<mod> w=1; for(int k=0;k<mid;k++,w=w*wn){ modint<mod> x=a[j+k],y=w*a[j+k+mid]; a[j+k]=x+y,a[j+k+mid]=x-y; } } } if(typ<0){ modint<mod> inv=modint<mod>(1)/n; for(int i=0;i<n-1;i++)a[i]*=inv; } }
多項式加、減、乘
poly operator +(poly a,poly b){ int n=max(a.size(),b.size());a.resize(n),b.resize(n); for(int i=0;i<n;i++)a[i]+=b[i];return a; } poly operator -(poly a,poly b){ int n=max(a.size(),b.size());a.resize(n),b.resize(n); for(int i=0;i<n;i++)a[i]-=b[i];return a; } inline poly operator*(poly a,poly b){ int n=a.size()+b.size()-1,k=ext(n); a.resize(1<<k),b.resize(1<<k),init(k); ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i]; ntt(a,k,-1),a.resize(n);return a; } inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; } inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; } inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }
多項式求逆
如果多項式\(F\)只有一項,那麼顯然\(G_0\)就是\(F_0\)的逆元。
若有\(n\)項,遞迴求解。
假如我們已知\(F(x)H(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
又顯然\(F(x)G(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
那麼\(F(x)(G(x)-H(x))\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
即\(G(x)-H(x)\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
兩邊同時平方。
\(G(x)^2+H(x)^2-2G(x)H(x) \equiv 0 \pmod{x^n}\)
兩邊同時乘\(F(x)\),再由\(F(x)G(x) \equiv 1 \pmod{x^n}\),得出:
\(G(x) \equiv 2H(x)-F(x)H(x)^2 \pmod{x^n}\)
poly inv(poly F,int k){
int n=1<<k;F.resize(n);
if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
poly G,H=inv(F,k-1);
G.resize(n),H.resize(n<<1),F.resize(n<<1);
for(int i=0;i<n/2;i++)G[i]=H[i]*2;
init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
ntt(H,k+1,-1),H.resize(n);
for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
int n=a.size();
a=inv(a,ext(n)),a.resize(n);return a;;
}
求導數、原函式
\((x^a)^\prime=ax^{a-1}\)
\(\int x^adx=\frac{1}{a+1}x^{a+1}\)
inline poly deriv(poly a){//求導
int n=a.size()-1;
for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
a.resize(n);return a;
}
inline poly inter(poly a){//求原
int n=a.size()+1;a.resize(n);
for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
a[0]=0;return a;
}
多項式\(ln\)
\(G(x)=F(A(x)),F(x)=ln(x)\)
兩邊同時求導,得到:
\(G(x)^\prime=F^\prime(A(x))A^\prime(x)=\frac{A^\prime(x)}{A(x)}\)
在求一遍逆即可。
inline poly ln(poly a){
int n=a.size();
a=inter(deriv(a)*inv(a));
a.resize(n);return a;
}
多項式\(exp\)
計算\(F(x)=e^{A(x)}\)
變形得\(\ln F(x)-A(x)=0\)
設\(G(F(x))=lnF(x)-A(x)\),那麼就是要求這一個函式的零點。
那麼我們把\(F(x)\)看做變數,\(A(x)\)看做常數,對這個進行求導,得\(G^\prime(F(x))=\frac{1}{F(x)}\)
那麼代入牛頓迭代的公式得\(F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod{x^n}\)
\(F(x)\equiv F_0(x)(1-\ln F_0(x)+A(x))\pmod{x^n}\)
(這裡的\(F_0(x)\)是指在\(\bmod x^{\frac n 2}\)下的答案)
然後因為\(A(0)=0\),所以\(F(x)\)的常數項為\(1\)
poly exp(poly a,int k){
int n=1<<k;a.resize(n);
if(n==1)return one();
poly f0=exp(a,k-1);f0.resize(n);
return f0*(one()+a-ln(f0));
}
poly exp(poly a){
int n=a.size();
a=exp(a,ext(n));a.resize(n);return a;
}
多項式開根
設\(H^2(x)\equiv F(x)(\pmod\ x^{\frac n2})\)
\(G(x)\equiv H(x)(\pmod\ x^{\frac n2})\)
\(G(x)-H(x)\equiv 0(\pmod\ x^{\frac n2})\)
\((G(x)-H(x))^2\equiv 0(\pmod\ x^{\frac n2})\)
\(G^2(x)-2H(x)*G(x)+H^2(x)\equiv 0(\pmod\ x^n)\)
\(F(x)-2H(x)*G(x)+H^2(x)\equiv 0(\pmod\ x^n)\)
\(G(x)=\frac {F(x)+H^2(x)}{2H(x)}\)
\(a_0=1\)
poly sqrt(poly F,int k){
int n=1<<k;F.resize(n);
if(n==1){F[0]=1;return F;}
poly G,H=sqrt(F,k-1);
H.resize(n);init(k);poly invH=inv(H);
G.resize(n),H.resize(n<<1),F.resize(n<<1);
init(k+1),ntt(H,k+1,1);
for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i];
ntt(H,k+1,-1),H.resize(n);
for(int i=0;i<n;i++)G[i]=H[i]+F[i];
G=G*invH;G.resize(n);for(int i=0;i<n;i++)G[i]/=2;
return G;
}
inline poly sqrt(poly a){
int n=a.size();
a=sqrt(a,ext(n)),a.resize(n);return a;;
}
\(a_0≠1\)留坑
多項式快速冪
\(a_0=1\)
\(G(x)=(F(x))^k\pmod{x^n}\)
\(\ln G(x)=kF(x)\pmod{x^n}\)
\(G(x)=e^{kF(x)}\pmod{x^n}\)
inline poly pow(poly a,modint<mod> k){//保證a[0]=1
a=ln(a);for(int i=0;i<a.size();i++)a[i]*=k.x;
return exp(a);
}
\(a_0≠1\)
把最低項改為\(1\)
\(F(x)^k=\left(\frac{F(x)}{ax^t} \right)^k(ax)^{tk}\)
實現了一個modpair
,因為\(k\)既要對\(mod\)取模,也要對\(\varphi(mod)\)取模
struct modpair{
//用於快速冪中的次數
modint<mod>k1;modint<mod-1>k2;
struct trueint{
double lg;int x;
trueint &operator=(int o){return x=o,lg=log10(o),*this;}
trueint &operator+=(int o){return lg<=8&&(x+=o,lg=log10(x)),*this;}
trueint &operator*=(int o){return x*=o,lg+=log10(o),*this;}
}k3;
modpair(modint<mod>_1,modint<mod-1>_2,trueint _3):k1(_1),k2(_2),k3(_3){}
modpair(int o=0){k1=o,k2=o,k3=o;}
modpair &operator = (int o){return k1=o,k2=o,k3=o,*this;}
modpair &operator +=(int o){return k1+=o,k2+=o,k3+=o,*this;}
modpair &operator *=(int o){return k1*=o,k2*=o,k3*=o,*this;}
friend modpair operator +(modpair a,int o){return a+=o;}
friend modpair operator *(modpair a,int o){return a*=o;}
modpair operator-(){return modpair(k1,k2,k3);}
};
inline poly pow2(poly a,modpair m){//不保證a[0]=1
int k=0;modint<mod> val;
while(a[k]==0&&k<a.size())k++;
if(k==a.size()||k!=0&&m.k3.lg>8||1ll*m.k3.x*k>=a.size()){
for(int i=0;i<a.size();i++)a[i]=0;return a;}//bye~
val=a[k];poly b;b.resize(a.size()-k);
for(int i=0;i<b.size();i++)b[i]=a[i+k]/val;
b=pow(b,m.k1);for(int i=0;i<a.size();i++)a[i]=0;
for(int i=0;i<b.size()&&i+k*m.k1.x<a.size();i++)
a[i+k*m.k1.x]=b[i]*(val^m.k2.x);
return a;
}
多項式三角函式
由尤拉公式得到:
\[\begin{cases}e^{iF(x)}=\cos(F(x))+i\sin(F(x))\\ e^{-iF(x)}=\cos(F(x))-i\sin(F(x)) \end{cases}\]解方程得到
\[\cos F(x)=\frac{e^{iF(x)}+e^{-iF(x)}}{2} \]\[\sin F(x)=\frac{e^{iF(x)}-e^{-iF(x)}}{2i} \]在取模條件下,\(i=\sqrt{-1}=\sqrt{mod-1}\),即\(mod-1\)的二次剩餘
inline poly sin(poly a){return (exp(a*I)-exp(-a*I))/(I*2);}
inline poly cos(poly a){return (exp(a*I)+exp(-a*I))/2;}
多項式反三角函式
\(G(x)\equiv\sin^{-1}F(x)\space (\text{mod }x^n)\)
\(G'(x)\equiv \frac{F'(x)}{\sqrt{1-F(x)^2}}\space ({\text{mod }x^n})\)
\(G(x)\equiv \int \frac{F'(x)}{\sqrt{1-F(x)^2}}\text dx\space ({\text{mod }x^n})\)
\(H(x)\equiv\tan^{-1}F(x)\space (\text{mod }x^n)\)
\(H'(x)\equiv\frac{F'(x)}{1+F(x)^2}\space (\text{mod }x^n)\)
\(H(x)\equiv\int\frac{F'(x)}{1+F(x)^2} \text dx\space (\text{mod }x^n)\)
inline poly asin(poly a){
poly G=-a*a;G[0]++;
return inter(deriv(a)*inv(sqrt(G)));
}
inline poly atan(poly a){
poly G=a*a;G[0]++;
return inter(deriv(a)*inv(G));
}
多項式黑科技
任意模數NTT(MTT)
把兩個多項式 \(A(x),B(x)\) 分別拆成 \(A_0(x),A_1(x),B_0(x),B_1(x)\)後,求出它們的點值表示。
構造兩個多項式:
\[P(x)=A(x)+iB(x) \]\[Q(x)=A(x)+iB(x) \]由於\(A,B\)的虛部都為\(0\),\(P,Q\)的每一項係數都互為共軛。對\(P\)做一遍DFT,可以std::conj
在\(O(n)\)求出\(Q\)的點值表示。然後\(A_0,A_1,B_0,B_1\)的點值就是有手就行了
接下求\(A,B\)之間的兩兩乘積,乘出來對\(A_0B_0,A_1B_0,A_0B_0,A_0B_1\)做IDFT。
由於虛部不為0,但仍然可以借用上述的方法。構造兩個多項式:
\[P(x)=A_0(x)B_0(x)+iA_1(x)B_0(x) \]\[Q(X)=A_0(x)B_1(x)+iA_1(x)B_1(x) \]分別做IDFT,由於 \(A_0B_0,A_0B_1,A_1B_0,A_1B_1\)這四個多項式捲起來後的係數表示中虛部一定為\(0\),那麼\(P\)的實部與虛部就分別為\(A_0(x)B_0(x)\)和\(A_1(x)B_0(x)\),而\(Q\)就位\(A_0(x)B_1(x)\)與\(A_1(x)B_1(x)\)
給出部分程式碼:
typedef std::complex<double>complex;
const int N=4e6+10;const double PI=acos(-1);const complex I=complex(0,1);
int rev[N];complex Wn[N];int M,mod;
int ksm(int x,int y) {
int re=1;
for(;y;y>>=1,x=1LL*x*x%mod)if(y&1)re=1LL*re*x%mod;
return re;
}
inline long long num(complex x){double d=x.real();return d<0?(long long)(d-0.5)%mod:(long long)(d+0.5)%mod;}
struct poly{
std::vector<complex>a0,a1;
int size(){return a0.size();}
void resize(int n){a0.resize(n);a1.resize(n);}
void set(int x,long long y){
y%=mod;
a0[x]=y/M;
a1[x]=y%M;
}
long long get(int x){return (M*M*num(a0[x].real())%mod +
M*(num(a0[x].imag())+num(a1[x].real()))%mod+num(a1[x].imag()))%mod;}
long long val(int x){return (long long)(M*a0[x].real()+a1[x].real()+mod)%mod;}
};
poly operator+(poly a,poly b){
int n=std::max(a.size(),b.size());a.resize(n),b.resize(n);
for(int i=0;i<n;i++)a.set(i,a.val(i)+b.val(i));return a;
}
poly operator-(poly a,poly b){
int n=std::max(a.size(),b.size());a.resize(n),b.resize(n);
for(int i=0;i<n;i++)a.set(i,a.val(i)-b.val(i));return a;
}
inline poly one(){poly a;a.resize(1);a.set(0,1);return a;}
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){
int n=1<<k;
for(int i=0;i<n;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
for(int i=0;i<n;i++)
Wn[i]={cos(PI/n*i),sin(PI/n*i)};
}
void FFT(std::vector<complex>&A,int n,int t){
if(t<0)for(int i=1;i<n;i++)if(i<(n-i))std::swap(A[i],A[n-i]);
for(int i=0;i<n;i++)
if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(int m=1;m<n;m<<=1)
for(int i=0;i<n;i+=m<<1)
for(int k=i;k<i+m;k++){
complex W=Wn[1ll*(k-i)*n/m];
complex a0=A[k],a1=A[k+m]*W;
A[k]=a0+a1;A[k+m]=a0-a1;
}
if(t<0)for(int i=0;i<n;i++)A[i]/=n;
}
void MTT(poly &A,int n,int t){
for(int i=0;i<n;i++)A.a0[i]=A.a0[i]+I*A.a1[i];
FFT(A.a0,n,t);
for(int i=0;i<n;i++)A.a1[i]=std::conj(A.a0[i?n-i:0]);
for(int i=0;i<n;i++){
complex p=A.a0[i],q=A.a1[i];
A.a0[i]=(p+q)*0.5;A.a1[i]=(q-p)*0.5*I;
}
}
inline poly operator*(poly a,poly b){
int n=a.size()+b.size()-1,k=ext(n);
a.resize(1<<k),b.resize(1<<k),init(k);
MTT(a,1<<k,1);MTT(b,1<<k,1);
for(int i=0;i<(1<<k);i++){
complex p=a.a0[i]*b.a0[i]+I*a.a1[i]*b.a0[i];
complex q=a.a0[i]*b.a1[i]+I*a.a1[i]*b.a1[i];
a.a0[i]=p,a.a1[i]=q;
}
FFT(a.a0,1<<k,-1);FFT(a.a1,1<<k,-1);a.resize(n);
long long tmp;for(int i=0;i<n;i++)
tmp=a.get(i),a.set(i,tmp);
return a;
}
inline poly deriv(poly a){//求導
int n=a.size()-1;
for(int i=0;i<n;i++)a.set(i,a.val(i+1)*(i+1));
a.resize(n);return a;
}
inline poly inter(poly a){//求原
int n=a.size()+1;a.resize(n);
for(int i=n;i>=1;i--)a.set(i,a.val(i-1)*ksm(i,mod-2));
a.set(0,0);return a;
}
poly inv(poly F,int k){
int n=1<<k;F.resize(n);
if(n==1){F.set(0,ksm(F.val(0),mod-2));return F;}
poly G,H=inv(F,k-1);
G.resize(n),H.resize(n<<1),F.resize(n<<1);
for(int i=0;i<n/2;i++)G.set(i,H.val(i)*2);
H=H*H;H.resize(n);H=H*F;H.resize(n);
G=G-H;return G;
}
inline poly inv(poly a){
int n=a.size();
a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly ln(poly a){
int n=a.size();
a=inter(deriv(a)*inv(a));
a.resize(n);return a;
}
別問我vector
太慢怎麼辦下次一定改