【总结】多项式基础操作
多项式乘法
\(A * B=\sum_{k=0}^{n}(\sum_{k=i+j}a_ib_j)x^k\)
快速傅里叶变换(FFT)
#include<bits/stdc++.h>
#define ui unsigned int
#define ll long long
#define db double
#define ld long double
#define ull unsigned long long
#define ft first
#define sd second
#define PII std::pair<int,int>
#define PLL std::pair<ll,ll>
#define pb(a) push_back(a)
#define mp(a,b) std::make_pair(a,b)
#define ITR(a,b) for(auto a:b)
#define REP(a,b,c) for(register int a=(b),a##end=(c);a<=a##end;++a)
#define DEP(a,b,c) for(register int a=(b),a##end=(c);a>=a##end;--a)
const int MAXN=1<<21;
const db Pi=acos(-1);
int n1,n2,n,m,cnt,rev[MAXN];
struct Complex{
db real,imag;
inline Complex operator + (const Complex &A) const {
return (Complex){real+A.real,imag+A.imag};
};
inline Complex operator - (const Complex &A) const {
return (Complex){real-A.real,imag-A.imag};
};
inline Complex operator * (const Complex &A) const {
return (Complex){real*A.real-imag*A.imag,real*A.imag+imag*A.real};
};
};
Complex a[MAXN],b[MAXN];
template<typename T> inline void read(T &x)
{
T data=0,w=1;
char ch=0;
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
x=data*w;
}
template<typename T> inline void write(T x,char ch='\0')
{
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
if(ch!='\0')putchar(ch);
}
template<typename T> inline bool chkmin(T &x,T y){return y<x?(x=y,true):false;}
template<typename T> inline bool chkmax(T &x,T y){return y>x?(x=y,true):false;}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline void FFT(Complex *A,int tp)
{
REP(i,0,n-1)if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(register int l=2;l<=n;l<<=1)
{
Complex wn=(Complex){cos(2*Pi/l),sin(tp*2*Pi/l)};
for(register int i=0;i<n;i+=l)
{
Complex w=(Complex){1,0};
for(register int j=0;j<(l>>1);++j)
{
Complex A1=A[i+j],A2=A[i+j+(l>>1)]*w;
A[i+j]=A1+A2;A[i+j+(l>>1)]=A1-A2;
w=w*wn;
}
}
}
}
int main()
{
read(n1);read(n2);n1++;n2++;
REP(i,0,n1-1)scanf("%lf",&a[i].real);
REP(i,0,n2-1)scanf("%lf",&b[i].real);
m=n1+n2-1;
for(n=1;n<m;n<<=1,++cnt);
REP(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
FFT(a,1);FFT(b,1);
REP(i,0,n-1)a[i]=a[i]*b[i];
FFT(a,-1);
REP(i,0,m-1)printf("%d ",(int)(a[i].real/n+0.5));
return 0;
}
快速数论变换(NTT)
用原根代替单位根
#include<bits/stdc++.h>
#define ui unsigned int
#define ll long long
#define db double
#define ld long double
#define ull unsigned long long
#define ft first
#define sd second
#define PII std::pair<int,int>
#define PLL std::pair<ll,ll>
#define pb(a) push_back(a)
#define mp(a,b) std::make_pair(a,b)
#define ITR(a,b) for(auto a:b)
#define REP(a,b,c) for(register int a=(b),a##end=(c);a<=a##end;++a)
#define DEP(a,b,c) for(register int a=(b),a##end=(c);a>=a##end;--a)
const int MAXN=1<<21,Mod=998244353;
int n1,n2,n,m,rev[MAXN],cnt;
ll pw[MAXN],ipw[MAXN],a[MAXN],b[MAXN];
template<typename T> inline void read(T &x)
{
T data=0,w=1;
char ch=0;
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
x=data*w;
}
template<typename T> inline void write(T x,char ch='\0')
{
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
if(ch!='\0')putchar(ch);
}
template<typename T> inline bool chkmin(T &x,T y){return y<x?(x=y,true):false;}
template<typename T> inline bool chkmax(T &x,T y){return y>x?(x=y,true):false;}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline ll qexp(ll a,ll b)
{
ll res=1;
while(b)
{
if(b&1)res=res*a%Mod;
a=a*a%Mod;
b>>=1;
}
return res;
}
inline void NTT(ll *A,int tp)
{
REP(i,0,n-1)if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(register int l=2;l<=n;l<<=1)
{
ll wn=tp>0?pw[l]:ipw[l];
for(register int i=0;i<n;i+=l)
{
ll w=1;
for(register int j=0;j<(l>>1);++j)
{
ll A1=A[i+j],A2=1ll*A[i+j+(l>>1)]*w%Mod;
A[i+j]=(A1+A2)%Mod;A[i+j+(l>>1)]=(A1-A2+Mod)%Mod;
w=1ll*w*wn%Mod;
}
}
}
}
int main()
{
read(n1);read(n2);n1++;n2++;
REP(i,0,n1-1)read(a[i]);
REP(i,0,n2-1)read(b[i]);
m=n1+n2-1;
for(n=1;n<m;n<<=1,++cnt);
ll inv3=qexp(3,Mod-2),invn=qexp(n,Mod-2);
for(register int l=1;l<=n;l<<=1)pw[l]=qexp(3,(Mod-1)/l),ipw[l]=qexp(inv3,(Mod-1)/l);
REP(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
NTT(a,1);NTT(b,1);
REP(i,0,n-1)a[i]*=b[i];
NTT(a,-1);
REP(i,0,m-1)write(1ll*a[i]*invn%Mod,' ');
return 0;
}
任意模数情况下的NTT(MTT)
将每个系数分解为 \(a_i=Ap_i*t+Ak_i~(b_i=Bp_i*t+Bk_i)\),\(t\) 自己定,相当于转成 \(t\) 进制数,但要保证 \(t\) 范围FFT不会爆精度。
\(\begin{aligned} A * B & = \sum_{k=0}^n(\sum_{k=i+j}a_ib_j)x^k \\ & = \sum_{k=0}^n(\sum_{k=i+j}(Ap_{i}t+Ak_{i})(Bp_jt+Bk_j))x^k \\ & = \sum_{k=0}^n(\sum_{k=i+j}Ap_iBp_jt^2+Ap_iBk_jt+Ak_iBp_jt+Ak_iBk_j)x^k \end{aligned}\)
然后做四次FFT,最后转回 \(10\) 进制累加起来
(图省事,没写共轭优化,并且还多了一次IDFT)
#include<bits/stdc++.h>
#define ui unsigned int
#define ll long long
#define db double
#define ld long double
#define ull unsigned long long
#define ft first
#define sd second
#define PII std::pair<int,int>
#define PLL std::pair<ll,ll>
#define pb(a) push_back(a)
#define mp(a,b) std::make_pair(a,b)
#define ITR(a,b) for(auto a:b)
#define REP(a,b,c) for(register int a=(b),a##end=(c);a<=a##end;++a)
#define DEP(a,b,c) for(register int a=(b),a##end=(c);a>=a##end;--a)
const int MAXN=1<<21,M=32767;
const ld Pi=acos(-1);
int n1,n2,n,m,cnt,p,rev[MAXN];
ll ans[MAXN];
struct Complex{
ld real,imag;
inline Complex operator + (const Complex &A) const {
return (Complex){real+A.real,imag+A.imag};
};
inline Complex operator - (const Complex &A) const {
return (Complex){real-A.real,imag-A.imag};
};
inline Complex operator * (const Complex &A) const {
return (Complex){real*A.real-imag*A.imag,real*A.imag+imag*A.real};
};
};
Complex a1[MAXN],a2[MAXN],b1[MAXN],b2[MAXN],c[MAXN],w[MAXN];
template<typename T> inline void read(T &x)
{
T data=0,w=1;
char ch=0;
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
x=data*w;
}
template<typename T> inline void write(T x,char ch='\0')
{
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
if(ch!='\0')putchar(ch);
}
template<typename T> inline bool chkmin(T &x,T y){return y<x?(x=y,true):false;}
template<typename T> inline bool chkmax(T &x,T y){return y>x?(x=y,true):false;}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline void FFT(Complex *A,int tp)
{
REP(i,0,n-1)if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(register int l=2;l<=n;l<<=1)
{
for(register int i=0;i<(l>>1);++i)w[i]=(Complex){cos(2*Pi*i/l),sin(tp*2*Pi*i/l)};
for(register int i=0;i<n;i+=l)
for(register int j=0;j<(l>>1);++j)
{
Complex A1=A[i+j],A2=A[i+j+(l>>1)]*w[j];
A[i+j]=A1+A2;A[i+j+(l>>1)]=A1-A2;
}
}
}
inline void MTT(Complex *a,Complex *b,ll Bas)
{
REP(i,0,n-1)c[i]=a[i]*b[i];FFT(c,-1);
REP(i,0,m-1)(ans[i]+=1ll*(ll)(c[i].real/n+0.5)%p*Bas%p)%=p;
}
int main()
{
read(n1);read(n2);read(p);n1++;n2++;
int x;
REP(i,0,n1-1)read(x),a1[i].real=(ld)((int)x/M),a2[i].real=(ld)((int)x%M);
REP(i,0,n2-1)read(x),b1[i].real=(ld)((int)x/M),b2[i].real=(ld)((int)x%M);
m=n1+n2-1;
for(n=1;n<m;n<<=1,++cnt);
REP(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
FFT(a1,1);FFT(a2,1);FFT(b1,1);FFT(b2,1);
MTT(a1,b1,1ll*M*M%p);
MTT(a1,b2,M%p);
MTT(a2,b1,M%p);
MTT(a2,b2,1);
REP(i,0,m-1)write(ans[i],' ');
return 0;
}
多项式求导
多项式积分
多项式开根
已知 \(B(x)\) ,求 \(A(x)\) 满足\(A(x)^2 \equiv B(x) \pmod{x^n}\)
假设已知 \(F(x)'\) ,使得 \(F(x)'^2\equiv B(x) \pmod{x^{\lceil\frac{p}{2}\rceil}}\),想要得到 \(F(x)\) ,使得 \(F(x)^2 \equiv B(x) \pmod{x^p}\)
是有两个解的,这里只推其中一个解的式子
递归、倍增均可求解
多项式逆元
已知 \(B(x)\) ,求 \(A(x)\) 满足\(A(x)B(x) \equiv 1 \pmod{x^n}\)
假设已知 \(F(x)'\) ,使得 \(F(x)'B(x)\equiv 1 \pmod{x^{\lceil\frac{p}{2}\rceil}}\),想要得到 \(F(x)\) ,使得 \(F(x)B(x) \equiv 1 \pmod{x^p}\)
递归、倍增均可求解
多项式带余除法
已知 \(F(x)\) 和 \(G(x)\) ,求 \(Q(x)\) 和 \(R(x)\) 满足\(F(x)=G(x)Q(x)+R(x)\)
前置:对于 \(n\) 次多项式 \(F(x)\),\(F_R(x)=x^nF(\frac{1}{x})\) ,发现即是 \(reverse\) ;\(F(x)\) 是 \(n\) 次多项式,\(G(x)\) 是 \(m\) 次多项式,那么 \(Q(x)\) 次数为 \(n-m\) ,\(R(x)\) 次数小于 \(m\)
因为 \(Q(x)\) 的次数为 \(n-m\) ,所以可以求出。求出 \(Q(x)\) 后,\(R(x)=F(x)-G(x)Q(x)\)
多项式对数
已知 \(B(x)\) ,求 \(A(x)\) 满足\(A(x) \equiv \ln B(x) \pmod{x^n}\)
同时求导
多项式指数
已知 \(B(x)\) ,求 \(A(x)\) 满足\(A(x) \equiv \exp B(x) \pmod{x^n}\)
带入牛顿迭代公式 \(F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^n}\) ,\(G(F(x))=\ln F(x)-B(x) \pmod{x^n}\) 可以得到:
递归、倍增均可求解
封装
板子要满足的几个条件:
-
模数是 \(998244353\)
-
开根,对数时输入多项式的 \(0\) 次项是 \(1\)
-
指数时输入多项式的 \(0\) 次项是 \(0\)
namespace poly{
const int Mod=998244353;
int rev[MAXN];
ll b[MAXN],c[MAXN],tmp[MAXN],inv2,inv3,invm,pw[MAXN],ipw[MAXN],iv[MAXN];
inline ll qexp(ll a,ll b)
{
ll res=1;
while(b)
{
if(b&1)res=res*a%Mod;
a=a*a%Mod;
b>>=1;
}
return res;
}
inline void init(int n)
{
inv2=qexp(2,Mod-2);inv3=qexp(3,Mod-2);
for(register int l=1;l<=(n<<2);l<<=1)
pw[l]=qexp(3,(Mod-1)/l),ipw[l]=qexp(inv3,(Mod-1)/l);
REP(i,1,n+n+1)iv[i]=qexp(i,Mod-2);
}
inline void NTT(ll *A,int n,int tp)
{
REP(i,0,n-1)if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(register int l=2;l<=n;l<<=1)
{
ll wn=tp>0?pw[l]:ipw[l];
for(register int i=0;i<n;i+=l)
{
ll w=1;
for(register int j=0;j<(l>>1);++j)
{
ll A1=A[i+j],A2=1ll*A[i+j+(l>>1)]*w%Mod;
A[i+j]=(A1+A2)%Mod,A[i+j+(l>>1)]=(A1-A2+Mod)%Mod;
w=1ll*w*wn%Mod;
}
}
}
}
inline void Get_Mul(int n,ll *A,ll *B,ll *C)
{
int m,cnt;
for(m=1,cnt=0;m<n;m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
NTT(A,m,1);NTT(B,m,1);
REP(i,0,m-1)A[i]=1ll*A[i]*B[i]%Mod;
NTT(A,m,-1);
invm=qexp(m,Mod-2);
REP(i,0,n-1)C[i]=1ll*A[i]*invm%Mod;
}
inline void Get_Inv(int len,ll *A,ll *B)
{
if(len==1)
{
B[0]=qexp(A[0],Mod-2);
return ;
}
Get_Inv((len+1)>>1,A,B);
int m,cnt;
for(m=1,cnt=0;m<(len<<1);m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
REP(i,0,len-1)tmp[i]=A[i];
REP(i,len,m-1)tmp[i]=0;
NTT(tmp,m,1);NTT(B,m,1);
REP(i,0,m-1)B[i]=(2ll*B[i]%Mod-1ll*B[i]*B[i]%Mod*tmp[i]%Mod+Mod)%Mod;
NTT(B,m,-1);
invm=qexp(m,Mod-2);
REP(i,0,len-1)B[i]=1ll*B[i]*invm%Mod;
REP(i,len,m-1)B[i]=0;
}
inline void Get_Sqr(int len,ll *A,ll *B)
{
if(len==1)
{
B[0]=1;
return ;
}
Get_Sqr((len+1)>>1,A,B);
memset(b,0,sizeof(b));Get_Inv(len,B,b);
int m,cnt;
for(m=1,cnt=0;m<(len<<1);m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
REP(i,0,len-1)tmp[i]=A[i];
REP(i,len,m-1)tmp[i]=0;
NTT(B,m,1);NTT(tmp,m,1);NTT(b,m,1);
REP(i,0,m-1)B[i]=1ll*(tmp[i]+1ll*B[i]*B[i]%Mod)%Mod*b[i]%Mod*inv2%Mod;
NTT(B,m,-1);
invm=qexp(m,Mod-2);
REP(i,0,len-1)B[i]=1ll*B[i]*invm%Mod;
REP(i,len,m-1)B[i]=0;
}
inline void Get_Div(int n1,ll *A,int n2,ll *B,ll *C,ll *D)
{
memset(c,0,sizeof(c));
REP(i,0,n2-1)c[i]=B[n2-1-i];
Get_Inv(n1-n2+1,c,b);
int m,cnt;
for(m=1,cnt=0;m<((n1-n2+1)<<1);m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
REP(i,0,n1-n2)tmp[i]=A[n1-1-i];
REP(i,n1-n2+1,m-1)tmp[i]=0;
NTT(tmp,m,1);NTT(b,m,1);
REP(i,0,m-1)tmp[i]=1ll*tmp[i]*b[i]%Mod;
NTT(tmp,m,-1);
invm=qexp(m,Mod-2);
memset(c,0,sizeof(c));
REP(i,0,n1-n2)C[n1-n2-i]=c[n1-n2-i]=1ll*tmp[i]*invm%Mod;
for(m=1,cnt=0;m<n1;m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
NTT(B,m,1);NTT(c,m,1);
REP(i,0,m-1)B[i]=1ll*B[i]*c[i]%Mod;
NTT(B,m,-1);
invm=qexp(m,Mod-2);
REP(i,0,n2-2)D[i]=(A[i]-1ll*B[i]*invm%Mod+Mod)%Mod;
}
inline void Get_Der(int n,ll *A,ll *B)
{
REP(i,0,n-2)B[i]=1ll*A[i+1]*(i+1)%Mod;B[n-1]=0;
}
inline void Get_Int(int n,ll *A,ll *B)
{
REP(i,1,n-1)B[i]=1ll*A[i-1]*iv[i]%Mod;B[0]=0;
}
inline void Get_Ln(int n,ll *A,ll *B)
{
memset(b,0,sizeof(b));
Get_Inv(n,A,b);Get_Der(n,A,c);
int m,cnt;
for(m=1,cnt=0;m<(n<<1);m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
NTT(b,m,1);NTT(c,m,1);
REP(i,0,m-1)b[i]=1ll*b[i]*c[i]%Mod;
NTT(b,m,-1);
invm=qexp(m,Mod-2);
REP(i,0,n-1)b[i]=1ll*b[i]*invm%Mod;
Get_Int(n,b,B);
}
inline void Get_Exp(int len,ll *A,ll *B)
{
if(len==1)
{
B[0]=1;
return ;
}
Get_Exp((len+1)>>1,A,B);Get_Ln(len,B,c);
int m,cnt;
for(m=1,cnt=0;m<(len<<1);m<<=1,++cnt);
REP(i,0,m-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
REP(i,0,len-1)tmp[i]=(A[i]-c[i]+(i==0)+Mod)%Mod;
REP(i,len,m-1)tmp[i]=0;
NTT(tmp,m,1);NTT(B,m,1);
REP(i,0,m-1)B[i]=1ll*tmp[i]*B[i]%Mod;
NTT(B,m,-1);
invm=qexp(m,Mod-2);
REP(i,0,len-1)B[i]=1ll*B[i]*invm%Mod;
REP(i,len,m-1)B[i]=0;
}
}