2^x 的不动点近似值，输入范围为 s5.26
Fixed point approximation of 2^x, with input range of s5.26
如何使用 exp2() 的极小极大多项式逼近实现 2^x 定点算术 s5.26 并且输入值在 [-31.9, 31.9] 范围内
如何使用下面提到的Sollya Tool生成多项式link
由于定点运算通常不包括表示溢出结果的 "infinity" 编码，因此 exp2() 对 s5.26 格式的任何实现都将限于区间内的输入（ -32, 5), 导致输出为 [0, 32).
超越函数的计算通常包括参数约简、核心逼近、最终结果构造。在 exp2(a) 的情况下，合理的参数缩减方案是将 a 拆分为整数部分 i 和小数部分 f，使得 a == i + f，其中 f 在 [-0.5, 0.5] 中。然后计算 exp2(f)，并将结果缩放 2i，这对应于定点算术中的移位：exp2(a) = exp2(f) * exp2(i).
计算 exp2(f) 的常见设计选择是 exp2() 的表格值插值或多项式逼近。由于我们需要 31 个结果位作为最大参数，精确插值可能需要使用二次插值来保持 table 大小合理。由于许多现代处理器（包括嵌入式系统中使用的处理器）提供了快速整数乘法器，因此我将在这里重点介绍多项式逼近。为此，我们需要一个具有 minimax 属性的多项式，即与参考值相比最大误差最小的多项式。
商业和免费工具都提供了生成极小极大近似值的内置功能，例如Mathematica 的 MiniMaxApproximation 命令、Maple 的 minimax 命令和 Sollya 的 fpminimax 命令。人们也可能选择基于 Remez algorithm 构建自己的基础设施，这就是我使用的方法。与通常使用最接近或偶数舍入的浮点运算相反，定点运算通常仅限于截断中间结果。这会在表达式评估期间增加额外的错误。因此，尝试基于启发式搜索对生成的近似值的系数进行小幅调整以部分平衡那些累积的单边误差通常是个好主意。
因为我们在结果中需要最多 31 位，并且因为核心近似中的系数通常在幅度上小于单位，所以我们不能使用本机定点精度，这里 s5.26，用于多项式评估.相反，我们希望通过动态调整我们正在使用的定点格式来扩大中间计算中的操作数，以充分利用 32 位整数的可用范围。出于效率的原因，安排计算似乎是可取的乘法使用重新归一化右移 32 位。这通常允许消除 32 位处理器上的显式移位。
由于中间计算使用带符号的数据，因此将发生带符号的负操作数的右移。我们希望这些右移映射到算术右移指令，这是 C 标准 而不是 保证的。但是在最常用的平台上，C 编译器会做我们想要的事情。否则，可能需要求助于内部函数或内联汇编。我在 x64 平台上使用 Microsoft 编译器开发了以下代码。
在对 exp2(f) 的多项式逼近的评估中，原始浮点系数、动态缩放和启发式调整都清晰可见。下面的代码不能完全准确地处理大参数。最大的绝对误差是 1.10233e-7，对于参数 0x12de9c5b = 4.71739332: fixed_exp2() returns 0x693ab6a3 而准确的结果是 0x693ab69c。大概可以通过将多项式核心近似的次数增加一个来实现完全准确。
#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;
}
基于 table 的替代方案将权衡 table 存储以减少执行的计算量。根据 L1 数据缓存的大小，这可能会也可能不会提高性能。一种可能的方法是将 [0, 1) 中的 f 的 2f-1 制表。将函数参数拆分为一个整数 i 和一个分数 f，使得 f 在 [0, 1) 中。为了使 table 保持合理的小，使用二次插值，从三个连续的 table 项中动态计算多项式的系数。结果通过启发式确定的偏移量进行了轻微调整，以在某种程度上补偿定点算法的截断性质。
table 由小数 f 的前导位索引。使用七位索引（导致 table 的 128+2 个条目），准确性比之前的极小极大多项式近似略差。最大绝对误差为 1.74935e-7。它发生在 0x11580000 = 4.33593750 的参数中，其中 fixed_exp2() returns 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;
}
如何使用 exp2() 的极小极大多项式逼近实现 2^x 定点算术 s5.26 并且输入值在 [-31.9, 31.9] 范围内
如何使用下面提到的Sollya Tool生成多项式link
由于定点运算通常不包括表示溢出结果的 "infinity" 编码，因此 exp2() 对 s5.26 格式的任何实现都将限于区间内的输入（ -32, 5), 导致输出为 [0, 32).
超越函数的计算通常包括参数约简、核心逼近、最终结果构造。在 exp2(a) 的情况下，合理的参数缩减方案是将 a 拆分为整数部分 i 和小数部分 f，使得 a == i + f，其中 f 在 [-0.5, 0.5] 中。然后计算 exp2(f)，并将结果缩放 2i，这对应于定点算术中的移位：exp2(a) = exp2(f) * exp2(i).
计算 exp2(f) 的常见设计选择是 exp2() 的表格值插值或多项式逼近。由于我们需要 31 个结果位作为最大参数，精确插值可能需要使用二次插值来保持 table 大小合理。由于许多现代处理器（包括嵌入式系统中使用的处理器）提供了快速整数乘法器，因此我将在这里重点介绍多项式逼近。为此，我们需要一个具有 minimax 属性的多项式，即与参考值相比最大误差最小的多项式。
商业和免费工具都提供了生成极小极大近似值的内置功能，例如Mathematica 的 MiniMaxApproximation 命令、Maple 的 minimax 命令和 Sollya 的 fpminimax 命令。人们也可能选择基于 Remez algorithm 构建自己的基础设施，这就是我使用的方法。与通常使用最接近或偶数舍入的浮点运算相反，定点运算通常仅限于截断中间结果。这会在表达式评估期间增加额外的错误。因此，尝试基于启发式搜索对生成的近似值的系数进行小幅调整以部分平衡那些累积的单边误差通常是个好主意。
因为我们在结果中需要最多 31 位，并且因为核心近似中的系数通常在幅度上小于单位，所以我们不能使用本机定点精度，这里 s5.26，用于多项式评估.相反，我们希望通过动态调整我们正在使用的定点格式来扩大中间计算中的操作数，以充分利用 32 位整数的可用范围。出于效率的原因，安排计算似乎是可取的乘法使用重新归一化右移 32 位。这通常允许消除 32 位处理器上的显式移位。
由于中间计算使用带符号的数据，因此将发生带符号的负操作数的右移。我们希望这些右移映射到算术右移指令，这是 C 标准 而不是 保证的。但是在最常用的平台上，C 编译器会做我们想要的事情。否则，可能需要求助于内部函数或内联汇编。我在 x64 平台上使用 Microsoft 编译器开发了以下代码。
在对 exp2(f) 的多项式逼近的评估中，原始浮点系数、动态缩放和启发式调整都清晰可见。下面的代码不能完全准确地处理大参数。最大的绝对误差是 1.10233e-7，对于参数 0x12de9c5b = 4.71739332: fixed_exp2() returns 0x693ab6a3 而准确的结果是 0x693ab69c。大概可以通过将多项式核心近似的次数增加一个来实现完全准确。
#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;
}
基于 table 的替代方案将权衡 table 存储以减少执行的计算量。根据 L1 数据缓存的大小，这可能会也可能不会提高性能。一种可能的方法是将 [0, 1) 中的 f 的 2f-1 制表。将函数参数拆分为一个整数 i 和一个分数 f，使得 f 在 [0, 1) 中。为了使 table 保持合理的小，使用二次插值，从三个连续的 table 项中动态计算多项式的系数。结果通过启发式确定的偏移量进行了轻微调整，以在某种程度上补偿定点算法的截断性质。
table 由小数 f 的前导位索引。使用七位索引（导致 table 的 128+2 个条目），准确性比之前的极小极大多项式近似略差。最大绝对误差为 1.74935e-7。它发生在 0x11580000 = 4.33593750 的参数中，其中 fixed_exp2() returns 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;
}