2^x 的定点逼近，输入范围为 s5.26答案

【问题标题】：Fixed point approximation of 2^x, with input range of s5.262^x 的定点逼近，输入范围为 s5.26
【发布时间】：2018-12-12 05:52:01
【问题描述】：

如何使用 exp2() 的极小多项式逼近来实现 2^x 定点算术 s5.26 并且输入值在 [-31.9, 31.9] 范围内如何使用以下链接中提到的 Sollya 工具生成多项式 Power of 2 approximation in fixed point

【问题讨论】：

标签： fixed-point approximation

【解决方案1】：

由于定点算术通常不包括表示溢出结果的“无穷大”编码，因此 exp2() 对 s5.26 格式的任何实现都将限于区间 (-32, 5) 中的输入，从而导致在 [0, 32) 中输出。

超越函数的计算通常包括参数缩减、核心逼近、最终结果构造。对于exp2(a)，合理的参数缩减方案是将a拆分为整数部分i和小数部分f，这样a == i + f，f在[-0.5, 0.5]中。然后计算exp2(f)，并将结果缩放2ⁱ，这对应于定点算术中的移位：exp2(a) = exp2(f) * exp2(i)。

计算exp2(f) 的常见设计选择是在exp2() 的列表值中插值，或多项式逼近。由于最大参数需要 31 个结果位，因此精确插值可能需要使用二次插值来保持表大小合理。由于许多现代处理器（包括嵌入式系统中使用的处理器）都提供了快速整数乘法器，因此我将在这里重点介绍多项式逼近。为此，我们需要一个具有极小极大属性的多项式，即与参考相比，最大误差最小化的多项式。

商业和免费工具都提供了生成极小极大近似值的内置功能，例如Mathematica 的MiniMaxApproximation 命令、Maple 的minimax 命令和Sollya 的fpminimax 命令。也可以选择基于Remez algorithm 构建自己的基础架构，这是我使用的方法。与通常使用最近舍入或舍入舍入的浮点算术相反，定点算术通常仅限于截断中间结果。这在表达式评估期间增加了额外的错误。因此，尝试基于启发式搜索对生成的近似值的系数进行小的调整以部分平衡那些累积的单边误差通常是一个好主意。

因为我们在结果中需要多达 31 位，并且因为核心近似中的系数通常小于 1，所以我们不能使用原生定点精度（此处为 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 数据缓存的大小，这可能会或可能不会提高性能。一种可能的方法是将 f 的 2^f-1 列在 [0, 1) 中。将函数参数拆分为整数i 和分数f，这样f 在[0, 1) 中。为了使表格保持合理的小，请使用二次插值，并根据三个连续的表格条目动态计算多项式的系数。通过启发式确定的偏移量对结果进行了轻微调整，以在一定程度上补偿定点算术的截断性质。

该表由分数f 的前导位索引。使用 7 位作为索引（产生一个包含 128+2 个条目的表），准确度比之前的 minimax 多项式逼近略差。最大绝对误差为 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;
}

【讨论】：

我只是想知道您是否将类似的概念应用于定点 log2 函数（结果为 s15.16 格式）。
@Ali：我过去有过（大约 15 年前，在 ARM 处理器上的浮点单元普及之前）。
@Ali 你总是可以问一个问题，这就是这个网站的用途:-)
有时我会惊讶于人们如此接近于找出多项式 LUT 却又停下来。您可以使用预先计算的多项式系数的 LUT，而不是使用原始值的 LUT，然后从 3 个这样的连续值推导二次多项式，因此核心近似将只是计算索引，然后进行 2 次乘法和 2 次加法。只需决定每个 LUT 条目的范围代表什么，然后在该范围内拟合一个多项式，并将系数存储在其中。我多年来一直在使用这种技术，有时甚至使用五分表，它很棒。
如果您寻找“分段多项式近似”，您应该能够找到相关文献。我见过的实现没有使用三次多项式（样条）以上的任何东西。我根据this question 的线性近似对一个进行了逆向工程。正如我在这里展示的那样，与动态系数计算的权衡似乎是总表存储量的增加，这在我的无 FP 嵌入式系统（~200 MHz ARM32 ) 没有那么吸引人。