How do I compute a fixed point binary logarithm?

Question

The goal is this: given an unsigned (64 bit) integer, return the binary logarithm of it as a Q64.64 fixed point decimal. This means the lower 64 bits represent the fractional part of the logarithm when taken as a ratio with 2^64. I only need 5, maybe 6 bits of precision here.

It's easy to use the "count leading zeroes" trick to compute the integer part. However, I'm having trouble with the fractional part.

Rust is the recommended language as it has an unsigned 128 bit integer.

When I write the following equation, in which f is the fractional part and w is the whole part (where b is easily computable):

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

I get the following simplification, which doesn't really help:

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

Since the number we're taking the logarithm of would be getting larger, not smaller. There is no recursive base case.

So how would you compute this fractional part of the logarithm without using floating point numbers or other libraries?

This code is kind of doing what I want, that is calculating the log base 2 of in that case a fixed point number: https://github.com/Uniswap/v3-core/blob/d8b1c635c275d2a9450bd6a78f3fa2484fef73eb/contracts/libraries/TickMath.sol#L112

I'm not sure how it works though. Is anyone able to understand?

Answer 1

If you really only need 5 to 6 bits of precision, you can get away with approximating the logarithm by a second-order polynomial in the interval [1, 2). This gives at least seven bits of precision, and is trivial to implement:

pub fn log2_1p_approx(x: u64) -> u64 {
    let y = x as u128 * x.wrapping_neg() as u128 >> 64;
    x + (y * 6269126026662125184 >> 64) as u64
}

Both the argument and the return value are the fractional part of a fixed-point binary number, and the function approximates log2(1 + x). The function computes the interpolation polynomial that returns exact values for x = 0.0, x = 0.5 and x = 1.0. I obtained the integer constant in the function by computing the exact value of log2_1p(0.5) and simplifying the polynomial.

By comparing the results to the exact values, I confirmed that this gives 7 bits of precision at worst.

Note that this can easily be further optimized. We don't need the 128 bit multiplications, since the precision is rather low anyway.

Answer 2

I previously showed and explained how to compute the binary logarithm of a pure fraction (that is, 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). One therefore splits the integer input x like so: x = a * 2**expo, and computes expo via clz and log2(a) using the algorithm I demonstrated previously.

I do not know Rust, so the code below was written in ISO-C99, which I hope can be translated to Rust in a straightforward manner. The clz functionality is implemented manually here. Obviously one would want to avail oneself of appropriate built-ins for this where available.

The specification given by asker stated that the result be returned in UQ64.64 format which seemed unusual to me, so the ilog2() function here returns a UQ16.16 result. If a result in UQ64.64 format is definitely required, a conversion is easily accomplished by converting the result of ilog2() to a __uint128_t and performing a left shift by 48 on he result.

I have included some light-weight test scaffolding to demonstrate the correct operation of the code. Using built-ins for clz should result in about 50 to 60 machine instructions being generated on common CPU architectures like x86-64 and ARM64 when using the latest compilers, so performance should be sufficient for various use cases (the asker did not mention specific performance goals).

#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;
}