tho*_*han 1 fixed-point approximation
如何使用 exp2() 的极小极大多项式近似来实现 2^x 定点算术 s5.26 并且输入值在 [-31.9, 31.9] 范围内 如何使用以下链接中提到的 Sollya 工具生成多项式 Power of 2 定点逼近
由于定点算术通常不包括“无限”编码代表溢出的结果，任何实施exp2()为一个s5.26格式将被限制为在区间（-32，5）的输入，从而导致输出在[0，32）。
超越函数的计算通常包括参数约简、核逼近、最终结果构造。在 的情况下exp2(a)，合理的参数缩减方案是将 拆分a为整数部分i和小数部分f，使得a == i + f，f在 [-0.5, 0.5] 中。然后计算exp2(f)，并将结果按比例缩放 2 i，这对应于定点算术中的移位：exp2(a) = exp2(f) * exp2(i)。
计算 的常见设计选择是 的exp2(f)列表值中的插值exp2()或多项式近似。由于最大参数需要 31 个结果位，因此精确插值可能需要使用二次插值来保持表大小合理。由于许多现代处理器（包括嵌入式系统中使用的处理器）提供了快速整数乘法器，因此我将在此重点介绍多项式逼近。为此，我们需要一个具有极小极大属性的多项式，也就是说，与参考相比，可以最小化最大误差的多项式。
商业和免费工具都提供了生成极小极大近似的内置功能，例如 Mathematica 的MiniMaxApproximation命令、Maple 的minimax命令和 Sollya 的fpminimax命令。也可能选择基于Remez 算法构建自己的基础设施，这是我使用的方法。与通常使用最接近或偶数舍入的浮点运算相反，定点运算通常仅限于截断中间结果。这会在表达式计算过程中增加额外的错误。因此，尝试基于启发式搜索对生成的近似值的系数进行小幅调整以部分平衡那些累积的单边误差通常是一个好主意。
因为我们在结果中需要多达 31 位，并且因为核心近似中的系数通常在幅度上小于统一，所以我们不能使用本机定点精度（此处s5.26为 ）进行多项式评估。相反，我们希望通过动态调整我们正在使用的定点格式来扩大中间计算中的操作数以充分利用 32 位整数的可用范围。乘法使用重新归一化右移 32 位。这通常可以消除 32 位处理器上的显式移位。
由于中间计算使用有符号数据，因此会发生有符号负操作数的右移。我们希望这些右移映射到算术右移指令，这是 C 标准不能保证的。但是在大多数常用的平台上，C 编译器会做我们想要的事情。否则，可能需要求助于内部函数或内联汇编。我在 x64 平台上使用 Microsoft 编译器开发了以下代码。
在评估exp2(f)原始浮点系数的多项式近似时，动态缩放和启发式调整都清晰可见。下面的代码对于大参数并不能完全达到准确度。最大绝对误差为 1.10233e-7，对于0x12de9c5b= 4.71739332的参数：fixed_exp2()返回0x693ab6a3而准确结果为0x693ab69c。据推测，可以通过将多项式核心近似的次数增加 1 来实现完全准确。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
/* on 32-bit architectures, there is often an instruction/intrinsic for this */
int32_t mulhi (int32_t a, int32_t b)
{
return (int32_t)(((int64_t)a * (int64_t)b) >> 32);
}
/* compute exp2(a) in s5.26 fixed-point arithmetic */
int32_t fixed_exp2 (int32_t a)
{
int32_t i, f, r, s;
/* split a = i + f, such that f in [-0.5, 0.5] */
i = (a + 0x2000000) & ~0x3ffffff; // 0.5
f = a - i;
s = ((5 << 26) - i) >> 26;
f = f << 5; /* scale up for maximum accuracy in intermediate computation */
/* approximate exp2(f)-1 for f in [-0.5, 0.5] */
r = (int32_t)(1.53303146e-4 * (1LL << 36) + 996);
r = mulhi (r, f) + (int32_t)(1.33887795e-3 * (1LL << 35) + 99);
r = mulhi (r, f) + (int32_t)(9.61833261e-3 * (1LL << 34) + 121);
r = mulhi (r, f) + (int32_t)(5.55036329e-2 * (1LL << 33) + 51);
r = mulhi (r, f) + (int32_t)(2.40226507e-1 * (1LL << 32) + 8);
r = mulhi (r, f) + (int32_t)(6.93147182e-1 * (1LL << 31) + 5);
r = mulhi (r, f);
/* add 1, scale based on integral portion of argument, round the result */
r = ((((uint32_t)r * 2) + (uint32_t)(1.0*(1LL << 31)) + ((1U << s) / 2) + 1) >> s);
/* when argument < -26.5, result underflows to zero */
if (a < -0x6a000000) r = 0;
return r;
}
/* convert from s5.26 fixed point to double-precision floating point */
double fixed_to_float (int32_t a)
{
return a / 67108864.0;
}
int main (void)
{
double a, res, ref, err, maxerr = 0.0;
int32_t x, start, end;
start = -0x7fffffff; // -31.999999985
end = 0x14000000; // 5.000000000
printf ("testing fixed_exp2 with inputs in [%.9f, %.9f)\n",
fixed_to_float (start), fixed_to_float (end));
for (x = start; x < end; x++) {
a = fixed_to_float (x);
ref = exp2 (a);
res = fixed_to_float (fixed_exp2 (x));
err = fabs (res - ref);
if (err > maxerr) {
maxerr = err;
}
}
printf ("max. abs. err = %g\n", maxerr);
return EXIT_SUCCESS;
}
基于表的替代方案将权衡表存储以减少执行的计算量。根据 L1 数据缓存的大小，这可能会也可能不会提高性能。一种可能的方法是将[0, 1) 中的2 f -1制成表格f。将函数参数拆分为一个整数i和一个分数f，使得f在 [0, 1) 中。为了使表格保持合理小，请使用二次插值，从三个连续的表格条目中动态计算多项式的系数。结果通过启发式确定的偏移量稍微调整，以在一定程度上补偿定点算术的截断性质。
该表由分数的前导位索引f。使用 7 位作为索引（产生 128+2 个条目的表），精度比之前的极小极大多项式近似略差。最大绝对误差为 1.74935e-7。它发生在参数0x11580000= 4.33593750 时，其中fixed_exp2()返回0x50c7d771，而准确的结果是0x50c7d765。
/* For i in [0,129]: (exp2 (i/128.0) - 1.0) * (1 << 31) */
static const uint32_t expTab [130] =
{
0x00000000, 0x00b1ed50, 0x0164d1f4, 0x0218af43,
0x02cd8699, 0x0383594f, 0x043a28c4, 0x04f1f656,
0x05aac368, 0x0664915c, 0x071f6197, 0x07db3580,
0x08980e81, 0x0955ee03, 0x0a14d575, 0x0ad4c645,
0x0b95c1e4, 0x0c57c9c4, 0x0d1adf5b, 0x0ddf0420,
0x0ea4398b, 0x0f6a8118, 0x1031dc43, 0x10fa4c8c,
0x11c3d374, 0x128e727e, 0x135a2b2f, 0x1426ff10,
0x14f4efa9, 0x15c3fe87, 0x16942d37, 0x17657d4a,
0x1837f052, 0x190b87e2, 0x19e04593, 0x1ab62afd,
0x1b8d39ba, 0x1c657368, 0x1d3ed9a7, 0x1e196e19,
0x1ef53261, 0x1fd22825, 0x20b05110, 0x218faecb,
0x22704303, 0x23520f69, 0x243515ae, 0x25195787,
0x25fed6aa, 0x26e594d0, 0x27cd93b5, 0x28b6d516,
0x29a15ab5, 0x2a8d2653, 0x2b7a39b6, 0x2c6896a5,
0x2d583eea, 0x2e493453, 0x2f3b78ad, 0x302f0dcc,
0x3123f582, 0x321a31a6, 0x3311c413, 0x340aaea2,
0x3504f334, 0x360093a8, 0x36fd91e3, 0x37fbefcb,
0x38fbaf47, 0x39fcd245, 0x3aff5ab2, 0x3c034a7f,
0x3d08a39f, 0x3e0f680a, 0x3f1799b6, 0x40213aa2,
0x412c4cca, 0x4238d231, 0x4346ccda, 0x44563ecc,
0x45672a11, 0x467990b6, 0x478d74c9, 0x48a2d85d,
0x49b9bd86, 0x4ad2265e, 0x4bec14ff, 0x4d078b86,
0x4e248c15, 0x4f4318cf, 0x506333db, 0x5184df62,
0x52a81d92, 0x53ccf09a, 0x54f35aac, 0x561b5dff,
0x5744fccb, 0x5870394c, 0x599d15c2, 0x5acb946f,
0x5bfbb798, 0x5d2d8185, 0x5e60f482, 0x5f9612df,
0x60ccdeec, 0x62055b00, 0x633f8973, 0x647b6ca0,
0x65b906e7, 0x66f85aab, 0x68396a50, 0x697c3840,
0x6ac0c6e8, 0x6c0718b6, 0x6d4f301f, 0x6e990f98,
0x6fe4b99c, 0x713230a8, 0x7281773c, 0x73d28fde,
0x75257d15, 0x767a416c, 0x77d0df73, 0x792959bb,
0x7a83b2db, 0x7bdfed6d, 0x7d3e0c0d, 0x7e9e115c,
0x80000000, 0x8163daa0
};
int32_t fixed_exp2 (int32_t x)
{
int32_t f1, f2, dx, a, b, approx, idx, i, f;
/* extract integer portion; 2**i is realized as a shift at the end */
i = (x >> 26);
/* extract fraction f so we can compute 2^f, 0 <= f < 1 */
f = x & 0x3ffffff;
/* index table of exp2 values using 7 most significant bits of fraction */
idx = (uint32_t)f >> (26 - 7);
/* difference between argument and next smaller sampling point */
dx = f - (idx << (26 - 7));
/* fit parabola through closest 3 sampling points; find coefficients a,b */
f1 = (expTab[idx+1] - expTab[idx]);
f2 = (expTab[idx+2] - expTab[idx]);
a = f2 - (f1 << 1);
b = (f1 << 1) - a;
/* find function value offset for argument x by computing ((a*dx+b)*dx) */
approx = a;
approx = (int32_t)((((int64_t)approx)*dx) >> (26 - 7)) + b;
approx = (int32_t)((((int64_t)approx)*dx) >> (26 - 7 + 1));
/* combine integer and fractional parts of result, round result */
approx = (((expTab[idx] + (uint32_t)approx + (uint32_t)(1.0*(1LL << 31)) + 22U) >> (30 - 26 - i)) + 1) >> 1;
/* flush underflow to 0 */
if (i < -27) approx = 0;
return approx;
}