「オイラープロジェクト321~330」の編集履歴(バックアップ)一覧に戻る

オイラープロジェクト321~330 - (2012/09/20 (木) 17:05:34) の編集履歴(バックアップ)


問い324

3x3xn の塔を 2x1x1 のブロックで埋める場合の数を f(n) とする。
ブロックは好きな方向に回転することができる;しかし、塔自身を回転・反転等させたものは区別して数える。
たとえば(q=100000007 として):
f(2) = 229,
f(4) = 117805,
f(10) mod q = 96149360,
f(10^3) mod q = 24806056,
f(10^6) mod q = 30808124.
f(10^10000) mod 100000007 を求めよ。



何か漸化式があるのだろうか?
数式を思いつかないのでとりあえず動的計画法で解こうとしてみる。

2段を作り、2段を2つ結合で4段、4段を2つ結合して8段、8段を2つ結合して16段、、として10^10000段まで到達し後は10^10000を2進数に分解2進数で1のついた段だけを掛け算していけば理論上計算できるはず。
多分相当遅いが自分はこれしか思いつかない。
10億段くらいの小さい数だったら簡単なのにな。

とりあえず小さな数で正しい答えを出せるところまで持ってきた。
最初とても遅い動的計画法から始め何とか高速化して10^3や10^6段で正しい答えを出せるところまで前進した。
だがまだ遅いとても遅い。
10^6段を計算するのに役10分。
10^50段程度なら80分、10^10000段、2進数で3万数千ケタを取り扱うと16000分も計算時間がいることになる。
時間にすると約270時間だ。
これでは話にならない。

色々考えて数値をピーピングしてみた結果16段以上からは組み合わせの集合が5種類に分割されているような感じがする。

5つのまとまりに分かれていて、まとまり内部で組み合わせの和をとる感じである。
一度5つにわかれるとまとまりはまとまり内部のものとしか掛け算の関係を持たない?



#include<stdio.h>
#include<map>
#include<set>
#include<string.h>
#include<iostream>
#include<vector>
const __int64 base=100000007;
struct S{
int blocks[2][3][3];//上面、下面のみをもってブロックを判別する
//立てた棒は11、横か縦の棒は1
//上面は隙間を許さない
bool operator<(const S& s)const{
	for(int h=0;h<2;h++){
		for(int r=0;r<3;r++){
			for(int c=0;c<3;c++){
				if(blocks[h][r][c]!=s.blocks[h][r][c])return blocks[h][r][c]<s.blocks[h][r][c];
			}
		}
	}
	return false;
}
};
std::map<S,__int64> memo,nextMemo,ansMemo,ansNext;
std::map<S,std::vector<S> > okPerms,nextOKPerms;
std::vector<S> vec;
bool isConnectOK(S& s1,S& s2){
//s1が下のブロック、s2が上のブロック
int upSide,downSide;
for(int r=0;r<3;r++){
	for(int c=0;c<3;c++){
		upSide=s1.blocks[1][r][c];//接合部下面
		downSide=s2.blocks[0][r][c];//接合部上面
		if(downSide==0&&upSide!=11)return false;//結合部上面に隙間があった
		if(downSide!=0&&upSide==11)return false;//結合部上面と下面に出っ張りがあった
	}
}
//ぴったり結合出来た
return true;
}
S connectBlock(S& s1,S& s2){
S s3;
//結合可能なので結合し上面と下面をもって種類を判別する
memcpy(s3.blocks[1],s2.blocks[1],sizeof(s2.blocks[1]));
memcpy(s3.blocks[0],s1.blocks[0],sizeof(s1.blocks[0]));
return s3;
}
void first(int deep,int perm[2][3][3],int sukima,int hBar){
//ブロックの組み合わせ2段を作る
//上面に隙間は許されない
if(deep>=18){
	S s;
	memcpy(s.blocks,perm,sizeof(s.blocks));
	if(memo.find(s)==memo.end()){
		memo[s]=1;
	}else{
		memo[s]++;
	}
	return ;
}else{
	int h=deep/9;
	int r=(deep%9)/3;
	int c=deep%3;
	//この部分にはもうブロックがあった
	if(perm[h][r][c]!=0){
		first(deep+1,perm,sukima,hBar);
		return ;
	}
	//ここを隙間にする、隙間が許されるのは下面のみ
	if(h==0){
		first(deep+1,perm,sukima+1,hBar);
	}
	
	//横
	perm[h][r][c]=1;
	if(c+1<3&&perm[h][r][c+1]==0){
		perm[h][r][c+1]=1;
		first(deep+1,perm,sukima,hBar);
		perm[h][r][c+1]=0;
	}
	//縦
	
	if(r+1<3&&perm[h][r+1][c]==0){
		perm[h][r+1][c]=1;
		first(deep+1,perm,sukima,hBar);
		perm[h][r+1][c]=0;
	}
	
	//高さ
	perm[h][r][c]=11;
	if(h==0){
		perm[h][r][c]=1;
		perm[h+1][r][c]=1;
	}
	first(deep+1,perm,sukima,hBar+1);
	
	if(h==0)perm[h+1][r][c]=0;
	perm[h][r][c]=0;
}

}
bool downCheck(S& s1,S& s2){
for(int r=0;r<3;r++){
	for(int c=0;c<3;c++){
		if(s1.blocks[0][r][c]!=s2.blocks[0][r][c])return false;
	}
}
return true;
}
S createNext(S& s1){
S s2;
memset(s2.blocks,0,sizeof(s2.blocks));
for(int r=0;r<3;r++){
	for(int c=0;c<3;c++){
		if(s1.blocks[1][r][c]==11)s2.blocks[0][r][c]=0;
		else s2.blocks[0][r][c]=1;
	}
}
return s2;
}
void ansCalc(char c){
if(c=='0')return;
if(ansMemo.empty()){
	ansMemo.insert(memo.begin(),memo.end());
}else{
	ansNext.clear();
	std::map<S,__int64>::iterator it1,it2;
	S s1,s2,s3,s4;
	for(it1=ansMemo.begin();it1!=ansMemo.end();it1++){
		s1=(*it1).first;
		s2=createNext(s1);
		for(it2=memo.lower_bound(s2);it2!=memo.end()&&downCheck((*it2).first,s2);it2++){
			s4=connectBlock(s1,(*it2).first);
			if(ansNext.find(s4)==ansNext.end()){
				ansNext[s4]=((*it1).second*(*it2).second)%base;
			}else{
				ansNext[s4]=(ansNext[s4]+((*it1).second*(*it2).second)%base)%base;;
			}
		}
	}
	ansMemo.clear();
	ansMemo.insert(ansNext.begin(),ansNext.end());
}
}
bool lastCheck(S& s1){
int t;
for(int r=0;r<3;r++){
	for(int c=0;c<3;c++){
		if(s1.blocks[0][r][c]==0)return false;
		if(s1.blocks[1][r][c]!=1)return false;
	}
}
return true;
}
int main(){
int perm[2][3][3];
memset(perm,0,sizeof(perm));
first(0,perm,0,0);
std::map<S,__int64>::iterator it1,it2,itD,itU;
S s1,s2,s3,s4;
//char text[30]="01111100010000101111";//10^6の逆転
char text[30]="0001011111";//10^3の逆転
int len=strlen(text);
std::cout<<len<<"\n";
int keta;
for(keta=1;keta<len;keta++){
	ansCalc(text[keta]);
	nextMemo.clear();
	for(it1=memo.begin();it1!=memo.end();it1++){
		s1=(*it1).first;
		s2=createNext(s1);
		for(it2=memo.lower_bound(s2);it2!=memo.end()&&downCheck((*it2).first,s2);it2++){
			s4=connectBlock(s1,(*it2).first);
			if(nextMemo.find(s4)==nextMemo.end()){
				nextMemo[s4]=((*it1).second*(*it2).second)%base;
			}else{
				nextMemo[s4]=(nextMemo[s4]+((*it1).second*(*it2).second)%base)%base;;
			}
		}
	}
	memo.clear();
	memo.insert(nextMemo.begin(),nextMemo.end());
	std::cout<<keta<<" "<<ansMemo.size()<<" "<<memo.size()<<"\n";
}

//蓋はいらなくなった

__int64 ans=0,t;
for(it1=ansMemo.begin();it1!=ansMemo.end();it1++){
	s1=(*it1).first;
	if(lastCheck(s1)==true)ans=(ans+(*it1).second)%base;
}
std::cout<<"\nans="<<ans;
}