我假设您感兴趣的是求幂函数的实现,这些函数可以在HLL的标准数学库中找到,尤其是C / C
++。这些包括的功能
exp(),
exp2(),
exp10(),和
pow(),以及单精度对应
expf(),
exp2f(),
exp10f(),和
powf()。
您提到的幂计算方法(例如k进制,滑动窗口)通常用于基于密码的加密算法(例如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))],主要近似区间中。在第二步中,我们现在近似Ë
ř用多项式。这样的近似值可根据各种设计标准,如最小化绝对或相对误差进行设计。多项式可以通过各种方式,包括霍纳方案和埃斯特林的方案进行评估。
下面的代码通过采用minimax逼近来使用一种非常常见的方法,该方法将整个逼近间隔内的最大误差最小化。用于计算这种近似值的标准算法是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> float quick_and_dirty_rintf (float a){ float cvt_magic = 0x1.800000p+23f; return (a + cvt_magic) - cvt_magic;}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;}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;}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;}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;}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;}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;}


