【loj3058】【hnoi2019】白兔之舞-优快云博客

本文详细解析了一个涉及单位根反演和快速傅里叶变换（FFT）的组合计数问题，通过将问题转化为矩阵幂运算，利用单位根反演技巧和快速傅里叶变换进行高效计算，解决了大规模网格上白兔移动路径计数的问题。

题意

有一个$(L+1)*n$ 的网格图，初始时白兔在$(0,X)$ , 每次可以向横坐标递增，纵坐标随意的位置移动，两个位置之间的路径条数只取决于纵坐标，用$w(i,j)$ 表示，如果要求白兔停下的点纵坐标为$Y$ 依次输出移动的步数对$k$ 取模为 $0 - k -1 $的方案数；

$p$为质数且$10^8 \lt p \lt 2^{30} , 1 \le n \le 3 , 1 \le x , y \le n , 0 \le w(i,j) \lt p , 1 \le k \le 65536 , \le L \le 10^8 $，

满足$k$是$p-1$的一个约数；

题解

即求：$ f_r = \sum_{i=0}^{L}(^L_i) (A^i)_{X,Y} [ i mod k = r] $
$Part 1 $
单位根反演：

\[根据求和引理：\\ \frac{1}{n} \sum_{j=0}^{n-1}\omega_n^{ij} = \begin{cases} 1 \ \ \ i \ mod \ n = 0 \\ 0 \ \ \ i \ mod \ n \neq 0 \\ \end{cases} \\这个等比数列求和即可证明； \\用这个来换掉式子中的取模：\\ \begin{align} &=\sum_{i=0}^{L}(^L_i)A^i_{x,y}\frac{1}{k}\sum_{j=0}^{k-1} \omega_k^{(i-r)j}\\ &=\frac{1}{k}\sum_{i=0}^{k-1}\omega_k^{-ir}\sum_{j=0}^{L}(^L_j)(A^j)_{x,y}\times \omega^{ij}_k\\ &=\frac{1}{k}\sum_{i=0}^{k-1}\omega_k^{-ir}(w^i_kA+I)^L_{x,y}\\ &=\frac{1}{k}\sum_{i=0}^{k-1}a_i\omega_k^{-ir} \end{align}\]

当$k = 2^x$次方时，可以用裸的$fft$求出答案;
$Part 2$
但是没有必要，当 $k$ 为任意数，由于$k \ | \ mod - 1$，所以$\omega_k$可以用原根的$(mod-1)/k$来代替;
同时注意到： $ ab = (^{a+b}_2) - (^a_2) - (^b_2) $
$ f_r = \frac{1}{k}w_k^{(^r_2)} \sum_{i=0}^{k-1}a_i \omega _k^{(^i_2)} \times \omega_k^{-(^{i+r}_2)} $

直接reverse+mtt即可;

//注意mtt的精度;
#include<bits/stdc++.h>
#define ll long long  
#define ld double
using namespace std;
const int N=1<<20;
int n,k,l,X,Y,mod,a[N],b[N],A[N],B[N],C[N],fac[N],tot;
char gc(){
  static char*p1,*p2,s[1000000];
  if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin);
  return(p1==p2)?EOF:*p1++;
}
int rd(){
  int x=0;char c=gc();
  while(c<'0'||c>'9')c=gc();
  while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+c-'0',c=gc();
  return x;
}
int pw(int x,int y){
  int re=1;
  while(y){
      if(y&1)re=(ll)re*x%mod;
      y>>=1;x=(ll)x*x%mod;
  }
  return re;
}
int find_rt(int p){
  int tmp=p-1;
  for(int i=2;i*i<=tmp;++i)if(tmp%i==0){
      fac[++tot]=i;
      while(tmp%i==0)tmp/=i;
  }
  if(tmp!=1)fac[++tot]=tmp;
  for(int i=2;;++i){
      int fg=0;
      for(int j=1;j<=tot;++j)if(pw(i,(p-1)/fac[j])==1){fg=1;break;}
      if(!fg)return i;
  }
  return 0;
}
struct Mat{
  int a[3][3];
  void init(){
      memset(a,0,sizeof(a));
      for(int i=0;i<n;++i)a[i][i]=1;
  }
  Mat operator *(const int&A){
      Mat re;
      for(int i=0;i<n;++i)
      for(int j=0;j<n;++j)re.a[i][j]=(ll)a[i][j]*A%mod;
      return re;
  }
  Mat operator *(const Mat&A){
      Mat re;
      for(int i=0;i<n;++i)
      for(int j=0;j<n;++j){
          re.a[i][j]=0;
          for(int k=0;k<n;++k){
              re.a[i][j]+=(ll)a[i][k]*A.a[k][j]%mod;
              if(re.a[i][j]>=mod)re.a[i][j]-=mod;
          }
      }
      return re;
  }
  Mat operator +(const Mat&A){
      Mat re;
      for(int i=0;i<n;++i)
      for(int j=0;j<n;++j){
          re.a[i][j]=a[i][j]+A.a[i][j];
          if(re.a[i][j]>=mod)re.a[i][j]-=mod;
      }
      return re;
  }
}mp,I;
int pw(Mat x,int y){
  Mat re;re.init();
  while(y){
      if(y&1)re=re*x;
      y>>=1;x=x*x;
  }
  return re.a[X][Y];
}
struct com{
  ld x,y;
  com(ld _x=0,ld _y=0):x(_x),y(_y){};
  com operator +(const com&A)const{return com(x+A.x,y+A.y);}
  com operator -(const com&A)const{return com(x-A.x,y-A.y);}
  com operator *(const com&A)const{return com(x*A.x-y*A.y,x*A.y+y*A.x);}
  com operator /(const ld&A)const{return com(x/A,y/A);}
  com operator !()const{return com(x,-y);}
};
namespace poly{
  const ld pi=acos(-1);
  int L,len,rev[N];com Wn[N];
  void init(int l){
      for(L=0,len=1;len<=l<<1;len<<=1,++L);
      for(int i=0;i<len;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
      //Wn[0]=com(1,0);Wn[1]=com(cos(pi/len),sin(pi/len));
      //for(int i=2;i<len;++i)Wn[i]=Wn[i-1]*Wn[1];
      for(int i=0;i<len;++i)Wn[i]=com(cos(pi*i/len),sin(pi*i/len));
  }
  void fft(com*A,int f){
      for(int i=0;i<len;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
      for(int i=1;i<len;i<<=1){
          for(int j=0;j<len;j+=i<<1){
              for(int k=0;k<i;++k){
                  com w=~f?Wn[len/i*k]:!Wn[len/i*k];//
                  com x=A[j+k],y=w*A[j+k+i];
                  A[j+k]=x+y;A[j+k+i]=x-y;
              }
          }
      }
      if(!~f)for(int i=0;i<len;++i)A[i]=A[i]/len;
  }
  void mtt(int*A,int*B,int*C){
      static com t1[N],t2[N],t3[N],t4[N];
      int all=(1<<15)-1;
      for(int i=0;i<len;++i){
          t1[i]=com(A[i]>>15,A[i]&all); 
          t2[i]=com(B[i]>>15,B[i]&all);
      }
      fft(t1,1);fft(t2,1);
      for(int i=0;i<len;++i){
          com x1=t1[i],y1=!t1[len-i&len-1];
          com x2=t2[i],y2=!t2[len-i&len-1];
          t3[i] = (x1+y1)*x2*com(0.5,0);//(x1+y1)/2(x2+y2)/2 + (x1+y1)/2(x2-y2)/2i * i ; 
          t4[i] = (x1-y1)*x2*com(0,-0.5);//(x1-y1)/2i(x2+y2)/2 + (x1-y1)/2i(x2-y2)/2i * i ; 
      }
      fft(t3,-1);fft(t4,-1);
      for(int i=0;i<len;++i){
          C[i] = (((ll)(t3[i].x+0.5)%mod <<30)%mod + 
                  ((ll)(t3[i].y+0.5)%mod <<15)%mod + 
                  ((ll)(t4[i].x+0.5)%mod <<15)%mod + 
                   (ll)(t4[i].y+0.5)%mod ) %mod ;
      }
  }
}
int main(){
//    freopen("dance.in","r",stdin);
//    freopen("dance.out","w",stdout);

  n=rd();k=rd();l=rd();X=rd()-1;Y=rd()-1;mod=rd();I.init();
  for(int i=0;i<n;++i)for(int j=0;j<n;++j)mp.a[i][j]=rd();
  int G=find_rt(mod),wk=pw(G,(mod-1)/k);

  for(int i=0,w=1;i<k;++i,w=(ll)w*wk%mod)a[i]=pw((mp*w+I),l);
  for(int i=0;i<=k-1<<1;++i)b[i]=pw(wk,(ll)i*(i-1)/2%(mod-1));
  for(int i=0;i<k;++i)A[i]=(ll)a[i]*b[i]%mod;
  for(int i=0;i<=k-1<<1;++i)B[i]=pw(b[(k-1<<1)-i],mod-2);

  poly::init(k-1<<1);
  poly::mtt(A,B,C);
  int iv=pw(k,mod-2);
  for(int i=0;i<k;++i){
      int now = (ll)iv*b[i]%mod*C[(k-1<<1)-i]%mod;
      printf("%d\n",now);
  }

  return 0;
}

转载于:https://www.cnblogs.com/Paul-Guderian/p/10787195.html