[营业日志]萌新的多项式入门
同步更新于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;
}
多项式三角函数
由欧拉公式得到：
解方程得到
在取模条件下，\(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)\)后，求出它们的点值表示。
构造两个多项式：
由于\(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，但仍然可以借用上述的方法。构造两个多项式：
分别做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太慢怎么办下次一定改
