これ、Hyukiさんをはじめ多くの方が疑問に思っていらっしゃるようなので、いまのうちに答えておきましょう。blogで書く以上、書く順番は順不同で構わないのですし。

アルゴリズム百選 - フィボナッチ数列にO()を学ぶ - www.textfile.org
フィボナッチ数列の一般項を求める式を使ったときってO(1)って言えるのだろうか?

まずは、論より証拠、というわけで実測値をご覧下さい。Cでフィボナッチ数をO(n)アルゴリズムと公式を使ってそれぞれ100万回計算した時にかかった時間をプロットしたものです。最適化をかけていないものと-O3で最適化をかけたものと双方を計測しています。fib(73)まで計算したのは、doubleで整数で保っておける限界がそこまでだったので。ソースは本entryの最後に。

graph

見ての通り、双方とも近似曲線も出しています。O(n)アルゴリズムの方は1次式で、公式利用の方は対数で。R2に注目してください。O(n)アルゴリズムの方は0.95を超えているのに対し、公式利用の方はかなり低い数字となっています。

どう見ても、O(log n)ではなくO(n)のようです。

しかし、これは確かに直感に反します。普通にベキ乗していたのでは、一つづつ掛け算してO(n)になるか、まとめて掛け算して(このアルゴリズムも面白いが今は省略)O(log n)になるかのどちらかのはずです。

一体libmのpow()は、どんな魔法をつかっているのでしょうか。

ここはやはり本連載のタネ本でもある、奥村先生の「C言語による最新アルゴリズム事典」に登場願いましょう。「最新」はさすがに今では問題ありますが、それ以外は今も使える名著です。ちなみに本連載は、初心者本と奥村本の橋渡しをするということが目的の一つ、Suffix Arrayなど奥村本にはやや新しすぎて登場しないアルゴリズムを紹介するのが今ひとつの目的でもあります。

P. 105

指数関数 exponential funciton

[中略]

収束を速くする工夫として、x= loge2 + t, |t| ≦ 1/2 loge2 を満たす t と整数 k を求め、etを上述のようにして計算し、ex = et2kとする。基数2の浮動小数点計算なら2k倍するのは指数部にkを加えるだけでよい。C言語にはそのためのライブラリ関数ldxp(k, k)がある。

見てのとおり、2kの計算はO(1)です。あとは etがO(1)であることを示せばよいわけです。奥村先生のソースに語ってもらいましょう。なお、奥村バージョンではlong doubleとなっているところをここではdoubleにしてあります。

P.106 より弾改変
#define LOG2 log(2)
#define N 22 /* 本文参照 */
double okumura_exp(double x){
  int i, k;
  double x2, w;
  k = (int)(x / LOG2 + (x >= 0 ? 0.5 : -0.5));
  x -= k * LOG2;
  x2 = x * x; w = x2 / N;
  for (i = N - 4; i >= 6; i -= 4) w = x2 / (w + i); /* ここに注目! */
  return ldexp((2 + w + x) / (2 + w - x), k);
}

このループの回数、見てのとおり常に一定です。上の22というのは、これだけ繰り返せば必ずdoubleの精度に収まるというマジックナンバーです。詳しくは奥村本をご覧ください。

指数関数がO(1)ということは、ベキ乗も自然とO(1)となります。もう少し厳密に言うと、精度に限りがあるのであれば、指数関数もO(1)でよい、ということです。

アルゴリズムって、ほんとすごいですねえ。

Dan the Inconstant

fibbench.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define sqrt5 sqrt(5)
#define phi_major ((1+sqrt5)/2)
#define phi_minor ((1-sqrt5)/2)

double time_hires(){
    struct timeval now;
    gettimeofday(&now, NULL);
    return now.tv_sec + now.tv_usec/1e6;
}

double fib_iter(int n){
    double f1=1, f2=1, t;
    if (n <= 2) return 1;
    while(n-- > 2){
        t = f1; f1 += f2; f2 =t;
    }
    return f1;
}

double fib_by_formula(int n){
    return (pow(phi_major, n) - pow(phi_minor, n))/sqrt5;
}

int main(int argc, char *argv[]){
    int i,j;
    double ff, fi, tf, ti;
    for (i = 1; i < 74; i++) {
         ti = time_hires();
         for (j = 0; j < 1e6; j++) fi = fib_iter(i);
         ti = time_hires() - ti;
         tf = time_hires();
         for (j = 0; j < 1e6; j++) ff = fib_by_formula(i);
         tf = time_hires() - tf;
         printf("fib(%d)\t =%.0f|%.0f\t%.0f\t%.0f\n", i, fi, ff, ti*1e6, tf*1e6);
    }
}