2016-06-13

ベルヌーイ数の説明は大変なので省きます。wikipediaとかに詳しく書いてあるのでそちらを参照してください。

とりあえず、計算する際には以下の漸化式が便利だそうです。
$B_0 = 1$
$B_n = -\frac{1}{n+1}\sum_{k=0}^{n-1} \left( \begin{array}{c} n+1 \\ k \end{array} \right) B_k$
ここで縦の括弧はコンビネーションです。
今回はこれをそのまま使おうと思います。

まずはテンプレートメタプログラミングではなく普通のプログラミングの場合

int combination(int n, int m)
{
	if (m == n || m == 0){
		return 1;
	}
	else{
		return n * combination(n - 1, m - 1) / m;
	}
}
fraction bernoulli(int n)
{
	if (n == 0){
		return fraction(1, 1);
	}else{
		fraction sum(0, 1);
		for (int k = 0; k < n; ++k){
			sum += bernoulli(k) * combination(n + 1, k);
		}
		return sum / -(n + 1);
	}
}

さて、ベルヌーイ数は有理数になるので出力は分数です。動的な変数の分数ライブラリは標準にはまだ存在しないので適当なライブラリをfractionという名前で定義させてください。

一つ目の関数がコンビネーションです。ちなみに計算には以下の公式を利用しています。
$k \left( \begin{array}{c} n \\ k \end{array} \right) = n \left( \begin{array}{c} n-1 \\ k-1 \end{array} \right)$
「なんでこんな公式をいちいち覚えているんだ？」と思われる人もいるかと思われるが以下のような方法で暗記している。
「あるｎ人のメンバーからｋ人の委員を選び、その中から１人委員長を選ぶ組み合わせ」
これだと公式の左辺になるが、「委員より委員長を先に一人選んでから残りk-1人の委員を選ぶ」という方法も考えられて、それが右辺になる。

さてこの式、20あたりまではまあ普通に計算できるが20を過ぎたあたりからだんだん重くなると思う。高速化するにはまず、ベルヌーイ数の性質である「１を除く奇数の場合は０になる」を使う方法がある。もう一つは、必ずｎが低い値を計算するのでテーブル化する方法がある。後者を使う場合は実装を大きく変える必要があるが

さて本題に移ってこれをテンプレートメタプログラミングで頑張ってみます。動的な分数ライブラリは存在しないのになぜか静的な分数ライブラリはC++11から導入されています。理由はよく分かりませんがなんにせよこれを使わない手はありません。

template< size_t n, size_t m >
struct combination{
	enum{ val = n * combination< n - 1, m - 1 >::val / m };
};

template< size_t n >
struct combination< n, 0 >{
	enum{ val = 1 };
};

template< size_t n >
struct combination< n, n >{
	enum{ val = 1 };
};

template< int n >
struct Bernoulli{
private:
	template< int nn >
	struct Sum{
		typedef
			std::ratio_add<
				std::ratio_multiply<
					std::ratio< combination< n + 1, nn >::val, 1 >,
					typename Bernoulli< nn >::val
				>,
				typename Sum< nn - 1 >::val
			> val;
	};
	template<>
	struct Sum< 0 >{
		typedef std::ratio< 1, 1 > val;
	};
public:
	typedef
		typename std::ratio_divide<
			typename Sum< n - 1 >::val,
			std::ratio< -(n + 1), 1 >
		> val;
};

template<>
struct Bernoulli< 0 >{
	typedef std::ratio< 1, 1 > val;
};

計算方法は以下の通り。結果はstd::ratioの型で出てきます。

typedef Bernoulli< 10 >::val num;
std::cout << num::num << " / " << num::den << std::endl; }

最初の3つはコンビネーション部分です。定番なテンプレートメタプログラミングなので特筆することはないかと思います。
さて面倒なのはベルヌーイ数本体です。何が面倒かというとforによるシグマ計算があります。テンプレートメタプログラミングでは変数が全て定数になるので変数を用意してそこに足していくという方法は使えませんので、ループを再起関数にして足し算していくのが定番ですね。privateに押し込めているsum構造体がそれに相当します。

さて、実際に計算したら意外とコンパイルが速かったです。なんでかな？ひょっとして静的だからテーブルがコンパイル中に出来てるのかな？まあ根拠のない想像ですが。

久々にテンプレートメタプログラミングでかなりてこずりました。いやどちらかというとratioでてこずったというか・・・。なれないライブラリは時間がかかります。でもまあ練習になったのでよしとします。

ちなみにこのプログラム、VS2012だとビルドが通りませんでした。理由はよく分かりませんがおそらくstd::ratio_divideにバグがあるのでは？だれか確認よろしくお願いします（ぇ

trinoの走り書き

久々にかなりピンチな状況になりました

ベルヌーイ数をテンプレートメタプログラミングで頑張ってみた。