Mon*_*nti 4 c compiler-construction algorithm math cryptography
我一直在学习更快的取幂算法（k-ary、滑动门等），并且想知道 CPU/编程语言中使用了哪些？（我不清楚这是在 CPU 中还是通过编译器发生的）
只是为了踢，哪个是最快的？
关于广泛性的编辑：这是故意广泛的，因为我知道有很多不同的技术可以做到这一点。检查的答案有我想要的。
我假设您的兴趣是实现可以在 HLL 的标准数学库中找到的幂函数，特别是 C/C++。这些包括的功能exp()，exp2()，exp10()，和pow()，以及单精度对应expf()，exp2f()，exp10f()，和powf()。
您提到的取幂方法（例如k-ary、滑动窗口）通常用于加密算法，例如基于取幂的RSA。它们通常不用于通过math.h或提供的幂函数cmath。标准数学函数的实现细节exp()不同，但通用方案遵循三步过程：
辅助步骤通常是处理特殊情况。这些可能与特殊的数学情况有关，例如log(0.0)，或特殊的浮点操作数，例如 NaN（非数字）。
下面的 C99 代码以expf(float)示例方式显示了这些步骤在具体示例中的样子。参数a首先被拆分，使得exp(a)= e r * 2 i，其中i是一个整数，r并且在 [log(sqrt(0.5), log(sqrt(2.0)]] 中，主要的近似区间。在第二步中，我们现在近似e r与多项式。可以根据各种设计标准来设计此类近似值，例如最小化绝对或相对误差。可以以各种方式评估多项式，包括霍纳方案和埃斯特林方案。
下面的代码使用了一种非常常见的方法，即采用极小极大逼近，从而将整个逼近区间的最大误差降至最低。计算这种近似值的标准算法是 Remez 算法。评估是通过霍纳的方案；使用 提高了该评估的数值精度fmaf()。
这个标准数学函数实现了所谓的融合乘加法或 FMA。这将在加法过程中a*b+c使用完整的乘积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接近统一并且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;
}