CPU/programming 语言使用哪些求幂算法？

Question

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

为了好玩，哪个最快？

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

Answer 1

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

您提到的求幂方法（例如 k 元、滑动 window）通常用于基于求幂的密码算法，例如 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()。在“工业强度”库代码中，通常在这里使用特定于机器的习语，它利用 float 的 IEEE-754 二进制算法。最后，代码清除了上溢和下溢的情况。

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

CPU/programming 语言使用哪些求幂算法？

Which exponentiation algorithms do CPU/programming languages use?

c

compiler-construction

algorithm

math

cryptography