多项式全家桶(纯模板)
已经完成的：
- 多项式乘法
- 多项式乘法逆
- 多项式\(\ln\)
- 多项式\(\exp\)
- 多项式开根(\(a_0\)为二次剩余)
- 多项式快速幂(\(a_0\)可以不为1)
- 多项式三角函数(\(\sin,\cos\))
- 多项式除法、取模
- 分治\(FFT\)
- 任意模数\(NTT\)(\(MTT\))
- 多项式多点求值
- 多项式快速插值
- 常系数齐次线性递推
- BM
还未完成的
- 多项式非齐次线性递推
- 普通幂转下降幂/下降幂转普通幂
- 多项式复合逆
总结
一般多项式题，都是通过一些奇妙的操作（DP、生成函数）后化成一个能用多项式快速完成的东西
一些常见套路：
\[f_n=\sum_{i=1}^{n}g_ih_{n-i}
\]
就是卷积的基础柿子
\[f_n=\sum_{i=1}^{n-1}f_ig_{n-i}
\]
分治 \(FFT\)：直接 \(cdq\) 分治，每次计算 \(f_{l\sim mid}\) 与 \(g_{1\sim r-l}\) 卷积对 \(f_{mid+1\sim r}\) 的贡献，复杂度 \(\mathcal O(n\log^2 n)\)。（当然可以多项式求逆，但分治 \(FFT\) 可用于更广泛的情况）
特殊情况：自己卷自己：
\[f_n=\sum_{i=1}^{n-1}g_ih_{n-i}
\]
其中 \(g,h\) 都是与 \(f\) 有关的函数，需要先求出 \(f_x\) 才能得到 \(g_x,h_x\)。
还是考虑之前的做法，但在 \(solve(l,r)\) 时，可能会出现 \(r-l<mid\) 的情况导致需要的 \(h\) 没有被算出来。注意到情况只会在 \(l=1\) 时出现，因此在 \(l=1\) 时，先只计算 \(g_{1\sim mid}\) 与 \(h_{1\sim mid}\) 卷积对 \(f_{mid+1\sim r}\) 的贡献。在 \(l\not=1\) 时，在上一条的卷积之外，额外增加 \(g_{1\sim r-l}\) 与 \(h_{l\sim mid}\) 卷积对 \(f_{mid+1\sim r}\) 的贡献，弥补之前缺失的部分即可。
减法卷积：
\[f_n=\sum_{i=1}^{d}g_ih_{n+i}
\]
将\(g(x)\)\(reverse\)一下，于是就变成了：
\[f_n=\sum_{i=1}^{d}g_{d-i}h_{n+i}
\]
于是又可以卷积了，设\(g*h=A\),\(f_n=A_{n+d}\)
\[f_n=\sum_{i=1}^{n}\binom nig_{i}h_{n-i}
\]
设\(F,G,H\)分别是\(f(x),g(x),h(x)\)的\(EGF\),于是有\(F=GH\)
5.指数公式定理：
如果存在两个EGF \(F(x)\) 与 \(G(x)\)，满足\(e^{F(x)}=G(x)\)，\(F(x)\) 是 \(f_i\) 的 EGF，那么\(G(x)\) 是
\[g_n=\sum_{\pi=\{S_1,S_2,\dots,S_n\}}\prod_{i=1}^n f_{|S_i|}
\]
的EGF，其中\(\pi\)是 \([n]\) 的划分。
Update on 2022.3.22 Berlekamp–Massey 算法
BM求解最短递推式
```cpp namespace BM{ vectorUpdate on 2022.2.26 常系数齐次线性递推
%x^n\bmod f(x)$
```cpp int main(){ // freopen("poly.in","r",stdin); // freopen("1.out","w",stdout); n=Rint;k=Rint; init_poly((k+1)<<1); for(int i=1;i<=k;++i) f[k-i]=-Rint,f[k-i]=(f[k-i]%mod+mod)%mod; f[k]=1;GINV=f.rev().inv(); poly x,ret; x[1]=ret[0]=1; for(;n;n>>=1,x=(x*x)%f) if(n&1) ret=(ret*x)%f; int ans=0; for(int i=0,x;i波斯坦茉莉算法
```cpp inline int solve(poly f,poly g,ll n){ for(;n;n>>=1){ poly h=g; for(int i=1;iUpdate on 2021.3.29 数组版本新鲜出炉，跑的比较快但难用的一匹。
Update on 2021.3.26 多点求值与快速插值
多点求值与快速插值
#define lc (p<<1)
#define rc (p<<1|1)
#define mid ((l+r)>>1)
namespace Multiask{
int a[N],n,m,ans[N];
poly T[100010<<2],F[100010<<2],f;
inline void build(int p,int l,int r){
if(l==r){T[p][1]=mod-a[l];T[p][0]=1;return ;}
build(lc,l,mid);build(rc,mid+1,r);
T[p]=T[lc]*T[rc];
}
inline poly mulT(const poly &F,const poly &G,int siz,int tp){
poly f=F,g=G;
if(f.size()<=100){
poly ret;
for(int i=0;i<f.size();++i){
int fr=0;
if(tp==1) fr=max(fr,i-(int)(f.size()-g.size()));
for(int j=fr;j<=i&&j<g.size();++j) inc(ret[i-j],1ll*f[i]*g[j]%mod);
}
return ret;
}
int len=lg[f.size()*tp],lim=1<<len+1,gg=g.size();
init(lim);
reverse(g.v.begin(),g.v.end());
NTT(f,lim,1);NTT(g,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
NTT(f,lim,0);
return vec(f.v.begin()+gg-1,f.v.begin()+gg+siz-1);
}
inline void ask(int p,int l,int r){
// cout<<l<<" "<<r<<" ";F[p].print(r-l+1);
if(l==r){
ans[l]=F[p][0];
return ;
}
F[lc]=mulT(F[p],T[rc],mid-l+1,1);F[rc]=mulT(F[p],T[lc],r-mid,1);
ask(lc,l,mid);
ask(rc,mid+1,r);
}
inline void query(int n){
build(1,1,n);
T[1]=T[1].inv();T[1].resize(n);f.resize(n);
F[1]=mulT(f,T[1],n,2);
ask(1,1,n);
}
}
namespace lagrange{
int x[N],y[N];
poly g[100010<<2];
inline void build(int p,int l,int r){
if(l==r){g[p][0]=dec(0,x[l]);g[p][1]=1;return ;}
build(lc,l,mid);build(rc,mid+1,r);
g[p]=g[lc]*g[rc];
}
inline poly query(int p,int l,int r){
if(l==r) return 1ll*y[l]*ksm(Multiask::ans[l],mod-2)%mod;
return (query(lc,l,mid)*g[rc])+(query(rc,mid+1,r)*g[lc]);
}
inline void work(int n){
build(1,1,n);
Multiask::f=g[1].der();
for(int i=1;i<=n;++i) Multiask::a[i]=x[i];
Multiask::query(n);
query(1,1,n).print(n);
}
}
#undef lc
#undef rc
#undef mid
Update on 2021.3.25 出现了第二代多项式板子，方便使用且常数较优，但考场上基本上写不出来
Update on 2021.1.13 大幅优化了下方模板的常数
Update on 2021.1.12
任意模数NTT
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2e7+10;
const int M=6;
int n,v,mod,q;
int k[M][M],b[M][M];
//--------iobuff--------
namespace iobuff{
const int LEN=1000000;
char in[LEN+5], out[LEN+5];
char *pin=in, *pout=out, *ed=in, *eout=out+LEN;
inline char gc(void)
{
#ifdef LOCAL
return getchar();
#endif
return pin==ed&&(ed=(pin=in)+fread(in, 1, LEN, stdin), ed==in)?EOF:*pin++;
}
inline void pc(char c)
{
pout==eout&&(fwrite(out, 1, LEN, stdout), pout=out);
(*pout++)=c;
}
inline void flush()
{ fwrite(out, 1, pout-out, stdout), pout=out; }
template<typename T> inline void read(T &x)
{
static int f;
static char c;
c=gc(), f=1, x=0;
while(c<'0'||c>'9') f=(c=='-'?-1:1), c=gc();
while(c>='0'&&c<='9') x=10*x+c-'0', c=gc();
x*=f;
}
template<typename T> inline void putint(T x, char div)
{
static char s[15];
static int top;
top=0;
x<0?pc('-'), x=-x:0;
while(x) s[top++]=x%10, x/=10;
!top?pc('0'), 0:0;
while(top--) pc(s[top]+'0');
pc(div);
}
}
using namespace iobuff;
//---------math----------
int base[N],inv[N],jc[N];
inline ll add(ll x,ll y,ll m=mod){return (x+y>=m)?x+y-m:x+y;}
inline ll dec(ll x,ll y,ll m=mod){return (x-y<0)?x-y+m:x-y;}
inline void init(){
base[0]=base[1]=inv[0]=inv[1]=jc[0]=jc[1]=1;
for(int i=2;i<N;++i){
base[i]=1ll*base[i-1]*i%mod;
inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
jc[i]=1ll*jc[i-1]*inv[i]%mod;
}
}
inline int binom(int x,int y){
return 1ll*base[x]*jc[y]%mod*jc[x-y]%mod;
}
inline int calc(int l,int r,int j,int u){
return dec(binom(r+j+1,j+u+1),binom(l+j,j+u+1));
}
inline ll ksm(ll x,ll y,ll mod){
ll ret=1ll;
for(;y;y>>=1,x=x*x%mod) (y&1)&&(ret=ret*x%mod);
return ret;
}
//--------poly----------
const ll mod1=998244353;
const ll mod2=1004535809;
const ll mod3=469762049;
const ll mod12=mod1*mod2;
const ll iv1=669690699;
const ll iv2=354521948;
const ll gg=3;
namespace Container{
ll Wn[2][N];
int lg[N];
struct poly{
vector<ll>v;
inline ll& operator[](int x){while(x>=v.size())v.push_back(0);return v[x];}
inline poly(int x=0):v(1){v[0]=x;}
inline int size(){return v.size();}
inline void resize(int x){v.resize(x);}
inline void mem(int l,int lim){
int t=min(lim,(int)v.size());
for(int i=l;i<t;++i) v[i]=0;
}
};
int p1,p2;
inline void init_poly(int n,ll mod){
p1=1;
while((p1<<1)<=n) p1<<=1;p2=p1<<1;
ll giv=(mod+1)/3;
ll wn1=ksm(3,(mod-1)/p2,mod),wn2=ksm(giv,(mod-1)/p2,mod);
Wn[0][p1]=Wn[1][p1]=1;
for(int i=p1+1;i<=p2;++i)
Wn[0][i]=Wn[0][i-1]*wn1%mod,Wn[1][i]=Wn[1][i-1]*wn2%mod;
for(int i=p1-1;i>=1;--i) Wn[0][i]=Wn[0][i<<1],Wn[1][i]=Wn[1][i<<1];
for(int i=2;i<n;++i)
lg[i]=lg[i>>1]+1;
}
}
using namespace Container;
namespace Poly{
int r[N];
ll fr[N];
inline void NTT(int lim,poly& f,int tp,ll mod){
int kk=tp==1?0:1;
for(int i=0;i<lim;++i) fr[i]=f[r[i]];
for(int mid=1;mid<lim;mid<<=1){
int len=mid<<1;
for(int l=0;l+len-1<lim;l+=len){
for(int k=l;k<=l+mid-1;++k){
ll w1=fr[k],w2=Wn[kk][k-l+mid]*fr[k+mid]%mod;
fr[k]=add(w1,w2,mod);fr[k+mid]=dec(w1,w2,mod);
}
}
}
if(tp==-1){
ll t=ksm(lim,mod-2,mod);
for(int i=0;i<lim;++i) fr[i]=fr[i]*t%mod;
}
for(int i=0;i<lim;++i) f[i]=fr[i];
}
inline ll merge(ll x1,ll x2,ll x3){
ll k1=(x2+mod2-x1)%mod2*iv1%mod2*mod1+x1;
return (((x3-k1%mod3+mod3)%mod3*iv2%mod3*(mod12%mod))%mod+k1)%mod;
}
inline poly poly_mul(int n,int m,poly f,poly g){
int lim=1,len=0;
while(lim<(n+m)) lim<<=1,len++;
for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(len-1));
poly f1=f,f2=f,f3=f,g1=g,g2=g,g3=g;
init_poly(lim,mod1);
NTT(lim,f1,1,mod1);NTT(lim,g1,1,mod1);
for(int i=0;i<lim;++i) f1[i]=f1[i]*g1[i]%mod1;
NTT(lim,f1,-1,mod1);
init_poly(lim,mod2);
NTT(lim,f2,1,mod2);NTT(lim,g2,1,mod2);
for(int i=0;i<lim;++i) f2[i]=f2[i]*g2[i]%mod2;
NTT(lim,f2,-1,mod2);
init_poly(lim,mod3);
NTT(lim,f3,1,mod3);NTT(lim,g3,1,mod3);
for(int i=0;i<lim;++i) f3[i]=f3[i]*g3[i]%mod3;
NTT(lim,f3,-1,mod3);
for(int i=0;i<lim;++i) f[i]=merge(f1[i],f2[i],f3[i]);
return f;
}
}
using namespace Poly;
//----------main------------
poly f,g;
int m;
int main(){
read(n);read(m);read(mod);++n;++m;
for(int i=0;i<n;++i) read(f[i]);
for(int i=0;i<m;++i) read(g[i]);
f=poly_mul(n,m,f,g);
for(int i=0;i<n+m-1;++i) putint(f[i],' ');flush();
return 0;
}
Update on ???
- 多项式求牛顿级数，也算是个比较有用的东西吧：
求\(g(x)\)的牛顿级数\(c_k\)：
\[c_k=\triangle^kg(0)=\sum_k\binom{n}{k}(-1)^{n-k}g(k)
\]
于是令\(f(x)=(-1)^x\)，故\(C=F*G\)，其中\(C,F,G\)都是\(EGF\)
多项式牛顿级数
inline int mi(int x){return x&1?mod-1:1;}
inline poly Newton(int n,poly g){
poly c;
for(int i=0;i<n;++i)
g[i]=1ll*g[i]*jc[i]%mod;
for(int i=0;i<n;++i)
c[i]=1ll*mi(i)*jc[i]%mod;
c=poly_mul(n,n,c,g);
for(int i=0;i<n;++i)
c[i]=1ll*c[i]*base[i]%mod;
return c;
}//这里我们假设已经求出了$g(i)$的值，存在$g$中
第二代模板
view code
#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
const int N=(1<<21)+20;
const int mod=998244353;
typedef vector<int> vec;
inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
inline void inc(int &x,int y){x=add(x,y);}
inline void rec(int &x,int y){x=dec(x,y);}
inline int ksm(int x,int y){
int ret=1;
for(;y;y>>=1,x=1ll*x*x%mod) if(y&1) ret=1ll*ret*x%mod;
return ret;
}
namespace IO {
inline char nc(){
// return getchar();
static char buf[500005],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,500000,stdin),p1==p2)?EOF:*p1++;
}
char out[500005],*pout=out,*eout=out+500000;
inline void flush() { fwrite(out,1,pout-out,stdout),pout=out; }
inline void pc(char c) { pout==eout&&(fwrite(out,1,500000,stdout),pout=out); (*pout++)=c; }
template<typename T> inline void put(T x,char suf) {
static char stk[65];int top=0; while(x) stk[top++]=x%10,x/=10;
!top?pc('0'),0:0; while(top--) pc(stk[top]+'0'); pc(suf);
}
inline int read(){
char ch=nc(); int sum=0; for(;ch<'0'||ch>'9';ch=nc());
while(ch>='0'&&ch<='9')sum=(sum<<1)+(sum<<3)+ch-48,ch=nc();
return sum;
}
}
#define Rint IO::read()
using IO::put;
using IO::nc;
namespace Prep{
int Wn[N<<1],lg[N],p2,mx=1,r[N],tot;
inline void init_poly(int n){
int p=1;while(p<=n)p<<=1;
for(int i=2;i<=p;++i) lg[i]=lg[i>>1]+1;
for(int i=1;i<p;i<<=1){
int wn=ksm(3,(mod-1)/(i<<1));
Wn[++tot]=1;
for(int j=1;j<i;++j) ++tot,Wn[tot]=1ll*Wn[tot-1]*wn%mod;
}
}
inline void init(int lim){
// cout<<"init"<<lim<<endl;
int len=lg[lim]-1;
for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<len);
}
int iv[N],tp;
inline void init_inv(int n){
if(!tp){tp=2;iv[0]=iv[1]=1;}
for(;tp<=n;++tp) iv[tp]=1ll*(mod-mod/tp)*iv[mod%tp]%mod;
}
}
using namespace Prep;
namespace Cipolla{
int I,fl=0;
mt19937 rnd(time(0));
struct pt{
int a,b;
pt(int _a=0,int _b=0){a=_a;b=_b;}
};
inline pt operator *(pt x,pt y){
pt ret;
ret.a=add(1ll*x.a*y.a%mod,1ll*x.b*y.b%mod*I%mod);
ret.b=add(1ll*x.a*y.b%mod,1ll*x.b*y.a%mod);
return ret;
}
inline bool check(int x){
return ksm(x,(mod-1)/2)==1;
}
inline int random(){return rnd()%mod;}
inline pt qpow(pt a,int b){
pt ret=pt(1,0);
for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
return ret;
}
inline int cipolla(int n){
if(!check(n)) return 0;
int a=random();
while(!a||check(dec(1ll*a*a%mod,n))) a=random();
I=dec(1ll*a*a%mod,n);
int ans=qpow(pt(a,1),(mod+1)/2).a;
return min(ans,(int)mod-ans);
}
}
using namespace Cipolla;
int ddddd=0;
bool flag=0;
namespace Poly{
struct poly{
vec v;
inline poly(int w=0):v(1){v[0]=w;}
inline poly(const vec&w):v(w){}
inline int operator [](int x)const{return x>=v.size()?0:v[x];}
inline int& operator [](int x){
if(x>=v.size()) v.resize(x+1);
return v[x];
}
inline int size(){return v.size();}
inline void resize(int x){v.resize(x);}
inline void read(int n){v.resize(n);for(int i=0;i<n;++i) v[i]=Rint;}
inline void print(int n)const{for(int i=0;i<n-1;++i) put(operator[](i),' ');put(operator[](n-1),'\n');}
inline poly slice(int len)const{
if(len<=v.size()) return vec(v.begin(),v.begin()+len);
vec ret(v);ret.resize(len);
return ret;
}
inline poly operator *(const int &x)const{
poly ret(v);
for(int i=0;i<v.size();++i) ret[i]=1ll*ret[i]*x%mod;
return ret;
}
inline poly operator -()const{
poly ret(v);
for(int i=0;i<v.size();++i) ret[i]=dec(0,ret[i]);
return ret;
}
inline poly operator*(const poly &g)const;
inline poly operator/(const poly &g)const;
inline poly operator%(const poly &g)const;
inline poly der()const{
vec ret(v);
for(int i=0;i<ret.size()-1;++i) ret[i]=1ll*ret[i+1]*(i+1)%mod;
ret.pop_back();
return ret;
}
inline poly jifen()const{
vec ret(v);
init_inv(ret.size());ret.push_back(0);
for(int i=ret.size()-1;i;--i) ret[i]=1ll*ret[i-1]*iv[i]%mod;ret[0]=0;
return ret;
}
inline poly rev()const{
vec ret(v);
reverse(ret.begin(),ret.end());
return ret;
}
inline poly inv()const;
inline poly div(const poly &FF)const;
inline poly ln()const;
inline poly exp()const;
inline poly pow(int k)const;
inline poly sqrt()const;
inline poly mulT(const poly &g,int siz,int tp)const;
};
inline poly operator +(const poly &x,const poly &y){
vec v(max(x.v.size(),y.v.size()));
for(int i=0;i<v.size();++i) v[i]=add(x[i],y[i]);
return v;
}
inline poly operator -(const poly &x,const poly &y){
vec v(max(x.v.size(),y.v.size()));
for(int i=0;i<v.size();++i) v[i]=dec(x[i],y[i]);
return v;
}
ull fr[N];
const ull Mod=998244353;
inline void NTT(poly& f,int lim,int tp){
// ddddd+=lim*__builtin_ctz(lim);
for(int i=0;i<lim;++i) fr[i]=f[r[i]];
for(int mid=1;mid<lim;mid<<=1){
for(int len=mid<<1,l=0;l+len-1<lim;l+=len){
for(int k=l;k<l+mid;++k){
ull w1=fr[k],w2=fr[k+mid]*Wn[mid+k-l]%Mod;
fr[k]=w1+w2;fr[k+mid]=w1+Mod-w2;
}
}
}
for(int i=0;i<lim;++i) fr[i]>=Mod?fr[i]%=Mod:0;
if(!tp){
reverse(fr+1,fr+lim);
int iv=ksm(lim,mod-2);
for(int i=0;i<lim;++i) fr[i]=fr[i]*iv%mod;
}
for(int i=0;i<lim;++i) f[i]=fr[i];
}
inline poly poly::operator *(const poly &G)const{
poly f(v),g=G;
int rec=f.size()+g.size()-1,d=max(f.size(),g.size());
int len=lg[rec],lim=1<<len+1;
init(lim);
NTT(f,lim,1);NTT(g,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
NTT(f,lim,0);
return f.slice(rec);
}
inline poly poly::inv()const{
poly g,g0,d;
g[0]=ksm(v[0],mod-2);
for(int lim=2;(lim>>1)<v.size();lim<<=1){
g0=g;d=slice(lim);
init(lim);
NTT(g0,lim,1);NTT(d,lim,1);
for(int i=0;i<lim;++i) d[i]=1ll*g0[i]*d[i]%mod;
NTT(d,lim,0);
fill(d.v.begin(),d.v.begin()+(lim>>1),0);
NTT(d,lim,1);
for(int i=0;i<lim;++i) d[i]=1ll*d[i]*g0[i]%mod;
NTT(d,lim,0);
for(int i=lim>>1;i<lim;++i) g[i]=dec(g[i],d[i]);
}
return g.slice(v.size());
}//10E(n)
inline poly poly::div(const poly &FF)const{
if(v.size()==1) return 1ll*v[0]*ksm(FF[0],mod-2)%mod;
int len=lg[v.size()],lim=1<<len+1,nlim=lim>>1;
poly F=FF,G0=FF.slice(nlim);
G0=G0.inv();
poly H0=slice(nlim),Q0;
init(lim);
NTT(G0,lim,1);NTT(H0,lim,1);
for(int i=0;i<lim;++i) Q0[i]=1ll*G0[i]*H0[i]%mod;
NTT(Q0,lim,0);Q0.resize(nlim);
poly ret=Q0;
NTT(Q0,lim,1);NTT(F,lim,1);
for(int i=0;i<lim;++i) Q0[i]=1ll*Q0[i]*F[i]%mod;
NTT(Q0,lim,0);
fill(Q0.v.begin(),Q0.v.begin()+nlim,0);
for(int i=nlim;i<lim&&i<v.size();++i) Q0[i]=dec(Q0[i],v[i]);
NTT(Q0,lim,1);
for(int i=0;i<lim;++i) Q0[i]=1ll*Q0[i]*G0[i]%mod;
NTT(Q0,lim,0);
for(int i=nlim;i<lim;++i) ret[i]=dec(ret[i],Q0[i]);
return ret.slice(v.size());
}
inline poly poly::ln()const{
return der().div(*this).jifen();
}//13E
namespace EXP{
const int logB=4;
const int B=16;
poly f,ret,g[30][B];
inline void exp(int lim,int l,int r){
// cout<<lim<<" "<<l<<" "<<r<<endl;
if(r-l<=64){
for(int i=l;i<r;++i){
ret[i]=(!i)?1:1ll*ret[i]*iv[i]%mod;
for(int j=i+1;j<r;++j) inc(ret[j],1ll*ret[i]*f[j-i]%mod);
}
return ;
}
int k=(r-l)/B;
poly bl[B];
for(int i=0;i<B;++i) bl[i].resize(k<<1);
int len=1<<lim-logB+1;
for(int i=0;i<B;++i){
if(i>0){
init(len);NTT(bl[i],len,0);
for(int j=0;j<k;++j)
inc(ret[l+i*k+j],bl[i][j+k]);
}
exp(lim-logB,l+i*k,l+(i+1)*k);
if(i<B-1){
poly H;H.resize(k<<1);
for(int j=0;j<k;++j) H[j]=ret[j+l+i*k];
init(len);NTT(H,len,1);
for(int j=i+1;j<B;++j)
for(int t=0;t<(k<<1);++t) inc(bl[j][t],1ll*H[t]*g[lim][j-i-1][t]%mod);
}
}
}
inline void init_exp(){
ret.resize(f.size());
for(int i=0;i<f.size();++i) f[i]=1ll*f[i]*i%mod,ret[i]=0;
int mx=lg[f.size()]+1;
init_inv(1<<mx);
for(int lim=mx;lim>=logB;lim-=logB){
int bl=1<<(lim-logB),tot=0,ll=1<<(lim-logB+1);
init(ll);
for(int i=0;i<B-1;++i){
g[lim][i].resize(bl<<1);
for(int j=0;j<(bl<<1);++j) g[lim][i][j]=f[j+bl*i];
NTT(g[lim][i],ll,1);
}
}
}
}
inline poly poly::exp()const{
EXP::f=*this;
EXP::init_exp();
EXP::exp(lg[v.size()]+1,0,1<<lg[v.size()]+1);
return EXP::ret.slice(v.size());
}
inline poly poly::pow(int k)const{
return ((*this).ln()*k).exp();
}
inline poly poly::operator /(const poly &Q)const{
if(v.size()<Q.v.size()) return 0;
int p=v.size()-Q.v.size()+1;
poly fr=rev(),qr=Q.rev();
fr.resize(p);qr.resize(p);
return fr.div(qr).rev();
}
inline poly poly::operator %(const poly &Q)const{
poly F(v);
return (F-(Q*(F/Q))).slice(Q.v.size()-1);
}
inline poly poly::sqrt()const{
poly g,h,gf,F1,F2,F3,f(v);
g[0]=cipolla(operator[](0));
h[0]=ksm(g[0],mod-2);
gf[0]=g[0];gf[1]=g[0];
int iv=(mod+1)/2;
init(1);
for(int lim=1;lim<v.size();lim<<=1){
for(int i=0;i<lim;++i) F1[i]=1ll*gf[i]*gf[i]%mod;
NTT(F1,lim,0);
for(int i=0;i<lim;++i) F1[i+lim]=dec(F1[i],f[i]),F1[i]=0;
int nlim=lim<<1;init(nlim);
for(int i=lim;i<nlim;++i) rec(F1[i],f[i]);
F2=h;F2.resize(lim);
NTT(F1,nlim,1);NTT(F2,nlim,1);
for(int i=0;i<nlim;++i) F1[i]=1ll*F1[i]*F2[i]%mod;
NTT(F1,nlim,0);
for(int i=lim;i<nlim;++i) g[i]=dec(0,1ll*F1[i]*iv%mod);
if(nlim<v.size()){
gf=g;
NTT(gf,nlim,1);
for(int i=0;i<nlim;++i) F3[i]=1ll*gf[i]*F2[i]%mod;
NTT(F3,nlim,0);
fill(F3.v.begin(),F3.v.begin()+lim,0);
NTT(F3,nlim,1);
for(int i=0;i<nlim;++i) F3[i]=1ll*F3[i]*F2[i]%mod;
NTT(F3,nlim,0);
for(int i=lim;i<nlim;++i) rec(h[i],F3[i]);
}
}
return g.slice(v.size());
}
inline poly poly::mulT(const poly &G,int siz,int tp)const{
poly f(v),g=G;
if(f.size()<=100){
poly ret;
for(int i=0;i<f.size();++i){
int fr=0;
if(tp==1) fr=max(fr,i-(int)(f.size()-g.size()));
for(int j=fr;j<=i&&j<g.size();++j) inc(ret[i-j],1ll*f[i]*g[j]%mod);
}
return ret;
}
int len=lg[f.size()*tp],lim=1<<len+1,gg=g.size();
init(lim);
reverse(g.v.begin(),g.v.end());
NTT(f,lim,1);NTT(g,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
NTT(f,lim,0);
return vec(f.v.begin()+gg-1,f.v.begin()+gg+siz-1);
}
}
using namespace Poly;
inline int sread(){
long long m=mod,x=0;char ch=nc();
while(!isdigit(ch)){ch=nc();}
while(isdigit(ch)){x=(x<<3)+(x<<1)+(ch^48);x%=m;ch=nc();}
return x;
}
//remember to init_poly(多项式最高可能次数)
int n,k;poly f;
/*
int main(){
// freopen("poly2.in","r",stdin);
// freopen("1.out","w",stdout);
n=Rint+1,k=Rint;
f.read(n);
init_poly(n<<1);
poly g=f+poly(2)-poly(f[0]);
g=g-(f.sqrt().inv().jifen().exp());
(g.ln()+poly(1)).pow(k).der().print(n-1);
IO::flush();
return 0;
}
*/
int main(){
n=Rint;
f.read(n);
init_poly(n);
f.exp().print(n);
IO::flush();
return 0;
}
多点求值与快速插值
view code
#define lc (p<<1)
#define rc (p<<1|1)
#define mid ((l+r)>>1)
namespace Multiask{
int a[N],n,m,ans[N];
poly T[100010<<2],F[100010<<2],f;
inline void build(int p,int l,int r){
if(l==r){T[p][1]=mod-a[l];T[p][0]=1;return ;}
build(lc,l,mid);build(rc,mid+1,r);
T[p]=T[lc]*T[rc];
}
inline poly mulT(const poly &F,const poly &G,int siz,int tp){
poly f=F,g=G;
if(f.size()<=100){
poly ret;
for(int i=0;i<f.size();++i){
int fr=0;
if(tp==1) fr=max(fr,i-(int)(f.size()-g.size()));
for(int j=fr;j<=i&&j<g.size();++j) inc(ret[i-j],1ll*f[i]*g[j]%mod);
}
return ret;
}
int len=lg[f.size()*tp],lim=1<<len+1,gg=g.size();
init(lim);
reverse(g.v.begin(),g.v.end());
NTT(f,lim,1);NTT(g,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
NTT(f,lim,0);
return vec(f.v.begin()+gg-1,f.v.begin()+gg+siz-1);
}
inline void ask(int p,int l,int r){
// cout<<l<<" "<<r<<" ";F[p].print(r-l+1);
if(l==r){
ans[l]=F[p][0];
return ;
}
F[lc]=mulT(F[p],T[rc],mid-l+1,1);F[rc]=mulT(F[p],T[lc],r-mid,1);
ask(lc,l,mid);
ask(rc,mid+1,r);
}
inline void query(int n){
build(1,1,n);
T[1]=T[1].inv();T[1].resize(n);f.resize(n);
F[1]=mulT(f,T[1],n,2);
ask(1,1,n);
}
}
namespace lagrange{
int x[N],y[N];
poly g[100010<<2];
inline void build(int p,int l,int r){
if(l==r){g[p][0]=dec(0,x[l]);g[p][1]=1;return ;}
build(lc,l,mid);build(rc,mid+1,r);
g[p]=g[lc]*g[rc];
}
inline poly query(int p,int l,int r){
if(l==r) return 1ll*y[l]*ksm(Multiask::ans[l],mod-2)%mod;
return (query(lc,l,mid)*g[rc])+(query(rc,mid+1,r)*g[lc]);
}
inline void work(int n){
build(1,1,n);
Multiask::f=g[1].der();
for(int i=1;i<=n;++i) Multiask::a[i]=x[i];
Multiask::query(n);
query(1,1,n).print(n);
}
}
数组版poly：略快但难以使用
view code
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<19)+20;
const int mod=998244353;
typedef vector<int> vec;
typedef unsigned long long ull;
namespace IO{
inline char nc(){
static char buf[500005],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,500000,stdin),p1==p2)?EOF:*p1++;
}
char out[500005],*pout=out,*eout=out+500000;
inline void flush() { fwrite(out,1,pout-out,stdout),pout=out; }
inline void pc(char c) { pout==eout&&(fwrite(out,1,500000,stdout),pout=out); (*pout++)=c; }
template<typename T> inline void put(T x,char suf) {
static char stk[40];int top=0; while(x) stk[top++]=x%10,x/=10;
!top?pc('0'),0:0; while(top--) pc(stk[top]+'0'); pc(suf);
}
inline int read(){
char ch=nc(); int sum=0; for(;ch<'0'||ch>'9';ch=nc());
while(ch>='0'&&ch<='9')sum=(sum<<1)+(sum<<3)+ch-48,ch=nc();
return sum;
}
}
#define Rint IO::read()
using IO::put;
using IO::nc;
//IObuff没判负数
namespace Math{
inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
inline void inc(int &x,int y){x=add(x,y);}
inline void rec(int &x,int y){x=dec(x,y);}
inline int ksm(int x,int y){
int ret=1;
for(;y;y>>=1,x=1ll*x*x%mod) if(y&1) ret=1ll*ret*x%mod;
return ret;
}
int iv[N],tp;
inline void init_inv(int n){
if(!tp){iv[0]=iv[1]=1;tp=2;}
for(;tp<=n;++tp) iv[tp]=1ll*(mod-mod/tp)*iv[mod%tp]%mod;
}
}
using namespace Math;
namespace Cipolla{
int I,fl=0;
mt19937 rnd(time(0));
struct pt{int a,b;pt(int _a=0,int _b=0){a=_a;b=_b;}};
inline pt operator *(pt x,pt y){
pt ret;
ret.a=add(1ll*x.a*y.a%mod,1ll*x.b*y.b%mod*I%mod);
ret.b=add(1ll*x.a*y.b%mod,1ll*x.b*y.a%mod);
return ret;
}
inline bool check(int x){return ksm(x,(mod-1)/2)==1;}
inline int random(){return rnd()%mod;}
inline pt qpow(pt a,int b){
pt ret=pt(1,0);
for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
return ret;
}
inline int cipolla(int n){
if(!check(n)) return 0;
int a=random();
while(!a||check(dec(1ll*a*a%mod,n))) a=random();
I=dec(1ll*a*a%mod,n);
int ans=qpow(pt(a,1),(mod+1)/2).a;
return min(ans,(int)mod-ans);
}
}
using namespace Cipolla;
struct poly{
vec v;
inline poly(int w=0):v(1){v[0]=w;}
inline poly(const vec&w):v(w){}
inline int operator [](int x)const{return x>=v.size()?0:v[x];}
inline int& operator [](int x){if(x>=v.size()) v.resize(x+1);return v[x];}
inline int size(){return v.size();}
inline void resize(int x){v.resize(x);}
inline void read(int n){v.resize(n);for(int i=0;i<n;++i) v[i]=Rint;}
inline void print(int n)const{for(int i=0;i<n-1;++i) put(operator[](i),' ');put(operator[](n-1),'\n');}
inline poly slice(int len)const{
if(len<=v.size()) return vec(v.begin(),v.begin()+len);
vec ret(v);ret.resize(len);
return ret;
}
inline poly operator *(const int &x)const{
poly ret(v);
for(int i=0;i<v.size();++i) ret[i]=1ll*ret[i]*x%mod;
return ret;
}
inline poly operator +(const poly &x)const{
vec ret(max(v.size(),x.v.size()));
for(int i=0;i<v.size();++i) ret[i]=add(v[i],x[i]);
return ret;
}
inline poly operator -(const poly &x)const{
vec ret(max(v.size(),x.v.size()));
for(int i=0;i<v.size();++i) ret[i]=dec(v[i],x[i]);
return ret;
}
};
int Wn[N<<1],lg[N],r[N],tot;
inline void init_poly(int n){
int p=1;while(p<=n)p<<=1;
for(int i=2;i<=p;++i) lg[i]=lg[i>>1]+1;
for(int i=1;i<p;i<<=1){
int wn=ksm(3,(mod-1)/(i<<1));
Wn[++tot]=1;
for(int j=1;j<i;++j) ++tot,Wn[tot]=1ll*Wn[tot-1]*wn%mod;
}
}
inline void init_pos(int lim){
int len=lg[lim]-1;
for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<len);
}
ull fr[N];
const ull Mod=998244353;
inline void NTT(int *f,int lim,int tp){
for(int i=0;i<lim;++i) fr[i]=f[r[i]];
for(int mid=1;mid<lim;mid<<=1){
for(int len=mid<<1,l=0;l+len-1<lim;l+=len){
for(int k=l;k<l+mid;++k){
ull w1=fr[k],w2=fr[k+mid]*Wn[mid+k-l]%Mod;
fr[k]=w1+w2;fr[k+mid]=w1+Mod-w2;
}
}
}
for(int i=0;i<lim;++i) f[i]=fr[i]%Mod;
if(!tp){
reverse(f+1,f+lim);
int iv=ksm(lim,mod-2);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*iv%mod;
}
}
inline poly to_poly(int *a,int n){
poly ret;
ret.resize(n);
memcpy(ret.v.data(),a,n<<2);
return ret;
}
inline poly der(poly F,int n){
vec ret;ret.resize(n-1);
for(int i=0;i<n-1;++i) ret[i]=1ll*F[i+1]*(i+1)%mod;
return ret;
}
inline poly jifen(poly F,int n){
vec ret;ret.resize(n);
init_inv(n);
for(int i=n-1;i;--i) ret[i]=1ll*F[i-1]*iv[i]%mod;
return ret;
}
namespace Basic{
int f[N],g[N];
inline poly mul(poly F,poly G,int n,int m){
int rec=n+m-1,d=max(n,m);F.resize(n);G.resize(m);
memcpy(f,F.v.data(),sizeof(int)*(n));
memcpy(g,G.v.data(),sizeof(int)*(m));
if(d<=150){
poly ret;ret.resize(rec);
for(int i=0;i<n;++i)
for(int j=0;j<m;++j) inc(ret[i+j],1ll*f[i]*g[j]%mod);
return ret;
}
int len=lg[rec],lim=1<<len+1;
init_pos(lim);
NTT(f,lim,1);NTT(g,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
NTT(f,lim,0);
return to_poly(f,rec);
}
int g0[N],f0[N];
inline poly Inv(poly F,int n){
int l=1<<lg[n]+1;
memset(g,0,sizeof(int)*(l));
memset(g0,0,sizeof(int)*(l));
memset(f0,0,sizeof(int)*(l));
memcpy(f0,F.v.data(),sizeof(int)*(n));
g[0]=ksm(f0[0],mod-2);
for(int lim=2,nlim=1;nlim<n;lim<<=1,nlim<<=1){
memcpy(g0,g,sizeof(int)*(nlim));
memcpy(f,f0,sizeof(int)*(lim));
init_pos(lim);
NTT(g0,lim,1);NTT(f,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*g0[i]*f[i]%mod;
NTT(f,lim,0);
fill(f,f+nlim,0);
NTT(f,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g0[i]%mod;
NTT(f,lim,0);
for(int i=nlim;i<lim;++i) rec(g[i],f[i]);
}
return to_poly(g,n);
}
inline void Inv(int *F,int *G,int n){
int l=1<<lg[n]+1;
memset(g,0,sizeof(int)*(l));
memset(g0,0,sizeof(int)*(l));
memset(f,0,sizeof(int)*(l));
memcpy(f0,F,sizeof(int)*(n));
memset(f0+n,0,sizeof(int)*(l-n));
g[0]=ksm(f0[0],mod-2);
for(int lim=2,nlim=1;nlim<n;lim<<=1,nlim<<=1){
memcpy(g0,g,sizeof(int)*(nlim));
memcpy(f,f0,sizeof(int)*(lim));
init_pos(lim);
NTT(g0,lim,1);NTT(f,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*g0[i]*f[i]%mod;
NTT(f,lim,0);
fill(f,f+nlim,0);
NTT(f,lim,1);
for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g0[i]%mod;
NTT(f,lim,0);
for(int i=nlim;i<lim;++i) rec(g[i],f[i]);
}
memcpy(G,g,sizeof(int)*(n));
}
}
namespace Ln{
int f[N],g0[N],h0[N],q0[N],ret[N],rc[N];
inline poly div(poly H,poly F,int n,int m){
if(n==1) return 1ll*H[0]*ksm(F[0],mod-2)%mod;
int len=lg[n],lim=1<<len+1,nlim=lim>>1;
memcpy(f,F.v.data(),sizeof(int)*(m));
memcpy(g0,F.v.data(),sizeof(int)*(nlim));
memcpy(h0,H.v.data(),sizeof(int)*(nlim));
memcpy(rc,H.v.data(),sizeof(int)*(min(lim,m)));
memset(f+n,0,sizeof(int)*(lim-n));
memset(g0+nlim,0,sizeof(int)*(nlim));
memset(h0+nlim,0,sizeof(int)*(nlim));
memset(ret+nlim,0,sizeof(int)*(nlim));
Basic::Inv(g0,g0,nlim);
init_pos(lim);
NTT(g0,lim,1);NTT(h0,lim,1);
for(int i=0;i<lim;++i) q0[i]=1ll*g0[i]*h0[i]%mod;
NTT(q0,lim,0);memset(q0+nlim,0,sizeof(int)*(nlim));
memcpy(ret,q0,sizeof(int)*(nlim));
NTT(q0,lim,1);NTT(f,lim,1);
for(int i=0;i<lim;++i) q0[i]=1ll*q0[i]*f[i]%mod;
NTT(q0,lim,0);
memset(q0,0,sizeof(int)*(nlim));
for(int i=nlim;i<lim&&i<m;++i) rec(q0[i],rc[i]);
NTT(q0,lim,1);
for(int i=0;i<lim;++i) q0[i]=1ll*q0[i]*g0[i]%mod;
NTT(q0,lim,0);
for(int i=nlim;i<lim;++i) rec(ret[i],q0[i]);
return to_poly(ret,n);
}
inline poly getln(poly F,int n){
return jifen(div(der(F,n),F,n,n),n);
}//11E
}
namespace Exp{
const int logB=4;
const int B=16;
int f[N],ret[N],H[N];
poly g[4][B];
inline void exp(int lim,int l,int r,int dep){
if(r-l<=128){
for(int i=l;i<r;++i){
ret[i]=(!i)?1:1ll*ret[i]*iv[i]%mod;
for(int j=i+1;j<r;++j)
inc(ret[j],1ll*ret[i]*f[j-i]%mod);
}
return ;
}
int k=(r-l)/B;
int len=1<<lim-logB+1;
vector<unsigned long long> bl[B];
for(int i=0;i<B;++i) bl[i].resize(k<<1);
for(int i=0;i<B;++i){
if(i>0){
init_pos(len);
for(int j=0;j<(k<<1);++j) H[j]=bl[i][j]%mod;
NTT(H,len,0);
for(int j=0;j<k;++j)
inc(ret[l+i*k+j],H[j+k]);
}
exp(lim-logB,l+i*k,l+(i+1)*k,dep+1);
if(i<B-1){
memcpy(H,ret+l+i*k,sizeof(int)*(k));
memset(H+k,0,sizeof(int)*(k));
init_pos(len);NTT(H,len,1);
for(int j=i+1;j<B;++j)
for(int t=0;t<(k<<1);++t)
bl[j][t]+=1ll*H[t]*g[dep][j-i-1][t];
}
}
}
inline poly getexp(poly F,int n){
memcpy(f,F.v.data(),sizeof(int)*(n));
int mx=lg[n]+1;init_inv(1<<mx);
for(int i=0;i<n;++i) f[i]=1ll*f[i]*i%mod;
memset(ret,0,sizeof(int)*(1<<mx));
for(int lim=mx,dep=0;lim>=8;lim-=logB,dep++){
int len=1<<(lim-logB+1);
init_pos(len);
for(int i=0;i<B-1;++i){
g[dep][i].resize(len);
memcpy(g[dep][i].v.data(),f+(len>>1)*i,sizeof(int)*(len));
NTT(g[dep][i].v.data(),len,1);
}
}
exp(mx,0,1<<mx,0);
return to_poly(ret,n);
}
inline poly getpow(poly F,int n,int k){
F=Ln::getln(F,n);
return getexp(F*k,n);
}
}
namespace Sqrt{
int f[N],g[N],h[N],gf[N],F1[N],F2[N],F3[N];
inline void allmem(int n){
memset(f,0,sizeof(int)*(n));
memset(g,0,sizeof(int)*(n));
memset(h,0,sizeof(int)*(n));
memset(gf,0,sizeof(int)*(n));
memset(F1,0,sizeof(int)*(n));
memset(F2,0,sizeof(int)*(n));
memset(F3,0,sizeof(int)*(n));
}
inline poly sqrt(poly F,int n){
// allmem(1<<lg[n]+1);
memcpy(f,F.v.data(),sizeof(int)*(n));
g[0]=cipolla(f[0]);
h[0]=ksm(g[0],mod-2);
gf[0]=g[0];gf[1]=g[0];
int iv=(mod+1)/2;
init_pos(1);
for(int lim=1;lim<n;lim<<=1){
for(int i=0;i<lim;++i) F1[i]=1ll*gf[i]*gf[i]%mod;
NTT(F1,lim,0);
for(int i=0;i<lim;++i) F1[i+lim]=dec(F1[i],f[i]),F1[i]=0;
int nlim=lim<<1;init_pos(nlim);
for(int i=lim;i<nlim;++i) rec(F1[i],f[i]);
memcpy(F2,h,sizeof(int)*(lim));
NTT(F1,nlim,1);NTT(F2,nlim,1);
for(int i=0;i<nlim;++i) F1[i]=1ll*F1[i]*F2[i]%mod;
NTT(F1,nlim,0);
for(int i=lim;i<nlim;++i) g[i]=dec(0,1ll*F1[i]*iv%mod);
if(nlim<n){
memcpy(gf,g,sizeof(int)*(nlim));
NTT(gf,nlim,1);
for(int i=0;i<nlim;++i) F3[i]=1ll*gf[i]*F2[i]%mod;
NTT(F3,nlim,0);
memset(F3,0,sizeof(int)*(lim));
NTT(F3,nlim,1);
for(int i=0;i<nlim;++i) F3[i]=1ll*F3[i]*F2[i]%mod;
NTT(F3,nlim,0);
for(int i=lim;i<nlim;++i) rec(h[i],F3[i]);
}
}
return to_poly(g,n);
}
}
int n,k;
poly f;
inline int sread(){
long long m=mod,x=0;char ch=nc();
while(!isdigit(ch)){ch=nc();}
while(isdigit(ch)){x=(x<<3)+(x<<1)+(ch^48);x%=m;ch=nc();}
return x;
}
int main(){
// freopen("poly.in","r",stdin);
// freopen("1.out","w",stdout);
n=Rint+1,k=Rint;
f.read(n);
init_poly(n);
poly g=f+poly(2)-poly(f[0]);
f=Basic::Inv(Sqrt::sqrt(f,n),n);
g=Ln::getln(g-Exp::getexp(jifen(f,n),n),n);
der(Exp::getpow((g+poly(1)),n,k),n).print(n-1);
IO::flush();
return 0;
}
