Writing/The Unbearable Lightness of Coding
Chapter 24 · Draft · 59 min read
Performance: is Python Really that Slow?
Is Python Really that Slow?
The second most common complaint about Python — after the indentation, of course — is that it's slow. You've probably heard it, too; but how many times have you actually measured it, or even cared to do so? In this chapter, we'll figure out if it's true, and it turns out that the answer is "yes, but only in a very precise sense, and mostly no in the sense that most people mean".
Before we do, however, I'd like to fulfill two promises. In the chapter about code as an object, I claimed that interpreted languages, like Python, can still achieve speeds comparable to, and sometimes surprisingly higher than, those of compiled languages, like C. In the chapter on type hinting, I said that the performance argument for types is something we'd cover later; well, this is later. And performance is not an afterthought, but part of design: a good design is environmentally friendly, in that it doesn't waste the machine's time any more than it does the reader's — and it's performant, or at least, it's designed in such a way that it can be optimized on demand.
So here's the plan: we'll start with measuring, because you can't really fix what you can't see, and anyway, you're probably wrong about what the problem is; we'll look at benchmarking and profiling, both for time and for memory. Then, we'll talk about optimizing orders of magnitude: algorithms and data structures, where real change happens. Then, we'll talk about optimizing constants: the small and unglamorous techniques of allocating less, storing more compactly and copying nothing. Then, we'll leave pure Python behind — calling C through foreign function interfaces, both ctypes and CFFI; writing Python-with-types with Cython; compiling on the fly with numba; and finally, PyPy and its strange cousin, RPython — all on the same running example: a Mandelbrot fractal, which we'll squeeze from 700 milliseconds to 16, and the fastest version won't be the one written in C.
But first, a word of caution about numbers: everything I'm going to say about performance is going to come with numbers, but those numbers are going to be measured on one particular machine — a recent laptop with CPython 3.13, to be precise — so yours will be different. It doesn't really matter, though: don't memorize the milliseconds, just appreciate the ratios between them; they're the story.
And before we jump into the story, another word of caution, about attitude: most of the speedups we'll see in this chapter are not due to cleverness, but due to looking. The best performance tool I know is my own suspicion of my own intuition — and underneath the technical question is a human one that we'll return to at the end: why do we believe the things we believe about performance, and how many of these things did we actually check?
Benchmarking
The term benchmarking comes from land surveying, where a bench mark would be chiselled into stone at a known elevation, into which a bench, or bracket, would be set to rest a levelling rod on, so subsequent measurements can be made with reference to a fixed point. I used to tell a story about how rifles were fixed to such a bench, so the shooter's skill could be measured, but it turns out to be folk etymology, so let's move on.
In software, it means to measure some aspect, or metric, of the system against a reference, so as to identify its bottleneck — by definition, the slowest part, and the only one worth optimizing, since speeding up anything else will not actually make the system faster.
There are different metrics to measure: running time, memory usage, I/O operations, and with the cloud, even money; we can even consider user experience — what do customers complain about the most? — but some metrics can be very hard to define at all: how exactly do you evaluate a neural network that detects chairs, short of labelling a billion images yourself?
The famous computer scientist Donald Knuth had a pithy quote about this: "premature optimization is the root of all evil". It's often quoted out of context, so I'd like to tell you the whole story. What he actually says is that we should forget about small efficiencies, say about 97% of the time — but also that we should not pass up our opportunities in that critical 3%.
So Knuth is not saying, "don't optimize your code" — he's saying, "don't prematurely optimize your code". It's an important nuance, and there are two very good reasons for it.
The first one is that optimizations usually mean adding more code — not less; and if your system doesn't actually need this, then it's just more code to maintain, with no gain to justify the additional complexity.
The second reason, which people tend to forget, is that these "premature" optimizations are almost always done in the wrong place, anyway — because before you measure your system, you don't know what the bottleneck is, now do you?
But Knuth's quote is not an excuse for bad code, either: if you're implementing a sort algorithm, and you know that there's a simple and sensible way to do it that's not $O(n^2)$ — then yeah, maybe don't do it like that. Don't do stupid things, obviously — but also, don't invest time and effort in doing clever things before it's clear that you need them.
Let's say we have a sorting algorithm, and we want to see how well it performs. To take a contrived example, suppose it's an insertion sort, where we go over the elements and insert them in their right place amongst the elements we already considered — much like we sort a hand of cards, going over them one by one, and inserting each card in its right place amongst those we've already placed in our hand (or, if it doesn't fit anywhere, leaving it at the end). Here's how it might look like:
>>> def sort(items):
... result = []
... for item in items:
... for i, anchor in enumerate(result):
... if item < anchor:
... result.insert(i, item)
... break
... else:
... result.append(item)
... return result
>>> sort([3, 1, 2])
[1, 2, 3]
Let's give it a try on 100,000 numbers in reverse order, and time it by hand; my algorithms course taught me that this is the worst case scenario for insertion sort, so let's see:
>>> import time
>>> items = range(100_000, 0, -1)
>>> started = time.time()
>>> result = sort(items)
>>> elapsed = time.time() - started
>>> print(f'{elapsed:.2f} seconds')
1.02 seconds
And it's about a second. Doing it by hand every time is pretty tedious, though, so let's write a context manager to do it for us (remember our chapter on the lifecycle and context management?).
A few words on implementation: we use time.perf_counter, not time.time: the latter is the wall clock, and it can be adjusted from under us — the former is a monotonic clock, and the highest resolution timer we have. The result is in milliseconds, because it's the most common resolution for timings, but feel free to use anything else. Here it is:
>>> class Timer:
... def __enter__(self):
... self.started = time.perf_counter()
... def __exit__(self, exception, error, traceback):
... elapsed = time.perf_counter() - self.started
... print(f'{elapsed * 1000:.3f} ms')
>>> with Timer():
... result = sort(items)
1028.875 ms
For smaller snippets of code, we have the timeit module in the standard library: it executes a piece of code multiple times (with garbage collection disabled, by default) to get a more stable measurement — after all, a single run might coincide with packets being downloaded from the internet, or some other noise.
One thing to note, though, is that it's intended for small snippets, not entire programs: if your program takes 20 seconds, then running it a thousand times and dividing by this number will take a while. Anyway, here's how it might look like:
>>> import timeit
>>> timeit.timeit("''.join(str(i) for i in range(100))", number=100_000)
0.5962906247004867
>>> timeit.timeit("''.join([str(i) for i in range(100)])", number=100_000)
0.48483220813795924
>>> timeit.timeit("''.join(map(str, range(100)))", number=100_000)
0.4935883330181241
And notice how, for 100 items, the list comprehension is actually faster than the generator — probably because creating the generator object and resuming it to generate each item costs more than creating a small list. But don't generalize this: measure it.
We can even use it from the command line, in which case it figures out how many loops to run on its own, and reports the best of five repeats:
$ python -m timeit "'-'.join(map(str, range(100)))"
50000 loops, best of 5: 5.29 usec per loop
As you can see, manual benchmarking is to performance what print-debugging is to bugs — it's great if you already have a suspicion as to what might be wrong, but if you don't, you might want a tool that can automatically watch your entire execution and tell you where it's spending its time.
Profiling
This is known as profiling: a dynamic analysis that builds a profile of a program as it runs — not just how long or how much, but where, when, and on whose behalf. In other words, it's like an FBI profiler, trying to create a psychological portrait of a wanted criminal from its patterns; we want to create a psychological portrait of our program, and hopefully not of whom it murders, but of where it spends most of its time.
Profiling requires the profiler to cooperate with its subject, and its logic actually slows the subject down, so you never quite measure the system as it runs in production; there's something a bit quantum in the act of measuring changing the result. There are several techniques: event-based profiling, where we subscribe to events like function calls and returns, which the runtime must support (Python lets us register a function to be called on every call, which is also how coverage knows which lines were executed); instrumentation, where we inject monitoring code into the subject, either by hand (like with our prints) or automatically, into its source code, binary code, compiler output or even dynamic runtime; or statistical profiling, where we sample the program every so often, which is less accurate (because we might hit a lucky or unlucky streak), but much less intrusive (the py-spy library, for example, samples a running Python process from the outside).
The standard library has its own event-based profiler, cProfile; let's use it on our sort and see what happens:
$ python -m cProfile sort.py
100004 function calls in 1.051 seconds
Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 1.051 1.051 {built-in method builtins.exec}
1 0.001 0.001 1.051 1.051 sort.py:1(<module>)
1 0.027 0.027 1.051 1.051 sort.py:1(sort)
99999 1.024 0.000 1.024 0.000 {method 'insert' of 'list' objects}
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
1 0.000 0.000 0.000 0.000 {method 'append' of 'list' objects}
And that's... something. The ncalls column indicates how many times each function was called, while tottime is how long it took — not including any functions it called, which is captured by cumtime. It's a bit confusing, but the gist is that out of a second, almost all of it was spent on 99,999 calls to list.insert; the algorithm itself, the part that my intuition was so sure about, barely registers.
What's happening? Well, a Python list is actually implemented as a dynamic array of pointers; and while it's fast to append an element at its end (amortized), inserting an element at its beginning is a linear operation, because all the existing pointers need to be shifted one slot to the right. When the input is reversed, we end up inserting every single item at the beginning, which is why it takes so long; and the deeper truth is that for this particular version of insertion sort, having a reversed input is actually the best case scenario as far as comparisons go — because while the textbook insertion sort scans backwards from the end of the sorted part, this one scans forward, so every new item is smaller than the first anchor and breaks out immediately. All this quadratic work actually happens inside insert — but it's written in C and tucked away from sight, so we don't see it in our code.
Luckily, there's a simple fix: we need a data structure into which it's cheap to insert an item at the beginning. Python provides such a data structure, collections.deque, which is a doubly-linked list of fixed-size blocks, and supports constant-time operations at both ends. All we have to do is change one line:
# sort2.py
import collections
def sort(items):
result = collections.deque()
for item in items:
for i, anchor in enumerate(result):
if item < anchor:
result.insert(i, item)
break
else:
result.append(item)
return result
if __name__ == '__main__':
sort(range(100_000, 0, -1))
And we get:
$ python -m cProfile sort2.py
100004 function calls in 0.031 seconds
That's more like it: from about a second to 31 milliseconds, some 35 times faster. If I hadn't bothered to measure it, I would have replaced the algorithm; all we needed to do was change one line. Not knowing how Python works under the hood actually made me use a data structure in exactly the wrong way.
There's another profiling tool that I'm a bit ambivalent about: line_profiler. The idea is nice — instead of reporting per function, it reports per line, because supposedly we think in lines more than we think in function calls. The interface, though, is pretty terrible: we need to decorate our function with @profile, but never actually import it — because kernprof will inject it into our builtins when it runs our file, which will crash with a NameError if we ever try to run it normally. Then, it will generate a .lprof file, which we need to parse with another command:
$ kernprof -l sort_lp.py
Wrote profile results to 'sort_lp.py.lprof'
Inspect results with:
python -m line_profiler -rmt sort_lp.py.lprof
$ python -m line_profiler sort_lp.py.lprof
Timer unit: 1e-06 s
Total time: 1.14704 s
File: sort_lp.py
Function: sort at line 1
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 @profile
2 def sort(items):
3 1 0.0 0.0 0.0 result = []
4 100001 22156.0 0.2 1.9 for item in items:
5 100000 29510.0 0.3 2.6 for i, anchor in enumerate(result):
6 99999 20688.0 0.2 1.8 if item < anchor:
7 99999 1050076.0 10.5 91.5 result.insert(i, item)
8 99999 24608.0 0.2 2.1 break
9 else:
10 1 0.0 0.0 0.0 result.append(item)
11 1 1.0 1.0 0.0 return result
But when it works, the result is actually pretty lovely; 91.5% of the time spent on the insert line. And notice the Hits column: the inner for line ran 100,000 times, and the comparison ran 99,999 times — one comparison per item, just like we figured out ourselves.
As an aside, this interface reminds me of medieval paintings of horses: it's clear that the artists were very skilled and well-intentioned, and yet nobody seems to have noticed that the horse looks like that. The human mind is a wonderful thing, isn't it?
Profiling Memory
Of course, we don't just care how long our code runs, but also how many resources it consumes; and in many cases, the most important resource is memory, especially on servers and other long-running processes. If you notice a process's memory growing, it might indicate a memory leak: that is, some resources that are being allocated and somehow not released, until the machine eventually runs out of memory and the process dies. It's not as common in Python as it is in C, but it can happen — or it might be intentional, like a cache, in which case you want to see the memory growing, because if it doesn't, then it's not being used, which means your predictions about the demand weren't correct. Either way, first things first: we need to measure.
The way we do this in Python is with sys.getsizeof. This returns the size of an object — or at least, the object's own size, in bytes:
>>> import sys
>>> sys.getsizeof(1)
28
>>> sys.getsizeof(2**100)
40
>>> class A:
... def __init__(self, x):
... self.x = x
>>> sys.getsizeof(A(1))
48
>>> sys.getsizeof(list(range(1000)))
8056
Note how small integers take 28 bytes — that's how big an object has to be, in Python, even if it only stores a single digit. A regular int in C would take 4 bytes, or even a single byte; but then, it wouldn't be an object — it wouldn't have a reference count, a type, or a size. Big integers, however, do grow; and that tiny object is actually 48 bytes. But what about the list? It contains a thousand integers, so it should be 28,000 bytes, right? In fact, it's only 8,056: that is, 8 bytes per element, plus some header — because the list actually stores pointers to the integers, and sys.getsizeof is a shallow measurement. The integers themselves are stored elsewhere, and their size isn't taken into account.
To get the full picture, we need tracemalloc. To use it, we simply start it, and it will hook Python's memory allocators, recording every allocation that happens along with the line that made it. When we're ready to see the results, we take a snapshot and ask for statistics, grouped by line:
>>> import tracemalloc
>>> tracemalloc.start()
>>> x = list(range(1_000_000))
>>> snapshot = tracemalloc.take_snapshot()
>>> tracemalloc.stop()
>>> print(snapshot.statistics('lineno')[0])
<stdin>:1: size=38.1 MiB, count=999745, average=40 B
That's 38.1 MiB — or 40 bytes per allocation, of which there are almost a million. This all makes sense: we're dealing with integers, so each one's 28 bytes, rounded up to 32 by the allocator's alignment, plus 8 bytes for the pointer to it, in the list. The number of allocations is a bit less than a million because for small integers, between 0 and 256, Python actually preallocates singletons, which it then reuses instead of allocating over and over again.
What's nice about tracemalloc is that we can also compare two snapshots, in which case it will only show us the difference between them — so it's natural to use it as a context manager:
>>> import contextlib
>>> @contextlib.contextmanager
... def profile_memory():
... tracemalloc.start()
... before = tracemalloc.take_snapshot()
... try:
... yield
... finally:
... after = tracemalloc.take_snapshot()
... tracemalloc.stop()
... stats = after.compare_to(before, 'filename')
... print(f'{sum(stat.size_diff for stat in stats) / 1024:,.1f} KiB')
>>> with profile_memory():
... x = list(range(1_000_000))
39,055.1 KiB
We'll be using this for the rest of the chapter; just note that it's not as reliable as you might expect, since Python's garbage collector runs when it wants, and I've seen it allocate a big object and report less memory than before, simply because it decided to free up some other object at this particular moment.
Optimizing Orders
There are two kinds of change: revolution and evolution. The former is about big, dramatic shifts, while the latter is about gradual, incremental improvements. In Zero to One, Peter Thiel calls them going from 0 to 1, and going from 1 to n. The first is about creating something new, while the second is about copying something that already works, and scaling it up. Algorithms, the course you hopefully took, are about revolutions; not about making your quadratic thing slightly faster, but about making it linear.
Big-O, little-o and Theta describe the behaviour in the grand scheme of things, asymptotically; but you should take them with a grain of salt, because our world is painfully finite. Fibonacci heaps are wonderful asymptotically, but their constant overhead is so large that a plain binary heap usually wins. It's often a trade-off with simplicity, too.
The example everyone uses to illustrate this, and with good reason, is Fibonacci. What's not so beautiful about it is that if you want to calculate it recursively, like it's defined, then we're in for some trouble:
>>> def fib(n):
... return n if n < 2 else fib(n - 1) + fib(n - 2)
>>> with Timer():
... fib(35)
9227465
976.073 ms
It took about a second, and incurred some 30 million calls — because each number recursively computes the two previous numbers, all the way down. This is, quite literally, the natural definition of the Fibonacci sequence, and it's exponential. Here's how it looks iteratively:
>>> def fib(n):
... a, b = 0, 1
... for _ in range(n):
... a, b = b, a + b
... return a
>>> with Timer():
... fib(35)
9227465
0.006 ms
>>> with Timer():
... x = fib(10_000)
1.721 ms
>>> len(str(x))
2090
Which can be done in microseconds. In fact, since Python supports integers of arbitrary precision, we can calculate the 10,000th term, which has 2,090 digits, in less than two milliseconds. This is the classic example for teaching people about loops, because it really drives the point home.
The surprise, then, is that there is a closed-form formula, known as Binet's formula. It uses the golden ratio and its conjugate, so it's not immediately obvious, but once you compute those constants, it's just a couple of powers:
>>> sq5 = 5 ** 0.5
>>> PHI = (1 + sq5) / 2
>>> phi = (1 - sq5) / 2
>>> def binet(n):
... return round((PHI**n - phi**n) / sq5)
>>> with Timer():
... binet(35)
9227465
0.004 ms
>>> binet(35) == fib(35)
True
>>> next(n for n in range(100) if binet(n) != fib(n))
71
>>> binet(1500)
OverflowError: (34, 'Result too large')
It's constant time, and actually quite beautiful — except, these are floating-point numbers, and as such, have only 53 bits of precision. This means that this method can only compute Fibonacci numbers up to fib(70) accurately; from fib(71) it's already wrong, and around fib(1500) it overflows. This is a mathematician's solution: beautiful, but limited by the constraints of the machine.
A student once suggested that there is in fact a way to calculate it even faster. It's a bit underhanded, but really gets the job done: just precompute all the numbers you care about, however long it takes, and then simply look them up:
>>> FIBS = [fib(n) for n in range(1000)]
>>> def fib_table(n):
... return FIBS[n]
>>> with Timer():
... fib_table(32)
2178309
0.002 ms
And that's it — constant time, and exact. It's an interesting solution, but it's also an engineer's solution, not a mathematician's. You see, to a mathematician, it doesn't make sense to define a function as a list of its values; to a computer scientist, it's a time/space trade-off, and it generalizes to any function over a small enough finite domain. We actually did this in the chapter on decorators, in the lazy version, called memoization (which is available in the standard library as functools.cache).
The other half of algorithms is data structures. For example, if we want to check if something is in a list, we have to scan it from the start, which is linear in its length; but if it's in a set (or a dict), which are implemented as hashtables, this can be done in constant time on average:
>>> items = list(range(1_000_000))
>>> targets = 0, 250_000, 500_000, 750_000, 1_000_000
>>> with Timer():
... found = [target for target in targets if target in items]
15.886 ms
>>> items = set(items)
>>> with Timer():
... found = [target for target in targets if target in items]
0.003 ms
>>> found
[0, 250000, 500000, 750000]
This takes 15 milliseconds, but if we put the same numbers in a set, the same checks take just 3 microseconds — thousands of times faster, and most of it is the timer's overhead. If you're repeatedly checking whether some item is in a collection, it's worth considering using a set or a dictionary. A student once pointed out that if we really need to check if some item is in a list, we can at least iterate over it once and check each item against our five targets, right? — and this is true, but this would still take a million steps, instead of five.
Optimizing Constants: Allocations
So we've optimized orders, using various algorithms and data structures. But what if we're already using the standard algorithm, maybe even the state of the art one — and it's still not fast enough? In this case, maybe we're doing something rather wasteful in the way we describe it in Python — and here, the most effective solutions are also the most anti-climactic.
First and foremost, there are allocations: it may seem negligible, but if we're accumulating a million characters by concatenating them to a string, we can actually do so by collecting them in a list, and concatenating everything at once:
>>> with Timer():
... s = ''
... for i in range(1_000_000):
... s += 'a'
9144.369 ms
>>> with Timer():
... s = []
... for i in range(1_000_000):
... s.append('a')
... s = ''.join(s)
73.939 ms
And it makes a difference: 9 seconds for a million s += 'a', 74 milliseconds for append and join. This is because strings are immutable, so every += operation allocates a new string and copies everything so far over to it, which means we're allocating a million strings and copying an ever-increasing amount of data to each — it's quadratic in total. join, on the other hand, calculates the total length and allocates everything in one go.
However, there's a twist: CPython actually has a neat trick, where it can grow the string in-place, as long as it's the only reference to this string that we're working with — but in CPython 3.13, this only works for local variables in a function, which is why it doesn't work at the top level (or in a notebook). If we put it in a function, it suddenly takes 48 milliseconds; and if we add just one more reference to the string, it grows to almost 10 seconds again, because it can no longer be modified behind anyone's back:
>>> def concat(n):
... s = ''
... for i in range(n):
... s += 'a'
... return s
>>> with Timer():
... s = concat(1_000_000)
48.523 ms
>>> def concat_alias(n):
... s = ''
... for i in range(n):
... s += 'a'
... alias = s
... return s
>>> with Timer():
... s = concat_alias(1_000_000)
9292.502 ms
It's not that this optimization is so important, but one that you can't see is one that you can lose without even realizing it. With bytes, this doesn't work at all — and this is actually where it's most important, when we're receiving data from a socket. In this case, accumulating the chunks in a list and joining them at the end can save you a lot of time:
>>> with Timer():
... data = b''
... for i in range(100_000):
... data += b'a'
184.281 ms
>>> with Timer():
... chunks = []
... for i in range(100_000):
... chunks.append(b'a')
... data = b''.join(chunks)
7.826 ms
184 milliseconds for a hundred thousand data += b'a', and 8 milliseconds for a list of chunks and b''.join; it gets even worse the more data you receive. We'll see this pattern — receive the size you still need, append the chunk, raise an error if the connection is closed early, join at the end — again when we talk about socket programming.
Another "data structure" optimization that comes to mind is that a list is cheap at the end, and expensive at the beginning. A deque is cheap on both ends:
>>> items = [0] * 1_000_000
>>> with Timer():
... items.insert(0, 1)
0.285 ms
>>> with Timer():
... items.pop(0)
1
0.188 ms
>>> import collections
>>> items = collections.deque(items)
>>> with Timer():
... items.appendleft(1)
0.002 ms
>>> with Timer():
... items.popleft()
1
0.006 ms
Hundreds of microseconds for every operation on the list, since every pointer has to be shifted; a few microseconds for the deque, which is mostly the timer measuring itself at this point. If you're treating something as a queue, just use a queue.
Finally, there are generators. In the generators chapter, we saw how we can actually compute things on demand, rather than accumulating them in advance. When we want to reverse a list, for example, we could create a new list in reverse order, or simply yield the elements from the end:
>>> items = list(range(1_000_000))
>>> def reverse(items):
... result = []
... for i in range(len(items)):
... result.append(items[len(items) - i - 1])
... return result
>>> with profile_memory():
... r = reverse(items)
8,251.3 KiB
>>> def reverse(items):
... for i in range(len(items)):
... yield items[len(items) - i - 1]
>>> with profile_memory():
... r = reverse(items)
0.8 KiB
>>> with profile_memory():
... for item in reverse(items):
... pass
0.6 KiB
The list version allocates some 8 megabytes of pointers (not the numbers themselves, which it reuses); the generator version allocates almost nothing, not just because nothing is happening when it's created (at which point, no memory is needed to store the result), but also when it's fully consumed, since it never held more than the next item at any given moment. When we iterate over it, it takes pretty much the same time as it did to construct the list; but it doesn't use any memory, as it simply traverses it in reverse order (the same as the built-in reversed function, which is exactly such a lazy iterator).
Optimizing Constants: Compression and Views
These are all fairly simple optimizations; but there are two less trivial ideas, too, and those are compression, or storing data more compactly, and references, or reusing the same data through pointers, so there's less of it to store altogether.
Let's start with compression. Suppose we have a list of a million numbers. Each number takes 28 to 32 bytes, and we need another 8 to point to it, so we end up with a list that's some 39 megabytes big. If we store it in C, though, we can use an unsigned int, which only takes 4 bytes; the array module lets us do just that. It has an annoyingly named array class, but once we wrap our head around this unfortunate decision, it's pretty cool: we pass in a type code, like 'I' for unsigned int, and it knows to store the raw values in a contiguous fashion, while providing list-like functionality:
>>> from array import array
>>> a = array('I', [1, 2, 3])
>>> a.append(4)
>>> a[0], len(a)
(1, 4)
>>> a.tobytes()
b'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00'
>>> with profile_memory():
... numbers = list(range(1_000_000))
39,055.6 KiB
>>> with profile_memory():
... numbers = array('I', range(1_000_000))
3,996.7 KiB
>>> a.append(-1)
OverflowError: can't convert negative value to unsigned int
It's 39 MB versus 4 MB, for the same interface. The price we pay is that it's constrained: it can only hold unsigned integers, so if we try to add a negative value we get an OverflowError. As you can see, we can even serialize it with tobytes, and load it back with frombytes — straight into a compact representation, without ever having to allocate a list. Notice that this is actually another trend we'll see again: the more we know about the type of data we're dealing with, the more we can optimize it. People often argue in favor of adding types to Python in the name of safety, but this is quite a different argument, and in my opinion, a much more interesting one.
But what if we're not dealing with numbers? Every Python object has a __dict__, a hash table mapping any attribute names to any values; __slots__ tells Python what attributes we're going to use, so it can lay them out as a fixed struct instead:
>>> class Point:
... def __init__(self, x, y):
... self.x = x
... self.y = y
>>> with profile_memory():
... points = [Point(i, i) for i in range(1_000_000)]
125,434.0 KiB
>>> class Point:
... __slots__ = 'x', 'y'
... def __init__(self, x, y):
... self.x = x
... self.y = y
>>> with profile_memory():
... points = [Point(i, i) for i in range(1_000_000)]
86,368.2 KiB
>>> p = Point(1, 2)
>>> p.z = 3
AttributeError: 'Point' object has no attribute 'z' and no __dict__ for setting new attributes
>>> p.__dict__
AttributeError: 'Point' object has no attribute '__dict__'. Did you mean: '__dir__'?
That's 125 megabytes down to 86 — and if you're familiar with CPython's internals, you might know that recent versions actually do store attributes in a compact fashion, as long as we don't touch __dict__ (in which case they need to switch to a full hash table), so the gap used to be even bigger. But as you can see, the price is having no new attributes, and no __dict__ either. After all, if you think about it, the whole reason Python is so dynamic is because every object has a __dict__, so in a way, we're giving up the dynamicity that makes Python so cool in exchange for memory savings. But it's not a bad tradeoff, really, if we're dealing with a lot of objects that have a fixed set of fields, like datetime does. In fact, if you look at the pure-Python implementation of that class, you'll see that this is exactly what they did. As for methods, they're not stored on the object either way, but as part of the class.
Now, let's talk about references. Suppose we want to read a 100 MB file into a mutable bytearray: if we do bytearray(reader.read()), this actually allocates the data twice — once as the bytes read returns, and again when bytearray copies them. If we preallocate a bytearray once, however, we can get a file to read directly into it, and avoid both the extra allocation and the copy — saving us about a sixth of the time here, just by knowing that this is even possible:
>>> path = 'file.bin'
>>> with open(path, 'wb') as writer:
... writer.write(os.urandom(100 * 2**20))
104857600
>>> os.path.getsize(path)
104857600
>>> with Timer():
... with open(path, 'rb') as reader:
... data = bytearray(reader.read())
31.795 ms
>>> with Timer():
... with open(path, 'rb') as reader:
... data = bytearray(os.path.getsize(path))
... reader.readinto(data)
104857600
26.809 ms
The more common sin here is slicing. Whenever we take a slice of a list, or a bytes or bytearray, Python allocates a copy of the relevant subsequence. If we just want to look at it, though, we can use a memoryview instead — a simple view on the original memory, which can be indexed and sliced all the same, but without actually copying anything:
>>> with profile_memory():
... tail = data[len(data) // 2:]
51,200.6 KiB
>>> view = memoryview(data)
>>> with profile_memory():
... tail = view[len(view) // 2:]
0.6 KiB
>>> tail[0] == data[len(data) // 2]
True
>>> tail[0] = 0
>>> data[len(data) // 2]
0
So, 50 megabytes down to less than a kilobyte. And notice how, because it's the same memory, if we do end up modifying it, we actually modify the original — which is the whole point, but also something we have to be careful about.
This idea can be taken even further: the numpy library, for example, is entirely built around such smart references. In fact, an array is actually one contiguous buffer, accompanied by a shape and a series of strides: how many bytes we need to move in order to jump one item forward or back along any axis.
>>> import numpy as np
>>> a = np.arange(12).reshape(3, 4)
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> a.strides
(32, 8)
>>> a.T
array([[ 0, 4, 8],
[ 1, 5, 9],
[ 2, 6, 10],
[ 3, 7, 11]])
>>> a.T.strides
(8, 32)
>>> b = a[:2, 1:3]
>>> b
array([[1, 2],
[5, 6]])
>>> b.strides
(32, 8)
>>> b[0, 0] = 100
>>> a[0, 1]
np.int64(100)
>>> b.copy().strides
(16, 8)
For a 3 by 4 array of 8-byte integers, then, the next column is 8 bytes away, but the next row is 32. Transposing the array doesn't actually move anything around — it just swaps the strides; and a slice keeps the original strides, just with an offset into the original buffer, so that writing to it writes to the original all the same. The only way to get compact strides of our own is to actually copy the array, so:
>>> big = np.zeros((10_000, 10_000))
>>> big.nbytes
800000000
>>> with Timer():
... t = big.T
0.001 ms
>>> with Timer():
... t = big.T.copy()
369.895 ms
That's a microsecond for the transpose, compared to more than a third of a second for the actual copy. The point of a smart layout is that we can save on resources without having to write any C.
And that's really all this has been about: speaking Python well. There's a funny exercise: if you ask a group of Israelis if the Hebrew word for a junction, tsomet, is masculine or feminine, a lot of them will get it wrong (it's masculine). If you're writing a dictionary, or reading the news on TV, maybe it's important; but in a regular conversation, everyone will understand you just fine — and that's really what language is about. There's something charming about both the grammar pedant and the poet; but really, it's alright to write non-performant Python code when you don't care about performance, and it's good to know how to write performant Python code when you do.
Python versus C
Sometimes, though, we write our system as well as we can; and it's still not fast enough. This is when we start counting clock cycles and shaving off milli-, micro- and nanoseconds, and it's a very different game: more evolution than revolution, but there's still value in evolution — especially since it's all we have until the next revolution. If you had to choose between a horse and a car, you'd probably pick the car; and if you had to choose between the cars, you'd probably pick the Lamborghini over the Toyota.
This is where the criticism lands, so let's try to address it as precisely as we can. A language can't really be slow asymptotically: complexity is a property of the algorithm, not the language in which it is implemented — the latter is just used to describe it. In other words, Python isn't slower than C in orders of magnitude: it's slower in its constants. And that's actually an interesting tradeoff, because it can also go the other way: Python's flexibility lets us implement smarter algorithms, which would otherwise be impractical. Implementing a Fibonacci heap in C is decidedly less fun than doing it in Python, and for large enough $n$, the better algorithm will win anyway, whatever the constants — which is why Instagram can run on Python with billions of users, but I'd still write a real-time packet parser in C.
Still, we work in Python for everything else it gives us, and we'd like it to be fast, too. As an example, let's take the computation of a Mandelbrot fractal, courtesy of a talk by Matti Picus (one of PyPy's maintainers) at PyCon Israel. It's a cool image that keeps repeating itself when you zoom in on it; and it's computed using some neat complex-number arithmetic that you don't actually have to understand, as long as you follow the flow. Here's how it looks in Python:
>>> class Image:
... def __init__(self, width, height):
... self.width = width
... self.height = height
... self.data = bytearray(width * height)
>>> def mandelbrot(x, y, max_iters):
... c = z = complex(x, y)
... for n in range(max_iters):
... z = z*z + c
... if abs(z) >= 2:
... return n
... return max_iters
>>> def fractalize(image, max_iters, func):
... x_size = 3 / image.width
... y_size = 2 / image.height
... for y in range(image.height):
... imag = y * y_size - 1
... yy = y * image.width
... for x in range(image.width):
... real = x * x_size - 2
... image.data[yy + x] = func(real, imag, max_iters)
>>> image = Image(1200, 800)
>>> with Timer():
... fractalize(image, 20, mandelbrot)
729.644 ms
>>> reference = bytes(image.data)
An image is a width, a height, and a bytearray of width times height bytes. mandelbrot takes a point and returns after how many iterations it escapes, and fractalize calls it on every pixel of the 960,000 pixels that comprise our image. It accepts the function as a parameter, so that we can swap it with different implementations later on; but for now, this is it, and it takes some 730 milliseconds to compute.
Let's do the same in C, this time in the most Sisyphean version we can write, as something to aspire to. First, a header:
// fractalize.h
typedef struct Image {
int width;
int height;
unsigned char* data;
} Image;
unsigned char mandelbrot(double x, double y, unsigned char max_iters);
void fractalize(Image image, unsigned char max_iters);
And the implementation itself:
// fractalize.c
#include <complex.h>
#include "fractalize.h"
unsigned char mandelbrot(double x, double y, unsigned char max_iters) {
double complex c = CMPLX(x, y);
double complex z = c;
for (unsigned char n = 0; n < max_iters; n++) {
z = z*z + c;
if (creal(z)*creal(z) + cimag(z)*cimag(z) >= 4) {
return n;
}
}
return max_iters;
}
void fractalize(Image image, unsigned char max_iters) {
double x_size = 3.0 / image.width;
double y_size = 2.0 / image.height;
for (int y = 0; y < image.height; y++) {
double imag = y * y_size - 1;
int yy = y * image.width;
for (int x = 0; x < image.width; x++) {
double real = x * x_size - 2;
image.data[yy + x] = mandelbrot(real, imag, max_iters);
}
}
}
Notice how much less pleasant it is, from complex.h and its CMPLX to the explicit types of every variable and the cumbersome loops. And that's even without objects — this is an easy case, just some arithmetic. As an aside, when benchmarking, it's only fair to make sure all implementations compute the same thing — and indeed, an earlier version of this example would escape at abs(z) >= 4 in Python, but at a squared magnitude of 4 (that is, abs(z) >= 2) in C, and would use floats in C but doubles in Python; here, I've made sure to use doubles everywhere, and that the escape radius is the same for every version — in fact, every version's output image is compared to the Python one, byte for byte, just to make sure they're exactly the same.
Next, we need a small shared library and a demo.c that allocates the image, times the call to fractalize with gettimeofday, prints out the milliseconds and writes the raw image to a file — as well as a makefile to compile them together:
# makefile
demo: demo.c fractalize.h libfractalize.so
gcc demo.c -o demo -L. -Wl,-rpath,. -lfractalize
libfractalize.so: fractalize.c fractalize.h
gcc -shared -fPIC fractalize.c -o libfractalize.so
clean:
@rm -f demo libfractalize.so image.raw
$ make
gcc -shared -fPIC fractalize.c -o libfractalize.so
gcc demo.c -o demo -L. -Wl,-rpath,. -lfractalize
$ ./demo
96.028 ms
So it's indeed faster: some 96 milliseconds, as opposed to our 730 — that's almost eight times as fast as our Python version. Keep these numbers in mind, though, and the fact that the makefile doesn't pass any optimization flags; we'll be back to it later.
Foreign Function Interfaces
A foreign function interface, or FFI, is a mechanism for calling code written in one language from another language, by actually running it as part of the same process. This is not to be confused with serializing data over the wire using some standard format, like HTTP or protobufs. Protobufs do generate code for each language, which is clever; but it's not FFI, because it doesn't run as part of the same process.
The gist of a foreign function interface is simple: you load a binary, translate the arguments to the right types, jump to the function's code and execute it, then translate its return value back. Of course, there are some complications, like garbage collection (the other language might collect an object that you still have a reference to), or runtime environments that don't have a counterpart on the other side. Luckily, CPython is written in C, so we get a C foreign function interface as part of the standard library: ctypes. Other implementations might come with their own FFIs: Jython, for example, can talk to Java, while IronPython can communicate with .NET — Python, after all, is a standard, like SQL, with multiple implementations.
Anyway, with ctypes, we can describe the C struct as a Structure subclass with _fields_, load the shared library, and declare each function's argument and return types — although, if we don't specify argtypes and restype, it will assume everything is an int:
>>> import ctypes
>>> class CImage(ctypes.Structure):
... _fields_ = [
... ('width', ctypes.c_int),
... ('height', ctypes.c_int),
... ('data', ctypes.POINTER(ctypes.c_ubyte)),
... ]
... def __init__(self, width, height):
... self.width = width
... self.height = height
... self.buffer = (ctypes.c_ubyte * (width * height))()
... self.data = self.buffer
>>> c_lib = ctypes.CDLL('./libfractalize.so')
>>> c_fractalize = c_lib.fractalize
>>> c_fractalize.argtypes = [CImage, ctypes.c_ubyte]
>>> c_fractalize.restype = None
>>> c_image = CImage(1200, 800)
>>> with Timer():
... c_fractalize(c_image, 20)
93.708 ms
>>> bytes(c_image.buffer) == bytes(image.data)
True
Let's go over it, because it's a bit confusing: (ctypes.c_ubyte * n)() is how we create an array type and instantiate it; we also keep a reference to it in the object (self.buffer), so that it doesn't get garbage collected while C is still pointing to it. And the timing? It's about 93 milliseconds — the same as in C, with just a sliver of overhead to manage the call, but also exactly the same image, byte for byte.
However, there's a catch: to do it, we had to implement everything in C and then figure out how to wrap it as well. What if we kept our nice, simple loop in Python, but used mandelbrot in C to do the heavy lifting?
>>> c_mandelbrot = c_lib.mandelbrot
>>> c_mandelbrot.argtypes = [ctypes.c_double, ctypes.c_double, ctypes.c_ubyte]
>>> c_mandelbrot.restype = ctypes.c_ubyte
>>> c_mandelbrot(0, 0, 20), c_mandelbrot(1, 1, 20)
(20, 0)
>>> with Timer():
... fractalize(image, 20, c_mandelbrot)
706.701 ms
This is… underwhelming. It takes about 700 milliseconds, which is barely better than the pure Python implementation; and no wonder, since every one of the 960,000 calls involves crossing the border: translating the arguments, jumping through the generic call machinery to get the result, and translating it back — which is just as expensive as the arithmetic we were doing. It makes sense, then, that FFIs become useful when we don't cross them so often, but when we do, we actually need to get something done.
That's what CFFI is for: it does most of the wrapping for us, by parsing the C declarations itself. In other words, we can use the same header file we already wrote:
>>> import cffi
>>> ffi = cffi.FFI()
>>> with open('fractalize.h') as reader:
... ffi.cdef(reader.read())
>>> ffi_lib = ffi.dlopen('./libfractalize.so')
>>> ffi_lib.fractalize
<cdata 'void(*)(Image, unsigned char)' 0x104c3c538>
>>> ffi_image = ffi.new('Image *')
>>> ffi_image.width, ffi_image.height = 1200, 800
>>> ffi_data = ffi.new('unsigned char[]', 1200 * 800)
>>> ffi_image.data = ffi_data
>>> with Timer():
... ffi_lib.fractalize(ffi_image[0], 20)
94.232 ms
>>> bytes(ffi.buffer(ffi_data)) == bytes(image.data)
True
Here, we use ffi.new to allocate memory that's owned by the returned object, which is then freed when it's garbage collected — so we need to keep a reference to ffi_data in a variable, rather than assigning a fresh allocation straight into the struct, which would leave it pointing to memory that was already freed. And once again, the timing is the same as ctypes, but with much less ceremony. This is CFFI's ABI mode — it also provides an API mode, which compiles a proper extension module, and can be even faster.
>>> with Timer():
... fractalize(image, 20, ffi_lib.mandelbrot)
320.013 ms
It's about 320 milliseconds — more than twice as fast as our ctypes version, because CFFI's calling path is considerably more efficient. So if you have a shared library that you want to use from Python, CFFI is usually a better choice.
Cython and numba
And what if we didn't have to write C at all? Cython is a whole new programming language, sort of: it's basically Python, but with static type declarations (and a few other features) that are used to compile it to highly optimized C, which also generates the FFI glue for us, so its interface remains Python. It can also be used to wrap existing C libraries — but that's something we've already seen.
This is where our tangent about the array returns: types are absolutely key to performance. I've heard people argue for statically typed languages on the basis that they are safer — but to me, it's always been the less interesting argument, for all the reasons we covered at length in the type hints chapter; the interesting point is that they are more efficient, because the compiler can use more efficient instructions to deal with them, and store them in a more compact format.
To illustrate it, consider this code, in cy_fractalize.pyx:
# cy_fractalize.pyx
cpdef unsigned char cy_mandelbrot(double x, double y, unsigned char max_iters):
cdef unsigned char n
cdef double complex c, z
c = z = complex(x, y)
for n in range(max_iters):
z = z*z + c
if z.real*z.real + z.imag*z.imag >= 4:
return n
return max_iters
cpdef cy_fractalize(image, unsigned char max_iters):
cdef int x, y, yy
cdef int width = image.width, height = image.height
cdef double x_size = 3.0 / width, y_size = 2.0 / height, imag, real
cdef unsigned char[:] data = image.data
for y in range(height):
imag = y * y_size - 1
yy = y * width
for x in range(width):
real = x * x_size - 2
data[yy + x] = cy_mandelbrot(real, imag, max_iters)
It might look like a strange, alien dialect — but it's actually quite clever. cdef is used to define variables at the C level, which are typed accordingly (double complex and unsigned char), and there are three types of functions: def is a regular Python function, cdef is a pure C one (so it's fast, but not callable from Python), and cpdef is a combination of both — a C definition, but with a Python-callable wrapper. unsigned char[:] is a typed memoryview — a slice of memory of a particular type, which in this case is used to efficiently access our image object's memory, no copy. The main overhead (besides the weird syntax) is that we have to annotate everything — but we don't actually have to do that, just the bits we want to optimize; the rest of the code will remain slower. After a few days with the tutorial, it becomes pretty natural (if still annoying). The quickest way to try it out is pyximport: it compiles .pyx files when they're imported, but in a real project you'd probably want to use cythonize as part of your packaging:
>>> import pyximport
>>> pyximport.install(language_level=3)
(None, <pyximport.pyximport.PyxImportMetaFinder object at 0x105768440>)
>>> from cy_fractalize import cy_fractalize, cy_mandelbrot
>>> image = Image(1200, 800)
>>> with Timer():
... cy_fractalize(image, 20)
93.116 ms
>>> bytes(image.data) == reference
True
>>> with Timer():
... fractalize(image, 20, cy_mandelbrot)
178.932 ms
Which is, indeed, about 93 milliseconds — as fast as our C; and the hybrid version is about 180, because calling a cpdef function from Python is cheaper than going through ctypes or CFFI.
Numba, on the other hand, requires even less work: instead of us having to declare types, it attempts to infer them, and does an amazing job at that. All we need to do is decorate our function with @numba.jit, and the first time we call it, it will look at the arguments' types and compile the function's bytecode into machine code (using LLVM); subsequent calls will reuse the same binary, so this first call will be slower — but we don't end up with an entire new language.
And what's more, this optimization is polymorphic: if we call it with ints, it will generate a version of the function specialized for int — and if we call it with floats, it will compile another version, dispatching based on the argument types. It's an absolutely brilliant solution, and the interface is so elegant it almost feels like magic.
First, here's what I had in my old notebooks: it jitted fractalize and applied it to our Image object:
>>> import numba
>>> @numba.jit
... def numba_mandelbrot(x, y, max_iters):
... c = z = complex(x, y)
... for n in range(max_iters):
... z = z*z + c
... if abs(z) >= 2:
... return n
... return max_iters
>>> @numba.jit
... def numba_fractalize(image, max_iters):
... x_size = 3 / image.width
... y_size = 2 / image.height
... for y in range(image.height):
... imag = y * y_size - 1
... yy = y * image.width
... for x in range(image.width):
... real = x * x_size - 2
... image.data[yy + x] = numba_mandelbrot(real, imag, max_iters)
>>> numba_fractalize(Image(1200, 800), 20)
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
...
- argument 0: Cannot determine Numba type of <class '__main__.Image'>
And it used to work... kind of. Numba couldn't really type our Python object, so it silently gave up, and reverted to object mode, in which it's basically the same as Python, only with some extra steps involved (and a nice warning message that nobody read). Since version 0.59, it defaults to nopython mode, which refuses to work with such an object; but it's actually quite easy to fix. All we need to do is give it some data that it can understand: the width, the height, and a numpy array over the image's bytearray, which np.frombuffer creates as a view (not a copy); and since a bytearray is writable, so is the array, and numba happily writes the pixels straight into our image:
>>> import numpy as np
>>> @numba.jit
... def numba_fractalize(data, width, height, max_iters):
... x_size = 3 / width
... y_size = 2 / height
... for y in range(height):
... imag = y * y_size - 1
... yy = y * width
... for x in range(width):
... real = x * x_size - 2
... data[yy + x] = numba_mandelbrot(real, imag, max_iters)
>>> image = Image(1200, 800)
>>> data = np.frombuffer(image.data, dtype=np.uint8)
>>> with Timer():
... numba_fractalize(data, image.width, image.height, 20)
253.503 ms
>>> with Timer():
... numba_fractalize(data, image.width, image.height, 20)
31.458 ms
>>> bytes(image.data) == reference
True
>>> with Timer():
... fractalize(image, 20, numba_mandelbrot)
197.684 ms
So, it takes around 250 milliseconds the first time, mostly to compile, and then some 31 milliseconds: that's three times faster than our C code. Not because numba is smarter than C, mind you, but because it actually compiled with optimizations on — something that, as we shall see, our C code didn't. And the hybrid version, by the way, is around 200 milliseconds.
So next time something mathematical is taking too long — and as long as it's not a database query, which no amount of clever optimizations will help — try slapping a @numba.jit decoration on top and see what happens; it's really cool.
PyPy and RPython
There's one more technology I'd like to discuss, and that's PyPy. It's one of the most interesting software projects I've ever heard about, and well worth knowing about even if you have no intention of using it; just appreciate it for the mind-bending, cool idea that it is.
A group of programming-language researchers decided to write the Python interpreter in Python. The first, obvious question is, why? — and the second, obvious question is, how? To address the first question, they wanted to see how much easier it is to write an interpreter in Python than in C; to address the second question, it actually isn't that impossible, but we'll get there. What I found amazing about the project, besides its audacity, is that it actually is better — a lot better, in fact.
At first glance, it doesn't really make sense: in order to run an interpreter written in Python, you'd have to have some other Python to run it on — so it can't be as fast as that other Python, by definition. And at first, that's exactly what happened.
But at the same time, there's something profoundly cool in the idea of writing an interpreter in Python. Like we saw before, an interpreter includes a lexer, to break the code into tokens, a parser, to group the tokens in a syntax tree, a compiler, to translate the syntax tree into some bytecode, and a virtual machine, to execute that bytecode by performing the corresponding operations, one opcode at a time (we saw these steps in action in our chapters on code objects and modules, and saw the resulting .pyc files in __pycache__):
>>> import dis
>>> def f():
... x = 1
... y = 2
... return x + y
>>> dis.dis(f)
1 RESUME 0
2 LOAD_CONST 1 (1)
STORE_FAST 0 (x)
3 LOAD_CONST 2 (2)
STORE_FAST 1 (y)
4 LOAD_FAST_LOAD_FAST 1 (x, y)
BINARY_OP 0 (+)
RETURN_VALUE
>>> len(dis.opmap)
150
However high-level our Python code is, it eventually boils down to a fairly small set of opcodes: in CPython 3.13, dis.opmap contains about 150 different ones. If you implement all of them, you basically implemented Python — anything else, like metaclasses and generators, are really just code that runs on top of this opcode "kernel".
And if you do implement it in Python, it suddenly becomes so much easier to experiment with: you can just edit a function and see if your idea works. In fact, it becomes so much easier to experiment that they found a lot of ideas that work a lot better than what CPython does. One such idea was the compact dictionary layout that Raymond Hettinger proposed, which is more memory-efficient, and happens to be ordered as a side effect — and in fact, it was first implemented in PyPy, and only then included in CPython (in version 3.6), so that's the real reason dictionaries are ordered nowadays. It just goes to show how a higher-level language can drastically reduce the gulf of execution, that mental chasm between an idea and its implementation.
Anyway, the reason this is interesting, performance-wise, is that the way they actually did it is by implementing a translator. What it does is take Python code and translate it into C; the resulting C is then compiled into a binary that can be executed. In other words, they fed their Python interpreter into their Python translator, which translated it into C and compiled it into a binary.
This, in turn, presented them with a problem: in some cases, it was impossible to translate the code into C; so they had to restrict the language in which they wrote the Python interpreter so that it is possible to translate. There's one small lie in the previous paragraph, then: PyPy isn't actually written in Python, but in RPython, which is a restricted subset of Python that they can analyze and translate into C. So when you install the pypy binary, it's not actually Python running on Python: it's an interpreter written in RPython and translated into a binary, pretty much like CPython is written in C and compiled into a binary. On top of that, the translator is able to generate a just-in-time compiler for the interpreter, which can further optimize and compile hot loops into machine code on the fly.
It's a lot of big words, but using it is actually quite simple: if we go back to our fractal script and run it on PyPy, we get:
$ time python fractalize.py
real 0m0.773s
user 0m0.745s
sys 0m0.025s
$ time pypy3 fractalize.py
real 0m0.175s
user 0m0.161s
sys 0m0.012s
That's the exact same script, mind you; not a single line is different, and the only thing I changed was the command used to run it. It went from about 770 milliseconds to just 175 — and that's PyPy in action. There's nothing else to show, really, but I was curious to see what RPython is about. It's quite poorly documented, so it took some trial and error: no bytearray, just a list of ints; no complex numbers, so I had to track the real and imaginary parts separately as two floats; a main function that accepts an argv and returns an exit code, and a target function that returns it; and writing out the file by joining all the characters together:
# rfractalize.py
import time
class Image:
def __init__(self, width, height):
self.width = width
self.height = height
self.data = [0] * (width * height)
def mandelbrot(x, y, max_iters):
a, b = x, y
for n in range(max_iters):
a, b = a*a - b*b + x, 2*a*b + y
if a*a + b*b >= 4:
return n
return max_iters
def fractalize(image, max_iters):
x_size = 3.0 / image.width
y_size = 2.0 / image.height
for y in range(image.height):
imag = y * y_size - 1
yy = y * image.width
for x in range(image.width):
real = x * x_size - 2
image.data[yy + x] = mandelbrot(real, imag, max_iters)
def main(argv):
image = Image(1200, 800)
started = time.time()
fractalize(image, 20)
elapsed = time.time() - started
print('%f ms' % (elapsed * 1000))
with open('image.raw', 'wb') as writer:
writer.write(''.join([chr(b) for b in image.data]))
return 0
def target(*args):
return main, None
So, translating it with RPython (which took about half a minute, and printed some colourful ASCII art along the way — they're a rather eccentric bunch) and running the resulting binary:
$ pypy2 ../pypy/rpython/bin/rpython rfractalize.py 2>&1 | tail -4
[Timer] source_c --- 0.9 s
[Timer] compile_c --- 0.8 s
[Timer] =========================================
[Timer] Total: --- 27.0 s
$ ./rfractalize-c
19.242048 ms
$ cmp image.raw ../c/image.raw && echo same image
same image
That is, it runs in 16 to 30 milliseconds, whereas my hand-written C takes about 95, which I'd like to believe is decent enough C; when I first saw these results, I was sure I had some bug, because there's no way — but the image is identical, byte for byte, so there's no room for a bug. So for years, I taught it exactly like that: thirty lines of restricted Python, translated automatically, beat hand-written C — and I'd let my students marvel at how mind-blowing it is.
To be fair, you should look at what it generates: some 40,000 lines of C, and nobody in their right mind would call it readable — it's full of gotos and blocks. And if you look at our mandelbrot function, you'll see it became a loop of float macros, with the complex math spelled out in explicit terms; and if you look closer, you'll see that max_iters is replaced with the constant 20, because the translator was smart enough to analyze the entire program and realize this is the only value ever passed to it, so it might as well be inlined:
// implement.c (generated by RPython, excerpt)
Signed pypy_g_mandelbrot(double l_x_1, double l_y_1, Signed l_max_iters_2) {
...
/* block0: (inlined) */
l_x_3 = l_x_1;
l_y_3 = l_y_1;
l_v460 = 1L;
l_v461 = 20L;
l_next_3 = 0L;
goto block1;
block1:
while (1) {
OP_FLOAT_MUL(l_x_3, l_x_3, l_v462);
OP_FLOAT_MUL(l_y_3, l_y_3, l_v463);
OP_FLOAT_SUB(l_v462, l_v463, l_v464);
OP_FLOAT_ADD(l_v464, l_x_1, l_x_2);
OP_FLOAT_MUL(2.0, l_x_3, l_v465);
OP_FLOAT_MUL(l_v465, l_y_3, l_v466);
OP_FLOAT_ADD(l_v466, l_y_1, l_y_2);
OP_FLOAT_MUL(l_x_2, l_x_2, l_v467);
OP_FLOAT_MUL(l_y_2, l_y_2, l_v468);
OP_FLOAT_ADD(l_v467, l_v468, l_v469);
OP_FLOAT_GE(l_v469, 4.0, l_v470);
if (l_v470) break;
...
}
...
}
Still, it's quite something. Then one day, I was looking at the makefile RPython generates — and noticed it was compiling the C with the -O3 flag, while mine was using no optimizations at all, which means gcc's default behavior of -O0. Years of lectures, and nobody noticed this — least of all myself. So I did have a bug after all: in my benchmark. If I compare the same C code, compiled with -O2:
$ ./demo
95.408 ms
$ gcc -O2 demo.c fractalize.c -o demo-O2
$ ./demo-O2
18.262 ms
We get 18 milliseconds, which is pretty much neck and neck with RPython. It's still pretty astonishing: that thirty lines of restricted Python, translated automatically into C, can rival hand-written and optimized C, but it's also a lesson in miniature of what we've been talking about this whole time: a benchmark is only as fair as its setup. Here's the full table of results from python shootout.py:
$ python shootout.py 2>/dev/null
Method Time
Python 722 ms
C (gcc -O0) 92 ms
C (gcc -O2) 18 ms
ctypes 92 ms
Python + ctypes 698 ms
CFFI 94 ms
Python + CFFI 313 ms
Cython 93 ms
Python + Cython 174 ms
numba (first call 536 ms) 31 ms
Python + numba 192 ms
PyPy 127 ms
RPython 16 ms
And it's quite interesting: first, CFFI is as fast as ctypes, and much, much simpler to use; second, Cython, numba and PyPy are all close to C, without actually writing any C at all; third, numba is as fast as Cython (in this particular case, even faster), and much, much simpler to use; fourth, PyPy is absolutely awesome; and fifth, RPython is ridiculously fast — but also nearly impossible to work with, since it's not actually Python, but an undocumented subset, intended as an implementation detail of PyPy, so they reserve the right to change it whenever they please.
Why Doesn't Everyone Use PyPy?
So if all it takes is changing the command, why isn't everyone using PyPy?
First, it's weird: a lot of people just don't know about it, and some who do don't trust it. Some programmers are wary of anything that smells too much of magic; but if you're a wizard, it's not really fair to refuse casting fireballs just because they're dangerous — and it's hard to deny that plenty of people will tell you that voodoo is bad, and that they only want to engage in boring, straightforward code.
Second, there's inertia: a really good reason. It's a big project, it works, and it's in CPython; why take the risk and do something? You know, cars were invented, too; but you've already got a horse, and a stable, and lots of hay. So you can understand why some people prefer to keep riding. Similarly, a lot of companies have huge, complex monoliths that work; and while they might want to switch, it's just not worth the hassle.
Third, my example is biased. It's a tight mathematical loop — the kind of thing that a just-in-time compiler loves to sink its teeth into. But if we were running a webserver, most of its time would be spent in recv and send — waiting, really; and in that case, whether the interpreter is faster or not barely matters. And if we're talking about small scripts, they might not even run long enough for the JIT to warm up.
Fourth, a lot of numeric code is already in numpy, anyway — which is C underneath, remember?
$ python squares.py
661 ms
$ pypy3 squares.py
20 ms
$ python squares_numpy.py
15 ms
So in CPython it takes 660 milliseconds; in PyPy it takes 20; and in numpy, which runs on CPython, it takes 15. If you're already vectorized, there's little left for PyPy to offer; and anyway, numpy, like any other library that's built as a CPython C extension, is available in PyPy only through a compatibility layer that's slower and less complete. If even one of the libraries your project depends on is such a case, switching becomes painful. Also, PyPy tends to trail CPython's versions (at the time of writing, it supports Python 3.11).
Finally, most people choose Python for its flexibility, not its performance; and rewriting things just to go faster might not be the right thing to do in their situation. For those who do care about performance, though, PyPy is often an easy win.
On a practical level, then, it's a neat trick: PyPy is literally better than you (or at least, me). On a philosophical level, though, the question I'm interested in is why do we have beliefs to begin with.
Tim Urban, the author of the Wait But Why blog, wrote about this in his series The Story of Us. He mentions brain-imaging research, which showed that when people were presented with evidence that contradicted their non-political beliefs, their brain's decision-making machinery was engaged, as you would expect. But when evidence contradicted their political beliefs, their default mode network was engaged, too — the part of the brain associated with our sense of self. In other words, some opinions we have are about something; but some opinions we have are about us. We reason about the former; but we identify with the latter.
I suspect that for most of you, before you started reading this book, you thought Python was slow. But why did you think that? Did you measure it? Or did you hear it from someone? My guess is that most people never bothered to measure, and just repeated what they heard, until it became part of who they are as programmers.
Tim's suggestion is to invent new, neutral terminology; solve the problem in this abstract realm, where we have no preconceptions, no anchors or connotations; and then project the conclusions back onto reality. It's how math works, after all: you start university and think that 1 is an apple, so -1 must be an apple you owe; and then along come complex numbers, so what the fuck is i apples? Nothing, really — it's an abstraction, but the things that happen in this abstraction have very real consequences.
I don't know why, but whenever I had to solve a performance problem, my first intuition was always wrong; it was always better the moment I measured, or the moment I thought about the problem in terms that I didn't have an opinion about yet.
Let me finish with an army story: I wrote an important system there, and it started crawling, so I was called up for reserve duty to fix it. The team tried running it on PyPy ("like in your lectures"), but it didn't work out of the box, because some libraries were C extensions; so they defined my task as "make it work on PyPy". It got pretty complicated, and at some point I asked why; and they said because it's slow. But what's slow? Did you benchmark it? They didn't; PyPy was kind of the no-brainer, the silver bullet, just change the command and everything will be faster.
So I spent an hour profiling the system, and discovered that out of every ten seconds, it spent nine and a half waiting for one database query: the tables were missing an index. With the right index, the system ran about 30 times faster; and if we had managed to run it on PyPy, it would execute the Python part a bit faster, and then wait for the same database anyway. So first, you benchmark; then, you don't optimize what's not the bottleneck.
Anyway, that system was slow not because of computation, but because of waiting; and it turns out that most programs spend most of their time waiting — for databases, for disks, but mostly for networks. To understand this waiting, let's take a look at how the internet actually works.