How to get the calculate the binary logarithm of a fixed-point number with 18 decimals?

Question

Suppose that:

1. You're using a programming language that doesn't support fixed-point or floating-point numbers.
2. The said programming language rounds down divisions, e.g. 12 / 5 = 2.
3. You're adding 18 trailing decimals to represent fixed-point numbers in your algorithm. for instance, you'd represent 0.5 as 500000000000000000, equivalent to 5e17 in scientific notation.

Now, say you want to calculate the integer part of the binary logarithm of a number x that adheres to the format above.

If x is greater than or equal to 1e18, then the calculation is easy: you just divide x by 1e18, and calculate the most significant bit of the rounded down result. For example, 4e18 / 1e18 is 4, of which the msb is 2. The integer part is thus 2.

But what if x is between 0 and 1e18 (both exclusive)? For example, 1e17? The binary logarithm of a number between 0 and 1 is negative ..

Answer 1

I figured it out, eventually, by making the following observations:

1. $\log_2{(x)} = -log_2{(\frac{1}{x})}$
2. To accurately calculate $\frac{1}{x}$, multiply $1e18$ by itself then divide by $x$.

My implementation looks something like this (written in Solidity):

pragma solidity ^0.8.0;

function log2(uint256 x) internal pure returns (int256 result) {
    require(x > 0);
    unchecked {
        // This works because log2(x) = -log2(1/x).
        int256 sign;
        if (x >= 1e18) {
            sign = 1;
        } else {
            sign = -1;
            unchecked {
                x = 1e36 / x;
            }
        }

        // Calculate the integer part of the logarithm and add it to the result and finally calculate y = x * 2^(-n).
        uint256 n = mostSignificantBit(uint256(x / 1e18));

        // The integer part of the logarithm as a signed 59.18-decimal fixed-point number. The operation can't overflow
        // beacuse n is maximum 255, SCALE is 1e18 and sign is either 1 or -1.
        result = int256(n) * 1e18;

        // This is y = x * 2^(-n).
        int256 y = int256(x) >> n;

        // If y = 1, the fractional part is zero.
        if (y == 1e18) {
            return result * sign;
        }

        // Calculate the fractional part via the iterative approximation.
        // The "delta >>= 1" part is equivalent to "delta /= 2", but shifting bits is faster.
        for (int256 delta = 5e17; delta > 0; delta >>= 1) {
            y = (y * y) / 1e18;

            // Is y^2 > 2 and so in the range [2,4)?
            if (y >= 2 * 1e18) {
                // Add the 2^(-m) factor to the logarithm.
                result += delta;

                // Corresponds to z/2 on Wikipedia.
                y >>= 1;
            }
        }
        result *= sign;
    }
}

Where mostSignificantBit is a function that calculates the most significant bit of a number that can go as high as $2^{256} -1 $.