フィボナッチ数（1）

フィボナッチ数を求めるのにO(1)とか言ってた人がいたが、実際のところはどうなのだろう。高校でも習う（今はどうか知らないが）公式を使えば結局べき乗になる。べき乗の計算はlog取って、掛け算して、expの計算をする。logってどうやって計算するのだろう。テーラー展開？テーラー展開だと、F_nを求めるのに、n桁のオーダーの掛け算をO(n)繰り返さなければならない。ややあやしいが、掛け算がO(nlogn)だとすると、logの計算をするのにO(n²logn)となる。単純に足し算していくと、n桁の足し算をn回だから、O(n)で済む。これを実証してみよう。

下のようにコードを書いてみた。大きな整数をintの配列にして、最後に基数変換するのがめんどうなので、一つのintで10⁹まで表すことにして、加算を繰り返した。これでFib₁₀₀〜Fib_1000000まで（抜粋して）10回求めて、それにかかった時間の最大と最小を除いた8つの平均を取った。

グラフを見ると、O(n²)であることがわかる。

#include 

#include 

#include 

#include 
using namespace std;
typedef unsigned int    uint;
const int   INT_LIMIT = (int)1e9;
void add(uint *a, const uint *b, int n);

void print_num(const uint *a, int n);
int main(int argc, char **argv) {

    int     n;

    try {

        if(argc != 2) {

            throw(1);

        }

        n = atoi(argv[1]);

        if(n < 0)

            throw(1);

    }

    catch(int e) {

        cerr << "usage : fib n." << endl;

        cerr << "    n >= 0." << endl;

        exit(1);

    }

    

    LARGE_INTEGER nFreq, nTime[2];

    QueryPerformanceFrequency(&nFreq);

    QueryPerformanceCounter(&nTime[0]);

    

    const double    phi = (1 + sqrt(5.)) / 2;

    int     num_items = (int)(n * (log(phi) / log( (double)INT_LIMIT)) + 1.5);

    

    uint    *a1 = (uint *)calloc(num_items, sizeof(uint));

    uint    *a2 = (uint *)calloc(num_items, sizeof(uint));

    if(a1 == NULL || a2 == NULL)

        exit(2);

    

    a1[0] = 0;

    a2[0] = 1;

    

    for(int k = 2; k <= n; k++) {

        uint    *cur, *prev;

        if(k & 1) {         // odd

            add(a2, a1, num_items);

        }

        else {

            add(a1, a2, num_items);

        }

    }

    

    print_num(n & 1 ? a2 : a1, num_items);

    

    free(a1);

    free(a2);

    

    QueryPerformanceCounter(&nTime[1]);

    int     t1 = (nTime[1].QuadPart - nTime[0].QuadPart) * 1e6 / nFreq.QuadPart;

    cerr << n << " " << t1 << " μs" << endl;

}
void add(uint *a, const uint *b, int n) {

    for(int i = 0; i < n; i++) {

        a[i] += b[i];

        if(a[i] >= INT_LIMIT) {

            a[i] -= INT_LIMIT;

            a[i+1]++;

        }

    }

}
void print_num(const uint *a, int n) {

    stringstream    ss;

    int     mode = 0;

    for(int i = n - 1; i >= 0; i--) {

        if(mode == 0) {

            if(a[i]) {

                ss << a[i];

                mode = 1;

            }

        }

        else {

            ss.width(9);

            ss.fill('0');

            ss << a[i];

        }

    }

    cout << ss.str() << endl;

}