CPU /編程語言使用哪種指數運算算法？

Question

我一直在學習更快的指數運算算法（K-ary，滑動門等），並想知道哪些用於CPU /編程語言？ （我在此事件是否發生在CPU或編譯通過模糊）

而只是踢，這是最快的？

關於寬度的編輯：這是有意廣泛的，因爲我知道有一堆不同的技術來做到這一點。檢查的答案有我正在尋找的。

Answer 1

我假設你的興趣是在執行的，可以在標準的數學庫中找到高級語言，尤其是C/C++的冪函數。這些包括功能exp(),exp2(),exp10()和pow()以及單精度對應件expf(),exp2f(),exp10f()和powf()。

你提到乘方方法（例如k元，滑動窗口）中的加密算法，例如RSA，它是基於冪通常採用。它們通常不用於通過math.h或cmath提供的求冪函數。實施細節爲標準數學函數等exp()不同，但一種常見方案遵循三個步驟：

1. 還原函數參數的一個初級近似於一個合適的基函數的 間隔
2. 近似初級近似間隔
3. 映射回用於主間隔的函數

的整個範圍內的結果的輔助步驟常常是特別CAS的處理ES。這些可以涉及特殊的數學情況，如log(0.0)，或特殊的浮點操作數，如NaN（非數字）。

爲expf(float)的C99碼下面示出了在示例性方式是什麼的那些步驟看起來爲一個具體的例子。的參數a是第一分割使得exp(a) = E ^{- [R} * 2 ^我，其中i是整數並且r是在[日誌（SQRT（0.5），日誌（SQRT（2.0）]，初級近似時間間隔。在第二步驟中，我們現在近似Ë^{- [R}用多項式。這樣的近似值可根據各種設計標準，如最小化絕對或相對誤差進行設計。多項式可以通過各種方式，包括霍納方案和埃斯特林的方案進行評估。

下面的代碼使用了一個非常常見的方法，它使用了一個極小極大值近似法，它可以使整個近似區間內的最大誤差最小化。用於計算這種近似值的rd算法是Remez算法。通過霍納的計劃進行評估;通過使用fmaf()可以提高此評估的數值準確性。

該標準數學函數實現了所謂的融合乘加或FMA。此計算在添加期間使用完整產品a*b計算a*b+c，並在最後應用單個舍入。在大多數現代硬件上，例如GPU，IBM Power CPU，最近的x86處理器（例如，Haswell），最近的ARM處理器（作爲可選的擴展），它直接映射到硬件指令。在缺乏這種指令的平臺上，fmaf()將映射到相當慢的仿真代碼，在這種情況下，如果我們對性能感興趣，我們不希望使用它。

最後的計算是乘以2 ⁱ，對此，C和C++提供函數ldexp()。在「工業強度」庫代碼中，人們通常在這裏使用特定於計算機的習慣用法，該工具利用IEEE-754二進制算術的優勢，用於float。最後，代碼清理溢出和下溢的情況。

x86處理器內部的x87 FPU有一條指令F2XM1，它計算[012,1]上的2 ^x -1。這可用於計算exp()和exp2()的第二步。在第三步中有一條指令FSCALE，用於乘以ⁱ。實現F2XM1本身的常見方式是使用有理或多項式近似的微碼。請注意，x87 FPU現在主要用於傳統支持。在現代x86平臺上，庫通常使用基於SSE和算法的純軟件實現，這些實現與下面所示的類似。有些將小表與多項式近似結合在一起。

pow(x,y)從概念上可以實現爲exp(y*log(x))，但是從精度顯著損失這個遭受當x是在C/C++標準規定的近團結y大的幅度，以及衆多的特殊情況下，不正確的處理。解決準確度問題的一種方法是以某種擴展精度的形式計算log(x)和產品y*log(x))。細節將填充整個冗長的單獨答案，而我沒有代碼方便地演示它。在各種C/C++數學庫中，pow(double,int)和powf(float, int)由單獨的代碼路徑計算，該代碼路徑應用「平方和乘法」方法以及整數指數的二進制表示的按位掃描。

#include <math.h> /* import fmaf(), ldexpf() */ 

/* Like rintf(), but -0.0f -> +0.0f, and |a| must be < 2**22 */ 
float quick_and_dirty_rintf (float a) 
{ 
    float cvt_magic = 0x1.800000p+23f; 
    return (a + cvt_magic) - cvt_magic; 
} 

/* Approximate exp(a) on the interval [log(sqrt(0.5)), log(sqrt(2.0))]. 
    Maximum ulp error = 0.70951 
*/ 
float expf_poly (float a) 
{ 
    float r; 

    r =    0x1.6ab95ep-10f; 
    r = fmaf (r, a, 0x1.126d28p-07f); 
    r = fmaf (r, a, 0x1.5558f8p-05f); 
    r = fmaf (r, a, 0x1.55543ap-03f); 
    r = fmaf (r, a, 0x1.fffffap-02f); 
    r = fmaf (r, a, 0x1.000000p+00f); 
    r = fmaf (r, a, 0x1.000000p+00f); 
    return r; 
} 

/* Approximate exp2() on interval [-0.5,+0.5] 
    Maximum ulp error = 0.79961 
*/ 
float exp2f_poly (float a) 
{ 
    float r; 

    r =    0x1.4308f2p-13f; 
    r = fmaf (r, a, 0x1.5f0722p-10f); 
    r = fmaf (r, a, 0x1.3b2bd4p-07f); 
    r = fmaf (r, a, 0x1.c6af80p-05f); 
    r = fmaf (r, a, 0x1.ebfbdep-03f); 
    r = fmaf (r, a, 0x1.62e430p-01f); 
    r = fmaf (r, a, 0x1.000000p+00f); 
    return r; 
} 

/* Approximate exp10(a) on [log(sqrt(0.5))/log(10), log(sqrt(2.0))/log(10)]. 
    Maximum ulp error = 0.77879 
*/ 
float exp10f_poly (float a) 
{ 
    float r; 

    r =    0x1.a65694p-3f; 
    r = fmaf (r, a, 0x1.158a00p-1f); 
    r = fmaf (r, a, 0x1.2bda9ap+0f); 
    r = fmaf (r, a, 0x1.046f72p+1f); 
    r = fmaf (r, a, 0x1.535248p+1f); 
    r = fmaf (r, a, 0x1.26bb1cp+1f); 
    r = fmaf (r, a, 0x1.000000p+0f); 
    return r; 
} 

/* Compute exponential base e. Maximum ulp error = 0.88790 */ 
float my_expf (float a) 
{ 
    float t, r; 
    int i; 

    t = a * 0x1.715476p+0f;   // 1/log(2) 
    t = quick_and_dirty_rintf (t); 
    i = (int)t; 
    r = fmaf (t, -0x1.62e430p-1f, a); // log_2_hi 
    r = fmaf (t, 0x1.05c610p-29f, r); // log_2_lo 
    t = expf_poly (r); 
    r = ldexpf (t, i); 
    if (a < -105.0f) r = 0.0f; 
    if (a > 105.0f) r = 1.0f/0.0f; // +INF 
    return r; 
} 

/* Compute exponential base 2. Maximum ulp error = 0.87922 */ 
float my_exp2f (float a) 
{ 
    float t, r; 
    int i; 

    t = quick_and_dirty_rintf (a); 
    i = (int)t; 
    r = a - t; 
    t = exp2f_poly (r); 
    r = ldexpf (t, i); 
    if (a < -152.0f) r = 0.0f; 
    if (a > 152.0f) r = 1.0f/0.0f; // +INF 
    return r; 
} 

/* Compute exponential base 10. Maximum ulp error = 0.95588 */ 
float my_exp10f (float a) 
{ 
    float r, t; 
    int i; 

    t = a * 0x1.a934f0p+1f; 
    t = quick_and_dirty_rintf (t); 
    i = (int)t; 
    r = fmaf (t, -0x1.344136p-2f, a); // log10(2)_hi 
    r = fmaf (t, 0x1.ec10c0p-27f, r); // log10(2)_lo 
    t = exp10f_poly (r); 
    r = ldexpf (t, i); 
    if (a < -46.0f) r = 0.0f; 
    if (a > 46.0f) r = 1.0f/0.0f; // +INF 
    return r; 
}