2016-12-05 84 views
2

對於我的項目,我需要做一個非常快速的計算pow(x,y)。希望它有點侷限於一個精確的領域,但如果速度不夠快,它也需要足夠高的內存效率。以優化的方式計算0(x,y)對於0 <= x <= 1且1 <y <2的C++

就像我說的那樣,它在0到1之間的x的短範圍內,y在1到2之間。因此,在整個範圍內必須足夠精確,以便遞歸地調用時會有輕微的縮小而不是拖延對數)

如果你們遇到了這樣的事情或有建議...

+1

@Mehrdad編輯 –

+1

'exp(y * log(x))'太慢了怎麼辦? - 你認爲什麼精度足夠了? (y-1)*(x-1)+(y-1)*((x-1) y-2)/ 2 *(x-1)^ 2 + O(x-1)^ 3'作爲近似值 – LutzL

+0

您是否考慮過嘗試查找表格? –

回答

3

通過

exp(y*log(x)) 

更換pow(x,y)這也是有時C庫實現,但給出稍微扭曲的結果對於許多整數輸入。因此,pow(x,y)實現通常有一個開銷,以捕獲指數爲0和1的平凡冪,整數冪,並進行一些其他轉換以獲得最精確的結果。削減這個開銷可能已經足夠加速。

+0

您認爲使用查找表優化日誌可以帶來更多的好處,因爲x作爲定義的範圍?順便說一句,這是我想要的那種答案 –

+0

這完全取決於處理器,必須進行測試。 'exp(y * log(x))'需要在FPU堆棧中加載兩個參數,執行兩個FPU內部操作並存儲結果。查找表需要FPU,CPU和處理器緩存的協調,這似乎更復雜。但是這是標準的PC架構。在微處理器上,其他限制可能導致其他偏好。 – LutzL

+0

根據我的經驗,這裏提出的替換通常會使用現代數學庫將指數的性能提高兩倍(可能不如舊版庫以x87 FPU爲準)。 – njuffa

1

超過兩個x & p最優化方法,例如有界區間是運行的exp2()log2()兩種優化近似值。由於來自exp2()極小極小多項式的誤差將呈指數形式複合,但對於0 < x ≤ 1,輸入始終爲p*log2(x) ≤ 0,因此該誤差將表現良好。注:在x = 0, log2(0) = -∞。以下是powf()例程的示例代碼以及和log2()exp2的多次最小最大近似值的錯誤圖,它們有望滿足您的絕對容錯性。

image

float aPowf(float x, float p){ 
    union { float f; uint32_t u; } L2x, e2e; 
    float pL2, pL2r, pL2i, E2; 

    if(x == 0) return(0.0f); 

    /* Calculate log2(x) Approximation */ 
    L2x.f = x; 
    pL2 = (uint8_t)(L2x.u >> 23) - 127;     // log2(m*2^e) = log2(m) + e 
    L2x.u = (L2x.u & ~(0xFF << 23)) | (0x7F << 23); 

    // Approximate log2(x) over 1 <= x < 2, use fma() fused multiply accumulate function for efficient evaluation 
    // pL2 += -0xf.4fb9dp-4 + L2x.f; // 4.303568e-2 
    // pL2 += -0x1.acc4cap0 + L2x.f * (0x2.065084p0 + L2x.f * (-0x5.847fe8p-4)); // 4.9397e-3 
    // pL2 += -0x2.2753ccp0 + L2x.f * (0x3.0c426p0 + L2x.f * (-0x1.0d47dap0 + L2x.f * 0x2.88306cp-4)); // 6.3717e-4 
    pL2 += -0x2.834a9p0 + L2x.f * (0x4.11f1d8p0 + L2x.f * (-0x2.1ee4fcp0 + L2x.f * (0xa.52841p-4 + L2x.f * (-0x1.4e4cf6p-4)))); // 8.761e-5 
    // pL2 += -0x2.cce408p0 + L2x.f * (0x5.177808p0 + L2x.f * (-0x3.8cfd5cp0 + L2x.f * (0x1.a19084p0 + L2x.f * (-0x6.aa30dp-4 + L2x.f * 0xb.7cb75p-8)))); // 1.2542058e-5 
    // pL2 += -0x3.0a3514p0 + L2x.f * (0x6.1cbb88p0 + L2x.f * (-0x5.5737a8p0 + L2x.f * (0x3.490a04p0 + L2x.f * (-0x1.442ae8p0 + L2x.f * (0x4.66497p-4 + L2x.f * (-0x6.925fe8p-8)))))); // 1.8533e-6 
    // pL2 += -0x3.3eb71cp0 + L2x.f * (0x7.2194ep0 + L2x.f * (-0x7.7cf968p0 + L2x.f * (0x5.c642f8p0 + L2x.f * (-0x2.faeb44p0 + L2x.f * (0xf.9e012p-4 + L2x.f * (-0x2.ef86f8p-4 + L2x.f * 0x3.dc524p-8)))))); // 2.831e-7 
    // pL2 += -0x3.6c382p0 + L2x.f * (0x8.23b47p0 + L2x.f * (-0x9.f803dp0 + L2x.f * (0x9.3b4f3p0 + L2x.f * (-0x5.f739ep0 + L2x.f * (0x2.9cb704p0 + L2x.f * (-0xb.d395dp-4 + L2x.f * (0x1.f3e2p-4 + L2x.f * (-0x2.49964p-8)))))))); // 4.7028674e-8 

    pL2 *= p; 
// if(pL2 <= -128)  return(0.0f); 

    /* Calculate exp2(p*log2(x)) */ 
    pL2i = floorf(pL2); 
    pL2r = pL2 - pL2i; 
    e2e.u = ((int)pL2i + 127) << 23; 

    // Approximate exp2(x) over 0 <= x < 1, use fma() fused multiply accumulate function for efficient evaluation. 
    // E2 = 0xf.4fb9dp-4 + pL2r; // 4.303568e-2 
    // E2 = 0x1.00a246p0 + pL2r * (0xa.6aafdp-4 + pL2r * 0x5.81078p-4); // 2.4761e-3 
    // E2 = 0xf.ff8fcp-4 + pL2r * (0xb.24b0ap-4 + pL2r * (0x3.96e39cp-4 + pL2r * 0x1.446bc8p-4)); // 1.0705e-4 
    E2 = 0x1.00003ep0 + pL2r * (0xb.1663cp-4 + pL2r * (0x3.ddbffp-4 + pL2r * (0xd.3b9afp-8 + pL2r * 0x3.81ade8p-8))); // 3.706393e-6 
    // E2 = 0xf.ffffep-4 + pL2r * (0xb.1729bp-4 + pL2r * (0x3.d79b5cp-4 + pL2r * (0xe.4d721p-8 + pL2r * (0x2.49e21p-8 + pL2r * 0x7.c5b598p-12)))); // 1.192e-7 
    // E2 = 0x1.p0 + pL2r * (0xb.17215p-4 + pL2r * (0x3.d7fb5p-4 + pL2r * (0xe.34192p-8 + pL2r * (0x2.7a7828p-8 + pL2r * (0x5.15bd08p-12 + pL2r * 0xe.48db2p-16))))); // 2.9105833e-9 
    // E2 = 0x1.p0 + pL2r * (0xb.17218p-4 + pL2r * (0x3.d7f7acp-4 + pL2r * (0xe.35916p-8 + pL2r * (0x2.761acp-8 + pL2r * (0x5.7e9f9p-12 + pL2r * (0x9.70c6ap-16 + pL2r * 0x1.666008p-16)))))); // 8.10693e-11 
    // E2 = 0x1.p0 + pL2r * (0xb.17218p-4 + pL2r * (0x3.d7f7b8p-4 + pL2r * (0xe.35874p-8 + pL2r * (0x2.764dccp-8 + pL2r * (0x5.76b95p-12 + pL2r * (0xa.15ca6p-16 + pL2r * (0xf.94e0dp-20 + pL2r * 0x1.cc690cp-20))))))); // 3.9714478e-11 

    return(E2 * e2e.f); 
} 

一旦選擇合適的最小最大近似確保實現與稠合的霍納多項式求乘法累加運算FMA()[其爲單週期指令。

爲了提高對數函數的精度,可以通過引入用於範圍縮減類型的LUT來進一步提高精度。

相關問題