この話題、以下の答えとしても適度なのでそのまま。

アルゴリズム百選 - フィボナッチ数列にO()を学ぶ - www.textfile.org
「O()が小さいからといって速いとは限らない」が抜けている。

ベキ乗アルゴリズム再考

ベキ乗のやり方として、すぐに思いつくのは以下の方法です。

function power(b, n){
  var result = 1;
  while(n--) result *= b; // b を n 回掛け算
  return result;
}

これがO(n)であることは、直感的にわかります。

ところが、これをO(log n)でやる方法も比較的すぐに思いつきます。

例えばbを21乗したいとします。21=16+4+1なので、b21はb(16 + 4 + 1)とも書けます。b16は((b2)2)2、b4は(b2)2です。ということは、nを二のベキで表現しておき--すなわち二進数そのままで表記しておき、桁ごとにbを二乗して、その桁が1の時だけbをかけるという方法をとれば、この演算の手間はnではなくnの二進法による桁数、すなわちlog2nに比例するはずです。JavaScriptで書くとこうなります。

function power(b, n){
  var result = 1;
  while(n != 0){
    if (n & 1) result *= b; // 最下位ビットが1ならbを掛ける
    b *= b                  // bを自乗
    n >>= 1                 // nを右シフト 
  }
  return result;
}

O(1)だから速いとは限らない

chart0

それでも404 Blog Not Found:アルゴリズム百選 - ベキ乗はO(1)でOK?で紹介したように、libmのpow()関数は、nの大きさに関わらずほぼ一定時間で答えを出しているように見えます。しかし、問題はこの一定時間の中身です。それほど大きなnを扱わないのに、常に1時間かかるというのでは、むしろO(n)やO(log n)のアルゴリズムを使った方が速く終る可能性だってあるでしょう。

実際に私のMacBook Proで調べてみた結果が、右のグラフです。ここでは3のベキ乗を64bit整数で0から39まで、それぞれ100万回計算してベンチマークを取っています。青がlibmのpow()関数を使ったもの、緑がO(n)アルゴリズム、そして赤がO(log n)アルゴリズムです。

なんということでしょう。一定時間で結果を出すはずのlibmが一番遅いという結果が出ました。これにはさまざまな理由が考えられます。手製の関数では全て整数演算しているのに対し、pow()では整数を浮動小数点に直す手間もそれを整数に戻す手間もかかっています。n=0で例外的に速いのはpow()がそれを特別なケースとして扱っていることが暗示されています。


When in doubt, Benchmark.
[疑わしければ、ベンチマーク。]
-- Mastering Algorithms with Perl

Benchmarkの罠の傾向と対策

実は草稿段階でベンチマークを取った時、私は一つの罠に落ちていました。「最適化した場合、この結果が変わる場合もある」というのが論旨だったのですが、過剰最適化を見落としていたのです。

詳しくは

を見ていただくとして、最適化の結果、プログラマーが「こうしてくれ」と命じたことと実際のコードがやることが食い違うことがあるのです。この場合、「100万回計算しろというが、最後の一回しか計算結果を利用していないなら繰り返さなくてもOK」とコンパイラーが判断したのが落とし穴でした。

C/C++でプログラムする場合、これを見つけるのは結構困難なのですが、これに対して効果的なのが、「プロファイルを取る」というものです。上記の記事では、プロファイラーを使って実際に関数呼び出しを追うことで、私が落ちた穴を見事に見つけています。

When the Benchmark is in doubt, Profile.
[ベンチマークも疑わしければ、プロファイルせよ]
-- Damian Conway

Dan the Benchmarker

追記: 精度が充分だったというのがわかったのでs/powl/pow/g。他の環境におけるベンチマーク結果も考慮

追々記:過剰最適化を防ぐようコード変更。id:odz++。その結果を受けて改稿。

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

#ifndef REPEAT
#define REPEAT 1e6
#endif

typedef unsigned int U32;
typedef unsigned long long U64;

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

static U64 count;

U64 uipower_o_logn(U64 b, U64 n){
    U64 result = 1;
    for (; n != 0 ; b *= b, n >>= 1) if (n & 1) result *= b;
    count++;
    return result;
}

U64 uipower_o_n(U64 b, U64 n){
    U64 result = 1;
    for (; n != 0 ; n--) result *= b;
    count++;
    return result;
}

U64 uipower_libm(U64 b, U64 n){
    U64 result = pow(b, n);
    count++;
    return result;
}

int main(int argc, char *argv[]){
    U32 i, j;
    U64 rl, rn, rm;
    double tl, tn, tm;
    for (i = 0; i < 40; i++) {
        tl = time_hires(); 
        count = 0;
        for (j = 0; j < REPEAT; j++) rl = uipower_o_logn(3, i);
        tl = time_hires() - tl;
 
        tn = time_hires(); 
        count = 0;
        for (j = 0; j < REPEAT; j++) rn = uipower_o_n(3, i);
        tn = time_hires() - tn;

        tm = time_hires(); 
        count = 0;
        for (j = 0; j < REPEAT; j++) rm = uipower_libm(3, i);
        tm = time_hires() - tm;
#if DEBUG
        printf("uipower(3,%u)=%llu/%llu/%llu\t%.0f\t%.0f\t%0.f\n",
               i, rl, rn, rm, tl*1e6, tn*1e6, tm*1e6);
#else
        printf("uipower(3,%u)=%llu\t%u\t%.0f\t%.0f\t%0.f\n",
               i, rl-rn+rm, i, tl*1e6, tn*1e6, tm*1e6);
#endif
     }
}