一開始也沒想到這個問題,但之後認真研究排列組合中的子集合時,卻發現有必要使用到 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 進行測試何者運算較快。結果如下所示:

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log(double) 0.05 0.51406 5.09812
log(long double) 0.0906 0.91436 9.05688
log(float) 0.0937 0.945 9.45124

由上表可知,實際上 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. 公式轉換法

對數中有所謂的換底公式,即 logab = logcb / logca;故,log2x 可進行對換為 log2x = logx / log2,或 log2x = log10x / log102,程式碼大致如下

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 即可,但此改善仍有限,這部份實測結果如下表所示

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log(double)*LOG_INVERSE 0.086 0.87876 8.65938
log10(double)*LOG10_INVERSE 0.0875 0.8825 8.71124
log(double)/log(2) 0.0859 0.88594 8.71938
log10(double)/log10(2) 0.0875 0.87968 8.72436

另使用換底公式時「有可能」會有另一個問題。正常而言,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 結果如下

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
shift 0.0813 0.9125 10.09594

測試結果 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;

實際結果如下所示

 

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log2_bitwise_magic_1_roll 0.0328 0.3228 3.23282
log2_bitwise_magic_1 0.0625 0.65 6.49

有沒有展開,效能差了快一倍!

 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;
}

 

結果如下所示

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log2_bitwise_magic_2_roll 0.0313 0.32094 3.23906
log2_bitwise_magic_2 0.0453 0.4578 4.61344

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];
}

 

測試結果

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log2_debruijn_int 0.025 0.26312 2.63374
log2_debruijn 0.0391 0.38844 3.92688

 

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;
}

 

 結果如下

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log2_table 0.0312 0.31406 3.30218
log2_table_macro 0.0329 0.32594 3.43656

 

 

吾人不建議用 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. 比較整理

  測試次數 / 每次測試迴圈次數
Method 10/10^6 50/10^7 50/10^8
log2_debruijn_int 0.025 0.26312 2.63374
log2_bitwise_magic_1_roll 0.0328 0.3228 3.23282
log2_bitwise_magic_2_roll 0.0313 0.32094 3.23906
log2_table 0.0312 0.31406 3.30218
log2_table_macro 0.0329 0.32594 3.43656
log2_debruijn 0.0391 0.38844 3.92688
log2_bitwise_magic_2 0.0453 0.4578 4.61344
log(double) 0.05 0.51406 5.09812
log2_bitwise_magic_1 0.0625 0.65 6.49
log(double)*LOG_INVERSE 0.086 0.87876 8.65938
log10(double)*LOG10_INVERSE 0.0875 0.8825 8.71124
log(double)/log(2) 0.0859 0.88594 8.71938
log10(double)/log10(2) 0.0875 0.87968 8.72436
log(long double) 0.0906 0.91436 9.05688
log(float) 0.0937 0.945 9.45124
shift 0.0813 0.9125 10.09594

10. 參考資料

bit twiddling hacks

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

MSDN

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()