CPU/编程语言使用哪些指数算法？

Question

我一直在学习更快速的指数算法（如k-ary、滑动门等），不知道哪些算法被用于CPU/编程语言中？（我不确定这是处理器还是编译器实现的）

只是出于好奇，哪个是最快的？

关于广度的编辑：我有意使问题比较宽泛，因为我知道有许多不同的技术可以做到这一点。选定答案就是我想要的。

- Monti

---

当你询问不同的CPU或编程语言时，你不能合理地询问“哪个算法”。也许不同的系统使用不同的算法。你需要更具体一些。 - Greg Hewgill

@GregHewgill 抱歉，语法错误。我只是在寻找一个例子，比如“x86使用滑动门”之类的东西。 - Monti

在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)的代码示例。 - njuffa

Answer 1

我假设你对实现标准数学库中的指数函数有兴趣，特别是C/C++中的函数，包括exp()、exp2()、exp10()和pow()，以及单精度版本的expf()、exp2f()、exp10f()和powf()。

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

1.将函数参数缩小到主逼近区间 2.在主逼近区间上逼近合适的基本函数 3.将主区间结果映射回整个函数范围

辅助步骤通常是处理特殊情况。这些可能涉及特殊的数学情况，如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。这样的逼近可以根据各种设计标准进行设计，例如最小化绝对误差或相对误差。多项式可以以各种方式进行评估，包括Horner方案和Estrin方案。

下面的代码使用了一种非常常见的方法，采用极小值逼近，该方法在整个逼近区间上最小化最大误差。计算这样的逼近的标准算法是Remez算法。评估通过Horner方案实现；这种评估的数值精度通过使用fmaf()得到增强。

这个标准数学函数实现了所谓的融合乘加或FMA。这使用完整乘积a*b在加法期间进行计算，并在最后应用单个舍入。在大多数现代硬件上，如GPU、IBM Power CPU、最新的x86处理器（例如Haswell）、最新的ARM处理器（作为可选扩展），这直接映射到硬件指令。在缺少这样的指令的平台上，fmaf()将映射到相当慢的仿真代码，这种情况下，如果我们关心性能，就不想使用它。

最终计算的是2的i次幂乘以系数，C和C++提供了函数ldexp()。在“工业级”库代码中，通常使用一种机器特定的构造方法，在此利用IEEE-754二进制算术来处理float类型。最后，该代码清除了溢出和下溢的情况。

x86处理器内的x87 FPU具有指令F2XM1，可以计算[-1, 1]上的2的x次幂减1。这可用于计算exp()和exp2()的第二步。有一个指令FSCALE，用于在第三步中乘以2的i次幂。实现F2XM1本身的常见方式是将其作为利用有理或多项式逼近的微码。请注意，现在x87 FPU主要用于遗留支持。在现代x86平台库中，通常使用基于SSE和类似于下面显示的算法的纯软件实现。其中一些结合了小表和多项式逼近。 pow(x，y)可以在概念上实现为exp(y*log(x))，但当x接近1且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;
}

- njuffa