What was missing is a 3D escape-time fractal that could be used like the Mandelbrot set. The constraints were that they should behave like complex numbers when the third coordinate was 0. In addition, only a squaring operation was needed, since the recursion function is u' = u ^ 2 + u_0. I tried by just expressing the 3D vertices in polar coordinates (r, phi, theta) and expressing the square with (r, phi, theta) -> (r^2, 2*phi, 2*theta). All points where the module of the series remains below a threshold (say 2000) after some iterations (15) are considered inside the fractal. The surface could be explored using a marching cubes algorithm. The result did not look so good: Mandelbrot was like a flat spaceship without many interesting features, so I abandoned the project.
The Mandelbulb function after a given number of iterations and its derivative w.r.t u_0 (useful to compute the surface normal), can be expressed as polynoms of high degrees. I used Maple to generate the C code that corresponds to them, which can be used as OpenCL almost as is. I ran it on a nVidia GTX 680 GPU to render the fractal in interactive time (ie. not-too-nervous people can navigate in the fractal). Here are few images of the results:
