zoukankan      html  css  js  c++  java
  • 【loj3058】【hnoi2019】白兔之舞

    题意

    有一个((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} = egin{cases} 1 i mod n = 0 \ 0 i mod n eq 0 \ end{cases} \这个等比数列求和即可证明; \用这个来换掉式子中的取模:\ egin{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} imes 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_iomega_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)} imes 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
      ",now);
      	}
      	
      	return 0;
      }
      
  • 相关阅读:
    Oracle中有大量的sniped会话
    Error 1130: Host '127.0.0.1' is not allowed to connect to this MySQL server
    汉字转换为拼音以及缩写(javascript)
    高效率随机删除数据(不重复)
    vs2010 舒服背景 优雅字体 配置
    mvc中的ViewData用到webfrom中去
    jquery ajax return值 没有返回 的解决方法
    zShowBox (图片放大展示jquery版 兼容性好)
    动感效果的TAB选项卡 jquery 插件
    loading 加载提示······
  • 原文地址:https://www.cnblogs.com/Paul-Guderian/p/10787195.html
Copyright © 2011-2022 走看看