Python logarithm speed

A few weeks ago a friend of mine brought up something to do with logarithms in a group chat. One thing led to another, and out of curiosity I timed a couple of Python’s built-in log functions:

>>> timeit.timeit('[math.log10(rand) for rand in r]',
... setup='import math;import random;r = [random.random() for _ in range(10000000)]',
... number=1)
2.0042254191357642
>>> timeit.timeit('[math.log(rand) for rand in r]',
... setup='import math;import random;r = [random.random() for _ in range(10000000)]', 
... number=1)
2.345342932967469

Surprisingly log base 10 is 14.5% faster than natural log. Why is this? I did a quick Google for ‘python log speed’ and got a bunch of unrelated articles on logging, so decided to take a look¹.

As some background, the majority of Python installations are backed by an interpreter called CPython². This shouldn’t be confused with Cython, which is a C-like extension that facilitates more performant code; CPython is the C that basically makes everything written in Python happen. This is where we’ll be digging for answers.

There are a couple of hoops to jump through to get to the code that actually calculates the logs. We start with the math module definition in cmathmodule.c. Module definitions are the entry point into CPython for each of their methods – each lives as a #define containing its name in Python (e.g. log, log10, tanh) and a pointer to the C function implementing it. Both of these are scattered throughout cmathmodule.c.h³, although most of the logic in this file is error handling around other functions back in cmathmodule.c which do the heavy lifting. Those are the ones we’re interested in, and I’ve reproduced the code for log below:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
static PyObject *
cmath_log_impl(PyObject *module,
               Py_complex x,
               PyObject *y_obj)
{
   Py_complex y;

   errno = 0;
   x = c_log(x);
   if (y_obj != NULL) {
       y = PyComplex_AsCComplex(y_obj);
       if (PyErr_Occurred()) {
           return NULL;
       }
       y = c_log(y);
       x = _Py_c_quot(x, y);
   }
   if (errno != 0)
       return math_error();
   return PyComplex_FromCComplex(x);
}

We can form a hypothesis by reading through this code. Even though c_log (line 7) returns the natural logarithm, this code also handles the general case of calculating the log in any base. As such, even if we’re just looking to calculate \(ln\) the conditional on line 8 must be checked, slowing things down regardless of whether it triggers.

This setup allows the Python math.log function to optionally take a base. If one is given then c_log is called again to do a base conversion (lines 13-14). We can see what this does to run time:

>>> timeit.timeit('[math.log(rand, 5) for rand in r]',
... setup='import math;import random;r = [random.random() for _ in range(10000000)]',
... number=1)
2.868814719840884

As expected, we get an increase. With this all established, how does log10 manage to be faster? Let’s take a look:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
static Py_complex
cmath_log10_impl(PyObject *module,
                 Py_complex z)
{
   Py_complex r;
   int errno_save;

   r = c_log(z);
   errno_save = errno;
   r.real = r.real / M_LN10;
   r.imag = r.imag / M_LN10;
   errno = errno_save;
   return r;
}

This code still uses c_log, but does the base conversion using a constant value of M_LN10 = \(log_e(10)\). It’s not great to come to conclusions on performance without profiling⁴, but it looks like doing away with the conditional and extra c_log call causes the difference.

What’s the lesson here? If performance really matters to you, then using math.log10 and doing the base conversion yourself will be fastest: the log base 5 example took 2.20 seconds (23% faster) when calculated as math.log10(rand) / log_10_5 rather than math.log(rand, 5). Surprisingly, using NumPy’s log function as a drop-in replacement took 10.81 seconds, but the real lesson is that if you really need the log of 10 million numbers at once then vectorisation is your friend:

>>> timeit.timeit('numpy.log(r)', 
... setup='import numpy;import random;r = [random.random() for _ in range(10000000)]',
... number=1)
0.6032462348230183