[倍增矩乘 FFT] LOJ#6275. 棋盘

考虑dp

$f_{i,s,j}$ 表示前 $i$ 列，最后一列的状态是 $s$ ，有 $j$ 个联通块的方案数

如果令多项式 $g_{i,s}(x)=\sum_{i=0}^n f_{i,s,i}x^i$，那么第二维的转移就是一个矩乘的形式

把多项式DFT一遍，每一位就可以单独做，直接倍增矩乘就好了

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#define copy(a,b) memcpy(a,b,sizeof(b))
#define set(a,b) memset(a,b,sizeof(a))

using namespace std;

const int N=500010,P=998244353;

typedef int poly[N];

int n,m,k,num,L,rev[N];
int w[2][N];

inline int Pow(int x,int y){
  int ret=1;
  for(;y;y>>=1,x=1LL*x*x%P) if(y&1) ret=1LL*ret*x%P;
  return ret;
}

inline void Pre(){
  int g=Pow(3,(P-1)/num);
  w[0][0]=w[1][0]=1;
  for(int i=1;i<num;i++) w[1][i]=1LL*w[1][i-1]*g%P;
  for(int i=1;i<num;i++) w[0][i]=w[1][num-i];
}

inline void FFT(int *a,int r){
  for(int i=1;i<num;i++) if(i>rev[i]) swap(a[i],a[rev[i]]);
  for(int i=1;i<num;i<<=1)
    for(int j=0;j<num;j+=(i<<1))
      for(int k=0;k<i;k++){
    int x=a[j+k],y=1LL*w[r][num/(i<<1)*k]*a[j+k+i]%P;
    a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
      }
  if(!r) for(int i=0,inv=Pow(num,P-2);i<num;i++) a[i]=1LL*a[i]*inv%P;
}

namespace Case1{
  typedef int mat[2][2];

  inline void mul(mat a,mat b,mat c){
    static mat tmp; set(tmp,0);
    for(int i=0;i<2;i++)
      for(int k=0;k<2;k++)
    for(int j=0;j<2;j++)
      tmp[i][j]=(tmp[i][j]+1LL*a[i][k]*b[k][j])%P;
    copy(c,tmp);
  }

  void solve(){
    static poly tmp0,tmp1,res; static mat cur,ret;
    tmp0[0]=1; tmp1[1]=1; FFT(tmp0,1); FFT(tmp1,1);
    set(res,0);
    for(int k=0;k<num;k++){
      ret[0][0]=ret[1][1]=1; ret[1][0]=ret[0][1]=0;
      cur[0][1]=tmp1[k]; cur[0][0]=cur[1][1]=cur[1][0]=tmp0[k];
      int c=m;
      for(;c;c>>=1,mul(cur,cur,cur)) if(c&1) mul(ret,cur,ret);
      res[k]=(ret[0][0]+ret[0][1])%P;
    }
    FFT(res,0); printf("%d\n",res[k]);
  }
}

namespace Case2{
  typedef int mat[3][3];

  inline void mul(mat a,mat b,mat c){
    static mat tmp; set(tmp,0);
    for(int i=0;i<3;i++)
      for(int k=0;k<3;k++)
    for(int j=0;j<3;j++)
      tmp[i][j]=(tmp[i][j]+1LL*a[i][k]*b[k][j])%P;
    copy(c,tmp);
  }

  void solve(){
    static poly tmp0,tmp1,tmp2,tmp3,tmp4,res; static mat cur,ret;
    tmp0[0]=1; FFT(tmp0,1);
    tmp1[1]=1; FFT(tmp1,1);
    for(int i=0;i<num;i++)
      tmp2[i]=2LL*tmp1[i]%P,tmp3[i]=2LL*tmp0[i]%P;
    tmp4[0]=1; tmp4[1]=1; FFT(tmp4,1);
    set(res,0);
    for(int k=0;k<num;k++){
      for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
      ret[i][j]=(i==j);
      cur[0][0]=tmp0[k]; cur[0][1]=tmp1[k]; cur[0][2]=tmp2[k];
      cur[1][0]=tmp0[k]; cur[1][1]=tmp0[k]; cur[1][2]=tmp3[k];
      cur[2][0]=tmp0[k]; cur[2][1]=tmp0[k]; cur[2][2]=tmp4[k];
      int c=m;
      for(;c;c>>=1,mul(cur,cur,cur)) if(c&1) mul(ret,cur,ret);
      for(int i=0;i<3;i++) res[k]=(res[k]+ret[0][i])%P;
    }
    FFT(res,0);
    printf("%d\n",res[k]);
  }
}

namespace Case3{           
  typedef int mat[7][7];

  inline void mul(mat a,mat b,mat c){
    static mat tmp; set(tmp,0);
    for(int i=0;i<7;i++)
      for(int k=0;k<7;k++)
    for(int j=0;j<7;j++)
      tmp[i][j]=(tmp[i][j]+1LL*a[i][k]*b[k][j])%P;
    copy(c,tmp);
  }

  #define t0 tmp[0][k]
  #define t1 tmp[1][k]
  #define t2 tmp[2][k]
  #define t3 tmp[3][k]
  #define t4 tmp[4][k]
  #define t5 tmp[5][k]
  #define t6 tmp[6][k]
  #define t7 tmp[7][k]

  void solve(){
    static poly tmp[8],res; mat cur,ret;
    tmp[0][num-1]=1;
    tmp[1][0]=1; tmp[2][1]=1; tmp[3][1]=2; tmp[4][2]=1;
    tmp[6][0]=2; tmp[7][0]=1; tmp[7][1]=1;
    set(res,0);
    FFT(tmp[0],1); FFT(tmp[1],1); FFT(tmp[2],1); FFT(tmp[4],1);
    for(int i=0;i<num;i++){
      tmp[3][i]=2LL*tmp[2][i]%P; tmp[6][i]=2LL*tmp[1][i]%P;
      tmp[7][i]=(tmp[1][i]+tmp[2][i])%P;
    }
    //for(int i=0;i<8;i++) FFT(tmp[i],1);
    for(int k=0;k<num;k++){
      for(int i=0;i<7;i++)
    for(int j=0;j<7;j++)
      ret[i][j]=(i==j);
      cur[0][0]=t1; cur[0][1]=t2; cur[0][2]=t3; cur[0][3]=t2; cur[0][4]=t3; cur[0][5]=t4; cur[0][6]=t5;
      cur[1][0]=t1; cur[1][1]=t1; cur[1][2]=t6; cur[1][3]=t1; cur[1][4]=t6; cur[1][5]=t5; cur[1][6]=t1;
      cur[2][0]=t1; cur[2][1]=t1; cur[2][2]=t7; cur[2][3]=t2; cur[2][4]=t7; cur[2][5]=t2; cur[2][6]=t5;
      cur[3][0]=t1; cur[3][1]=t1; cur[3][2]=t3; cur[3][3]=t1; cur[3][4]=t6; cur[3][5]=t4; cur[3][6]=t5;
      cur[4][0]=t1; cur[4][1]=t1; cur[4][2]=t7; cur[4][3]=t1; cur[4][4]=t6; cur[4][5]=t2; cur[4][6]=t5;
      cur[5][0]=t1; cur[5][1]=t0; cur[5][2]=t6; cur[5][3]=t2; cur[5][4]=t6; cur[5][5]=t1; cur[5][6]=t5;
      cur[6][0]=t1; cur[6][1]=t1; cur[6][2]=t6; cur[6][3]=t2; cur[6][4]=t6; cur[6][5]=t5; cur[6][6]=t1;

      int c=m;
      for(;c;c>>=1,mul(cur,cur,cur)) if(c&1) mul(ret,cur,ret);
      for(int i=0;i<7;i++) res[k]=(res[k]+ret[0][i])%P;
    }
    FFT(res,0); printf("%d\n",res[k]);
  }
}

int main(){
  scanf("%d%d%d",&n,&m,&k);
  if(k>(n*m+1)/2) return puts("0"),0;
  for(num=1;num<n*m;num<<=1,L++); num<<=1; Pre();
  for(int i=1;i<num;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
  if(n==1)
    Case1::solve();
  else if(n==2)
    Case2::solve();
  else
    Case3::solve();
  return 0;
}