固定点的 2 近似值的幂

Question

目前，我正在使用一个小型查找表和线性插值，它非常快且足够准确（最大误差小于 0.001）。但是我想知道是否有更快的近似值。

由于可以通过位移提取和计算指数的整数部分，因此近似值只需要在[-1,1]范围内工作 我试图找到一个切比雪夫多项式，但对于低阶多项式无法达到良好的精度。我猜我可以忍受大约 0.01 的最大误差，但我没有接近这个数字。高阶多项式不是一种选择，因为它们的效率远低于我目前基于查找表的解决方案。

Answer 1

由于没有说明特定的定点格式，我将演示使用相当常用的s15.16 定点算法来替代表查找的一种可能方法。基本思想是将输入a拆分为整数部分i和小数部分f，使得f在[-0.5,0.5]中，然后对exp2(f)使用极小极大多项式逼近[-0.5, 0.5] 并根据i进行最终缩放。

可以使用 Mathematica、Maple 或 Sollya 等工具生成 Minimax 近似值。如果这些工具均不可用，则可以使用 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;
}