Why are log2 and log1p so much faster than log and log10, in numpy?

Question

Whilst playing around with this question I noticed something I couldn't explain regarding the relative performance of np.log2, np.log and np.log10:

In [1]: %%timeit x = np.random.rand(100000)    ....: np.log2(x)    ....:  1000 loops, best of 3: 1.31 ms per loop  In [2]: %%timeit x = np.random.rand(100000) np.log(x)    ....:  100 loops, best of 3: 3.64 ms per loop  In [3]: %%timeit x = np.random.rand(100000) np.log10(x)    ....:  100 loops, best of 3: 3.93 ms per loop 

np.log2 is about 3x faster than np.log and np.log10. Perhaps even more counter-intuitively, np.log1p(x), which computes ln(x + 1), is on par with np.log2:

In [4]: %%timeit x = np.random.rand(100000) np.log1p(x)    ....:  1000 loops, best of 3: 1.46 ms per loop 

I obtained almost identical timings in numpy v1.10.1 and v1.8.2.

Is there an intuitive explanation for these discrepancies in runtime performance?

Answer 1

I found your question extremely interesting, so I spent a few hours doing a bit of reserach; I think I have found an explanation for the performance difference as it applies to integers (thank you Matteo Italia for your note) - It is unclear how that reasoning can be extended to floats though:

Computers use base 2 - According to the articles linked in reference, the computation of log2 is a 4 processor cycles process - that of log10 requires to multiply log2(val) by 1/log2(10) which adds another 5 cycles.

Finding log2 is a matter of finding the index of the least significant bit of a value. (video at about the 23rd minute mark onward).

bit hacks: Find integer log base 10 of an integer

bit hacks: Find the log base 2 of an N-bit integer in O(lg(N))

The integer log base 10 is computed by first using one of the techniques above for finding the log base 2. By the relationship log10(v) = log2(v) / log2(10), we need to multiply it by 1/log2(10), which is approximately 1233/4096, or 1233 followed by a right shift of 12. Adding one is needed because the IntegerLogBase2 rounds down. Finally, since the value t is only an approximation that may be off by one, the exact value is found by subtracting the result of v < PowersOf10[t].

This method takes 6 more operations than IntegerLogBase2. It may be sped up (on machines with fast memory access) by modifying the log base 2 table-lookup method above so that the entries hold what is computed for t (that is, pre-add, -mulitply, and -shift). Doing so would require a total of only 9 operations to find the log base 10, assuming 4 tables were used (one for each byte of v).

Of note: the use of DeBruijn sequences lookup & bit shifting techniques to calculate log2 in this MIT video: Lec 2 | MIT 6.172 Performance Engineering of Software Systems, Fall 2010 (video at the 36th minute mark onward).

Of note this StackOverflow post which demonstrates a method to make efficient log2 calculations truly cross platform with C++

Caveat: I have not checked numpy source code to verify that it indeed implements similar techniques, but it would be surprising that it did not. In fact, from the comments under the OP's post, Fermion Portal have checked:

Actually numpy uses math.h from glibc, you'll see the same difference in C/C++ if you use math.h/cmath.h. You can find the nicely commented source code for the two functions e.g. ic.unicamp.br/~islene/2s2008-mo806/libc/sysdeps/ieee754/dbl-64/… and ic.unicamp.br/~islene/2s2008-mo806/libc/sysdeps/ieee754/dbl-64/… – Fermion Portal[9]