我想你的兴趣在于实现指数函数,这些函数可以在HLL的标准数学库中找到,特别是C/C++。这些功能包括exp()
、exp2()
、exp10()
和pow()
,以及单精度对应项expf()
、exp2f()
、exp10f()
和powf()
。
您提到的指数化方法(如k元、滑动窗口)通常用于密码算法,如基于指数化的rsa。它们通常不用于通过math.h
或cmath
提供的指数函数。标准数学函数如exp()
的实现细节有所不同,但通用方案遵循三步过程:
函数参数约化到一次近似
间隔
一次近似区间上适当基函数的逼近
将主间隔的结果映射回函数的整个范围
辅助步骤通常是处理特殊情况的方法。它们可以适用于特殊的数学情况,如log(0.0)
,或特殊的浮点操作数,如nan(不是数字)。
下面的expf(float)
的C99代码以示例性的方式显示了这些步骤对于一个具体的示例是什么样子的参数a
首先拆分为exp(a)
=Er*2i,其中i
为整数,r
位于[log(qRT(0.5)),log(qRT(2)]中,主近似区间。在第二步中,我们用多项式逼近ER。这种近似可以根据各种设计准则设计,例如最小化绝对误差或相对误差。多项式可以用多种方法计算,包括horner方案和estrin方案。
下面的代码使用一个非常普遍的方法,采用极小极大近似,它最小化整个近似区间上的最大误差。计算这种近似的一个标准算法是ReMez算法。通过horner格式进行计算;使用fmaf()
可提高计算精度。
这个标准的数学函数实现了所谓的融合乘法加法或fma。这将在加法过程中使用全积a*b+c
计算a*b
,并在最后应用单个舍入。在大多数现代硬件上,如GPU、IBM Power CPU、最新的x86处理器(如Haswell)、最新的ARM处理器(作为可选扩展),这直接映射到硬件指令。在缺少这种指令的平台上,fmaf()
将映射到相当慢的仿真代码,在这种情况下,如果我们对性能感兴趣,我们就不想使用它。
最后的计算是2i的乘法,其中C和C++提供了函数ldexp()
。在“industrial strength”库代码中,这里通常使用特定于机器的习惯用法,该习惯用法利用了IEEE-754二进制算法最后,代码将清除溢出和下溢的情况。
x86处理器中的x87 FPU有一条指令float
在[-1,1]上计算2x-1。这可用于计算F2XM1
和exp()
的第二步。有一条指令exp2()
用于在第三步中乘以2i。实现FSCALE
本身的一种常见方式是使用有理或多项式逼近的微码。请注意,目前维护x87 FPU的主要目的是为了旧版支持。在现代x86平台上,库通常使用基于sse的纯软件实现和类似于下面所示的算法。一些结合小表与多项式近似。
F2XM1
可以在概念上被实现为pow(x,y)
,但是当exp(y*log(x))
接近统一和大的大小,以及在C/C++标准中指定的许多特殊情况的错误处理时,这会受到显著的精度损失。解决精度问题的一种方法是以某种扩展精度的形式计算x
和乘积y
详细信息将填充一个完整的、冗长的单独答案,我没有现成的代码来演示它。在各种C/C++数学库中,log(x)
和y*log(x))
是由一个单独的代码路径计算的,该代码路径应用“平方乘”的方法,对整数指数的二进制表示进行逐位扫描。
#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;
}