The Mandelbrot Project

I’ve published the source code of the program I wrote for my tech talk at the 2011 PyDay conference. It’s a Python script and a companion C library that calculates and draws the Mandelbrot set.

The objective of the tech talk was to show how to speed up Python programs using the power of native code.


A render of the Mandelbrot set as performed by the mandelbrot.py script. Computations were performed in C.

A render of the Mandelbrot set as performed by the mandelbrot.py script. Computations were performed in C.


What’s interesting about this program is that, although the core was written completely in Python, I wrote two compute backends for it: one in Python and one in C. The C code is interfaced with using the ctypes module.

The results of running the program are shown in the screenshot above. If you are interested in trying it, the full source code is hosted in GitHub, here: https://github.com/alesegovia/mandelbrot. I’ve licensed it under the GPLv3, so you can download, run it, test it and modify it.

As one would anticipate, the C implementation runs much faster than the Python one, even when taking into account the marshaling of objects from Python to C and back. Here’s the chart I prepared for the conference showing the specific numbers from my tests.

These tests were performed to compare the run times at different numer of iterations, note this is a logarithmic scale.


Comparison of the Python + C implementation vs a pure Python one. Scale is Logarithmic.

Comparison of the Python + C implementation vs a pure Python one. Scale is Logarithmic.


As you can see, Python programs can be significantly sped up using ctypes, especially when we are dealing with compute-intensive operations.

It might be possible to speed up the Python implementation to improve its performance to some extent, and now that the source code is available under the GPL, you are encouraged to! I would always expect well-written C code to outperform the Python implementation, but I would like to learn about your results if you happen to give it a go.

Happy hacking!

Implementing sqrt

I was reading some post about interview questions of 2011 and came across one that stated “find the square root of a number”.

Assuming we can’t use the sqrtf function of the standard C math library, let’s see how we can calculate the square root of a number x.

Given n, we know that its square root is a number x that holds:

\sqrt{n} = x

Let’s work on this equation a little. Raise both sides to the second power:

n = x^2

Move to the left of the equality:

0 = x^2 - n

If we found the roots of this last equation somehow, we would have found the square root of n. We can do this by using the Newton-Raphson iteration.

The Newton-Raphson iteration states that we can find the root of an equation using the following formula iteratively:

x_{n+1} = x_n - \frac{f(x)}{f'(x)}

Where f'(x) is the derivative of function f(x). We will approximate the derivative using the definition of derivative at a point (we could also note that the derivative could be trivially calculated; this method is more general).

f'(x) = \frac{f(x+h) - f(x)}{h}

The error of the Newton-Raphson iteration is given by:

|x_{n+1} - x_n|

Starting with a hardcoded seed value, we will perform this iteration in a loop until the error is less than a given value. I have chosen to iterate until the error is less than 1×10^(-10): 0.00000000001.

Let us see what a tentative “pythonesque” pseudocode for this loop could be:

def sqrt(n):
    f = function(x*x - n)
    x = 1 # seed
    xant = 0
    do:
        f1 = (f.eval(x+h) - f.eval(x)) / h
        xant = x
        x = x - f.eval(x)/f1
    while abs(x - xant) > err;
    return x

Assuming we have a symbolic function type, that loop does not seem too difficult. In order to code this in C, since the equation is always the same, I will hardcode it as a plain function.

typedef double real; // change to float for single precision

real f(real x, real n)
{
    return x*x - n;
}

real sqrt(real n)
{
    real err = 0.00000000001f;
    real h = 0.01f;

    real x = 1.0f; // seed
    real xant = 0.0f;

    do
    {
        xant = x;
        real df = (f(x+h, n) - f(x, n))/h;
        x = x - f(x, n)/df;
    }
    while (abs(x - xant) > err);

    return x;
}

Here are the results of running our custom square root function, compared to the standard version provided with the C programming language:

[ale@syaoran sqrt]$ ./sqrt 1.0
Custom sqrt of: 1 = 1
libm sqrt of: 1 = 1

[ale@syaoran sqrt]$ ./sqrt 2.0
Custom sqrt of: 2 = 1.41421
libm sqrt of: 2 = 1.41421

[ale@syaoran sqrt]$ ./sqrt 4.0
Custom sqrt of: 4 = 2
libm sqrt of: 4 = 2

[ale@syaoran sqrt]$ ./sqrt 16.0
Custom sqrt of: 16 = 4
libm sqrt of: 16 = 4

[ale@syaoran sqrt]$ ./sqrt 32.0
Custom sqrt of: 32 = 5.65685
libm sqrt of: 32 = 5.65685

[ale@syaoran sqrt]$ ./sqrt 100.0
Custom sqrt of: 100 = 10
libm sqrt of: 100 = 10

[ale@syaoran sqrt]$ ./sqrt 1000000.0
Custom sqrt of: 1e+06 = 1000
libm sqrt of: 1e+06 = 1000

Exploring the Android NDK

I have been testing the Android development tools. From what I have learned, the tools are separated into two main products: the Android Software Development Kit (SDK) and the Android Native Development Kit (NDK).

The SDK was the first development toolkit for Android, and it only allowed applications to be written using the Java programming language. The NDK was released some time later as a toolchain to enable developers to write parts of their applications using native programming languages (C and C++).

One of the first programs I developed to get a feeling of what the Android SDK looked like was an OpenGL ES App that painted the background using a color degrade. I wrote it a couple of months ago, but today was the first time I ran it on a real Android device.

The resulting image can be seen in the following picture:

Android OpenGL ES Test

Other than trying the SDK, what I really wanted to do was to experience how hard it would be to rewrite part of the App in C and then having both integrated. It turns out adding components built using the NDK is not very hard (for pure C code), so I decided to try moving all the rendering code to a plain C function.

I started a new project and coded all the rendering logic in a C function that I called “render“. Then, the NDK was used to compile the C code into a JNI-compliant shared library and, finally, I wrote a simple Java wrapper that calls into the shared library’s render function to do all the drawing.

The wrapper is responsible for creating the Android “Activity”, setting up an OpenGL ES context, and calling the native C function. The native C function clears the background and does all the drawing.

Getting all the JNI requirements in place in order to have the Java runtime recognize the native library and call a native method was not too hard, but it was not trivial either. It is definitely much more complicated that calling native libraries from Objective-C or even Python. After a few tries, the bond was made:

Android OpenGL ES rendering from C code.

Clearly this is a very simple example where the C code could be tailored to fit JNI’s requirements from the start. I expect porting an existing C++ codebase to be much more difficult. However, I am looking forward to continue delving into Android’s development tools.

Opaque Type Oriented Programming in C

Non-C / Non-C++ programmers usually tend to confuse both languages with one another. Even though C++ is (in general) a superset of C, they are at their root very different programming languages.

One of the main reasons C is so different from C++ is that C does not include classes and thus would be -in principle- not object oriented.

The fact that the language does not include provisions for creating classes does not preclude adopting some sort of object-oriented programming style, however. In this post I’m going to introduce Opaque-Type Oriented Programming, a technique that allows our C programs to be structured as if they were populated with objects.

Opaque-Type oriented programming is a great programming paradigm that fits perfectly on structured programming. Ever since the first time I discovered this programming style, I keep finding in every piece of C code I come across.

As usual, I’m going to introduce this concept using examples. In this case, let’s suppose we want to build a dynamic array type for the C programming language.

Opaque Type Declaration and Interface

Here are some operations we might want our dynamic array to include:

  1. Initialization (with a given capacity).
  2. Adding an object (we’ll stick to real numbers for this example).
  3. Random access.
  4. Destruction.

Structured programming would dictate that we have some sort of static array allocated on the heap and that we have functions that work on said array.

In Opaque-Type oriented programming, the idea is fundamentally the same, but instead of dealing with the array directly, we are going to hide it in a struct together with all the attributes that are related to it, namely, its current size and capacity.

typedef float real; //change this to double for double precision

typedef struct dynarray
{
    real* pdata;
    size_t size;
    size_t capacity;
} dynarray;

Size and capacity are different things. Size denotes the current number of elements stored in our dynamic array, whereas capacity denotes the maximum number of elements that could be stored in this array.

Note: if you don’t know what size_t is, you can think of it as being an “unsigned long int”. It is 32 bits wide when compiled for a 32-bit OS, 64-bits wide when compiled for a 64-bit OS and it doesn’t have a sign bit.

Operations

Now that we have a struct to hold our dynamic array’s state and data, we need a set of functions that operate on this array. The idea behind Opaque-Type oriented programming is that we are able to use this dynamic array without having to know anything about its inner state or data.

The “contract” for operating on the dynamic array is given by the functions that we develop that know how to work on it.

Let’s start with initialization. In this example, we are going to have the initialization function allocate all inner state but the struct must be allocated by the user. This gives the programmer the option to create the array in the stack or in the heap.

void dynarray_init(dynarray* parray, size_t capacity) 
{
     parray->pdata = (real *)malloc(capacity*sizeof(real));
     parray->capacity = capacity;
     parray->size = 0;
}

It is usually a good idea to prefix all the functions with the name of the opaque type. This helps the code be auto-documentative and hels prevent name clashes (remember that C doesn’t have namespaces).

Now that the init function has initialized all the state and prepared the array to be used by the other functions, let’s have a look at the rest of them.

One of the hardest issues to solve when dealing with dynamic arrays is handling the case where we need reallocate the inner data because we are out of space (i.e. the capacity is exhausted).

While in structured programming this logic may very well pollute all of our algorithms, in Opaque-Type programming we can just hide this mechanism in our “add” function and have the programmer completely isolated from it.

void dynarray_add(dynarray* parray, real elem)
{
    if ((parray->size + 1) == parray->capacity)
    {
        // out of space, expand the array 
        // to hold 10 more elements:

        parray->pdata = (real *)realloc(parray->pdata, (parray->capacity + 10)*sizeof(real));
        assert(parray->pdata);
        parray->capacity += 10;

    }

    parray->pdata[parray->size] = elem;
    parray->size += 1;

}

There it is. From the outside world, no one needs to know how our dynamic array expands to fit as many elements as we want to store.

The only constraint that we have is if we run out of memory, and a real-world implementation of a dynamic array should handle this case. Another possible improvement would be not to allocate a fixed 10 more reals but a variable number that depends on how fast the array is actually growing (just to prevent lots of calls to realloc, which is kind of expensive).

Up to this point you can imagine how the “at” function would work. You would be surprised to learn how many of my students actually get this function wrong the first time just because they don’t pause and consider that the dynamic array doesn’t have the actual data but a pointer to it.

I’m going to present this function here just for the sake of completeness:

real dynarray_at(dynarray parray*, size_t position)
{
    assert(position < parray->size); //sanity check

    return parray->pdata[position]; // of course parray[position] compiles too, but it's a gross logic error!

}

Finally, I’m going to present how to deallocate our dynamic array. Remember we didn’t allocate the struct, so we are going to free the inner data only. I completely subscribe to Objective-C’s trail of though of “whoever allocates something should also release it”.

void dynarray_free(dynarray* parray)
{
    free(parray->pdata);
    parray->capacity = 0;
    parray->size = 0;
}

Putting it all together

Now that our opaque-type dynamic array is ready, I’m going to show you how we can use it in practice.

An exercise for the reader would be to compare how much easier this is compared to trying to handle the same data structure without any wrapping (consider this homework ; ))

#include <stdio.h>
int main()
{
    dynarray darray; //alloc'd in the stack, no free needed
    dynarray_init(&darray, 10); // initialize our array
    for (size_t i = 0; i < 100; i++)
    {
        dynarray_add(&darray, i); // it will automatically grow
    }
    for (size_t i = 0; i < 100; i++)
    {
        printf("%f ", dynarray_at(&darray, i));
    }
    printf("\n");
    dynarray_free(&darray); // remember to free inner data!
    return 0;
}

There we go. We now have a dynamic array implementation that automatically expands to accomodate new elements that would not fit in the current allocated storage.

Just as we wrapped the dynamic array, you might wrap any concept: matrices, hashmaps, semaphores, threads, you name it. You can actually find many opaque types in the UNIX standard C library (the pthread and BSD sockets libraries are great examples).

I hope this post gave you a good first impression on Opaque-Type oriented programming and that you find it useful in your structured software designs.

Posted in C

How to Create Shared Objects under Mac OS X

Some time ago a friend from Colombia asked me how to create a Shared Object under Mac OS X in order to be able to use it from a Python application (via ctypes).

Although the compiler most commonly used on Mac OS X is the GCC compiler, the version supplied by Apple has a handful of modifications which are specific for OS X. In particular, the -shared compiler flag, used to create Shared Objects under GNU/Linux based systems, has no effect.

In spite of this, Apple did introduce its own mechanism for creating Shared Objects along its GCC extensions. Assuming we have a file called mylibrary.c, which contains the source code for our library, we should invoke GCC like so:

gcc -dynamiclib mylibrary.c -o mylibrary.dylib

If the compiling process succeeds, this command should create a file called mylibrary.dylib, which can now be dynamically linked from other applications.

Unlike Linux based Operating Systems, Mac OS X seems to automatically add the current working directory to the DYLD_LIBRARY_PATH environment variable (an equivalent to Linux’s own LD_LIBRARY_PATH), thus it should not be required to modify this variable in order to have our applications link against the newly created library at runtime, unless, of course, we wanted to move the Shared Object to a directory different from where our process will be executing.

Finally, it should be taken into account that this mechanism is just a simplified example. A correct way to manage the library creation process would be by means of make or equivalent software construction applications, which provide several advantages to software developers.