【矩阵乘法】CDOJ1610 黑红梅方

#include<cstdio>
#include<vector>
#include<iostream>
using namespace std;
typedef long long ll;
#define MOD 1000000009ll
typedef vector<ll> vec;
typedef vector<vec> mat;
mat I;
mat operator * (const mat &A,const mat &B){
	mat C(A.size(),vec(B[0].size()));
	for(int i=0;i<A.size();++i){
		for(int k=0;k<B.size();++k){
			for(int j=0;j<B[0].size();++j){
				C[i][j]=(C[i][j]+A[i][k]*B[k][j]%MOD)%MOD;
			}
		}
	}
	return C;
}
mat Quick_Pow(mat a,ll p)
{
    if(!p){
    	return I;
    }
    mat res=Quick_Pow(a,p>>1ll);
    res=res*res;
    if((p&1ll)==1ll){
    	res=res*a;
    }
    return res;
}
ll Quick_Pow(ll a,ll p)
{
    if(!p){
    	return 1;
    }
    ll res=Quick_Pow(a,p>>1ll);
    res=res*res%MOD;
    if((p&1ll)==1ll){
    	res=(a%MOD*res)%MOD;
    }
    return res;
}
ll n;
int main(){
	cin>>n;
	I.assign(64,vec(64));
	for(int i=0;i<64;++i){
		for(int j=0;j<64;++j){
			I[i][j]=0;
		}
	}
	for(int i=0;i<64;++i){
		I[i][i]=1;
	}
	mat A(64,vec(1));
	for(int i=0;i<64;++i){
		A[i][0]=1;
	}
	mat b(64,vec(64));
	for(int i=0;i<64;++i){
		for(int j=0;j<64;++j){
			b[i][j]=0;
		}
	}
	for(int i=0;i<4;++i){
		for(int j=0;j<4;++j){
			for(int k=0;k<4;++k){
				int hang=i*16+j*4+k;
				if(i==j || j==k || i==k){
					for(int l=0;l<4;++l){
						int lie=l*16+i*4+j;
						b[hang][lie]=1;
					}
				}
				else{
					for(int l=0;l<4;++l){
						if(l==i || l==j || l==k){
							int lie=l*16+i*4+j;
							b[hang][lie]=1;
						}
					}
				}
			}
		}
	}
//	for(int i=0;i<64;++i){
//		for(int j=0;j<64;++j){
//			printf("%I64d ",b[i][j]);
//		}
//		puts("");
//	}
	mat c=Quick_Pow(b,n-3ll)*A;
	ll n4=Quick_Pow(4ll,n);
	ll tmp=0;
	for(int i=0;i<64;++i){
		tmp=(tmp+c[i][0])%MOD;
	}
	cout<<(n4-tmp+MOD)%MOD<<endl;
	return 0;
}