数学の命題示しました

主に組合せ論について，読んだ本で出てきたことや，考えたことを書きます．

平面分割の列挙アルゴリズム

プログラミング

10月末に綾@ (@ayaHSYKZ)さんが，以下のような問題をツイートしました．

消えたときのために文章でも書いておくと，

以下の二条件を満たす 5\times 5行列 A=(a_{i,j})_{0\leq i,j\leq 4}の個数を求めよという問題です．

  • Aの各要素は 0か 1か 2である．
  • 全ての整数 0\leq i,j\leq 4に対して a_{i+1,j}\leq a_{i,j}かつ a_{i,j+1}\leq a_{i,j}が成立．つまり Aは各行各列について単調非増加である．


この問題は実は，「平面分割 (plane partition)」という名前でより一般に研究されています．
Plane partition - Wikipedia



今回は，一般に行列サイズが a\times bで，各要素が 0から n-1までの整数だとしたとき，上記の条件を満たす行列を列挙するアルゴリズムを設計し，c++で実装しました．

詳細は書くのがめんどくさいので省きますが，コードを載せておきます．
ちなみに，このアルゴリズムを使って求めた結果，綾@さんの問題の答えは「 19404通り」でした．

#include <iostream>
#include <vector>
#include <utility>

#define MIN(x,y) ((x<y)? x: y)

using namespace std;

void vector_vanish(vector<int> &vec);
void vector_inspect(vector<int> vec);
int vector_increment(vector<int> sup, vector<int> &data);
int level_increment(pair< vector<int>, vector<int> > sup,
					pair< vector<int>, vector<int> > &data);
int PP_increment(vector< pair< vector<int>, vector<int> > > &pp);
void PP_inspect(vector< pair< vector<int>, vector<int> > > pp);

int main(){
	const int r=5;
	const int c=5;
	//r<=cを仮定する．
	const int n=2;

	vector< pair< vector<int>,vector<int> > > pp;
	for(int i=r-1; i>=0; i--){
		vector<int> hor(c-i,0);
		vector<int> var(r-1-i,0);
		pair< vector<int>,vector<int> > adj(hor, var);
		pp.push_back(adj);
	}
	//平面分割を入れる箱を準備
	
	vector<int> hor(c+1,n);
	vector<int> var(r,n);
	pair< vector<int>, vector<int> > adj(hor, var);
	pp.push_back(adj);
	//平面分割の高さの上限を設定
	
	int count=0;
	do{
		PP_inspect(pp);
		count++;
		cout << endl;
	}while(PP_increment(pp));

	cout << "count = " << count << endl;

	double prod = 1;
	for(int i=0; i<=r-1; i++){
		for(int j=0; j<=c-1; j++){
			for(int k=0; k<=n-1; k++){
				prod *= ((double)(i+j+k+2))/((double)(i+j+k+1));
			}
		}
	}
	cout << "prod = " << prod << endl;

	return 0;
}

void PP_inspect(vector< pair< vector<int>, vector<int> > > pp){
	int r = pp.size()-1;
	int c = pp[r-1].first.size();

	vector< vector<int> > kekka;
	for(int i=0; i<=r-1; i++){
		vector<int> a(c,0);
		kekka.push_back(a);
	}

	for(int i=0; i<=r-1; i++){
		for(int j=0; j<=c-1; j++){
			if(i<=j){
				kekka[i][j] = (pp[r-1-i].first)[j-i];
			}else{
				kekka[i][j] = (pp[r-1-j].second)[i-j-1];
			}
		}
	}


	for(int i=0; i<=r-1; i++){
		for(int j=0; j<=c-1; j++){
			cout << kekka[i][j] << " ";
		}
		cout << endl;
	}

	//平面分割になってるかチェックする
	for(int i=0; i<=r-2; i++){
		for(int j=0; j<=c-2; j++){
			if(kekka[i][j] < kekka[i+1][j] ||
					kekka[i][j] < kekka[i][j+1]){
				cout << "Error: 平面分割になってません." << i << j << endl;
				exit(1);
			}
		}
	}
	return ;
}

int PP_increment(vector< pair< vector<int>, vector<int> > > &pp){
	int r = pp.size()-1;
	for(int i=0; i<=r-1; i++){
		if(level_increment(pp[i+1],pp[i])){
			//レベルiがインクリメントできたとき
			//レベルiより前を全部消す
			for(int j=0; j<=i-1; j++){
				vector_vanish(pp[j].first);
				vector_vanish(pp[j].second);
			}
			//mainに戻る
			return 1;
		}
	}
	//これ以上インクリメントできないとき
	return 0;
}


int level_increment(pair< vector<int>, vector<int> > sup,
					pair< vector<int>, vector<int> > &data){
	//verがインクリメントできるならverをインクリメント
	//horをインクリメントできるならhorをインクリメントしvarを消す
	//インクリメントできたなら1を返す
	
	//まず，関数vector_incrementに突っ込むためのベクトルを作る
	//verについて
	int ver_size = data.second.size();
	vector<int> vec_ver_sup(ver_size);
	for(int i=0; i<=ver_size-1; i++){
		vec_ver_sup[i] = MIN((sup.second)[i+1], (data.first)[0]);
	}
	//horについて
	int hor_size = data.first.size();
	vector<int> vec_hor_sup(hor_size);
	for(int i=0; i<=hor_size-1; i++){
		vec_hor_sup[i] = MIN((sup.first)[i+1], (sup.second)[0]);
	}
	
	//どちらもインクリメントできないなら0を返す
	if(vec_hor_sup == data.first && vec_ver_sup == data.second){
		return 0;
	}

	//varをインクリメント
	if(vector_increment(vec_ver_sup, data.second)){
		return 1;
	}
	//varがインクリメントできなかったとき
	//horをインクリメント
	vector_vanish(data.second);
	if(vector_increment(vec_hor_sup, data.first)){
		return 1;
	}

	//ここまで達するのはおかしい
	cout << "Error" << endl;
	exit(1);
}


int vector_increment(vector<int> sup, vector<int> &data){
	if(sup == data){
		return 0;
	}
	int r = sup.size();

	for(int j=r-1; 1<=j; j--){
		if(data[j-1] > data[j]){
			if(sup[j] > data[j]){
				data[j]++;
				for(int i=r-1; j<i; i--){
					data[i] = 0;
				}
				return 1;
			}
		}
	}
	if(sup[0] > data[0]){
		data[0]++;
		for(int i=r-1; 0<i; i--){
			data[i] = 0;
		}
		return 1;
	}else{
		cout << "Error" << endl;
		exit(1);
	}
}


void vector_vanish(vector<int> &vec){
	int size = vec.size();
	for(int i=0; i<=size-1; i++){
		vec[i] = 0;
	}
	return;
}

void vector_inspect(vector<int> vec){
	int size = vec.size();
	for(int i=0; i<=size-1; i++){
		cout << vec[i]  << " ";
	}
	cout << endl;
	return;
}