zoukankan      html  css  js  c++  java
  • #3058. 「HNOI2019」白兔之舞(单位根反演)

    题目描述

    https://loj.ac/problem/3058

    单位根反演

    因为ω太难写了所以用w代替

    ([n|k]=frac{1}{n}sum_{i=0}^{n-1} w_n^{ik})

    证明:

    当n|k时显然是1,否则(frac{w_n^{nk}-1}{w^n-1}=0)

    题解

    一开始想矩乘存多项式然后快速幂循环卷积,然后多乘了一个log直接炸飞

    先暴力找p-1的原根,然后即可求出k次单位根w,A是读入的矩阵

    (ans_t=sum_{i=0}^L [k|(i-t)]A_{x,y}^i inom{L}{i})

    (=frac{1}{k}sum_{i=0}^L sum_{j=0}^{k-1}w^{j(i-t)} A_{x,y}^i inom{L}{i})

    (=frac{1}{k}sum_{j=0}^{k-1}w^{-jt} sum_{i=0}^L w^{ij} A_{x,y}^i inom{L}{i})

    (=frac{1}{k}sum_{j=0}^{k-1}w^{-jt} (w^jA+1)^L_{x,y})

    (=frac{1}{k}sum_{j=0}^{k-1}w^{inom{j}{2}}w^{inom{t}{2}}w^{-inom{j+t}{2}} (w^jA+1)^L_{x,y})

    把j变成k-1-j即可卷积,要写mtt

    注意mtt求a-bi的点值时的关系是共轭而不是相等

    code

    #include <bits/stdc++.h>
    #define fo(a,b,c) for (a=b; a<=c; a++)
    #define fd(a,b,c) for (a=b; a>=c; a--)
    #define add(a,b) a=((a)+(b))%mod
    #define ll long long
    #define D 10000
    //#define file
    using namespace std;
    
    int N,len,n,K,L,x,y,mod,Mod,i,j,k,l;
    ll a[262144],b[262144],g,w,s;
    
    struct type{double x,y;} AA[262144];
    type operator + (type a,type b) {return {a.x+b.x,a.y+b.y};}
    type operator - (type a,type b) {return {a.x-b.x,a.y-b.y};}
    type operator * (type a,type b) {return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
    struct mat{ll a[3][3]; void clear() {memset(a,0,sizeof(a));}} I,A,B;
    mat operator * (mat a,mat b)
    {
    	static mat c;
    	int i,j,k;
    	
    	fo(i,0,n-1)
    	{
    		fo(j,0,n-1)
    		{
    			c.a[i][j]=0;
    			fo(k,0,n-1)
    			add(c.a[i][j],a.a[i][k]*b.a[k][j]);
    		}
    	}
    	return c;
    }
    mat Qpower(mat a,ll b) {mat ans=I; while (b) {if (b&1) ans=ans*a;a=a*a;b>>=1;} return ans;}
    
    ll qpower(ll a,ll b) {ll ans=1; while (b) {if (b&1) ans=ans*a%mod;a=a*a%mod;b>>=1;} return ans;}
    void dft(type *a,int tp)
    {
    	int i,j,k,l,S=N,s1=2,s2=1;
    	type w,W,u,v;
    	
    	fo(i,0,N-1)
    	{
    		j=i,k=0;
    		fo(l,1,len)
    		k=k*2+(j&1),j>>=1;
    		AA[i]=a[k];
    	}
    	memcpy(a,AA,N*sizeof(type));
    	
    	fo(i,1,len)
    	{
    		S>>=1;
    		fo(j,0,S-1)
    		{
    			fo(k,0,s2-1)
    			{
    				W={cos(2*M_PI*k/s1),sin(2*M_PI*k/s1)*tp};
    				u=a[j*s1+k],v=a[j*s1+k+s2]*W;
    				a[j*s1+k]=u+v;
    				a[j*s1+k+s2]=u-v;
    			}
    		}
    		s1<<=1,s2<<=1;
    	}
    }
    void mul(ll *a,ll *b)
    {
    	static type A[262144],B[262144],P[262144],Q[262144];
    	static ll ans[262144];
    	int i,j,k,l;
    	
    	fo(i,0,N-1) A[i]={a[i]/D,a[i]%D},B[i]={b[i]/D,b[i]%D};
    	dft(A,1),dft(B,1);
    	fo(i,0,N-1) P[i]=A[i]*B[i],Q[i]=type{A[(N-i)%N].x,-A[(N-i)%N].y}*B[i];
    	dft(P,-1),dft(Q,-1);
    	fo(i,0,N-1) P[i].x/=N,Q[i].x/=N,P[i].y/=N,Q[i].y/=N;
    	
    	fo(i,0,N-1)
    	ans[i]=(((ll)floor((P[i].x+Q[i].x)/2+0.5)%mod*D%mod*D)+(ll)floor(P[i].y+0.5)*D+(ll)floor((Q[i].x-P[i].x)/2+0.5))%mod;
    	memcpy(a,ans,N*8);
    }
    
    bool pd(ll g)
    {
    	int i,j,k,l,s=floor(sqrt(mod-1));
    	fo(i,1,s)
    	if (!((mod-1)%i))
    	{
    		if (qpower(g,i)==1) return 0;
    		if (i>1 && qpower(g,(mod-1)/i)==1) return 0;
    	}
    	return 1;
    }
    
    int main()
    {
    	#ifdef file
    	freopen("loj3058.in","r",stdin);
    	#endif
    	
    	scanf("%d%d%d%d%d%d",&n,&K,&L,&x,&y,&mod),--x,--y,len=ceil(log2(K*3)),N=pow(2,len),Mod=mod-2;
    	I.a[0][0]=I.a[1][1]=I.a[2][2]=1;
    	fo(i,0,n-1)
    	{
    		fo(j,0,n-1)
    		scanf("%lld",&A.a[i][j]);
    	}
    	for (g=1; !pd(g); ++g);
    	w=qpower(g,(mod-1)/K);
    	
    	fo(i,0,K*2-1) a[i]=qpower(qpower(w,1ll*i*(i-1)/2),Mod);
    	fo(i,0,K-1)
    	{
    		B=A,s=qpower(w,i);
    		fo(j,0,n-1)
    		{
    			fo(k,0,n-1)
    			B.a[j][k]=B.a[j][k]*s%mod+(j==k);
    		}
    		B=Qpower(B,L);
    		b[(K-1)-i]=qpower(w,1ll*i*(i-1)/2)*B.a[x][y]%mod;
    	}
    	mul(a,b);
    	fo(i,0,K-1)
    	printf("%lld
    ",(a[i+(K-1)]*qpower(K,Mod)%mod*qpower(w,1ll*i*(i-1)/2)%mod+mod)%mod);
    	
    	fclose(stdin);
    	fclose(stdout);
    	return 0;
    }
    
  • 相关阅读:
    PAT 1010. 一元多项式求导 (25)
    PAT 1009. 说反话 (20) JAVA
    PAT 1009. 说反话 (20)
    PAT 1007. 素数对猜想 (20)
    POJ 2752 Seek the Name, Seek the Fame KMP
    POJ 2406 Power Strings KMP
    ZOJ3811 Untrusted Patrol
    Codeforces Round #265 (Div. 2) 题解
    Topcoder SRM632 DIV2 解题报告
    Topcoder SRM631 DIV2 解题报告
  • 原文地址:https://www.cnblogs.com/gmh77/p/13819895.html
Copyright © 2011-2022 走看看