CPU /编程语言使用哪种幂运算算法？

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

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

1. 将函数自变量减少到一次近似间隔
2. 在主近似区间上近似合适的基函数
3. 将主要间隔的结果映射回函数的整个范围

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

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

下面的代码通过使用minimax逼近来使用一种非常常见的方法，该方法将整个逼近间隔内的最大误差最小化。用于计算这种近似值的标准算法是Remez算法。通过霍纳方案进行评估；通过使用可以提高此评估的数值准确性fmaf()。

这个标准的数学函数实现了所谓的融合乘加或FMA。这将在加法期间a*b+c使用完整乘积a*b进行计算，并在最后应用一次舍入。在大多数现代硬件上，例如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() */

/* 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))]. */
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.a5a000p-3f;  // 0.20587158
    r = fmaf (r, a, 0x1.155dcap-1f); // 0.54173118
    r = fmaf (r, a, 0x1.2bda68p+0f); // 1.17130136
    r = fmaf (r, a, 0x1.046fa8p+1f); // 2.03465748
    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 = 1.0f/0.0f;    // +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 = 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;            // log2(10); 3.321928
    t = quick_and_dirty_rintf (t);
    i = (int)t;
    r = fmaf (t, -0x1.344136p-2f, a);  // log10(2)_hi; 3.01030010e-1
    r = fmaf (t, 0x1.ec10c0p-27f, r);  // log10(2)_lo; 1.43209888e-8
    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;
}