CHAPTER
CYTHON WITH NUMPY ARRAYS
Below we have a similar Cython file, the original version for this approach was subbmited by Didrik Pinte of Enthought (thanks Didrik!). The main difference is the annotation of numpy arrays, see the tutorial for a great walkthrough: http://docs.cython.org/src/tutorial/numpy.html (and there’s a bit more detail in the wiki: http://wiki.cython.org/tutorials/numpy).
Using the numpy approach Python is able to address the underlying C data structures that are wrapped by numpy without the Python call overheads. This version of the Mandelbrot solver runs almost at the same speed as the Shed Skin solution (shown in the next section), making it the second fastest single-CPU implementation in this tutorial.
IAN_TODO I ought to remove Didrik’s local declaration of z = 0+0j to make it a fairer comparision with the rest of the code (though my gut says that this will have little effect on the runtime)
# calculate_z.pyx
# see ./cython_numpy_loop/cython_numpy_loop.py
from numpy import empty, zeros
cimport numpy as np
def calculate_z(np.ndarray[double, ndim=1] xs, np.ndarray[double, ndim=1] ys, int maxiter):
""" Generate a mandelbrot set """
cdef unsigned int i,j
cdef unsigned int N = len(xs)
cdef unsigned int M = len(ys)
cdef double complex q
cdef double complex z
cdef int iteration
cdef np.ndarray[int, ndim=2] d = empty(dtype=’i’, shape=(M, N))
for j in range(M):
for i in range(N):
# create q without intermediate object (faster)
q = xs[i] + ys[j] * 1j
z = 0+0j
for iteration in range(maxiter):
z = z * z + q
if z.real * z.real + z.imag* z.imag > 4.0:
break
else:
iteration = 0
d[j,i] = iteration
return d