zoukankan      html  css  js  c++  java
  • uoj#272. 【清华集训2016】石家庄的工人阶级队伍比较坚强(矩阵+三维FWT)

    传送门

    题解

    抄代码(20)分钟,搞懂题解在干嘛仨小时→_→

    到今天才算真正搞明白(FWT)在干吗了

    本题

    首先转移关系都是恒定的,设它为一个矩阵(B),那么要求的就是(f_n=f_0B^n)

    定义三进制不退位减法(ominus),三进制不退位加法(oplus),这两个互为逆运算(可以类比一下二进制的异或)

    根据(B)矩阵的定义以及剪刀石头布的性质易知(forall k lt 3^m,B_{ioplus k,joplus k}=B_{i,j})

    那么易知(B_{i,j}=B_{iominus j,0})。不难发现这个结论也可以推广到(B^n),即有(forall k lt 3^m,B^n_{ioplus k,joplus k}=B^n_{i,j}),可以用归纳法证,也可以感性理解一下

    那么对于(f_n=f_0B^n),它的第(i)项就为$$f_{n,i}=sum_k f_{0,k}B^n_{k,i}=sum_k f_{0,k}B^n_{0,iominus k}=sum_{xoplus y=i}f_{0,x}B^n_{0,y}$$
    于是这玩意儿可以写成一个卷积的形式

    我们发现上面的式子只需要(B)的第一行就够了,那么因为$$B^n_{0,i}=sum_k B^n_{0,k}B_{k,i}=sum_k B^n_{0,k}B_{0,iominus k}=sum_{xoplus y=i} B^n_{0,x}B_{0,y}$$
    那么(B)的第一行做(n)次卷积,然后(f)再和它做一次卷积就是答案了

    卷积

    我们现在需要找到一个三进制不退位加法的卷积变换。据LargestJN大佬说,如果下标运算是模意义下加法,这个卷积就叫做循环卷积

    然而就算它叫鸡卷也没用我们还是不会求

    首先一个卷积就是要满足(T(a)T(b)=T(aoplus b))(这里(oplus)为任意运算),设(T)为这个卷积的矩阵,(x_{i,j})为这个矩阵中的某个元素,那么上面的形式就可以表现为$$(sum_{j=0}^{n-1}x_{i,j}a_j)(sum_{k=0}^{n-1}x_{i,k}a_k)=sum_{j=0}^{n-1}x_{i,j}sum_{koplus t=j}a_kb_t$$
    于是考虑(a_k)(b_t)的贡献,得有$$x_{i,k}x_{i,t}a_kb_t=x_{i,koplus t}a_k,b_t$$
    那么就是对于(T)的每一行都得有$$x_{k}x_{t}=x_{koplus t}$$
    于是(T)的每一行都是方程组(x_{k}x_{t}=x_{koplus t})的一组解

    然而只有卷积万一没有卷积的逆变换(管它叫鸡卷好了)就gg了

    所以我们的(T)还得有逆矩阵,那么根据线性代数的芝士,(T)的行列式不为(0),那么方程组(x_{k}x_{t}=x_{koplus t})至少要有(n)组不同的解

    先来看看循环卷积的一般形式,长这样$$A imes B=sum_{i}sum_{(j+k)mod n=i}A_jB_k$$
    (T)满足(x_ix_j=x_{(i+j)mod n})

    然后我们发现(n)次单位根(omega^i)就是一组可行的解

    不知道(n)次单位根的可以回去看看(FFT)的原理

    然而需要(n)组解,于是矩阵长这样

    逆矩阵为$$frac{1}{n} egin{bmatrix}
    1& 1 & 1& ... & 1
    1& w_n^{-1}& w_n^{-2}& ... & w_n^{-(n - 1)}
    1& w_n^{-2} & w_n^{-4}& ... & w_n^{-2(n- 1)}
    ...& ...& ...& ...& ...
    1& w_n^{-(n - 1)}& w_n^{-2(n - 1)} & ... & w_n^{-(n - 1)(n - 1)}
    end{bmatrix}$$

    回到本题

    因为是三进制不进位加法,所以就选三次单位根(omega)。这里可以把所有的复数都表示为(a+bomega)的形式,那么因为三次单位根有(omega^2+omega+1=0),所以有(omega^2=-omega-1),所以所有的系数都是(0/1/-1),运算过程中就不用取模了。两个复数的乘法就是$$(a+bomega)(c+domega)=ac+(ad+bc)omega+bd(-omega-1)=(ac-bd)+(ad+bc-bd)omega$$
    然后根据上面,$$T= left[ egin{matrix} 1 & 1 & 1 1 & omega & omega^2 1 & omega^2 & omega end{matrix} ight]$$
    以及我们知道$$left[ egin{matrix} 1 & 1 & 1 1 & omega & omega^2 1 & omega^2 & omega end{matrix} ight] imes left[ egin{matrix} 1 & 1 & 1 1 & omega^2 & omega 1 & omega & omega^2 end{matrix} ight]=3E$$
    其中(E)是单位矩阵

    那么我们就可以把后面那个当做逆矩阵,带进去算

    因为(DFT)(IDFT)本质都是分治做向量对矩阵的乘法,它(IDFT)的时候每一次都多乘了个(3),总共有(log_3 n)层,那么总共多乘了(3^{log_3 n}=n)次,那么只要最后把所有元素都除以一个(n)就吼啦!

    奇怪的模数问题……

    最后是一个比较奇怪的模数问题,因为我们最后要除以(3^m),要求(3^m)有逆元,也就是(3 mid p),然后因为它这个奇怪的模数(p),假设(k=frac{p}{3})是个正整数,那么有$$frac{1}{k + 1} + frac{1}{k(k + 1)} = frac{1}{k} = frac{3}{p}$$
    矛盾,于是(3 mid p)

    然后……就真的没然后了……

    //minamoto
    #include<cstdio>
    #include<cstring>
    #include<map>
    #define R register
    #define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
    #define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
    #define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
    using namespace std;
    char buf[1<<21],*p1=buf,*p2=buf;
    inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
    int read(){
        R int res,f=1;R char ch;
        while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
        for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
        return res*f;
    }
    char sr[1<<21],z[20];int C=-1,Z=0;
    inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
    void print(R int x){
        if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
        while(z[++Z]=x%10+48,x/=10);
        while(sr[++C]=z[Z],--Z);sr[++C]='
    ';
    }
    const int N=6e5+5;
    int lim,m,P,t,b[25][25],a[N],inv2,invn;
    inline int add(R int x){return x>=P?x-P:x;}
    inline int dec(R int x){return x<0?x+P:x;}
    inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
    void exgcd(int &x,int &y,int a,int b){
    	if(!b)return (void)(x=1,y=0);
    	exgcd(y,x,b,a%b),y-=a/b*x;
    }
    inline int inv(int a){
    	int x,y;exgcd(x,y,a,P);
    	return (x%P+P)%P;
    }
    struct complex{
    	int x,y;
    	complex(int X=0,int Y=0):x(X),y(Y){}
    	inline complex operator +(const complex &b){return complex(add(x+b.x),add(y+b.y));}
    	inline complex operator -(const complex &b){return complex(dec(x-b.x),dec(y-b.y));}
    	inline complex operator *(const complex &b){
    		return complex(dec(mul(x,b.x)-mul(y,b.y)),dec(add(mul(x,b.y)+mul(y,b.x))-mul(y,b.y)));
    	}
    	inline bool operator <(const complex &b)const{
    		return x==b.x?y<b.y:x<b.x;
    	}
    }f[N],g[N];map<complex,complex>mp;
    complex ksm(complex x,R int y){
    	complex res(1,0);
    	for(;y;y>>=1,x=x*x)if(y&1)res=res*x;
    	return res;
    }
    complex calc1(complex b){return complex(dec(-b.y),dec(b.x-b.y));}
    complex calc2(complex b){return complex(dec(b.y-b.x),dec(-b.x));}
    void DFT(complex *A){
    	for(R int mid=1;mid<lim;mid*=3)
    		for(R int j=0;j<lim;j+=mid*3)
    			for(R int k=0;k<mid;++k){
    				complex x=A[j+k],y=A[j+k+mid],z=A[j+k+(mid<<1)];
    				A[j+k]=x+y+z;
    				A[j+k+mid]=x+calc1(y)+calc2(z);
    				A[j+k+(mid<<1)]=x+calc2(y)+calc1(z);
    			}
    }
    void IDFT(complex *A){
    	for(R int mid=1;mid<lim;mid*=3)
    		for(R int j=0;j<lim;j+=mid*3)
    			for(R int k=0;k<mid;++k){
    				complex x=A[j+k],y=A[j+k+mid],z=A[j+k+(mid<<1)];
    				A[j+k]=x+y+z;
    				A[j+k+mid]=x+calc2(y)+calc1(z);
    				A[j+k+(mid<<1)]=x+calc1(y)+calc2(z);
    			}
    	for(R int i=0;i<lim;++i)A[i].x=mul(A[i].x,invn);
    }
    int main(){
    //	freopen("testdata.in","r",stdin);
    	m=read(),t=read(),P=read();
    	lim=1;fp(i,1,m)lim*=3;
    	if(P==1){
    		fp(i,0,lim-1)print(0);
    		return Ot(),0;
    	}
    	fp(i,0,lim-1)a[i]=read();
    	fp(i,0,m)fp(j,0,m-i)b[i][j]=read();
    	fp(i,0,lim-1){
    		int tmp=i,cntw=0,cntl=0;
    		while(tmp){
    			int k=tmp%3;
    			k==1?++cntw:k==2?++cntl:0;
    			tmp/=3;
    		}g[i]=b[cntw][cntl],f[i]=a[i];
    	}
    	inv2=inv(2),invn=inv(lim);
    	DFT(f),DFT(g);
    	fp(i,0,lim-1)f[i]=f[i]*(mp.count(g[i])?mp[g[i]]:(mp[g[i]]=ksm(g[i],t)));
    	IDFT(f);
    	fp(i,0,lim-1)print(f[i].x);
    	return Ot(),0;
    }
    
  • 相关阅读:
    自我介绍+软工5问
    团队展示&选题
    团队展示&选题 (白衣天使队)
    结对编程(-java-实现)
    个人项目wc(Java)
    自我介绍+软工五问
    软件工程结课作业
    第四次博客作业-结对项目
    软件工程第三次作业——关于软件质量保障初探
    20194613 自动生成四则运算第一版报告
  • 原文地址:https://www.cnblogs.com/bztMinamoto/p/10239812.html
Copyright © 2011-2022 走看看