Yet Faster Fixed Point Mandelbrot Set

Interpreted languages such as python are slow. One thing which makes python particularly slow is its dynamic typing. Every time it meets a line such as x=x+y it has to determine what type of variables x and y are. Is the addition required integer addition, floating point addition, or string concatenation? Or is it an error condition because we are trying to add apples and pears, or some other incompatible types? Perhaps, given our example code, it is an GMP MPZ addition. Python's typing is so dynamic that the types may have changed from the last time that the statement was executed. The result is the the interpreter may take longer to determine what sort of addition is required than it takes to perform the addition!

A traditional compiled language such as C does not have this issue. So we will try rewritting the core loop of our code in C, and then attempt to use it from python. We cannot use numba in this case, for numba does not understand python's inbuilt arbitrary precision integers, and insists that they get converted to fixed-length C integers, and numba also does not understand the MPZ objects fromthe gmpy2 module.

Using GMP with C is not very pleasant. C does not permit the usual operators to be overloaded, so one is left calling functions to perform every operation, and the result looks rather too much like assembler.

#include <gmp.h>
  
int mandel (mpz_t x0, mpz_t y0, long bshift, int maxiter){
  mpz_t x,y,x2,y2,mod_test,tmp;
  int i;

  mpz_inits(x,y,x2,y2,mod_test,tmp,NULL);

  mpz_set(x0,MPZ(x0_py));
  mpz_set(y0,MPZ(y0_py));

  /* mod_test = 4<<(2*bshift) */
  mpz_set_ui(mod_test,(unsigned long)4);
  mpz_mul_2exp(mod_test,mod_test,2*bshift);
  
  mpz_set(x,x0);
  mpz_set(y,y0);
  mpz_mul(x2,x,x);
  mpz_mul(y2,y,y);

  for(i=0;i<maxiter;i++){
    /* tmp=x2+y2 */
    mpz_add(tmp,x2,y2);
    if (mpz_cmp(tmp,mod_test)>=0) break;
    /* y=((x*y)>>(bshift-1))+y0 */
    mpz_mul(y,x,y);
    mpz_tdiv_q_2exp(y,y,bshift-1);
    mpz_add(y,y,y0);
    /* x=((x2-y2)>>bshift)+x0 */
    mpz_sub(x,x2,y2);
    mpz_tdiv_q_2exp(x,x,bshift);
    mpz_add(x,x,x0);
    /* x2=x*x */
    mpz_mul(x2,x,x);
    /* y2=y*y */
    mpz_mul(y2,y,y);
  }
  
  mpz_clears(x,y,x2,y2,mod_test,tmp,NULL);
  return i;
}

But trying to interface this to python creates another layer of unpleasantness. Firstly python will pass its arguments as a tuple which needs unpacking, and will expect the return value to be a python object.

#define PY_SSIZE_T_CLEAN
#include <Python.h>
#include <gmp.h>
  
static PyObject *mandel(PyObject *self, PyObject *args){
  int maxiter;
  long bshift;
  PyObject *x0_py,*y0_py;
  mpz_t x,y,x0,y0,x2,y2,mod_test,tmp;
  int i;
  
  if (!PyArg_ParseTuple(args,"OOli",&x0_py,&y0_py,&bshift,&maxiter)){
    fprintf(stderr,"Argument error\n");
    exit(1);
  }

  mpz_inits(x,y,x0,y0,x2,y2,mod_test,tmp,NULL);
  mpz_set(x0,MPZ(x0_py));
  mpz_set(y0,MPZ(y0_py));

[...]

  return PyLong_FromLong((long)i);
}

And then there is the matter of adding an initialisation function too.

Finally one needs to compile it. The recommended route is via distutils.

A recipe

(This recipe confines iteself to a single directory, and does not clutter your collection of installed python modules. It is thus useful whilst developing. It needs the GMP development package to be installed, and some versions of Debian seem not to realise that that depends on two more Gnu arbitrary precision pacakges. One will probably need

  sudo apt-get install libgmp-dev libmpfr-dev libmpc-dev

as well as a standard C compilation environment.)

Place in the same directory

• setup.py
• pure_c.c
• mandel_tk_gmpc_par.py

Then build the module.

  $ python3 setup.py build_ext --inplace

which will place the module (shared library) in the current directory. It will have a very long name, something like
gmp_mandel.cpython-39-aarch64-linux-gnu.so
Then one should be able to run it.

  $ python3 mandel_tk_gmpc_par.py

Was it worth it? Continuing the timings presented earlier, the default set took 5.2s to display, and with the bit shift set to 160 it took just under 8s. It is still a lot slower than the standard double precision version using numba, which takes under 0.5s, but it is about 3.5x as fast as our previous best attempt at writing our own fixed point code, and it has shown how to write a C function which is callable directly from python without using cython.

Improvements?

One might wonder if trying to convert to C the whole of the loop generating a column would produce a further improvement. But what we have done has converted to C the majority of the code in terms of where the execution time is spent, and, once the number of iterations per pixel starts to rise, as it will at high zoom levels, the amount of work outside of the inner loop for a pixel becomes negligible.

But as an exercise to learn a few more techniques, perhaps moving the whole of the mandel_column function to C is worthwhile. To make it more interesting, we could pass the C function standard python integers, and have the C part convert them to the GMP MPZ type, leaving the python completely free of GMP and hence of the gmpy2 module too.

• setup.py
• pure_c2.c
• mandel_tk_gmpc_par_col.py

Done. It is about 0.7s faster on the two tests above, and the time saved is expected to depend on the number of pixels, but not the value of maxiter or the bit shift. The function PyLong2GMP shows how python stores integers. A count of "digits", with negative counts indicating a negative number, and the number of "digits" being the absolute value of the count. Each "digit" contains the next PyLong_SHIFT bits of the number stored. On 64-bit platforms a "digit" is a 32-bit unsigned integer, and PyLong_SHIFT is 30.

The code also shows how to return a newly-created Numpy array from C to python.