M_sea

洛谷5293 [HNOI2019]白兔之舞
LuoguLOJ分析orz laofu先看 $n=1$ 的情况。为了方便,设 $a=w_{1,1}$ 。设 $f_...
扫描右侧二维码阅读全文
03
2019/06

洛谷5293 [HNOI2019]白兔之舞

Luogu

LOJ

分析

orz laofu

先看 $n=1$ 的情况。为了方便,设 $a=w_{1,1}$ 。

设 $f_i$ 表示跳 $i$ 步的答案,显然有 $f_i=a^i{n\choose i}$ 。

那么有

$$\displaystyle ans_t=\sum_{i=0}^L[i\bmod k=t]{L\choose i}a^i$$

单位根反演,得到

$$\displaystyle\begin{aligned}ans_t&=\frac{1}{k}\sum_{i=0}^L\sum_{j=0}^{k-1}\omega_k^{j(i-t)}{L\choose i}a^i\\&=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-jt}\sum_{i=0}^L{L\choose i}\omega_k^{ij}a^i\\&=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-jt}(\omega_k^ja+1)^L\end{aligned}$$

设 $\displaystyle c_i=(\omega_k^ia+1)^L$ ,那么

$$\displaystyle ans_t=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-jt}c_j$$

考虑怎么处理这个 $\omega_k^{-jt}$ 。一个做法是将 $jt$ 拆成 $\frac{(j+t)^2}{2}-\frac{j^2}{2}-\frac{t^2}{2}$ ,但是 $\frac{x^2}{2}$ 可能不是整数。

于是我们换一种方法:$\displaystyle jt={j+t\choose 2}-{j\choose 2}-{t\choose 2}$ 。那么

$$\displaystyle\begin{aligned}ans_t&=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{{j-t\choose 2}-{j\choose 2}-{-t\choose 2}}c_j\\&=\frac{\omega_k^{-{-t\choose 2}}}{k}\sum_{j=0}^{k-1}\omega_k^{j-t\choose 2}\left(\omega_k^{-{j\choose 2}}c_j\right)\end{aligned}$$

显然是一个卷积的形式,那么直接 $\mathrm{MTT}$ 实现卷积即可。

$n\leq 3$ 的情况实际上就是扩展到了矩阵上,推式子的过程是一样的。

具体实现及细节见代码。

代码

不知道为什么写得特别长

// =================================
//   author: M_sea
//   website: http://m-sea-blog.com/
// =================================
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define re register
using namespace std;
typedef long long ll;

inline int read() {
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=600000;
const double PI=acos(-1);

int n,k,L,x,y,mod;

inline int mul(int a,int b) { return 1ll*a*b%mod; }
inline void add(int& a,int b) { a=(a+b)%mod; }

inline int qpow(int a,int b) {
    int c=1;
    for (;b;b>>=1,a=1ll*a*a%mod)
        if (b&1) c=1ll*c*a%mod;
    return c;
}

//求原根
inline int getroot(int x) {
    static int fac[100]; int cnt=0;
    for (re int i=2,t=x-1;i<=t;++i) {
        if (t%i) continue;
        fac[++cnt]=i;
        while (t%i==0) t/=i;
    }
    for (re int g=2;;++g) {
        int flag=1;
        for (re int j=1;j<=cnt;++j)
            if (qpow(g,(x-1)/fac[j])==1) { flag=0; break; }
        if (flag) return g;
    }
}

//矩阵
struct Matrix {
    int s[3][3];
    Matrix() { memset(s,0,sizeof(s)); }
    int* operator [](int i) { return s[i]; }
} A,I,S;

Matrix operator +(Matrix a,Matrix b) {
    Matrix c;
    for (re int i=0;i<3;++i)
        for (re int j=0;j<3;++j)
            c[i][j]=(a[i][j]+b[i][j])%mod;
    return c;
}

Matrix operator *(Matrix a,int b) {
    Matrix c;
    for (re int i=0;i<3;++i)
        for (re int j=0;j<3;++j)
            c[i][j]=mul(a[i][j],b);
    return c;
}

Matrix operator *(Matrix a,Matrix b) {
    Matrix c;
    for (re int i=0;i<3;++i)
        for (re int k=0;k<3;++k)
            for (re int j=0;j<3;++j)
                add(c[i][j],mul(a[i][k],b[k][j]));
    return c;
}

inline Matrix qpow(Matrix a,int b) {
    Matrix c=I;
    for (;b;b>>=1,a=a*a)
        if (b&1) c=c*a;
    return c;
}

//复数
struct Complex { double r,i; };
Complex operator +(Complex a,Complex b) { return (Complex){a.r+b.r,a.i+b.i}; }
Complex operator -(Complex a,Complex b) { return (Complex){a.r-b.r,a.i-b.i}; }
Complex operator *(Complex a,Complex b) { return (Complex){a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r}; }

//MTT
namespace MTT {
    Complex w[N]; int r[N];
    
    inline void FFT(Complex* A,int n,int op) {
        for (re int i=0;i<n;++i) if (i<r[i]) swap(A[i],A[r[i]]);
        for (re int i=1;i<n;i<<=1)
            for (re int j=0;j<n;j+=(i<<1))
                for (re int k=0;k<i;++k) {
                    Complex x=A[j+k],y=w[n/(i<<1)*k]*A[j+k+i];
                    A[j+k]=x+y,A[j+k+i]=x-y;
                }
        if (op==-1) for (re int i=0;i<n;++i) A[i].r=A[i].r/n+0.5;
    }
    
    inline void MTT(int* P,int n,int* Q,int m,int* ans) {
        const int LIM=32768;
        static Complex A[N],B[N],C[N],D[N],E[N],F[N],G[N],H[N];
        int l=0,lim=1;
        for (;lim<=n+m;lim<<=1,++l);
        for (re int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        for (re int i=0;i<lim;++i) {
            A[i].r=P[i]/LIM,A[i].i=0,B[i].r=P[i]%LIM,B[i].i=0;
            C[i].r=Q[i]/LIM,C[i].i=0,D[i].r=Q[i]%LIM,D[i].i=0;
        }
        for (re int i=0;i<lim;++i)
            w[i]=(Complex){cos(2*PI*i/lim),sin(2*PI*i/lim)};
        FFT(A,lim,1),FFT(B,lim,1),FFT(C,lim,1),FFT(D,lim,1);
        for (re int i=0;i<lim;++i)
            E[i]=A[i]*C[i],F[i]=A[i]*D[i],G[i]=B[i]*C[i],H[i]=B[i]*D[i];
        for (re int i=0;i<lim;++i)
            w[i]=(Complex){cos(2*PI*i/lim),-sin(2*PI*i/lim)};
        FFT(E,lim,-1),FFT(F,lim,-1),FFT(G,lim,-1),FFT(H,lim,-1);
        for (re int i=0;i<lim;++i) {
            ans[i]=0;
            add(ans[i],mul((ll)E[i].r%mod,(ll)LIM*LIM%mod));
            add(ans[i],mul((ll)F[i].r%mod,LIM%mod));
            add(ans[i],mul((ll)G[i].r%mod,LIM%mod));
            add(ans[i],(ll)H[i].r%mod);
        }
    }
}

int F[N],G[N],ans[N],w[N],c[N];

int main() {
    I[0][0]=I[1][1]=I[2][2]=1;
    
    n=read(),k=read(),L=read(),x=read()-1,y=read()-1,mod=read();
    
    w[0]=1,w[1]=qpow(getroot(mod),(mod-1)/k);
    for (re int i=2;i<k;++i) w[i]=mul(w[i-1],w[1]);
    
    for (re int i=0;i<n;++i)
        for (re int j=0;j<n;++j)
            A[i][j]=read();
    S[0][x]=1;
    for (re int i=0;i<k;++i) c[i]=(S*qpow(A*w[i]+I,L))[0][y];
    for (re int i=0;i<k+k;++i) F[i]=w[(k-1ll*i*(i-1)/2%k)%k];
    for (re int i=0;i<k;++i) G[i]=mul(c[i],w[1ll*i*(i-1)/2%k]);
    reverse(G,G+k+1); MTT::MTT(F,k+k,G,k,ans);
    
    for (re int i=0,invk=qpow(k,mod-2);i<k;++i)
        printf("%d\n",mul(mul(ans[i+k],invk),w[1ll*i*(i-1)/2%k]));
    return 0;
}
最后修改:2019 年 06 月 03 日 03 : 08 PM

发表评论