Edison.X. Blog

一開始也沒想到這個問題，但之後認真研究排列組合中的子集合時，卻發現有必要使用到 log2 此函式。在 cmath (math.h) 裡面，提供 log 與 log10 函式，其中 log 即為 ln；但log2 函式卻視 compiler 決定是否支援，吾人手邊之 MSCV 不論哪個版本均無 log2 這個東西，故特撰此文。下述將先探討一些問題，再進行解法與說明。

本文所探討之 log2，以求得整數為主，求得小數之浮點數，日後再予以探討，文中提供約五種方法以說明。

0. 筆者環境

作業系統：Win XP SP3
CPU：AMD Athlon(tm) II X2 245 2.91G 雙核
RAM：2.75 GB 實體位置延伸
顯卡：內建 NVIDIA GeForce 7025 / NVIDIA nFource 630a
Compiler : VC 2008

1. 確認電腦系統對於浮點數效能

在 math.h 中之函式，大多都有三種引數資料型態：float、double、long double；早期的書在整數強制轉型為浮點數時，都以 float 為轉型對象，於此以 log 進行測試何者運算較快。結果如下所示：

由上表可知，實際上 double 是最快的，約可節省其它資料型態 80% 之時間。雖在 VC 中 sizeof(double) = sizeof(long double)，但實際上運算卻比 double 還慢！讓部份人訝異的是，float 運算竟是最慢的。列出此表並非絕對，只是提醒，確認電腦系統何種浮點數運算會較快！也由於吾人手邊實際為 double 較快，故以下若不得已必需進行轉型，則以 double 為轉型對象。

[Lemma]  另 CPU 型號 (種類，如 AMD /  Intel) 對於浮點數、整數之處理速度並不同，且 compiler 對於內建數學函式之處理速度亦不相同，吾人手邊之 CPU 為 AMD，所測得之結果不盡正確，唯供參考。同時 long double , double, float 之處理速度，以現今電腦而言，若 sizeof(long double) = 10，則 long double > double > float；若 sizeof(long double)=8，則 double > long double > float (感謝 Jacob Lee 補充指導)。

2. 公式轉換法

對數中有所謂的換底公式，即 log_ab = log_cb / log_ca；故，log2x 可進行對換為 log₂x = logx / log2，或 log₂x = log₁₀x / log₁₀2，程式碼大致如下

double log2_log(double x)
{
    return log(x)/log(2.0);
}

double log2_log10(double x)
{
    return log10(x)/log10(2.0);
}

那要用哪種比較好？事實上這二個在速度上並沒有太大差異，甚至有較不佳的地方 - 在於 log10(2.0)，log(2.0)，由於這二個數都常用，所以可以先定義這二個東西 

const double LOG_INVERSE = 1.0/log(2.0);
const double LOG10_INVERSE = 1.0/log10(2.0);

到時調用時直接進行 log(x) * LOG_INVERSE 或 log10(x) * LOG10_INVERSE 即可，但此改善仍有限，這部份實測結果如下表所示

另使用換底公式時「有可能」會有另一個問題。正常而言，log(8) / log(2) = 3，但部份 compiler 出現是 2.999999 (不是所有 compiler 都有這個問題，據悉 BCB 即有此問題)；出現這問題就算了，如果直接要轉成整數的話就變成2，造成 log(8) / log(2) =2 的錯誤答案。若只是要求整數的話可考慮這麼做

// method1, 取接近之整數
double x = log(8.0) / log(2.0);
double rst = x+0.5 ;

若直接用上面這段，則就取不到小數之資訊，只能取到整數之資訊，故這裡再提供一版較麻煩的版本，速度也較慢(有待優化)。

// method2, 判斷是否與整數相當接近,若是則直接取代
#define EPS 1E-9
double x = log(8.0) / log(2.0);
double xc = ceil(x), xf = floor(x);
double rst;

if((xc-x) <= EPS ) rst=xc;
else if((x-xf) <= EPS) rst=xf; 
else rst = x;

3. shift

事實上 log2 若以二進制方式為思考的話很容易想出怎麼做。假設 y = log2(x)，觀念大致如下

(1) 初值設 y=-1。 ----> y=-1;
(2) 若 x > 0 則繼續；否則退出。-----> while(x)
(3) x=x/2 , y=y+1 ----> x>>=1, ++y;
(4) 回到 step (2)

初版的 shift 程式碼

unsigned log2_shift(unsigned x)
{
    unsigned result; 
    for(result = -1; x > 0; x >>= 1, ++result) ; 
    return result; 
}

這版程式碼速度並不快，照著上面的流程並簡化應該長這樣 (感謝 Jacob Lee 指導提供)

int log2(int n)
{ 
  int i=0;
  while (n>>=1) ++i;
  return i;
}

這版程式碼比初版 shift 還來得快，且不用考慮整數長度之問題。實際結果，使用 shift 結果如下

測試結果 shift 比換底公式還慢。但別太吃驚，上述已說明， CPU 對於整數及浮點數處理上之表示不一，吾人手邊之 CPU 為 AMD，若為 Intel CPU，則 shift 遠比換底公式快！

4. bitwise magic_number- ver1

bitwise 進行 log2 整數求取有許多方式，這裡只取出其中幾種出來。現在說的方法為針對 32bits 整數，所附之 table 化為二進制後表示如下

b0=0000 0000 0000 0000 0000 0000 0000 0010 (s0=1)
b1=0000 0000 0000 0000 0000 0000 0000 1100 (s1=2)
b2=0000 0000 0000 0000 0000 0000 0000 1111 (s2=4)
b3=0000 0000 0000 0000 1111 1111 0000 0000 (s3=8)
b4=1111 1111 1111 1111 0000 0000 0000 0000 (s4=16)

接下來只要進行五次動作即可完成，程式碼如下所示

unsigned log2_bitwise_magic_1(unsigned x)
{
    const unsigned int b[] = {0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000};
    const unsigned int S[] = {1, 2, 4, 8, 16};
    int i;

    register unsigned int r = 0; // result of log2(v) will go here
    for (i = 4; i >= 0; --i) // unroll for speed...
    {
       if (x & b[i])
       {
          x >>= S[i];
          r |= S[i];
       } 
    }
    return r;
}

上述程式碼用到的是 for-loop，但實際上展開後卻可以跑得更快！展開後程式碼如下所示

unsigned log2_bitwise_magic_1_roll(unsigned x)
{
    register unsigned int r = 0; // result of log2(v) will go here
       if(x & 0xFFFF0000) x>>=16, r|=16;
       if(x & 0x0000FF00) x>>=8 , r|=8 ;
       if(x & 0x000000F0) x>>=4 , r|=4 ;
       if(x & 0x0000000C) x>>=2 , r|=2 ;
       if(x & 0x00000002) x>>=1 , r|=1 ;
    return r;
}

若是考慮為 64bits ，依此類推，在最前面加上

if(x & 0xFFFFFFFF00000000) x>>=32, r|=32;

實際結果如下所示

有沒有展開，效能差了快一倍！

 5. bitwise magic_number- ver2

這版的 magic number 展成2進位後長這樣

b0=1010 1010 1010 1010 1010 1010 1010 1010
b1=1100 1100 1100 1100 1100 1100 1100 1100 
b2=1111 0000 1111 0000 1111 0000 1111 0000
b3=1111 1111 0000 0000 1111 1111 0000 0000
b4=1111 1111 1111 1111 0000 0000 0000 0000

程式碼如下所示

unsigned log2_bitwise_magic_2(unsigned x)
{    
    unsigned int i;
    static const unsigned int b[] = {0xAAAAAAAA, 0xCCCCCCCC, 
       0xF0F0F0F0, 0xFF00FF00, 0xFFFF0000};
    register unsigned int r = (x & b[0]) != 0;
    for (i = 4; i > 0; i--) // unroll for speed...
    {
       r |= ((x & b[i]) != 0) << i;
    }
    return r;
}

同樣的，for loop 可以予以展開，變成下面這樣

unsigned log2_bitwise_magic_2_roll(unsigned x)
{    
    register unsigned int r = (x & b[0]) != 0;
    r |= ((x & 0xFFFF0000) != 0) << 4;
    r |= ((x & 0xFF00FF00) != 0) << 3;
    r |= ((x & 0xF0F0F0F0) != 0) << 2;
    r |= ((x & 0xCCCCCCCC) != 0) << 1;
    r |= ((x & 0xAAAAAAAA) != 0) << 0;
    return r;
}

結果如下所示

6. debruijn

接下來說的有點複雜，建議直接到這個網站去看完，再回來看「應」可知道為什麼這麼做，這裡不再贅述。程式碼如下

unsigned log2_debruijn(unsigned x)
{    
    static const unsigned MultiplyDeBruijnBitPosition[32] = 
    {
       0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
       8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
    };

    x |= x >> 1; // first round down to one less than a power of 2
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;

    return MultiplyDeBruijnBitPosition[(unsigned)(x * 0x07C4ACDDU) >> 27];
}

 若確定引數必為 2^n 時，可改成下面型式，速度會更快！

unsigned log2_debruijn_int(unsigned x)
{    
    static const unsigned MultiplyDeBruijnBitPosition2[32] = 
    {
       0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 
       31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
    };
    return MultiplyDeBruijnBitPosition2[(unsigned)(x * 0x077CB531U) >> 27];
}

測試結果

7. log2_table

這個方式是只建一部份 table 進行分次處理。

unsigned log2_lookup_bitwise(unsigned x)
{
    static const int log2_table[256] = {
       0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
       6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
       7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
       7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
       8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
       8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
       8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
       8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8
    };

    int l = -1;
    while (x >= 256) { l += 8; x >>= 8; }
    return l + log2_table[x];
}

有另類似寫法也有 macro 版本

unsigned log2_table_macro(unsigned int v) 
{
    static const char LogTable256[256] = 
    {
#define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
       -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
       LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
       LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
    };

    unsigned r;     // r will be lg(v)
    register unsigned int t, tt; // temporaries

    if (tt = v >> 16)
    {
       r = (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
    }
    else 
    {
       r = (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
    }
    return r;
}

 結果如下

吾人不建議用 macro 取代之，一方面不容易讓人看懂；一方面「目前實測」速度並無較快！

8. 其餘方式

網路有兩段 code 事實上速度會更快，但是內崁 asm ，吾人不知如何修改，附上以供參考！

/*
// ===================================================
// bsf
unsigned log2_BSF(unsigned n) { __asm bsf eax, n }

// ===================================================
// fsater bsf
unsigned __fastcall log2_FastBSF(unsigned n) { __asm bsf eax, ecx }
*/

9. 比較整理

10. 參考資料

bit twiddling hacks

http://www.flounder.com/log2.htm

MSDN