第一个自己实现的多项式全家桶。打比赛时终于有板子了!
代码是从之前学转置原理的那篇博客里升级来的。
但是功能很残?效率很逊?接口很怪?评价是能用就行。
目前封装了:乘、逆、除(取模)、开方(无二次剩余,不会 qwq)、对数、指数、多点求值、快速插值、常系数齐次线性递推、Berlekamp-Massey 算法。
为什么突然开始多项式?原因是在尝试改鸡的时候发现了 sjy 的究极神秘题解,感觉 BM 算法竟然有意想不到的用处,于是学了 BM 算法,然后递归到了常系数其次线性递推,然后递归到了除法,然后升级了多项式板子。
不过话说 BM 算法其实与多项式没有很大关系不是吗?
#include <cstdio>
#include <vector>
#include <cassert>
#include <algorithm>
#define lc (p<<1)
#define rc (p<<1|1)
#define ALL(p) p.begin(),p.end()
#define IL inline
// #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<22,stdin)),p1==p2?EOF:*p1++)
using namespace std;
char buf[1<<22],*p1=buf,*p2=buf;
int read(){
char c=getchar();int x=0;bool f=0;
while(c<48||c>57) f|=(c=='-'),c=getchar();
do x=(x<<1)+(x<<3)+(c^48),c=getchar();
while(c>=48&&c<=57);
if(f) return -x;
return x;
}
typedef vector<int> vi;
const int P=998244353;
typedef long long ll;
IL int qp(int a,int b=P-2){
int res=1;
while(b){
if(b&1) res=(ll)res*a%P;
a=(ll)a*a%P;b>>=1;
}
return res;
}
int len,ilen,bt;
int rev[1<<20],cw[1<<20|1];
IL void init(int _len){ // mod x^len
len=1;bt=-1;
while(len<_len) len<<=1,++bt;
int w=qp(3,(P-1)>>(bt+1));
cw[0]=cw[len]=1;
for(int i=1;i<len;++i){
cw[i]=(ll)cw[i-1]*w%P;
rev[i]=(rev[i>>1]>>1)|((i&1)<<bt);
}
ilen=qp(len);
}
struct poly{
vi f;
poly():f(){}
poly(int Len):f(Len){}
poly(initializer_list<int> Init):f(Init){}
IL int& operator[](const int &x){return f[x];}
IL const int& operator[](const int &x)const{return f[x];}
IL void NTT(){
f.resize(len,0);
for(int i=1;i<len;++i) if(rev[i]<i) swap(f[rev[i]],f[i]);
for(int i=1,tt=len>>1;i<len;i<<=1,tt>>=1)
for(int j=0;j<len;j+=(i<<1))
for(int k=j,t=0;k<(j|i);++k,t+=tt){
int x=f[k],y=(ll)f[k|i]*cw[t]%P;
if((f[k]=x+y)>=P) f[k]-=P;
if((f[k|i]=x-y)<0) f[k|i]+=P;
}
}
IL void INTT(){
for(int i=1;i<len;++i) if(rev[i]<i) swap(f[rev[i]],f[i]);
for(int i=1,tt=len>>1;i<len;i<<=1,tt>>=1)
for(int j=0;j<len;j+=(i<<1))
for(int k=j,t=len;k<(j|i);++k,t-=tt){
int x=f[k],y=(ll)f[k|i]*cw[t]%P;
if((f[k]=x+y)>=P) f[k]-=P;
if((f[k|i]=x-y)<0) f[k|i]+=P;
}
for(int i=0;i<len;++i) f[i]=(ll)f[i]*ilen%P;
}
IL int size(){return f.size();}
IL void reduce(){while(!f.empty()&&!f.back()) f.pop_back();}
IL void rever(){reverse(ALL(f));}
IL void show(){
for(int x:f) printf("%d ",x);
putchar('\n');
}
IL void trunc(int lim){
if(lim<int(f.size())) f.erase(f.begin()+lim,f.end());
}
IL poly inv(int lim){ // mod x^lim
assert(f[0]);
poly cur({qp(f[0])});
for(int t=1;t<lim;t<<=1){
poly ff(t<<2);
copy(f.begin(),f.begin()+min(t<<1,size()),ff.f.begin());
init(t<<2);ff.NTT();cur.NTT();
poly tmp(len);
for(int i=0;i<len;++i){
tmp[i]=(2ll-(ll)cur[i]*ff[i])%P*cur[i]%P;
if(tmp[i]<0) tmp[i]+=P;
}
tmp.INTT();
cur.f.swap(tmp.f);
cur.trunc(t<<1);
}
cur.trunc(lim);
return cur;
}
friend poly operator+(poly A,poly B){
int mx=max(A.size(),B.size());
A.f.resize(mx,0);B.f.resize(mx,0);
poly C(mx);
for(int i=0;i<mx;++i){
C[i]=A[i]+B[i];
if(C[i]>=P) C[i]-=P;
}
return C;
}
friend poly operator-(poly A,poly B){
int mx=max(A.size(),B.size());
A.f.resize(mx,0);B.f.resize(mx,0);
poly C(mx);
for(int i=0;i<mx;++i){
C[i]=A[i]-B[i];
if(C[i]<0) C[i]+=P;
}
C.reduce();
return C;
}
friend poly operator*(poly A,poly B){
init(A.size()+B.size()-1);A.NTT();B.NTT();
poly C(len);
for(int i=0;i<len;++i) C[i]=(ll)A[i]*B[i]%P;
C.INTT();C.reduce();
return C;
}
poly operator*(int X){
int n=size();
poly F(n);
for(int i=0;i<n;++i) F[i]=(ll)f[i]*X%P;
return F;
}
friend poly operator%(poly F,poly G){
int n=F.size()-1,m=G.size()-1;
if(n<m) return F;
F.rever();G.rever();
poly Q=G.inv(n-m+1);
Q.prod(Q,F);Q.f.resize(n-m+1,0);
F.rever();G.rever();Q.rever();
poly R=F-Q*G;
R.reduce();
return R;
}
IL void prod(poly A,poly B){
init(A.size()+B.size()-1);A.NTT();B.NTT();
f.resize(len);
for(int i=0;i<len;++i) f[i]=(ll)A[i]*B[i]%P;
INTT();reduce();
}
IL void prodT(poly A,poly B){
int an=A.size()-1,bn=B.size()-1;
reverse(ALL(B.f));prod(A,B);
for(int i=0;i<=an;++i) f[i]=f[i+bn];
trunc(an+1);
}
IL int calc(int t){
int pw=1,res=0;
for(int x:f){
res=(res+(ll)pw*x)%P;
pw=(ll)pw*t%P;
}
return res;
}
IL poly deriv(){
int n=f.size();
poly D(n-1);
for(int i=1;i<n;++i) D[i-1]=(ll)f[i]*i%P;
return D;
}
IL poly integ(){
int n=f.size();
poly D(n+1),inv(n+1);
for(int i=1;i<=n;++i){
if(i>1) inv[i]=(ll)inv[P%i]*(P-P/i)%P;
else inv[1]=1;
D[i]=(ll)f[i-1]*inv[i]%P;
}
return D;
}
IL poly getln(int lim){ // mod x^lim
assert(f[0]==1);
poly F=deriv()*inv(lim);
F.f.resize(lim-1,0);
return F.integ();
}
IL poly getexp(int lim){ // mod x^lim
assert(f[0]==0);
poly F({1});
for(int t=1;t<lim;t<<=1){
poly ff(t<<1);
copy(f.begin(),f.begin()+min(t<<1,size()),ff.f.begin());
F=F*(poly({1})-F.getln(t<<1)+ff);
F.f.resize(t<<1,0);
}
F.trunc(lim);
return F;
}
IL poly sqrt(int lim){ // without quadratic residue , mod x^lim
assert(f[0]==1);
poly F({1});
for(int t=1;t<lim;t<<=1){
poly ff(t<<1);
copy(f.begin(),f.begin()+min(t<<1,size()),ff.f.begin());
F=(F+F.inv(t<<1)*ff)*((P+1)>>1);
F.f.resize(t<<1,0);
}
F.trunc(lim);
return F;
}
}F;
namespace extend{
const int N=200003;
int n;
int px[N],py[N],res[N];
namespace multieval{
poly G[N<<2],H[N<<2];
void calc(int p,int l,int r){
if(l==r){G[p]=poly({1,(P-px[r])%P});return;}
int mid=(l+r)>>1;
calc(lc,l,mid);
calc(rc,mid+1,r);
G[p].prod(G[lc],G[rc]);
}
void solve(int p,int l,int r){
H[p].trunc(r-l+1);
if(l==r){res[r]=H[p][0];return;}
int mid=(l+r)>>1;
G[p].f.clear();G[p].f.shrink_to_fit();
H[lc].prodT(H[p],G[rc]);
H[rc].prodT(H[p],G[lc]);
H[p].f.clear();H[p].f.shrink_to_fit();
solve(lc,l,mid);
solve(rc,mid+1,r);
}
void sol(){
// init n, px, py first
calc(1,1,n);
H[1].prodT(F,G[1].inv(n+1));
solve(1,1,n);
}
}
namespace interplot{
poly T[N<<2];
void eval(int p,int l,int r){
if(l==r){T[p]=poly({(P-px[r])%P,1});return;}
int mid=(l+r)>>1;
eval(lc,l,mid);
eval(rc,mid+1,r);
T[p].prod(T[lc],T[rc]);
}
poly proc(int p,int l,int r){
if(l==r) return poly({py[r]});
int mid=(l+r)>>1;
poly L=proc(lc,l,mid);
poly R=proc(rc,mid+1,r);
L.prod(L,T[rc]);
R.prod(R,T[lc]);
return L+R;
}
poly solve(){
// init n, px, py first
eval(1,1,n);
F=T[1].deriv();
multieval::sol();
for(int i=1;i<=n;++i) py[i]=(ll)py[i]*qp(res[i])%P;
return proc(1,1,n);
}
}
}
namespace berlekamp_massey{
const int N=200003;
int n,k,kk,las,pos,ex;
int s[N],ss[N],a[N],tmp[N];
int solve(){
n=read();ex=read();
for(int i=0;i<n;++i) a[i]=read();
bool fl=1;
for(int i=0;i<n;++i){
int cur=a[i];
for(int j=1;j<=k;++j) cur=(cur-(ll)a[i-j]*s[j])%P;
if(cur<0) cur+=P;
if(!cur) continue;
if(fl){
k=i+1;kk=0;pos=i;las=cur;fl=0;
continue;
}
int nk=max(k,i-pos+kk);
int val=(ll)cur*qp(las)%P;
for(int j=1;j<i-pos;++j) tmp[j]=0;
tmp[i-pos]=val;
for(int j=1;j<=kk;++j) tmp[i-pos+j]=(ll)(P-val)*ss[j]%P;
if(k-i<kk-pos){
kk=k;pos=i;las=cur;
for(int j=1;j<=k;++j) ss[j]=s[j];
}
k=nk;
for(int j=1;j<=k;++j){
s[j]+=tmp[j];
if(s[j]>=P) s[j]-=P;
}
}
poly M(k+1),F(2),G({1});
M[k]=1;
for(int i=1;i<=k;++i)
if(s[i]) M[k-i]=P-s[i];
F[1]=1;
while(ex){
if(ex&1) G=G*F%M;
F=F*F%M;ex>>=1;
}
G.f.resize(k,-1);
int ans=0;
for(int i=0;i<k;++i) ans=(ans+(ll)G[i]*a[i])%P;
if(ans<0) ans+=P;
return ans;
}
}
int main(){
int n=read();
poly F(n);
for(int i=0;i<n;++i) F[i]=read();
F=F.sqrt(n);
F.show();
return 0;
}
标签:return,lc,int,多项式,全家,mid,poly,rc
From: https://www.cnblogs.com/yyyyxh/p/17979236/poly