CPU/编程语言使用哪些求幂算法？ [关闭]答案

【问题标题】：Which exponentiation algorithms do CPU/programming languages use? [closed]CPU/编程语言使用哪些求幂算法？ [关闭]
【发布时间】：2015-06-05 12:41:45
【问题描述】：

我一直在学习更快的求幂算法（k-ary、滑动门等），并且想知道 CPU/编程语言中使用了哪些算法？（我不清楚这是发生在 CPU 中还是通过编译器）

只是为了踢球，哪个最快？

关于宽泛性的编辑：它是故意宽泛的，因为我知道有很多不同的技术可以做到这一点。检查的答案有我想要的。

【问题讨论】：

当您询问不同的 CPU 或编程语言时，您无法合理地询问“哪种算法”。也许不同的系统使用不同的算法。你必须更具体。
@GregHewgill 抱歉，语法错误。我只是在找一个例子，比如“x86 使用滑动门”之类的。
在 C/C++ 数学库中，pow(double,double) 和 pow(double,int) 通常是不同的代码路径，而 exp(double) 是另一个单独的函数。您对其中哪些感兴趣？ pow(double,int) 通常采用基于指数位扫描的平方乘法变体。 x86 处理器中的旧 x87 FPU 有一个 F2XM1 指令，可用于执行 exp() 和 pow()。如今，基于 x86 的系统通常使用 SSE 指令，x87 FPU 主要用于遗留支持。如果有帮助，我可以显示 expf(float) 的代码作为示例。

标签： c algorithm math compiler-construction cryptography

【解决方案1】：

我假设您的兴趣是实现指数函数，这些函数可以在 HLL 的标准数学库中找到，特别是 C/C++。其中包括函数exp()、exp2()、exp10() 和pow()，以及单精度对应函数expf()、exp2f()、exp10f() 和powf()。

您提到的求幂方法（例如 k-ary、滑动窗口）通常用于加密算法中，例如基于求幂的 RSA。它们通常不用于通过math.h 或cmath 提供的求幂函数。 exp() 等标准数学函数的实现细节有所不同，但通用方案遵循三个步骤：

将函数参数简化为初级近似值间隔
在主逼近区间上逼近合适的基函数
将主区间的结果映射回函数的整个范围

辅助步骤通常是处理特殊情况。这些可能与特殊的数学情况有关，例如 log(0.0)，或特殊的浮点操作数，例如 NaN（非数字）。

下面expf(float) 的 C99 代码以示例方式显示了这些步骤对于具体示例的样子。参数a 首先被拆分为exp(a) = e^r * 2ⁱ，其中i 是一个整数，r 在 [log( sqrt(0.5), log(sqrt(2.0)]，主要近似区间。在第二步中，我们现在用多项式近似 e^r。可以根据各种设计标准设计这种近似值，例如作为最小化绝对或相对误差。多项式可以通过多种方式进行评估，包括霍纳方案和埃斯特林方案。

下面的代码使用了一种非常常见的方法，即采用了极小极大近似，该方法可以最大限度地减少整个近似区间内的最大误差。计算这种近似值的标准算法是 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，它计算 [-1,1] 上的 2^x-1。这可用于exp() 和exp2() 计算的第二步。在第三步中有一条指令FSCALE 用于乘以2ⁱ。实现F2XM1 本身的一种常见方式是使用有理或多项式近似的微码。请注意，现在维护 x87 FPU 主要是为了支持旧版。在现代 x86 平台库上，通常使用基于 SSE 的纯软件实现和类似于如下所示的算法。有些将小表与多项式近似值结合在一起。

pow(x,y) 可以在概念上实现为exp(y*log(x))，但是当x 接近于一且y 的幅度很大时，精度会显着降低，以及对许多特殊情况的错误处理C/C++ 标准。解决精度问题的一种方法是以某种形式的扩展精度计算log(x) 和乘积y*log(x))。详细信息将填写一个完整的、冗长的单独答案，我没有方便的代码来演示它。在各种 C/C++ 数学库中，pow(double,int) 和 powf(float, int) 由单独的代码路径计算，该路径应用“平方和乘法”方法，对整数指数的二进制表示进行按位扫描。

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

/* Like rintf(), but -0.0f -> +0.0f, and |a| must be < 2**22 */
float quick_and_dirty_rintf (float a)
{
    const 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))]. */
float expf_poly (float a)
{ 
    float r;

    r =             0x1.694000p-10f;  // 1.37805939e-3
    r = fmaf (r, a, 0x1.125edcp-07f); // 8.37312452e-3
    r = fmaf (r, a, 0x1.555b5ap-05f); // 4.16695364e-2
    r = fmaf (r, a, 0x1.555450p-03f); // 1.66664720e-1
    r = fmaf (r, a, 0x1.fffff6p-02f); // 4.99999851e-1
    r = fmaf (r, a, 0x1.000000p+00f); // 1.00000000e+0
    r = fmaf (r, a, 0x1.000000p+00f); // 1.00000000e+0
    return r;
}

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

    r =             0x1.418000p-13f;  // 1.53303146e-4
    r = fmaf (r, a, 0x1.5efa94p-10f); // 1.33887795e-3
    r = fmaf (r, a, 0x1.3b2c6cp-07f); // 9.61833261e-3
    r = fmaf (r, a, 0x1.c6af8ep-05f); // 5.55036329e-2
    r = fmaf (r, a, 0x1.ebfbe0p-03f); // 2.40226507e-1
    r = fmaf (r, a, 0x1.62e430p-01f); // 6.93147182e-1
    r = fmaf (r, a, 0x1.000000p+00f); // 1.00000000e+0
    return r;
}

/* Approximate exp10(a) on [log(sqrt(0.5))/log(10), log(sqrt(2.0))/log(10)] */
float exp10f_poly (float a)
{ 
    float r;

    r =             0x1.a56000p-3f;  // 0.20574951
    r = fmaf (r, a, 0x1.155aa8p-1f); // 0.54170728
    r = fmaf (r, a, 0x1.2bda96p+0f); // 1.17130411
    r = fmaf (r, a, 0x1.046facp+1f); // 2.03465796
    r = fmaf (r, a, 0x1.53524ap+1f); // 2.65094876
    r = fmaf (r, a, 0x1.26bb1cp+1f); // 2.30258512
    r = fmaf (r, a, 0x1.000000p+0f); // 1.00000000
    return r;
}

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

    t = a * 0x1.715476p+0f;            // 1/log(2); 1.442695
    t = quick_and_dirty_rintf (t);
    i = (int)t;
    r = fmaf (t, -0x1.62e400p-01f, a); // log_2_hi; -6.93145752e-1
    r = fmaf (t, -0x1.7f7d1cp-20f, r); // log_2_lo; -1.42860677e-6
    t = expf_poly (r);
    r = ldexpf (t, i);
    if (a < -105.0f) r = 0.0f;
    if (a >  105.0f) r = INFINITY;     // +INF
    return r;
}

/* Compute exponential base 2. Maximum ulp error = 0.86770 */
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 = INFINITY;     // +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;            // log2(10); 3.321928
    t = quick_and_dirty_rintf (t);
    i = (int)t;
    r = fmaf (t, -0x1.344140p-2f, a);  // log10(2)_hi // -3.01030159e-1
    r = fmaf (t, 0x1.5ec10cp-23f, r);  // log10(2)_lo //  1.63332601e-7
    t = exp10f_poly (r);
    r = ldexpf (t, i);
    if (a < -46.0f) r = 0.0f;
    if (a >  46.0f) r = INFINITY;      // +INF
    return r;
}

#include <string.h>
#include <stdint.h>

uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

uint64_t double_as_uint64 (double a)
{
    uint64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

double floatUlpErr (float res, double ref)
{
    uint64_t i, j, err, refi;
    int expoRef;
    
    /* ulp error cannot be computed if either operand is NaN, infinity, zero */
    if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
        (res == 0.0f) || (ref == 0.0f)) {
        return 0.0;
    }
    /* Convert the float result to an "extended float". This is like a float
       with 56 instead of 24 effective mantissa bits.
    */
    i = ((uint64_t)float_as_uint32(res)) << 32;
    /* Convert the double reference to an "extended float". If the reference is
       >= 2^129, we need to clamp to the maximum "extended float". If reference
       is < 2^-126, we need to denormalize because of the float types's limited
       exponent range.
    */
    refi = double_as_uint64(ref);
    expoRef = (int)(((refi >> 52) & 0x7ff) - 1023);
    if (expoRef >= 129) {
        j = 0x7fffffffffffffffULL;
    } else if (expoRef < -126) {
        j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
        j = j >> (-(expoRef + 126));
    } else {
        j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
        j = j | ((uint64_t)(expoRef + 127) << 55);
    }
    j = j | (refi & 0x8000000000000000ULL);
    err = (i < j) ? (j - i) : (i - j);
    return err / 4294967296.0;
}

#include <stdio.h>
#include <stdlib.h>

int main (void)
{
    double ref, ulp, maxulp;
    float arg, res, reff;
    uint32_t argi, resi, refi, diff, sumdiff;

    printf ("testing expf ...\n");
    argi = 0;
    sumdiff = 0;
    maxulp = 0;
    do {
        arg = uint32_as_float (argi);
        res = my_expf (arg);
        ref = exp ((double)arg);
        ulp = floatUlpErr (res, ref);
        if (ulp > maxulp) maxulp = ulp;
        reff = (float)ref;
        refi = float_as_uint32 (reff);
        resi = float_as_uint32 (res);
        diff = (resi < refi) ? (refi - resi) : (resi - refi);
        if (diff > 1) {
            printf ("!! expf: arg=%08x res=%08x ref=%08x\n", argi, resi, refi);
            return EXIT_FAILURE;
        } else {
            sumdiff += diff;
        }
        argi++;
    } while (argi);
    printf ("expf   maxulp=%.5f  sumdiff=%u\n", maxulp, sumdiff);
    
    printf ("testing exp2f ...\n");
    argi = 0;
    maxulp = 0;
    sumdiff = 0;
    do {
        arg = uint32_as_float (argi);
        res = my_exp2f (arg);
        ref = exp2 ((double)arg);
        ulp = floatUlpErr (res, ref);
        if (ulp > maxulp) maxulp = ulp;
        reff = (float)ref;
        refi = float_as_uint32 (reff);
        resi = float_as_uint32 (res);
        diff = (resi < refi) ? (refi - resi) : (resi - refi);
        if (diff > 1) {
            printf ("!! expf: arg=%08x res=%08x ref=%08x\n", argi, resi, refi);
            return EXIT_FAILURE;
        } else {
            sumdiff += diff;
        }
        argi++;
    } while (argi);
    printf ("exp2f  maxulp=%.5f  sumdiff=%u\n", maxulp, sumdiff);

    printf ("testing exp10f ...\n");
    argi = 0;
    maxulp = 0;
    sumdiff = 0;
    do {
        arg = uint32_as_float (argi);
        res = my_exp10f (arg);
        ref = exp10 ((double)arg);
        ulp = floatUlpErr (res, ref);
        if (ulp > maxulp) maxulp = ulp;
        reff = (float)ref;
        refi = float_as_uint32 (reff);
        resi = float_as_uint32 (res);
        diff = (resi < refi) ? (refi - resi) : (resi - refi);
        if (diff > 1) {
            printf ("!! expf: arg=%08x res=%08x ref=%08x\n", argi, resi, refi);
            return EXIT_FAILURE;
        } else {
            sumdiff += diff;
        }
        argi++;
    } while (argi);
    printf ("exp10f maxulp=%.5f  sumdiff=%u\n", maxulp, sumdiff);
    

    return EXIT_SUCCESS;
}

【讨论】：