LOJ#3058. 「HNOI2019」白兔之舞 單位根反演+矩陣乘法+MTT
阿新 • • 發佈:2020-07-28
毒瘤多項式.
找原根+矩陣乘法+拆係數FFT大雜燴.
code:
#include <vector> #include <cmath> #include <cstdio> #include <cstring> #include <algorithm> #define N 100009 #define ll long long #define pb push_back #define db long double #define setIO(s) freopen(s".in","r",stdin) using namespace std; const db pi=acos(-1.0); vector<int>g; int gn,n,L,S,T,P,K,wi[N],val[5][5],seq[N<<1],os[N<<1]; int qpow(int x,int y,int mod) { int tmp=1; for(;y;y>>=1,x=(ll)x*x%mod) { if(y&1) tmp=(ll)tmp*x%mod; } return tmp; } int getroot(int x) { int phi=x-1; for(int i=2;i*i<=x;++i) { if(phi%i==0) { g.pb(i); while(phi%i==0) phi/=i; } } if(phi>1) g.pb(phi); phi=x-1; for(int i=2;;++i) { int flag=1; for(int j=0;j<g.size()&&flag;++j) { if(qpow(i,phi/g[j],x)==1) flag=0; } if(flag) return i; } } ll sqr(ll x) { return x*(x-1)/2; } void init() { gn=qpow(getroot(P),(P-1)/K,P); wi[0]=1; for(int i=1;i<K;++i) { wi[i]=(ll)wi[i-1]*gn%P; } } struct Matrix { int c[3][3]; Matrix(int x=0) { memset(c,0,sizeof(c)); for(int i=0;i<3;++i) c[i][i]=x; } int *operator[](int x){ return c[x]; } Matrix operator*(const Matrix &b) const { Matrix an; for(int i=0;i<3;++i) for(int j=0;j<3;++j) { for(int k=0;k<3;++k) (an.c[i][j]+=(ll)c[i][k]*b.c[k][j]%P)%=P; } return an; } friend Matrix operator^(Matrix x,int y) { Matrix an(1); for(;y;y>>=1,x=x*x) if(y&1) an=an*x; return an; } }W; const int base=1<<15; struct Poly { struct cp { db x,y; cp(db a=0,db b=0) { x=a,y=b; } cp operator+(const cp &b) const { return cp(x+b.x,y+b.y); } cp operator-(const cp &b) const { return cp(x-b.x,y-b.y); } cp operator*(const cp &b) const { return cp(x*b.x-y*b.y,x*b.y+y*b.x); } }f[2][N<<2],g[2][N<<2],ans[3][N<<2]; void FFT(cp *a,int len,int op) { for(int i=0,k=0;i<len;++i) { if(i>k) swap(a[i],a[k]); for(int j=len>>1;(k^=j)<j;j>>=1); } for(int l=1;l<len;l<<=1) { cp wn(cos(pi/l),op*sin(pi/l)); for(int i=0;i<len;i+=l<<1) { cp w(1,0),x,y; for(int j=0;j<l;++j) { x=a[i+j],y=w*a[i+j+l]; a[i+j]=x+y,a[i+j+l]=x-y; w=w*wn; } } } if(op==-1) { for(int i=0;i<len;++i) a[i].x/=len; } } ll actual(db x,int mod) { return (ll)((ll)(x+0.5))%mod; } void MTT(int *a,int n,int *b,int m,int mod,int *c) { int lim; for(lim=1;lim<(n+m+1);lim<<=1); for(int i=0;i<lim;++i) { f[0][i]=f[1][i]=g[0][i]=g[1][i]=cp(); ans[0][i]=ans[1][i]=ans[2][i]=cp(); } for(int i=0;i<=n;++i) { f[0][i].x=a[i]/base; f[1][i].x=a[i]%base; } for(int i=0;i<=m;++i) { g[0][i].x=b[i]/base; g[1][i].x=b[i]%base; } FFT(f[0],lim,1),FFT(f[1],lim,1); FFT(g[0],lim,1),FFT(g[1],lim,1); for(int i=0;i<lim;++i) { ans[0][i]=f[0][i]*g[0][i]; ans[1][i]=f[0][i]*g[1][i]+f[1][i]*g[0][i]; ans[2][i]=f[1][i]*g[1][i]; } for(int i=0;i<3;++i) { FFT(ans[i],lim,-1); } for(int i=0;i<=n+m;++i) { ll x=(ll)(actual(ans[0][i].x,mod)<<30ll)%mod; ll y=(ll)(actual(ans[1][i].x,mod)<<15ll)%mod; ll z=actual(ans[2][i].x,mod); c[i]=(ll)((ll)(x+y)%mod+z)%mod; } } }poly; void get_matrix(int x) { for(int i=0;i<3;++i) { for(int j=0;j<3;++j) { W[i][j]=(ll)val[i+1][j+1]*wi[x]%P; } ++W[i][i]; } } int AA[N<<1],ANS[N*3]; int main() { // setIO("input"); scanf("%d%d%d%d%d%d",&n,&K,&L,&S,&T,&P); init(); for(int i=1;i<=n;++i) { for(int j=1;j<=n;++j) scanf("%d",&val[i][j]); } for(int i=0;i<K;++i) { get_matrix(i); W=W^L; seq[i]=(ll)W[S-1][T-1]*wi[(int)((ll)sqr(i)%K)]%P; } for(int i=0;i<=2*K-2;++i) { os[i]=(ll)wi[(int)((ll)(K-(ll)sqr(i)%K+K)%K)]; } for(int i=0;i<K;++i) AA[i]=seq[K-1-i]; poly.MTT(AA,K-1,os,2*K-2,P,ANS); int inv=qpow(K,P-2,P); for(int i=0;i<K;++i) { printf("%d\n",(ll)inv*wi[(int)((ll)sqr(i)%K)]%P*ANS[K-1+i]%P); } return 0; }