目录
同步更新于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} ceil}})
又显然(F(x)G(x) equiv 1 pmod{x^{lceil frac{n}{2} ceil}})
那么(F(x)(G(x)-H(x))equiv 0 pmod{x^{lceil frac{n}{2} ceil}})
即(G(x)-H(x)equiv 0 pmod{x^{lceil frac{n}{2} ceil}})
两边同时平方。
(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))是指在(mod 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))^2equiv 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))^kpmod{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} ight)^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;
}
多项式三角函数
由欧拉公式得到：
[egin{cases}e^{iF(x)}=cos(F(x))+isin(F(x))\
e^{-iF(x)}=cos(F(x))-isin(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)equivsin^{-1}F(x)space ( ext{mod }x^n))
(G'(x)equiv frac{F'(x)}{sqrt{1-F(x)^2}}space ({ ext{mod }x^n}))
(G(x)equiv int frac{F'(x)}{sqrt{1-F(x)^2}} ext dxspace ({ ext{mod }x^n}))
(H(x)equiv an^{-1}F(x)space ( ext{mod }x^n))
(H'(x)equivfrac{F'(x)}{1+F(x)^2}space ( ext{mod }x^n))
(H(x)equivintfrac{F'(x)}{1+F(x)^2} ext dxspace ( ext{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太慢怎么办下次一定改