Find log2 of a number with Newton’s method

approximationcalculusnewton raphson

I am learning C++ and I find pow, log, log2, exp functions provided by cmath extremely slow, so that I would like to implement faster versions of them.

I currently have a version using polynomial approximations, which is much faster than the ones provided by the library and very accurate, but not as accurate as the library ones.

I wrote another version using lookup tables, it takes less than a nanosecond to complete and is by far the fastest, and it is also the least accurate:

#include <cmath>
#include <vector>
using std::vector;

vector<double> log2_double() {
    int lim = 1 << 24;
    vector<double> table(lim);
    for (uint64_t i = 0; i < lim; i++) {
        table[i] = log2(i << 29) - 1075;
    }
    return table;
}
const vector<double> LOG2_DOUBLE = log2_double();

inline double fast_log2_double(double d) {
    uint64_t bits = std::bit_cast<uint64_t>(d);
    uint64_t e = (bits >> 52) & 0x7ff;
    uint64_t m = bits & 0xfffffffffffff;
    return (e == 0 ? LOG2_DOUBLE[m >> 28] : e + LOG2_DOUBLE[(m | 0x10000000000000) >> 29]);
}

However it is really fast, and accurate to 6 decimal digits, so I would like to refine the value using Newton's method:

fast_log2(3.1415927f) = 1.651496887207031 : 0.32501220703125 nanoseconds
log2f(3.1415927f) = 1.651496171951294 : 33.90207290649414 nanoseconds

Newton's method from Wikipedia:

$$
x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}
$$

Derivative of $\log_2(x)$:

$$
\frac{d}{dx}\log_2(x) = \frac{1}{x \cdot \ln(2)}
$$

$\ln(2)$ is a constant, and the derivative is simple to compute.

The equation I am trying to solve is $\log_2(x_0) – x = 0$, a simple function is $f(x) = 2^x – x_0$.

Putting it all together:

$$
x_{n + 1} = x_{n} – (2^{x_n} – x_0) \cdot x_n \cdot \ln(2)
$$

I tried it in Python and it diverges:

import numpy as np

LOG2_DOUBLE = np.arange(1 << 24, dtype=np.uint64)
LOG2_DOUBLE = np.log2(LOG2_DOUBLE << 29) - 52
LN2 = np.log(2)

def fast_log2(f: float, lim: int) -> float:
    m, e = f.hex().split("p")
    e = int(e)
    m = int(m[4:], 16)
    l = LOG2_DOUBLE[m >> 28] if not e else e + LOG2_DOUBLE[(m | 0x10000000000000) >> 29]
    for _ in range(lim):
        l = l - (2 ** l - f) * l * LN2
    return l

In [839]: n = np.random.random()

In [840]: n
Out[840]: 0.14112379908067973

In [841]: np.log2(n)
Out[841]: -2.8249667907196945

In [842]: fast_log2(n, 0)
Out[842]: -2.824966918629066

In [843]: fast_log2(n, 16)
Out[843]: -2.82496890252251

In [844]: fast_log2(n, 128)
Out[844]: -99.49155122185653

What did I do wrong?

Best Answer

In your iterative scheme:

$$ š‘„_{n+1} = š‘„_nāˆ’(2^{x_n}āˆ’š‘„_0)ā‹…(x_n\ln(2)), $$

you use a function $f(x) = 2^x-x_0$. Then, the derivative of $f(x)$ is

$$ f'(x) = 2^x\ln(2), $$

Which means that your iterative scheme should be:

$$ x_{n+1} = x_nāˆ’\frac{2^{x_n}āˆ’š‘„_0}{2^{x_n}\ln(2)}. $$

Related Question