連分数

　連分数関連の処理をC++で書いて遊んでみました。

　連分数の基本については、Continued Fraction By Kardi Teknomo, PhD が分かりやすかったです。
　冪級数展開が発散するときでも連分数展開は収束する場合があり、かつ、どちらも収束するようなケースでは連分数展開の方が収束が早い場合が多いらしく、実用的なトピックです。黄金比やネイピア数を連分数展開すると綺麗な形になるというのもまた興味を唆られます。

分数を連分数展開する

p/qを連分数展開する場合を考えます。

p/q = [p/q] + (p%q)/q = [p/q] + 1/(q/(p%q))

となります。
ただし、[]はfloor関数です。

a₀ = [p/q], a₁ = [q/(p%q)], ....となるのでこれを再帰的に考えていくとユークリッドの互除法と同じパターンになることが分かります。

#include <iostream>
#include <vector>

using namespace std;

/*
 * p/q を連分数展開する
 */
void expand(int p, int q, vector<int> &terms) {
    if (q == 0)
        return;

    terms.push_back(p / q);
    expand(q, p % q, terms);

}

template <typename T>
void dump(const vector<T> &vs) {
    for (size_t i = 0; i < vs.size(); i++)
        cout << vs[i] << " ";
    cout << endl;
}

int main() {
    vector<int> terms;
    expand(109, 48, terms);
    dump(terms);  // 2 3 1 2 4 

    return 0;
}

連分数を分数に変換する

ボトムアップに計算する方法（右結合則を使って項数を減らしていくようなイメージ）もありますが、トップダウンに計算する方法の方が面白そうだったのでそちらをやってみます。

p₀ = a₀,  q₀ = 1,
p₁ = a₁ a₀ + 1,  q₁ = a_1, 
p_k = a_kp_k-1 + p_k-2(k > 1),
q_k = a_kq_k-1 + q_k-2(k > 1).

とすると、

[a₀; a₁, a₂, ..., a_k] = p_k / q_k

が成り立ちます（帰納法で証明できます）。

これを使えば連分数を分数に変換出来ます。
因みにa_kがすべて1の場合は、p_kとq_kが隣り合うフィボナッチ数になるため、連分数[1; 1,1,1,1,1,1,...]は黄金比になります。

#include <iostream>
#include <vector>
#include <iomanip>

using namespace std;

/*
 * 連分数を分数に変換する
 */
pair<int, int> compute(const vector<int> &terms) {
    int p[] = {terms[0], 1};
    int q[] = {1, 0};

    int len = terms.size();
    for (int i = 1; i < len; i++) {
        p[i & 1] = terms[i] * p[(i-1) & 1] + p[i & 1]; 
        q[i & 1] = terms[i] * q[(i-1) & 1] + q[i & 1];
    }
    --len;

    return make_pair(p[len & 1], q[len & 1]);
}


int main() {
    
    vector<int> terms(30, 1);
    pair<int, int> frac = compute(terms);
    cout << fixed << setprecision(15);
    cout << 1. * frac.first / frac.second << endl;   // 1.618033988750541

    frac = compute(vector<int> {2, 3, 1, 2, 3, 1});
    cout << fixed << setprecision(15);
    cout << 1. * frac.first / frac.second << endl;   // 2.270833333333333

    return 0;
}

実数を分数に変換する

実は、これをやりたく連分数について勉強したのです。

実数xを分数p/qで近似せよ。

という問題を考えます。

1. xを連分数展開する。
2. 連分数を分数に変換する。

の手順でp/qを導出することが出来ます。かなりテキトーなコードですが、お遊びレベルではそれっぽく動いています。

#include <cstdio>
#include <cmath>

using namespace std;

/*
 * 実数を分数に変換する
 */
void compute(double x, double eps = 1e-9) {

    double y = x;
    int a = (int)y;

    if (abs(a - x) < eps) {
        printf("%d / %d = %.15lf (diff = %.15lf)\n", a, 1, 1.*a, abs(a - x));
        return;
    }

    y = 1 / (y - a);
    int p[] = {a, 1};
    int q[] = {1, 0};

    int i = 0;
    while (i < 30) {
        a = (int)y;
        y = 1 / (y - a);
        ++i;
        p[i & 1] = a * p[(i-1) & 1] + p[i & 1]; 
        q[i & 1] = a * q[(i-1) & 1] + q[i & 1];

        if (abs(1. * p[i & 1] / q[i & 1] - x) < eps)
            break;
    }

    printf("%d / %d = %.15lf (diff = %.15lf)\n", 
           p[i&1], q[i&1], 1. * p[i&1] / q[i&1], abs(1. * p[i & 1] / q[i & 1] - x));
}

int main() {
    
    compute(sqrt(2));  // 47321 / 33461 = 1.414213562057320 (diff = 0.000000000315775)     
    compute(sqrt(3));  // 51409 / 29681 = 1.732050806913514 (diff = 0.000000000655363)  
    compute(M_PI);     // 103993 / 33102 = 3.141592653011902 (diff = 0.000000000577891) 
    compute(1.222);    // 611 / 500 = 1.222000000000000 (diff = 0.000000000000000)      

    return 0;
}

実は昔同じ問題を考察していたのを思い出しました。このときは二分探索/しゃくとり法で解こうとしていたみたいです。
今回連分数を勉強して、より実用的なツールを手に入れたなという感じです。