C++でシュテルンの二項級数（fusc関数）のn番目の項を求める方法

この記事では、シュテルンの二項級数（Stern's Diatomic Series）のn番目の項を求める方法を、C++のコード例とともに解説します。

シュテルンの二項級数とは

シュテルンの二項級数は、次のような数列です。

0, 1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4, …

この数列は「fusc関数」としても知られており、正の整数を二進法で表したときに現れる連続する1の個数と深い関係があることで有名です。また、有理数を一意に列挙できることから、数学のさまざまな分野で応用されています。

漸化式による定義

この数列は、以下の漸化式で定義されます。

• nが偶数のとき：p(n) = p(n / 2)
• nが奇数のとき：p(n) = p((n − 1) / 2) + p((n + 1) / 2)
• 初期値：p(0) = 0、p(1) = 1

動的計画法によるアプローチ

単純な再帰呼び出しでこの漸化式を実装すると、同じ値を何度も計算することになり非効率です。そこで動的計画法（Dynamic Programming）を用いて計算量を削減します。

手順は以下の通りです。

1. p(0)とp(1)という基本ケース（ベースケース）を配列に保存する。
2. インデックスi = 2からnまで順にループし、偶奇に応じた漸化式でp(i)を計算して配列に格納する。
3. 最後にtable[n]を返す。

すでに計算済みの小さなインデックスの値を再利用するため、各項は一度だけ計算すればよくなります。

C++でのサンプルコード

#include<iostream>
using namespace std;

int findTerm(int n) {
    int table[n+1];
    table[0] = 0;   // 基本ケース: p(0) = 0
    table[1] = 1;   // 基本ケース: p(1) = 1

    for (int i = 2; i <= n; i++) {
        if (i % 2 == 0)
            table[i] = table[i / 2];
        else
            table[i] = table[(i - 1) / 2] + table[(i + 1) / 2];
    }
    return table[n];
}

int main() {
    cout << "3 番目の項: " << findTerm(3) << endl;
    cout << "15 番目の項: " << findTerm(15) << endl;
    cout << "20 番目の項: " << findTerm(20) << endl;
    return 0;
}

実行結果

3 番目の項: 2
15 番目の項: 4
20 番目の項: 3

計算量の評価

この実装では、各項を一度だけ計算するため、時間計算量はO(n)、結果を保存するための配列が必要となるため空間計算量もO(n)となります。素朴な再帰実装では指数関数的な計算量になりうる場面がありますが、動的計画法を使うことで線形時間で効率的に求められる点が大きなメリットです。

なお、さらに最適化したい場合は、nが偶数のときは右シフト演算（n >> 1）、奇数のときはビット操作を組み合わせることで、メモリを使用せずO(log n)程度で直接計算する方法もあります。