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

我一直在学习更快的取幂算法（k-ary，滑动门等），并想知道在CPU /编程语言中使用哪些算法？ （我对这是否发生在CPU或编译器中都很模糊）

只是为了踢，这是最快的？

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

Answer 1

我假设您对实现可以在HLL的标准数学库中找到的取幂函数感兴趣，特别是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算法。 评估是通过霍纳的计划; 使用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接近于1且y在大量上时，以及对许多特殊情况的error handling时，这会严重损失精度。在C / C ++标准中规定。 解决准确性问题的一种方法是以某种forms的扩展精度计算log(x)和乘积y*log(x)) 。 细节将填充整个冗长的单独答案，我没有代码方便地演示它。 在各种C / C ++数学库中， pow(double,int)和powf(float, int)由单独的代码路径计算，该路径应用“square-and-multiply”方法，逐位扫描二进制表示forms。整数指数。

 #include  /* 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; }