2007年12月04日 23:00 [Edit]

アルゴリズム百選 - ベキ乗はO(1)でOK?

これ、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);
    }
}

この記事へのトラックバックURL

この記事へのトラックバック
あいかわらずの体調不良です。1つ前のエントリーで体調不良と数学について書いたら、hyuki先生から「お大事に」とのコメントをいただきました。これはもう「数学ガール/フェルマーの最終定理」買っちゃいなよ!ということなんだと思います。しかし、財布の中には千円しか入
[数学塾]フィボナッチ数列を求める公式とか計算量(オーダー)の話とか【医者を志す妻を応援する夫の日記】at 2009年02月20日 22:16
…結局気になってこればっかりかんがえちゃうから仕事できないじゃないか (笑) まぁA君 (S君?) の舌鋒鋭さに少し引き気味ではあるけども. というわけで,dankogai氏の次の記事のことです.また. アルゴリズム百選 - ベキ乗はO(1)でOK? アルゴリズム百選 - 用語の定義、ま
誤解の現況は「計算量」ということば?【okamoto7の日記】at 2007年12月06日 14:30
「実際に起こりうる状況ではご機嫌に速い」という意味で「O(1)」と書くことがある...
O(1)というのはご機嫌に速いということ?【Inquisitor】at 2007年12月05日 01:32
この記事へのコメント
> 近似計算の正しさ
doubleの精度が15桁くらいだったと思います。
φ^70 ≒ 4x10^14 ですから、n>70 ではもう正確ではないですね。というか
既に int32 には収まりませんね。n>100だと int64 でも収まりません
Posted by emanon at 2007年12月07日 16:21
近似計算の正しさはどう確認するのだろう。70<nだと、間違ってない?
Posted by anonymous at 2007年12月07日 02:01
手動トラックバック
http://www.algo.ics.tut.ac.jp/~yusuke_abe/html/2007-12-05.html#2007-12-05-2
Posted by schrodin at 2007年12月05日 20:48
極限(ε-δ論法)から勉強し直した方がいいんじゃないのかと訝しがってしまう程の苦しい実験ですね
Posted by hoge at 2007年12月05日 20:28
dankogai の O(1) は俺流の O
Posted by O Perl at 2007年12月05日 09:22
手動トラックバック。
http://d.hatena.ne.jp/textfile/20071204/algo
Posted by hyuki at 2007年12月05日 09:03
いや弾さん、精度を区切っちゃったらO記法の議論は成り立たないっすよ。極端な話、テーブル索引にすればどんな関数だって定数時間になります。
確かに、現場では厳密なオーダーの算出よりも限られた精度内で性能が出る方法を考える必要に迫られることの方が多いですが、これから学ぶ人に混乱を与えないためにも、「doubleの範囲では定数時間で計算するコンパクトなアルゴリズムがある」とかそういう言い方にしておいてもらえると助かります。
Posted by shiro at 2007年12月05日 07:15
奥村本の引用部分、
x = log_e 2 + t の右辺第一項は係数 k が必要では?
Posted by anonymous at 2007年12月05日 00:54