多项式运算封装
update on 21.12.30:添加了 polyeva;修补了 polymod 处多测时可能产生的 bug。
update on 22.2.7:重 写(前一版太丑了),改为完全封装版本(使用 std::vector 存放多项式系数,运算在命名空间 polynomial:: 里)
update on 22.2.11 第二版增加了快速求值插值(但是使用的是利用多项式取模进行降次的做法,常数巨大,关于转置原理的做法见作者在 Luogu 5050 的提交记录),
实现的并不优秀,但应该很稳(雾)。全部提交测试过,部分由 Vincra 这个账号提交。
索引:
Ver 2(施工中)
使用 poly A; 来声明一个名为 \(A\) 的多项式,其原型为 vector,较 Ver 1 更为优美、精简,去掉了没啥用、并且很丑的常数项不为 \(1\) 时的多项式快速幂。
下面是 P5158 的代码:
const int djq=998244353;
inline int ksm(int base,int p){
int ret=1;
while(p){
if(p&1) ret=1ll*ret*base%djq;
base=1ll*base*base%djq,p>>=1;
}
return ret;
}
namespace polynomial{
typedef vector poly;
const int G=3,invG=ksm(G,djq-2),inv2=ksm(2,djq-2);
int rev[1200005];
inline int BSGS(int a,int b){
map mp; b%=djq;
int t=ceil(sqrt((double)djq)),epow=1;
for(rg int j=0;j=0) return 1ll*i*t-mp[epow];
return b;
}
inline int modsqrt(int x){ return (!x||x==1)?x:ksm(3,BSGS(3,x)>>1); }
inline int initrev(const int n){
int len=1,lgn=0;
while(len<=n) len<<=1,++lgn;
for(rg int i=0;i>1]>>1)|((i&1)<<(lgn-1));
return len;
}
inline void NTT(poly &A,const int len,const int opt){
A.resize(len);
for(rg int i=0;i>1);
int len=initrev(n<<1);
NTT(A,len,1),NTT(B,len,1);
for(rg int i=0;i>1); A.resize(n);
poly B=ln(A,n);
for(rg int i=0;i>1),C=inv(B,n);
const int len=initrev(n<<1);
NTT(A,len,1),NTT(C,len,1);
for(rg int i=0;i>1;
evabuild(x<<1,l,mid),evabuild(x<<1|1,mid+1,r),p[x]=mul(p[x<<1],p[x<<1|1]);
}
void evasolve(int x,int l,int r,poly A){
if(l==r) return ey[l]=A[0],void();
const int mid=l+r>>1;
poly tmp=mod(A,p[x<<1]); evasolve(x<<1,l,mid,tmp);
tmp=mod(A,p[x<<1|1]); evasolve(x<<1|1,mid+1,r,tmp);
}
poly g;
int n,x[100005],y[100005];
poly intsolve(int x,int l,int r){
if(l==r) return poly(1,1ll*y[l]*ksm(ey[l],djq-2)%djq);
const int mid=l+r>>1;
poly tmp1=mul(intsolve(x<<1,l,mid),p[x<<1|1]),tmp2=mul(intsolve(x<<1|1,mid+1,r),p[x<<1]);
tmp1.resize(r-l+1),tmp2.resize(r-l+1);
for(rg int i=0;i
Ver 1
变量/常量:
N:多项式次数。djq G invG:模数、该模数的一个原根、该原根的逆元。ar1~ar5:不知道起什么名字的辅助数组。ivG~imG:看名字大概能猜出来什么用的辅助数组。rev1 rev2:蝴蝶变换数组。
函数( pulic ):
若无特殊说明,整数、多项式系数均为 \(\pmod {p}\) 意义下。
-
ksm(a,b):求 \(a^b\);快速幂;\(O(\log b)\)。 -
modsqrt(a):求 \(a^{\frac{1}{2}}\);原根 + BSGS;\(O(\sqrt {p})\)。 -
initrev(n,lgn,l,rev):预处理\(\pmod{x^n}\) 意义下 NTT 所用的蝴蝶变换数组;暴力;\(O(n)\)。 -
NTT(A,n,opt,rev):以 \(rev\) 作为蝴蝶变换数组,对 \(A(x)\pmod{x^n}\) 做快速数论变换(\(n^{-1}\) 放在外面乘);NTT;\(O(n\log n)\)。 -
polymul(A,B,C,n):求 \(C(x)=A(x)B(x)\pmod{x^{2n}}\);NTT;\(O(n\log n)\)。 -
polyinv(A,B,n):求 \(B(x)=A^{-1}(x)\pmod{x^n}\);牛顿迭代 + NTT;\(O(n\log n)\)。 -
polydif(A,B,n):求 \(B(x)=\dfrac{\mathrm{d}A(x)}{\mathrm{d}x}\pmod{x^n}\);微分;\(O(n)\)。 -
polyint(A,B,n):求 \(B(x)=\int A(x)\mathrm{d}x\pmod{x^n}\);积分;\(O(n)\)。 -
polysqrt(A,B,n):求 \(B(x)=A^{\frac{1}{2}}(x)\pmod{x^n}\)(对\(A_0\) 无要求);牛顿迭代 + NTT;\(O(n\log n)\)。 -
polyln(A,B,n):求 \(B(x)=\ln A(x)\pmod{x^n}\)(由麦克劳林级数定义);复合函数微分 + 积分 + NTT;\(O(n\log n)\)。 -
polyexp(A,B,n):求 \(B(x)=\exp A(x)\pmod{x^n}\)(由麦克劳林级数定义);牛顿迭代 +polyln;\(O(n\log n)\)。 -
polymod(A,B,C,n,m,D):有 \(A(x)\pmod{x^n}\) 和 \(B(x)\pmod{x^m}\) (\(m),求 \(D(x)=A(x)\bmod{B(x)}\pmod{x^{m-1}}\),\(C(x)=B^{-1}(x)(A(x)-D(x))\pmod{x^{n-m+1}}\);暴力 + polyinv;\(O(n\log n)\) -
polyksm(A,B,k,n):求 \(B(x)=A^k(x)\pmod{x^n}\)(要求 \(A_0=1\));初中数学 +polyln+polyexp;\(O(n\log n)\)。 -
sppolyksm(A,B,k1,k2,lim,n):求 \(B(x)=A^k(x)\pmod{x^n}\)(对\(A_0\) 无要求),其中 \(k1=k\pmod{p}\) , \(k1=k\pmod{\varphi(p)}\) , \(lim=\min(k,p)\)(用于特判);暴力 +polyksm;\(O(n\log n)\)。 -
polycdq(A,B,n,val):求 \(B_0=val,B_i=\sum_{j=0}^{i-1}A_{i-j}B_j\pmod{x^n}\) (最基础的半在线卷积);cdq 分治 + NTT;\(O(n\log^2 n)\)。 -
polyeva(A,B,n,m):给出 \(A(x)\pmod{x^n}\) 和 \(x_1\ldots x_m\),输出 \(A(x_1)\ldots A(x_m)\)。(多点快速求值)
const int N=100005,djq=998244353;
class poly{
private:
const int G=3,invG=332748118;
ll ar1[4*N],ar2[4*N],ar3[4*N],ar4[4*N],ar5[4*N],ar6[4*N];
ll ivG[4*N],sqG[4*N],lnG[4*N],sqtmp[4*N];
ll expG[4*N],modG[4*N],kmG[4*N];
int rev1[4*N],rev2[4*N];
void POLYINV(ll* A,ll* B,const int n){
if(n==1) return B[0]=ksm(A[0],djq-2),void();
POLYINV(A,B,(n+1)>>1);
int l,lgn; initrev(n<<1,lgn,l,rev2);
copy(A,A+n,ivG);
fill(ivG+n,ivG+l,0);
NTT(B,l,1,rev2),NTT(ivG,l,1,rev2);
for(rg int i=0;i mp; b%=djq;
int t=ceil(sqrt((double)djq)),epow=1;
for(rg int j=0;j=0)
return 1ll*i*t-mp[epow];
return b;
}
void CDQ(ll* F,ll* G,int l,int r){
if(l==r) return;
const int mid=l+r>>1;
CDQ(F,G,l,mid);
copy(F+l,F+mid+1,ar4),copy(G,G+r-l+1,ar5);
fill(ar4+mid-l+1,ar4+(r-l+1)*4,0),fill(ar5+r-l+1,ar5+(r-l+1)*4,0);
polymul(ar4,ar5,ar4,r-l+1);
for(rg int i=mid+1;i<=r;++i) F[i]=(F[i]+ar4[i-l])%djq;
CDQ(F,G,mid+1,r);
}
ll *p[4*N]; int len[4*N];
void up(int x){
len[x]=len[x<<1]+len[x<<1|1],p[x]=new ll[len[x]+1];
int l,lgn; initrev(len[x]+1,lgn,l,rev1);
fill(ar1,ar1+l,0),fill(ar2,ar2+l,0);
copy(p[x<<1],p[x<<1]+len[x<<1]+1,ar1),copy(p[x<<1|1],p[x<<1|1]+len[x<<1|1]+1,ar2);
NTT(ar1,l,1,rev1),NTT(ar2,l,1,rev1);
for(rg int j=0;j>1;
build(x<<1,l,mid,A);
build(x<<1|1,mid+1,r,A);
up(x);
}
void solve(int x,int l,int r,ll* A){
if(l==r) return printf("%d\n",(int)A[0]),void();
const int mid=l+r>>1;
ll tmp[len[x]*2+5];
polymod(A,p[x<<1],tmp,len[x],len[x<<1]+1,ar6);
solve(x<<1,l,mid,tmp);
polymod(A,p[x<<1|1],tmp,len[x],len[x<<1|1]+1,ar6);
solve(x<<1|1,mid+1,r,tmp);
}
public:
inline int ksm(ll base,int p){
ll ret=1;
while(p){
if(p&1) ret*=base,ret%=djq;
base*=base,base%=djq; p>>=1;
} return ret;
}
inline int modsqrt(int x){
if(!x||x==1) return x;
return ksm(3,BSGS(3,x)>>1);
}
inline void initrev(const int n,int& lgn,int& l,int* rev){
lgn=0,l=0;
while((1<>1]>>1)|((j&1)<<(lgn-1));
}
void NTT(ll* A,const int l,const int opt,int* rev){
for(rg int i=0;i>1);
polyinv(B,sqG,n);
for(rg int j=0;j>1);
polyln(B,expG,n);
for(rg int i=0;i=1ll*n) return fill(B,B+n,0),void();
for(rg int i=0;i=tmp;--i) B[i]=B[i-tmp];
fill(B,B+tmp,0);
}
void polycdq(ll* A,ll* B,const int n,const int va){
fill(A+n,A+4*n,0),fill(B,B+4*n,0),B[0]=va;
CDQ(B,A,0,n-1);
}
void polyeva(ll* A,ll* x,const int n,const int m){
build(1,1,m,x);
if(n>m) polymod(A,p[1],A,n,m+1,ar6);
solve(1,1,m,A);
for(rg int i=1;i<=4*m;++i) if(p[i]!=NULL) delete p[i],p[i]=NULL;
}
}work;