フィボナッチ数の最初の 100 個を表示する

任意精度の 10 進数加算でフィボナッチ数を求めるプログラム は求まる数の桁数が増えると遅くなるので、任意精度の 10 億進数加算に書き直し、「Five programming problems every Software Engineer should be able to solve in less than 1 hour - svpino.com」の問題 3 の答えにしておきます。

Write a function that computes the list of the first 100 Fibonacci numbers. By definition, the first two numbers in the Fibonacci sequence are 0 and 1, and each subsequent number is the sum of the previous two. As an example, here are the first 10 Fibonnaci numbers: 0, 1, 1, 2, 3, 5, 8, 13, 21, and 34.

ちなみに、 100 番目のフィボナッチ数は 64 ビット符号なし整数の最大値より大きいため、C++11 の規格だけで解くには 64 ビットより大きな整数同士の加算の実装が必須です。 問題を解くだけなら、 68 ビットあれば十分です。 ですが、 汎用性を求めて任意精度整数加算を使えば、 例えば先頭から 1 万個までのフィボナッチ数を求めることもできるようになります。 ライブラリ Boost.Multiprecision を使っても良いですし、加算だけなら自力で実装してもたいした手間ではありません。

ここでは、任意精度加算を組み込んである ruby 言語の CRuby 処理系用の下のコードを C++11 の main 関数に直訳できるように uintbig クラスを作ります。

f0, f1 = 0, 1
p f0, f1
(2 ... 100).each {
  f2 = f1 + f0
  f0, f1 = f1, f2
  p f1
}

ところで、基数を 10 億に選んだ理由ですが、 任意精度整数加算では、基数をバイナリに選んでも、デシマルに選んでも、速度差がさほどないため、10 進数文字列を作るのが楽なデシマルとしたわけです。 ここで、 基数にバイナリを選ぶと、 10 進数文字列化するために任意精度整数除算が必要になってしまいます。 そこで、 デシマルの中から桁上がり含めて符号なし 32 ビット整数の範囲内に入る 10 億を基数に選択しました。以前のやり方の 10 進数文字列同士を直接足して 10 進数文字列を求めたときは、文字列の末尾から先頭へ向かって計算をしていきました。今回はベクタ中の桁の配置を計算していく順番に合わせて 1 の位を先頭要素にしているため、ベクタの先頭から末尾へ向かって計算を進めます。一方、文字列化はベクタの末尾から先頭へ向かってゼロパディングの 9 桁の 10 進数文字列へと逆順に変換していきます。 他に、uintbig クラスには、std::swap を使えるようにするため、ムーブ・コンストラクタとムーブ代入演算子も定義しています。

#include <string>
#include <vector>
#include <iostream>

class uintbig {
    std::vector<uint32_t> frac;
public:
    enum {BASE = 1000000000, DIGITS = 9};

    uintbig () : frac (1, 0) {}
    ~uintbig () {}
    uintbig (uint32_t x0) : frac (1, x0) {}
    uintbig (uintbig&& x) : frac () { std::swap (frac, x.frac); }

    uintbig& operator= (uintbig&& x)
    {
        if (this != &x)
            std::swap (frac, x.frac);
        return *this;
    }

    uintbig operator+ (uintbig const& m1) const
    {
        if (frac.size () < m1.frac.size ())
            return m1 + *this;
        uintbig m;
        m.frac.resize (frac.size(), 0);
        uint32_t carry = 0;
        for (std::size_t i = 0; i < m1.frac.size (); ++i) {
            uint32_t const x = frac[i] + m1.frac[i] + carry;
            if (x >= BASE)
                carry = 1, m.frac[i] = x - BASE;
            else
                carry = 0, m.frac[i] = x;
        }
        for (std::size_t i = m1.frac.size (); i < frac.size (); ++i) {
            uint32_t const x = frac[i] + carry;
            if (x >= BASE)
                carry = 1, m.frac[i] = x - BASE;
            else
                carry = 0, m.frac[i] = x;
        }
        if (carry > 0)
            m.frac.push_back (carry);
        return m;
    }

    std::string to_string () const
    {
        std::size_t i = frac.size ();
        std::string s = std::to_string (frac[i - 1]);
        --i;
        for (; i > 0; --i) {
            std::string d = std::to_string (frac[i - 1]);
            std::string z (DIGITS - d.size (), '0');
            s += z + d;
        }
        return s;
    }
};

std::ostream&
operator << (std::ostream& o, uintbig const& m)
{
    std::string s = m.to_string ();
    o << s;
    return o;
}

int
main ()
{
    enum {N = 100};
    uintbig f0 (0);
    uintbig f1 (1);
    std::cout << f0 << std::endl;
    std::cout << f1 << std::endl;
    for (std::size_t i = 2; i < N; ++i) {
        uintbig f2 = f1 + f0;
        std::swap (f0, f1);
        std::swap (f1, f2);
        std::cout << f1 << std::endl;
    }
    return EXIT_SUCCESS;
}