Very nice article! I've never heard of PDL, though I do computational math (not numerical) everyday. I'm definitely going to find a use for it now.

Last line in MATLAB and Octave code, i.e., 'exit' causes that Matlab exits; hence, it is not possible to look at this fractal image.
mWo

Let me preface this by saying that PDL looks neat. But,

You're iterating over the whole image with every iteration, doing gobs of work for pixels/locations that have diverged. I don't believe that there is any useful information retained by your clamping-to-5 step. (It is pretty, though!) With Mandelbrot, if the distance from the origin is > 2 then divergence to infinity is inevitable. The Mandelbrot set itself is defined as those starting values which never diverge to infinity; I think but am not sure that your clamping makes it impossible to tell even what is in/out of the set.

(( I believe the canonical/typical coloring is based on #-of-iterations-before-divergence for those pixels that diverge and black for those that don't. ))

Here is a patch (sorry about lost whitespace/indenting):

37,39c37,39
< for(k=0;k < niter;k++){
< for (i=0;i < NPTS;i++) {
< for (j=0;j < NPTS;j++) {
---
> for (i=0;i < NPTS;i++) {
> for (j=0;j < NPTS;j++) {
> for(k=0;k < niter;k++){

48,52c48,51
< if (zRe[i][j] < -5.) zRe[i][j]=-5.;
< if (zRe[i][j] > 5.) zRe[i][j]=5.;
< if (zIm[i][j] < -5.) zIm[i][j]=-5.;
< if (zIm[i][j] > 5.) zIm[i][j]=5.;
---
> if (zRe[i][j]* zRe[i][j] + zIm[i][j]*zIm[i][j] > 4.) // Divergence criteria
> break;

How much improvement? 2000x2000 array of points, 1000 iterations, my laptop/gcc 4.1.2 (Note: I had to add 'static' to the local arrays to make it work):
- Old: 4m46.654s
- New: 0m16.373s

C wins again! ;)

I don't know how to get this speedup effect in any of the "whole-array" approaches/languages. (Probably because I don't know them well at all.) It seems hard to get any effect other than foreach(iteration) { foreach(location in image) { /*..do something.. */ } } ; with this order of operations there's no way to bail out of iterations early at diverging locations.

I would speculate that, for this problem, plain old perl would be faster than your PDL approach... ;)

If you replace the lines like

qgtfive= find(qRe > 5.);
zRe(qgtfive)=5.;

with

zRe(qRe>5)=5

at least in Octave the code is about twice as fast. I believe the same indexing trick can be used in recent versions of matlab as well. You are basically getting your speed improvement from the PDL clip function.

D.

Dan,

Something like

function z = mandel
npts=1000;
niter=51;
z=zeros(npts,npts);
[kRe,kIm] = meshgrid(linspace(-1.5,.5,npts),linspace(-1,1,npts));
k = kRe + 1i*kIm;
idx = [1:npts*npts];
while (niter)
z(idx) = z(idx).^ 2 + k(idx);
idx(abs(z(idx))>2) = [];
niter = niter - 1;
end
end

will do it for Octave/Matlab. However, at least in Octave the above indexed assignment can be quite expensive, and so the above is actually slower for this example. It might be faster with matlab...

D.

Hi David and Dan,

Here too, you can avoid indexed assignment by using a logical mask which is much faster in octave.

Cheers,
Alex

octave:7> tic; z=mandel(); toc;
Elapsed time is 43.521179 seconds.
octave:8> tic; z2=mandel2(); toc;
Elapsed time is 24.081490 seconds.
octave:9> max(abs(z(:)-z2(:)))
ans = 0

function z = mandel2
npts=1000;
niter=51;
z=zeros(npts,npts);
[kRe,kIm] = meshgrid(linspace(-1.5,.5,npts),linspace(-1,1,npts));
k = kRe + 1i*kIm;

while (niter)
mask = abs(z) < 2;
z(mask) = z(mask).^ 2 + k(mask);

niter = niter - 1;
end

That looks like it might work. Do octave/matlab have a way to apply a function to each element of a 2d array? Then the solution becomes very nice and clear, just creating a per-pixel 'mandelbrot' function and then applying it to each element of the (complex) image matrix.

gcc 4.x now provides for Fortran 90/95, which has a syntax similar to Matlab/Octave. A good Fortran 95 compiler should provide the ease of development of Matlab, with the speed quoted for Fortran 77. While you have the claim that Fortran is not rich in IO & maths functions you omit that there is more maths code in Fortran than in any other language - Netlib and Statlib being two of the best examples. I agree that Fortran is not terribly good for graphical IO.

Personally, i think C wins. I have to be honest and say that as a first approach to explore the problem I use the Matlab that my university provides me, but when it comes to flexibility, control, speed of execution and the ability to share code with my peers at the university, then C is the choice. Needless to say that there are available free paralelization tools to Fortran and C, and even though they need high computer skills, let's face it: if you are coping with mathematical, numerical computation software, high chances are that you are computer skilled.
I know, there's a lot of people (biologists, finance people, etc.) But you should look at the present in science, sneak into the current research journals and note that there's a tendency to people get more and more skilled in performing in silico simulations, coming from very different fields, not only math or computer science.
Nah!, IMHO, C (Fortran, if you like it better) win.

Compile with: gcc -lm -lpng mandlebrot.c .I used gcc 4.2.0.
#include <stdio.h>
#include <math.h>
#include <png.h>
#define NPTS 1000 //=width=height
png_uint_16 img[NPTS][NPTS][3];
void genImg();
void addPixel();
void writePng();
int main(void){ genImg(); writePng(); return 0; }
void genImg() {
int i,j,k,niter;
double zRe, zIm, kRe, kIm, qRe, qIm;
niter=51; // Number of iterations
for (i=0;i<NPTS;i++) {
for (j=0;j<NPTS;j++) {
// Generating z = 0 (real and imaginary part)
zRe=0.;
zIm=0.;
// Generating the constant k (real and imaginary part)
kRe=(double)i*2.0/((double)NPTS-1.)-1.5;
kIm=(double)j*2.0/((double)NPTS-1.)-1.;
//Iterating
for (k=0;k<niter;k++){
// Calculating q = z*z + k in complex space
// q is a temporary variable to store the result
qRe=zRe*zRe-zIm*zIm+kRe;
qIm=2.*zRe*zIm+kIm;
// Assigning the q values to z constraining between
// -10 and 10 to avoid numerical divergences
zRe=qRe;
zIm=qIm;
if (zRe < -10.) zRe=-10.; if (zRe > 10.) zRe=10.; if (zIm < -10.) zIm=-10.; if (zIm > 10.) zIm=10.;
}
addPixel(i,j,zRe,zIm);
}
}
}
void addPixel(int i, int j, double zRe, double zIm){
img[i][j][0] = (png_uint_16)(sin(log( sqrt(zRe*zRe+zIm*zIm)))*65536);
img[i][j][1] = (png_uint_16)(cos(log( sqrt(zRe*zRe+zIm*zIm)))*65536);
img[i][j][2] = (png_uint_16)(log( sqrt(zRe*zRe+zIm*zIm))*65536);
}
void writePng(){
FILE *fp = fopen("mandle.png", "wb");
png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, (png_voidp)NULL, NULL, NULL);
png_infop info_ptr = png_create_info_struct(png_ptr);
png_set_IHDR(png_ptr, info_ptr, NPTS, NPTS, 16, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
png_init_io(png_ptr, fp);
png_write_info(png_ptr, info_ptr);
png_set_swap(png_ptr); //big -> little endian
int i; for (i=0;i<NPTS;i++){ png_write_row(png_ptr, (png_bytep)img[i]); }
png_write_end(png_ptr,info_ptr);
png_destroy_write_struct(&png_ptr, &info_ptr);
}

Well, i know nothing of all this math. How do i get up my brownie points? Here are some interesting Linux fractal resources:

1. A cool open source screensaver which also uses fractal math to do it's thing. (http://electricsheep.org/)

2. Gnofract. The author says he think's it's the best fractal package currently available for Linux. Test it and decide for yourself (http://gnofract4d.sourceforge.net/)

3. For the purists, fractals are about math. For the arty, here's some fun fractal stuff you can do in GIMP (http://gentoo-wiki.com/TIP_GIMP_Fractal_Backgrounds)

4. Fyre is a tool for producing computational artwork based on histograms of iterated chaotic functions. Source is included on the site. (http://fyre.navi.cx/)

Fortran has a complex data type which would be appropriate here.
Also there is a FOSS implementation of IDL called GNUDL.
An finally R is a very interesting language that has not been looked at in this
comparison.