如何计算定点二进制对数？

Question

目标是这样的：给定一个无符号（64 位）整数，将其二进制对数作为 Q64.64 定点十进制返回。这意味着低 64 位表示与 2^64 进行比率时对数的小数部分。我这里只需要 5 位，也许 6 位精度。

使用“计算前导零”技巧来计算整数部分很容易。但是，我在处理小数部分时遇到了麻烦。

Rust 是推荐的语言，因为它具有无符号 128 位整数。

当我写出以下等式时，其中

f

是小数部分，

w

是整数部分（其中

b

很容易计算）：

w + f/2^64 = log2(2^w + b/2^w)

我得到以下简化，但这并没有真正帮助：

f = 2^64 * (log2(b + 2^2w) - 2w)

因为我们取对数的数字会变得更大，而不是更小。没有递归基本情况。

那么，在不使用浮点数或其他库的情况下，如何计算对数的小数部分呢？

这段代码是在做我想做的事情，即计算在这种情况下定点数的对数底数2：https://github.com/Uniswap/v3-core/blob/d8b1c635c275d2a9450bd6a78f3fa2484fef73eb/contracts/libraries/ TickMath.sol#L112

但我不确定它是如何工作的。有谁能看懂吗

Answer 1

我之前展示了并解释了如何计算纯分数的二进制对数（即 0 < a < 1) in fixed-point arithmetic. As the asker notes, the integer portion of a binary logarithm can be computed using a count-leading-zero primitive (

clz

）。因此，可以像这样分割整数输入

x

：

x = a * 2**expo

，并使用我之前演示的算法通过

expo

和

clz

计算

log2(a)

。

我不了解 Rust，所以下面的代码是用 ISO-C99 编写的，我希望可以直接将其翻译为 Rust。

clz

功能在这里手动实现。显然，人们希望在可用的情况下利用适当的内置功能。

asker给出的规范指出结果以UQ64.64格式返回，这对我来说似乎很不寻常，所以这里的

ilog2()

函数返回UQ16.16结果。如果确实需要 UQ64.64 格式的结果，可以通过将

ilog2()

的结果转换为

__uint128_t

并对其结果左移 48 来轻松完成转换。

我添加了一些轻量级测试脚手架来演示代码的正确操作。使用

clz

的内置函数应该会在使用最新编译器时在 x86-64 和 ARM64 等常见 CPU 架构上生成大约 50 到 60 条机器指令，因此性能应该足以满足各种用例（提问者没有提到）具体绩效目标）。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <math.h>

#define TEST_LIMIT (1ULL << 25)

uint32_t clz32 (uint32_t x)
{
    uint32_t n = 32;
    uint32_t y;
    
    y = x >> 16; if (y != 0) { n = n - 16; x = y; }
    y = x >>  8; if (y != 0) { n = n -  8; x = y; }
    y = x >>  4; if (y != 0) { n = n -  4; x = y; }
    y = x >>  2; if (y != 0) { n = n -  2; x = y; }
    y = x >>  1; if (y != 0) return n - 2;
    return n - x;
}

uint32_t clz64 (uint64_t a)
{
    uint32_t ahi = (uint32_t)(a >> 32);
    uint32_t alo = (uint32_t)(a);
    uint32_t res;
    if (ahi) {
        res = 0;
    } else {
        res = 32;
        ahi = alo;
    }
    res = res + clz32 (ahi);
    return res;
}

/* compute log2(x) for x in [1,2**64-1]; result is a UQ16.16 accurate to 6 fractional bits */
uint32_t ilog2 (uint64_t x)
{
    uint32_t a, t, one = ((uint32_t)1) << 31; // UQ1.31
    uint32_t r = 0; // UQ16.16
    uint32_t lz, expo;

    /* split: x = a * 2**expo, 0 < a < 1 */
    lz = clz64 (x);
    expo = 64 - lz;
    a = (uint32_t)((x << lz) >> (64 - 31)); // UQ1.31

    /* try factors (1+2**(-i)) with i = 1, ..., 8 */
    if ((t = a + (a >> 1)) < one) { a = t; r += 0x095c0; }
    if ((t = a + (a >> 2)) < one) { a = t; r += 0x0526a; }
    if ((t = a + (a >> 3)) < one) { a = t; r += 0x02b80; }
    if ((t = a + (a >> 4)) < one) { a = t; r += 0x01664; }
    if ((t = a + (a >> 5)) < one) { a = t; r += 0x00b5d; }
    if ((t = a + (a >> 6)) < one) { a = t; r += 0x005ba; }
    if ((t = a + (a >> 7)) < one) { a = t; r += 0x002e0; }
    if ((t = a + (a >> 8)) < one) { a = t; r += 0x00171; }

    r = (expo << 16) - r; // adjust for split of argument
    r = ((r + (1 << 9)) >> (16 - 6)) << (16 - 6); // round to 6 fractional bits

    return r;
}

int main (void)
{
    uint32_t ir;
    uint64_t errloc = 0, ix = 1;
    const double two_to_16 = 65536.0;
    double res, ref;
    double err, maxerr = 0;

    do {
        ir = ilog2 (ix);
        res = (double)ir / two_to_16;
        ref = log2 ((double)ix);
        err = fabs (res - ref);
        if (err > maxerr) {
            maxerr =  err;
            errloc = ix;
        }
        ix += 1;
        if ((ix & 0xffffff) == 0) printf ("\r%" PRIx64, ix);
    } while (ix < TEST_LIMIT);
    printf ("\nmaxerr = %.8f @ %" PRIu64 "\n", maxerr, errloc);
    return EXIT_SUCCESS;
}