LOJ#3058. 「HNOI2019」白兔之舞

LOJ#3058. 「HNOI2019」白兔之舞

https://loj.ac/problem/3058

分析

• 不妨設轉移矩陣爲\(f\)，且用\(f^i\)來表示走\(i\)步從\(x\)到\(y\)的方案數。
• 跳了\(i\)步的方案數爲\(\sum\limits_{j=i}^L\binom{j-1}{i-1}=\binom{L}{i}\)
• 直接推式子\(ans_t=\sum\limits_{i=0}^L[i\bmod k=t]f^i\binom{L}{i}\)
• 發現這是個挺套路的單位根反演

• \(ans_t\times k=\sum\limits_{i=0}^{L}\sum\limits_{j=0}^{k-1}w^{j(i-t)}f^i\binom{L}{i}\)

  • ​ $ =\sum\limits_{j=0}^{k-1}w^{-jt}\sum\limits_{i=0}^Lw^{ij}f^i\binom{L}{i}$
  • ​ \(=\sum\limits_{j=0}^{k-1}w^{-jt}(w^jf+I)\)
• 把乘積化成和的形式，有\(ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}\)
• \(ans_t\times k=\sum\limits_{j=0}^{k-1}w^{-\binom{j+t}{2}+\binom{j}{2}+\binom{t}{2}}g_j\)
• \(ans_t\times k\times w^{-\binom{t}{2}}=\sum\limits_{j=0}^{k-1}w^{-\binom{j+t}{2}}w^{\binom{j}{2}}g_j​\)

• 是一個卷積的形式，直接用\(mtt\)作就好了，用https://cmxrynp.github.io/2019/01/07/fft-optimization/中的\(4\)次或\(3.5\)次dft再加上一點點卡常便可經過此題。

代碼

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;
typedef long double f2;
#define N 262190
#define db(x) cerr<<#x<<" = "<<x<<endl
namespace groot {
    int pri[100],cnt;
    ll qp(ll x,ll y,ll p) {
        ll re=1;for(;y;y>>=1,x=x*x%p)if(y&1)re=re*x%p; return re;
    }
    bool check(int x,int phi) {
        int i;
        for(i=1;i<=cnt;i++) if(qp(x,phi/pri[i],phi+1)==1) return 0;
        return 1;
    }
    int get_g(int p,int phi) {
        int x=phi,i;
        for(i=2;i*i<=x;i++) {
            if(x%i==0) {
                pri[++cnt]=i;
                while(x%i==0) x/=i;
            }
        }
        if(x!=1) pri[++cnt]=x;
        int g=2;while(!check(g,phi)) g++;
        return g;
    }
}
int mod,n,K,L,X,Y,G;
ll g[N],w[N],f[N];
ll qp(ll x,ll y) {
    ll re=1; for(;y;y>>=1,x=x*x%mod) if(y&1) re=re*x%mod; return re;
}
ll C2(ll x) {return x*(x-1)/2%K;}
namespace matrix {
    struct Mat {
        ll a[3][3];
        Mat() {memset(a,0,sizeof(a));}
        Mat operator * (const Mat &u) const {
            int i,j,k; Mat re;
            for(i=0;i<n;i++)for(k=0;k<n;k++)for(j=0;j<n;j++) {
                re.a[i][j]+=a[i][k]*u.a[k][j];
            }for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]%=mod;
            return re;
        }
        Mat operator * (const ll &x) const {
            Mat re=*this; int i,j;
            for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]=re.a[i][j]*x%mod;
            return re;
        }
        Mat operator + (const Mat &u) const {
            Mat re=*this; int i,j;
            for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]=(re.a[i][j]+u.a[i][j])%mod;
            return re;
        }
    }I,O;
    Mat qp(Mat x,int y) {
        Mat re=I;
        for(;y;y>>=1,x=x*x)if(y&1)re=re*x; return re;
    }
    void solve() {
        int i;
        for(i=0;i<n;i++) I.a[i][i]=1;
        for(i=0;i<K;i++) {
            g[i]=qp(O*w[i]+I,L).a[X-1][Y-1]*w[C2(i)]%mod;
        }
    }
}
ll tmp[N];
namespace fft {
    const f2 pi=acos(-1);
    struct cp {
        f2 x,y;
        cp() {}
        cp(f2 x_,f2 y_) {x=x_,y=y_;}
        cp operator + (const cp &u) const {return cp(x+u.x,y+u.y);}
        cp operator - (const cp &u) const {return cp(x-u.x,y-u.y);}
        cp operator * (const cp &u) const {return cp(x*u.x-y*u.y,x*u.y+y*u.x);}
    }A[N],B[N],C[N],D[N];
    void dft(cp *a,int len) {
        int i,j,k,t; cp tmp,w,wn;
        for(i=k=0;i<len;i++) {
            if(i>k) swap(a[i],a[k]);
            for(j=len>>1;(k^=j)<j;j>>=1);
        }
        for(k=2;k<=len;k<<=1) {
            wn=cp(cos(2*pi/k),sin(2*pi/k));
            for(t=k>>1,i=0;i<len;i+=k) {
                for(w=cp(1,0),j=i;j<i+t;j++) {
                    tmp=a[j+t]*w; a[j+t]=a[j]-tmp; a[j]=a[j]+tmp; w=w*wn;
                }
            }
        }
    }
    int solve() {
        int i;
        w[K]=1;
        for(i=0;i<K+K;i++) f[i]=w[K-C2(i)];
        reverse(g,g+K+1);
        int n=K+K-1,m=K;
        int l=1;
        while(l<=(n+m))l<<=1;
        //for(i=0;i<l;i++) printf("%lld %lld\n",f[i],g[i]);
        for(i=0;i<=n;i++) {
            A[i].x=f[i]>>15;
            A[i].y=f[i]&32767;
        }
        for(i=0;i<=m;i++) {
            B[i].x=g[i]>>15;
            B[i].y=g[i]&32767;
        }
        dft(A,l),dft(B,l);
        for(i=0;i<l;i++) {
            int pl=(l-1)&(l-i);
            const cp ca=cp(A[pl].x,-A[pl].y),cb=cp(B[pl].x,-B[pl].y);
            const cp a=(A[i]+ca)*cp(0.5,0),b=(A[i]-ca)*cp(0,-0.5),c=(B[i]+cb)*cp(0.5,0),d=(B[i]-cb)*cp(0,-0.5);
            C[pl]=a*c+a*d*cp(0,1);
            D[pl]=b*c+b*d*cp(0,1);
        }
        dft(C,l),dft(D,l);
        ll invk=qp(K,mod-2);
        for(i=0;i<l;i++) C[i].x/=l,C[i].y/=l,D[i].x/=l,D[i].y/=l;
        //for(i=0;i<=n;i++) for(int j=0;j<=m;j++) {
            //tmp[i+j]=(tmp[i+j]+f[i]*g[j])%mod;
        //}
        for(i=K;i<K+K;i++) {
            ll a=C[i].x+0.5,b=C[i].y+0.5,c=D[i].x+0.5,d=D[i].y+0.5;
            a=((a%mod)<<30)+(((b+c)%mod)<<15)+d;
            a=(a%mod+mod)%mod;
            //a=tmp[i];
            //db(a);
            ll t=w[C2(i-K)]*invk%mod;
            printf("%lld\n",a*t%mod);
        }
        return 0;
    }
}
int main() {
    scanf("%d%d%d%d%d%d",&n,&K,&L,&X,&Y,&mod);
    G=groot::get_g(mod,mod-1);
    w[0]=1;w[1]=qp(G,(mod-1)/K);
    int i,j;
    for(i=2;i<K;i++) w[i]=w[i-1]*w[1]%mod;
    //printf("%d\n",G);
    for(i=0;i<n;i++)for(j=0;j<n;j++)scanf("%lld",&matrix::O.a[i][j]);
    matrix::solve();
    return fft::solve();
}