定点 Log2 近似

Question

我已经使用查找表和低阶多项式近似实现了定点 log2 函数，但对整个 32 位定点范围 [-1,+1) 的精度不太满意。输入格式为s0.31，输出格式为s15.16。

我在这里发布这个问题，以便其他用户可以发布他的答案（一些评论在另一个线程中交换，但他们更喜欢在单独的线程中提供全面的答案）。欢迎任何其他答案，如果您能提供算法及其实现的一些速度与准确性的详细信息，我将不胜感激。

谢谢。

Answer 1

通过简单地计算定点数中的前导零位x，就可以确定log2(x)最接近的严格较小的整数。在许多处理器架构上，存在“计数前导零”机器指令或内在指令。如果不可用，clz()可以通过多种方式构建相当有效的实现，其中之一包含在下面的代码中。

\n\n

为了计算对数的小数部分，两个主要明显的竞争者是表中的插值和极小极大多项式近似。在这种特定情况下，在相当小的表中进行二次插值似乎是更具吸引力的选择。x = 2 ⁱ * (1+f)，其中 0 \xe2\x89\xa4 f < 1。我们i如上所述确定并使用 的前导位f对表进行索引。通过此表项和下面的两个表项拟合抛物线，动态计算抛物线的参数。结果被四舍五入，并应用启发式调整来部分补偿定点算术的截断性质。最后，将整数部分相加，得到最终结果。

\n\n

应该注意的是，计算涉及有符号整数的右移，该整数可能为负数。我们需要这些右移来映射到机器代码级别的算术右移，这是ISO-C 标准无法保证的。然而，实际上大多数编译器都会做所需的事情。在本例中，我在运行 Windows 的 x64 平台上使用了 Intel 编译器。

\n\n

利用 32 位字的 66 项表，最大绝对误差可降至 8.18251e-6，从而s15.16实现完全精度。

\n\n

#include <stdio.h>\n#include <stdlib.h>\n#include <stdint.h>\n#include <math.h>\n\n#define FRAC_BITS_OUT (16)\n#define INT_BITS_OUT  (15)\n#define FRAC_BITS_IN  (31)\n#define INT_BITS_IN   ( 0)\n\n/* count leading zeros: intrinsic or machine instruction on many architectures */\nint32_t clz (uint32_t x)\n{\n    uint32_t n, y;\n\n    n = 31 + (!x);\n    if ((y = (x & 0xffff0000U))) { n -= 16;  x = y; }\n    if ((y = (x & 0xff00ff00U))) { n -=  8;  x = y; }\n    if ((y = (x & 0xf0f0f0f0U))) { n -=  4;  x = y; }\n    if ((y = (x & 0xccccccccU))) { n -=  2;  x = y; }\n    if ((    (x & 0xaaaaaaaaU))) { n -=  1;         }\n    return n;\n}\n\n#define LOG2_TBL_SIZE (6)\n#define TBL_SIZE      ((1 << LOG2_TBL_SIZE) + 2)\n\n/* for i = [0,65]: log2(1 + i/64) * (1 << 31) */\nconst uint32_t log2Tab [TBL_SIZE] =\n{\n    0x00000000, 0x02dcf2d1, 0x05aeb4dd, 0x08759c50, \n    0x0b31fb7d, 0x0de42120, 0x108c588d, 0x132ae9e2, \n    0x15c01a3a, 0x184c2bd0, 0x1acf5e2e, 0x1d49ee4c, \n    0x1fbc16b9, 0x22260fb6, 0x24880f56, 0x26e2499d, \n    0x2934f098, 0x2b803474, 0x2dc4439b, 0x30014ac6, \n    0x32377512, 0x3466ec15, 0x368fd7ee, 0x38b25f5a, \n    0x3acea7c0, 0x3ce4d544, 0x3ef50ad2, 0x40ff6a2e, \n    0x43041403, 0x450327eb, 0x46fcc47a, 0x48f10751, \n    0x4ae00d1d, 0x4cc9f1ab, 0x4eaecfeb, 0x508ec1fa, \n    0x5269e12f, 0x5440461c, 0x5612089a, 0x57df3fd0, \n    0x59a80239, 0x5b6c65aa, 0x5d2c7f59, 0x5ee863e5, \n    0x60a02757, 0x6253dd2c, 0x64039858, 0x65af6b4b, \n    0x675767f5, 0x68fb9fce, 0x6a9c23d6, 0x6c39049b, \n    0x6dd2523d, 0x6f681c73, 0x70fa728c, 0x72896373, \n    0x7414fdb5, 0x759d4f81, 0x772266ad, 0x78a450b8, \n    0x7a231ace, 0x7b9ed1c7, 0x7d17822f, 0x7e8d3846, \n    0x80000000, 0x816fe50b\n};\n\n#define RND_SHIFT     (31 - FRAC_BITS_OUT)\n#define RND_CONST     ((1 << RND_SHIFT) / 2)\n#define RND_ADJUST    (0x10d) /* established heuristically */\n\n/* \n   compute log2(x) in s15.16 format, where x is in s0.31 format\n   maximum absolute error 8.18251e-6 @ 0x20352845 (0.251622232)\n*/   \nint32_t fixed_log2 (int32_t x)\n{\n    int32_t f1, f2, dx, a, b, approx, lz, i, idx;\n    uint32_t t;\n\n    /* x = 2**i * (1 + f), 0 <= f < 1. Find i */\n    lz = clz (x);\n    i = INT_BITS_IN - lz;\n    /* normalize f */\n    t = (uint32_t)x << (lz + 1);\n    /* index table of log2 values using LOG2_TBL_SIZE msbs of fraction */\n    idx = t >> (32 - LOG2_TBL_SIZE);\n    /* difference between argument and smallest sampling point */\n    dx = t - (idx << (32 - LOG2_TBL_SIZE));\n    /* fit parabola through closest three sampling points; find coeffs a, b */\n    f1 = (log2Tab[idx+1] - log2Tab[idx]);\n    f2 = (log2Tab[idx+2] - log2Tab[idx]);\n    a = f2 - (f1 << 1);\n    b = (f1 << 1) - a;\n    /* find function value for argument by computing ((a*dx+b)*dx) */\n    approx = (int32_t)((((int64_t)a)*dx) >> (32 - LOG2_TBL_SIZE)) + b;\n    approx = (int32_t)((((int64_t)approx)*dx) >> (32 - LOG2_TBL_SIZE + 1));\n    approx = log2Tab[idx] + approx;\n    /* round fractional part of result */\n    approx = (((uint32_t)approx) + RND_CONST + RND_ADJUST) >> RND_SHIFT;\n    /* combine integer and fractional parts of result */\n    return (i << FRAC_BITS_OUT) + approx;\n}\n\n/* convert from s15.16 fixed point to double-precision floating point */\ndouble fixed_to_float_s15_16 (int32_t a)\n{\n    return a / 65536.0;\n}\n\n/* convert from s0.31 fixed point to double-precision floating point */\ndouble fixed_to_float_s0_31 (int32_t a)\n{\n    return a / (65536.0 * 32768.0);\n}\n\nint main (void)\n{\n    double a, res, ref, err, maxerr = 0.0;\n    int32_t x, start, end;\n\n    start = 0x00000001;\n    end =   0x7fffffff;\n    printf ("testing fixed_log2 with inputs in [%17.10e, %17.10e)\\n",  \n            fixed_to_float_s0_31 (start), fixed_to_float_s0_31 (end));\n\n    for (x = start; x < end; x++) {\n        a = fixed_to_float_s0_31 (x);\n        ref = log2 (a);\n        res = fixed_to_float_s15_16 (fixed_log2 (x));\n        err = fabs (res - ref);\n        if (err > maxerr) {\n            maxerr = err;\n        }\n    }\n\n    printf ("max. err = %g\\n", maxerr);\n    return EXIT_SUCCESS;\n}\n

\n\n

为了完整起见，我在下面展示了极小极大多项式近似。这种近似的系数可以通过几种工具生成，例如 Maple、Mathematica、Sollya，或者使用使用 Remez 算法的自制代码生成，这就是我在这里使用的。下面的代码显示了原始浮点系数、用于最大限度提高中间计算精度的动态缩放，以及用于减轻非舍入定点算术影响的启发式调整。

\n\n

计算 的典型方法log2(x)是使用 x = 2 ⁱ * (1+f) 并使用 log2(1+f) 的近似值来计算 [\xe2\x88\x9a\xc2\xbd, \xe2 \x88\x9a2]，这意味着我们p(f)在主近似区间 [\xe2\x88\x9a\xc2\xbd-1, \xe2\x88\x9a2-1] 上使用多项式。

\n\n

在我们希望使用 32 位操作作为其基本构建块的限制下，中间计算尽可能扩大操作数以提高准确性，mulhi因为这是许多 32 位架构上的本机指令，可以通过内联机器访问代码或作为内在函数。与基于表的代码一样，有符号数据的右移可能为负，并且此类右移必须映射到算术右移，ISO-C 不保证这一点，但大多数 C 编译器都会这样做。

\n\n

我设法将此变体的最大绝对误差降至 1.11288e-5，因此几乎完全s15.16准确，但比基于表格的变体稍差。我怀疑我应该在多项式中添加一项附加项。

\n\n

/* on 32-bit architectures, there is often an instruction/intrinsic for this */\nint32_t mulhi (int32_t a, int32_t b)\n{\n    return (int32_t)(((int64_t)a * (int64_t)b) >> 32);\n}\n\n#define RND_SHIFT  (25 - FRAC_BITS_OUT)\n#define RND_CONST  ((1 << RND_SHIFT) / 2)\n#define RND_ADJUST (-2) /* established heuristically */\n\n/* \n    compute log2(x) in s15.16 format, where x is in s0.31 format\n    maximum absolute error 1.11288e-5 @ 0x5a82689f (0.707104757)\n*/   \nint32_t fixed_log2 (int32_t x)\n{\n    int32_t lz, i, f, p, approx;\n    uint32_t t;\n    /* x = 2**i * (1 + f), 0 <= f < 1. Find i */\n    lz = clz (x);\n    i = INT_BITS_IN - lz;\n    /* force (1+f) into range [sqrt(0.5), sqrt(2)] */\n    t = (uint32_t)x << lz;    \n    if (t > (uint32_t)(1.414213562 * (1U << 31))) {\n        i++;\n        t = t >> 1;\n    }\n    /* compute log2(1+f) for f in [-0.2929, 0.4142] */\n    f = t - (1U << 31);\n    p =              + (int32_t)(-0.206191055 * (1U << 31) -  1);\n    p = mulhi (p, f) + (int32_t)( 0.318199910 * (1U << 30) - 18);\n    p = mulhi (p, f) + (int32_t)(-0.366491705 * (1U << 29) + 22);\n    p = mulhi (p, f) + (int32_t)( 0.479811855 * (1U << 28) -  2);\n    p = mulhi (p, f) + (int32_t)(-0.721206390 * (1U << 27) + 37);\n    p = mulhi (p, f) + (int32_t)( 0.442701618 * (1U << 26) + 35);\n    p = mulhi (p, f) + (f >> (31 - 25));\n    /* round fractional part of the result */\n    approx = (p + RND_CONST + RND_ADJUST) >> RND_SHIFT;\n    /* combine integer and fractional parts of result */\n    return (i << FRAC_BITS_OUT) + approx;\n}\n

\n