読者です 読者をやめる 読者になる 読者になる

clock-up-blog

go-mi-tech

巨大数:アッカーマン関数のC++実装

C++ 巨大数 アッカーマン関数 数学

本記事の概要

巨大数と、実直な実装と、スタック使いつぶしへの対策。
C++というかほぼCライクなC++にて。

補足

少し誤解が生まれそうなので補足しておきますが、寿司の小林銅蟲と当ブログのkobakeは別人です。苗字被ってますが血縁関係もありません。交流はあります。

巨大数という謎の数学分野

数学マニア達が嗜む分野のひとつとして 巨大数 という謎の領域があります。
数式とか関数とかでひたすらデカイ数を表現するという、よく分からない分野です。

何を思ったのか漫画で解説を試みている例があるので紹介しておくと、裏サンデー寿司 虚空編 (作:小林銅蟲) という作品があります。分かりやすいかどうかは分かりません。個人的な感想としては1~2話はなんとなく理解できましたが3話はよく分かんないですね。

f:id:kobake:20131208221234j:plain

先日、小林銅蟲先生に肉を焼いてもらいつつ巨大数の解説をしていただく機会があり、なんとなく理解が進みました。少しだけ。

アッカーマン関数のシンプルな実装

巨大数のひとつ、アッカーマン関数 というモノについて考えてみます。寿司の2話に登場します。これは文字通り関数なのでプログラムで簡単に実装できます。Wikipediaの項 からの引用になりますが、アッカーマン関数の定義は以下のような感じです。

f:id:kobake:20131208214845p:plain

そのままプログラムの関数に落とし込めまそうですね。僕は C++ が好きなので C++ で実装してみました。

int Ack(int m, int n)
{
	if(m == 0)return n + 1;
	else if(n == 0)return Ack(m - 1, 1);
	else return Ack(m - 1, Ack(m, n - 1));
}

超簡単です。うーん、こんなモノが巨大数を生むのだろうか…???

実行してみる

作った関数を以下のように実行します。

int main()
{
	printf("Ack(2, 1) = %d\n", Ack(2, 1));
	return 0;
}

実行結果

Ack(2, 1) = 5

Ack(2, 1) の結果は 5 らしいです。正しいのかどうかよく分からないですね。
幸い、Wikipediaの項アッカーマン関数の値の一覧(一部)が載っています。それによると Ack(2, 1) は 5 で正しいらしい。つまり今回組んだプログラムは正常らしいです。やった!完成だ!

もっと大きな数を突っ込んでみる

Ack(4, 1) あたりを実行してみよう。これの正解は 65533 らしい。

int main()
{
	printf("Ack(4, 1) = %d\n", Ack(4, 1));
	return 0;
}

実行結果

f:id:kobake:20131208224838p:plain:w400

落ちた。悲しい。

原因

これが Visual Studio のすごく優秀なところなのですが、デバッグ実行すると、だいたい落ちた原因が分かります。
f:id:kobake:20131208225400p:plain:w700

Ack の呼び出し時に Stack overflow が発生しています。
要するに再帰が深すぎてスタック容量使いつぶした。

再帰の深さを記録する

int g_maxdepth = 0;
int Ack(int m, int n, int depth = 0)
{
	// 深さの記録・表示
	if(depth > g_maxdepth)printf("depth: %d\n", g_maxdepth = depth);

	if(m == 0)return n + 1;
	else if(n == 0)return Ack(m - 1, 1, depth + 1);
	else return Ack(m - 1, Ack(m, n - 1, depth + 1), depth + 1);
}

はい。

実行結果

f:id:kobake:20131208231004p:plain:w600

再帰の深さが4千強あたりで落ちました。
ちなみにこれはデバッグビルド版の限界であり、リリースビルドだと深さ3万2千強あたりまで耐える。それでも計算は終わらない。

再帰コールを取り除こう

f:id:kobake:20131208231409j:plain:w400

寿司2話でいうところの関数Bが、今回のAckに当たります。処理自体はいわば再帰処理なんだけど、ガーっと展開される引数部分だけに注目します。要するに数字だけに注目します。
寿司の B(3, 3) の例でいうと、

3, 3
2, 3, 2
2, 2, 3, 1
2, 2, 2, 3, 0
2, 2, 2, 2, 1
2, 2, 2, 1, 2, 0
2, 2, 2, 1, 1, 1
2, 2, 2, 1, 0, 1, 0
2, 2, 2, 1, 0, 0, 1
2, 2, 2, 1, 0, 2

っていう感じに数列の変化だけを考える。

そう考えるとアラ不思議、再帰処理使わなくても数列をいじり倒すだけでイケそう。
で、書いたコードが以下。

// 再帰呼び出し使わない版
int Ack(int m_input, int n_input)
{
	// バッファ確保
	std::vector<int> vtable(1000000); // 100万個 = 4MB
	int ntable = (int)vtable.size();
	int* table = &vtable[0];
	int* table_end = &table[ntable]; // バッファ終端

	// 初期値入れる
	table[0] = m_input;
	table[1] = n_input;

	// 終端位置
	int* end = &table[2];

	// ひたすら計算
	while(1){
		// 要素数が1以下になっていたら計算終了
		int size = end - table;
		if(size <= 1)break;

		// 一番後ろの2つに注目
		int m = *(end - 2);
		int n = *(end - 1);
		if(m == 0){
			// 要素が1個減る
			*(end - 2) = n + 1; // まとめる
			// 終端を変更
			end--;
		}
		else if(n == 0){
			// 要素数は変わらない
			*(end - 2) = m - 1;
			*(end - 1) = 1;
		}
		else{
			// 要素が1個増える
			*(end - 2) = m - 1;
			*(end - 1) = m;
			*(end - 0) = n - 1; // 増えた分
			// 終端を変更
			end++;
			if(end >= table_end){
				printf("ERROR: Table overflow\n");
				exit(1);
			}
		}
	}

	// 計算結果
	return table[0];
}

実行結果

Ack(4, 1) = 65533

できた!
Ack(4, 1) の計算完了まで数十秒かかります。

関数展開の様子を表示してみよう

実装方式を数列の操作にしたおかげで、計算の途中経過もけっこう簡単に表示できます。
以下のような感じ。

	…
	…
	// ひたすら計算
	while(1){
		// 要素数が1以下になっていたら計算終了
		int size = end - table;
		if(size <= 1)break;

		// NEW: 途中経過表示
		for(int i = 0; i < size - 2; i++){
			printf("Ack(%d, ", table[i]);
		}
		printf("Ack(%d, %d)", *(end - 2), *(end - 1));
		for(int i = 0; i < size - 2; i++){
			printf(")");
		}
		printf("\n");

		// 一番後ろの2つに注目
		int m = *(end - 2);
		int n = *(end - 1);
		…
		…
	}

	// NEW: 計算結果表示
	printf("%d\n", table[0]);

	// 計算結果
	return table[0];

実行結果

Ack(4, 1) はさすがに表示が多くなるので、Ack(2, 1) あたりを計算してみました。

Ack(2, 1)
Ack(1, Ack(2, 0))
Ack(1, Ack(1, 1))
Ack(1, Ack(0, Ack(1, 0)))
Ack(1, Ack(0, Ack(0, 1)))
Ack(1, Ack(0, 2))
Ack(1, 3)
Ack(0, Ack(1, 2))
Ack(0, Ack(0, Ack(1, 1)))
Ack(0, Ack(0, Ack(0, Ack(1, 0))))
Ack(0, Ack(0, Ack(0, Ack(0, 1))))
Ack(0, Ack(0, Ack(0, 2)))
Ack(0, Ack(0, 3))
Ack(0, 4)
5
--
Ack(2, 1) = 5

寿司っぽくなりました。
ちなみに Ack(4, 1) の途中経過を表示すると、約28億行の表示になります。死ぬ。

Gist

今回のソースコードを Gist に上げました。
Visual Studio または gccコンパイル・実行できます。

Gist: Ackermann Function

未来

Ack(4, 2) とか計算してみたいですが、これは上記プログラムでは無理です。
まだまだ他の工夫が必要になります。心が折れなければ今後試みようと思います。

他言語

Haskell とかでは末尾再帰の解決を処理系が自動的にやってくれたりするらしいですね。悔しい。

おしまい

f:id:kobake:20131209232029j:plain:w400

});