不动点的 2 次方近似
Power of 2 approximation in fixed point
目前,我正在使用小型查找 table 和线性插值,速度非常快,也足够准确(最大误差小于 0.001)。但是我想知道是否有一个更快的近似值。
由于指数的整数部分可以通过移位来提取和计算,因此近似值只需要在[-1,1]范围内工作
我试图找到切比雪夫多项式,但无法为低阶多项式实现良好的精度。我猜我可以忍受 0.01 左右的最大误差,但我没有接近那个数字。高阶多项式不是一种选择,因为它们比我当前基于查找 table 的解决方案效率低得多。
由于没有说明具体的定点格式,我将演示使用 s15.16
定点算法替代 table 查找的可能方法,这是相当常用的。基本思想是将输入 a
拆分为整数部分 i
和小数部分 f
,使得 f
在 [-0.5,0.5] 中,然后使用exp2(f)
在 [-0.5, 0.5] 上的极小极大多项式逼近,并根据 i
.
执行最终缩放
可以使用 Mathematica、Maple 或 Sollya 等工具生成极小极大近似值。如果这些工具中的 none 可用,则可以使用 Remez 算法的自定义实现来生成极小极大近似值。
应使用霍纳方案来计算多项式。由于使用了定点运算,因此多项式的计算应在中间步骤中将操作数缩放到最大可能范围(即无溢出)以优化计算精度。
下面的 C 代码假定右移应用于有符号整数数据类型会导致算术移位操作,因此负操作数会适当移位。这是 ISO C 标准 不 保证的,但根据我的经验,它可以很好地与各种工具链一起工作。在最坏的情况下,可以使用内联汇编来强制生成所需的算术右移指令。
下面 fixed_exp2()
实施中包含的测试输出应如下所示:
testing fixed_exp2 with inputs in [-5.96484, 15)
max. rel. err = 0.000999758
这表明区间 [-5.96484, 15] 中的输入满足所需的误差范围 0.001。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
/* compute exp2(a) in s15.16 fixed-point arithmetic, -16 < a < 15 */
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 + 0x8000) & ~0xffff; // 0.5
f = a - i;
s = ((15 << 16) - i) >> 16;
/* minimax approximation for exp2(f) on [-0.5, 0.5] */
r = 0x00000e20; // 5.5171669058037949e-2
r = (r * f + 0x3e1cc333) >> 17; // 2.4261112219321804e-1
r = (r * f + 0x58bd46a6) >> 16; // 6.9326098546062365e-1
r = r * f + 0x7ffde4a3; // 9.9992807353939517e-1
return (uint32_t)r >> s;
}
double fixed_to_float (int32_t a)
{
return a / 65536.0;
}
int main (void)
{
double a, res, ref, err, maxerr = 0.0;
int32_t x, start, end;
start = 0xfffa0900;
end = 0x000f0000;
printf ("testing fixed_exp2 with inputs in [%g, %g)\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) / ref;
if (err > maxerr) {
maxerr = err;
}
}
printf ("max. rel. err = %g\n", maxerr);
return EXIT_SUCCESS;
}
目前,我正在使用小型查找 table 和线性插值,速度非常快,也足够准确(最大误差小于 0.001)。但是我想知道是否有一个更快的近似值。
由于指数的整数部分可以通过移位来提取和计算,因此近似值只需要在[-1,1]范围内工作 我试图找到切比雪夫多项式,但无法为低阶多项式实现良好的精度。我猜我可以忍受 0.01 左右的最大误差,但我没有接近那个数字。高阶多项式不是一种选择,因为它们比我当前基于查找 table 的解决方案效率低得多。
由于没有说明具体的定点格式,我将演示使用 s15.16
定点算法替代 table 查找的可能方法,这是相当常用的。基本思想是将输入 a
拆分为整数部分 i
和小数部分 f
,使得 f
在 [-0.5,0.5] 中,然后使用exp2(f)
在 [-0.5, 0.5] 上的极小极大多项式逼近,并根据 i
.
可以使用 Mathematica、Maple 或 Sollya 等工具生成极小极大近似值。如果这些工具中的 none 可用,则可以使用 Remez 算法的自定义实现来生成极小极大近似值。
应使用霍纳方案来计算多项式。由于使用了定点运算,因此多项式的计算应在中间步骤中将操作数缩放到最大可能范围(即无溢出)以优化计算精度。
下面的 C 代码假定右移应用于有符号整数数据类型会导致算术移位操作,因此负操作数会适当移位。这是 ISO C 标准 不 保证的,但根据我的经验,它可以很好地与各种工具链一起工作。在最坏的情况下,可以使用内联汇编来强制生成所需的算术右移指令。
下面 fixed_exp2()
实施中包含的测试输出应如下所示:
testing fixed_exp2 with inputs in [-5.96484, 15)
max. rel. err = 0.000999758
这表明区间 [-5.96484, 15] 中的输入满足所需的误差范围 0.001。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
/* compute exp2(a) in s15.16 fixed-point arithmetic, -16 < a < 15 */
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 + 0x8000) & ~0xffff; // 0.5
f = a - i;
s = ((15 << 16) - i) >> 16;
/* minimax approximation for exp2(f) on [-0.5, 0.5] */
r = 0x00000e20; // 5.5171669058037949e-2
r = (r * f + 0x3e1cc333) >> 17; // 2.4261112219321804e-1
r = (r * f + 0x58bd46a6) >> 16; // 6.9326098546062365e-1
r = r * f + 0x7ffde4a3; // 9.9992807353939517e-1
return (uint32_t)r >> s;
}
double fixed_to_float (int32_t a)
{
return a / 65536.0;
}
int main (void)
{
double a, res, ref, err, maxerr = 0.0;
int32_t x, start, end;
start = 0xfffa0900;
end = 0x000f0000;
printf ("testing fixed_exp2 with inputs in [%g, %g)\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) / ref;
if (err > maxerr) {
maxerr = err;
}
}
printf ("max. rel. err = %g\n", maxerr);
return EXIT_SUCCESS;
}