[数学记录]P5293 [HNOI2019]白兔之舞

· · 个人记录

题意 : 有一个 (L+1)\times n 的矩阵,位置由(i,j),(0\leq i\leq L,1\leq j\leq n) 表示。

给出权值数组 w[i][j]。

如果有 ui<vi ,一个点 (ui,uj) 到 (vi,vj) 的路径条数是 w[uj][vj]。

定义一支舞曲为从给定点 (0,x) 出发的一条任意路径。定义舞曲长度为经过的边的条数。

给定整数 k,y ,对于 t∈[0,k) ,分别询问有多少种舞曲使得长度 m=t\pmod k , 且终点的第二维坐标为 y。

答案对给定的质数 p 取模。n\leq 3,L\leq 10^8,k\leq65535,k|(p-1) , 时限 \texttt{2s}。

设 f[i][j] 表示舞曲长度为 i ,终点横坐标 j 的路径条数。

注意我们的行走方式允许随意跳过某些行,这不便于DP,我们考虑把行压紧。

设 g[i][j] 表示每次只能到下一行,舞曲长度为 i ,终点横坐标 j 的路径条数。

容易得知如下DP :

g[i][vj]=\sum\limits_{uj=1}^nw[uj][vj]g[i-1][uj]

可以把转移写成矩阵的形式,设向量 G_i 为 g[i][1...n]。

易得边界为 G_0[x]=1, 那么有 G_k=G_0W^k。

求得 g 之后, f 就相当于在 L 行中选取 i 行落脚,即 f[i][j]=\dbinom{L}{i}g[i][k]

{\rm Ans[t]}=\sum\limits_{i=1}^L[i\bmod k=t]f[i][y]

使用单位根反演。

=\sum\limits_{i=1}^L\dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{(i-t)j}f[i][y] \dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}\sum\limits_{i=1}^Lw_k^{ij}f[i][y] \dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}\sum\limits_{i=1}^Lw_k^{ij}\dbinom{L}{i}g[i][y]

使用 G_k=G_0W^k , 这样会得到一排向量,我们需要取出第 y 个。

[y]\dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}\sum\limits_{i=1}^Lw_k^{ij}\dbinom{L}{i}G_0W^i [y]\dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}G_0\sum\limits_{i=1}^L\dbinom{L}{i}(W·w_k^j)^i

发现是个二项式定理 (I和任意矩阵满足交换律)

[y]\dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}G_0(W·w_k^j+I)^L \dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}[y]G_0(W·w_k^j+I)^L

对于后面的 [y]G_0(W·w_k^j+I)^L ,可以矩阵快速幂求出。

前面的求和其实就是任意长度的 IDFT 变换,考虑使用 Bluestein 算法计算。

尝试使用 \dfrac{k^2+i^2-(k-i)^2}{2}=ik 但是 /2 需要 2k 次单位根,不一定存在。

考虑 \dbinom{i+j}{2}-\dbinom{i}{2}-\dbinom{j}{2}=ij , 发现没问题了。

具体来讲,如果需要计算 {\rm Ans}[t]=\dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{-tj}G[j]

=\dfrac{1}{k}\sum\limits_{j=0}^{k-1}w_k^{\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}}G[j] =\dfrac{w_k^\binom{t}{2}}{k}\sum\limits_{j=0}^{k-1}+w_k^\binom{j}{2}G[j]w_k^{-\binom{t+j}{2}}

能够发现就是差卷积。注意上界。

任意模数卷积需要MTT,求单位根需要原根。

复杂度O(k\log k+kn^3\log L+\sqrt{mod})。

个人觉得是一道不错的练习题,可是赛场上逼着选手写MTT有点毒瘤。

#include<algorithm>
#include<cstdio>
#include<cmath>
#define ld double
#define ll long long
#define Maxn 135000
using namespace std;
int mod;
ll powM(ll a,int t=mod-2){
  ll ret=1;
  while(t){
    if (t&1)ret=ret*a%mod;
    a=a*a%mod;t>>=1;
  }return ret;
}
namespace Grt
{
  int gcd(int a,int b)
  {return !b ? a : gcd(b,a%b);}
  int getphi(int n){
    int ret=n,p=2;
    while(p*p<=n){
      if (n%p==0){
        while(n%p==0)n/=p;
        ret=ret/p*(p-1);
      }p++;
    }if (n>1)ret=ret/n*(n-1);
    return ret;
  }
  int d[25],tn;
  void getft(int n){
    int sav=n,p=2;
    while(p*p<=n){
      if (n%p==0){
        while(n%p==0)n/=p;
        d[++tn]=sav/p;
      }p++;
    }if (n>1)d[++tn]=sav/n;
  }
  bool check(int n){
    for (int i=1;i<=tn;i++)
      if (powM(n,d[i])==1)
        return 0;
    return 1;
  }
  int get(){
    getft(getphi(mod));
    for (int i=1;;i++)
      if (gcd(i,mod)==1&&check(i))
        return i;
  }
}
const ld Pi=acos(-1);
struct CP{
  ld x,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};}
  CP operator * (const CP& B) const
  {return (CP){x*B.x-y*B.y,x*B.y+y*B.x};}
};
int tr[Maxn<<1];
void FFT(CP *f,int op,int n)
{
  static CP w[Maxn]={(CP){1.0,0.0}};
  for (int i=0;i<n;i++)
    if (i<tr[i])swap(f[i],f[tr[i]]);
  for(int l=1;l<n;l<<=1){
    CP tG=(CP){cos(Pi/l),sin(Pi/l)*op};
    for (int i=l-2;i>=0;i-=2)w[i+1]=(w[i]=w[i>>1])*tG;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        CP sav=w[p]*f[k|l|p];
        f[k|l|p]=f[k|p]-sav;
        f[k|p]=f[k|p]+sav;
      }
  }
}
void times(int *f,int *g,int n,int lim)
{
  static CP P1[Maxn<<1],P2[Maxn<<1],Q[Maxn<<1];
  for (int i=0,sav;i<lim;i++){
    P1[i]=(CP){f[i]>>15,f[i]&32767};
    P2[i]=(CP){f[i]>>15,-(f[i]&32767)};
    Q[i]=(CP){g[i]>>15,g[i]&32767};
  }for (int i=1;i<n;i++)
    tr[i]=tr[i>>1]>>1|((i&1)?n>>1:0);
  FFT(P1,1,n);FFT(P2,1,n);FFT(Q,1,n);
  for (int i=0;i<n;i++){
    Q[i].x/=n;Q[i].y/=n;
    P1[i]=P1[i]*Q[i];
    P2[i]=P2[i]*Q[i];
  }FFT(P1,-1,n);FFT(P2,-1,n);
  for (int i=0;i<lim;i++){
    ll a1b1=0,a1b2=0,a2b1=0,a2b2;
    a1b1=(ll)floor((P1[i].x+P2[i].x)/2+0.4)%mod;
    a1b2=(ll)floor((P1[i].y+P2[i].y)/2+0.4)%mod;
    a2b1=((ll)floor(P1[i].y+0.4)-a1b2)%mod;
    a2b2=((ll)floor(P2[i].x+0.4)-a1b1)%mod;
    f[i]=((((a1b1<<15)+(a1b2+a2b1))<<15)+a2b2)%mod;
    if (f[i]<0)f[i]+=mod;
  }
}
int n;
struct Mat
{
  int w[3][3];
  void clr(){
    for (int i=0;i<n;i++)
      for (int j=0;j<n;j++)
        w[i][j]=0;
  }
  Mat operator  * (const Mat &B) const{
    Mat R;R.clr();
    for (int k=0;k<n;k++)
      for (int i=0;i<n;i++)
        for (int j=0;j<n;j++)
          R.w[i][j]=(R.w[i][j]+1ll*w[i][k]*B.w[k][j])%mod;
    return R;
  }
  Mat operator * (const int &x) const{
    Mat R;R.clr();
    for (int i=0;i<n;i++)
      for (int j=0;j<n;j++)
        R.w[i][j]=1ll*w[i][j]*x%mod;
    return R;
  }
  Mat operator + (const Mat &B) const{
    Mat R;
    for (int i=0;i<n;i++)
      for (int j=0;j<n;j++)
        R.w[i][j]=(w[i][j]+B.w[i][j])%mod;
    return R;
  }
}T,I;
Mat powM(Mat a,int t=mod-2){
  Mat ret=I;
  while(t){
    if (t&1)ret=ret*a;
    a=a*a;t>>=1;
  }return ret;
}
int G,k,L,x,y,f[Maxn<<1],g[Maxn<<1],c[Maxn],pw[Maxn];
int main()
{
  scanf("%d%d%d%d%d%d",&n,&k,&L,&x,&y,&mod);
  x--;y--;
  G=powM(Grt::get(),(mod-1)/k);
  for (int i=0;i<n;i++)I.w[i][i]=1;
  for (int i=0;i<n;i++)
    for (int j=0;j<n;j++)
      scanf("%d",&T.w[i][j]);
  pw[0]=1;
  for (int i=1;i<=k;i++)pw[i]=1ll*pw[i-1]*G%mod;
  for (int i=0;i<k+k;i++){
    c[i]=1ll*i*(i-1)/2%k;
    g[i]=pw[k-c[i]];
  }for (int i=0;i<k;i++)
    f[i]=1ll*powM(T*powM(G,i)+I,L).w[x][y]*pw[c[i]]%mod;
  int len=1;while(len<4*k)len<<=1;
  reverse(g,g+k+k);times(f,g,len,k+k);reverse(f,f+k+k);
  int invk=powM(k);
  for (int i=0;i<k;i++)
    printf("%d\n",1ll*f[i]*pw[c[i]]%mod*invk%mod);
  return 0;
}