6

I'm trying to implement the Mandelbrot set in C, but I'm having a weird problem. My code is as follows:

#include <stdio.h>
#include <math.h>
#include <complex.h>

int iterate_pt(complex c);

int main() {
FILE *fp;
fp = fopen("mand.ppm", "w+");


double crmin = -.75;
double crmax = -.74;
double cimin = -.138;
double cimax = -.75; //Changing this value to -.127 fixed my problem.

int ncols = 256;
int nrows = 256;
int mand[ncols][nrows];
int x, y, color;
double complex c;

double dx = (crmax-crmin)/ncols;
double dy = (cimax-cimin)/nrows;

for (x = 0; x < ncols; x++){
    for (y = 0; y < nrows; y++){
        double complex imaginary = 0+1.0i;
        c = crmin+(x*dx) + (cimin+(y*dy)) * imaginary;
        mand[x][y] = iterate_pt(c);
    }
}

printf("Printing ppm header.");
fprintf(fp, "P3\n");
fprintf(fp, "%d %d\n255\n\n", ncols, nrows);

for (x = 0; x < ncols; x++) {
    for (y = 0; y < nrows; y++){
        color = mand[x][y];
        fprintf(fp, "%d\n", color);
        fprintf(fp, "%d\n", color);
        fprintf(fp, "%d\n\n", color); //Extra new line added, telling the ppm to go to next pixel.
    }
}
fclose(fp);

return 0;
}

int iterate_pt(double complex c){
double complex z = 0+0.0i;
int iterations = 0;
int k;
for (k = 1; k <= 255; k++) {
    z = z*z + c;
    if (sqrt( z*conj(z) ) > 50){
        break;
    }
    else
        ++iterations;
}
return iterations;
}

However, the output of this program, which is stored as a ppm file looks like this:

Converted to a GIF using GIMP. I can confirm that the GIF and original PPM look exactly the same as a PPM and GIF

Thanks for your help!

Nathan Jones
  • 4,904
  • 9
  • 44
  • 70
  • Check your logic against this pseudo code: http://en.wikipedia.org/wiki/Mandelbrot_set#For_programmers – Blender Oct 20 '11 at 00:31
  • I don't know if this relevant but I get `warning: unused variable 'real'` when I compile your code. – Sinan Ünür Oct 20 '11 at 00:36
  • I was experimenting with different techniques to represent the real and imaginary numbers, and forgot to remove that. – Nathan Jones Oct 20 '11 at 00:41
  • 1
    I don't see why you're doing this: `if (sqrt( z*conj(z) ) > 50){`. If you're trying to get the magnitude, use `cabs(z)`. – Tom Zych Oct 20 '11 at 00:44
  • I'm not exactly sure what the thinking behind this was. It was a replacement that my professor suggested. He claimed that it was more accurate than using the absolute value function directly. – Nathan Jones Oct 20 '11 at 00:46
  • Random +1 for posting a coding question with some geometrical content. – Peter Oct 20 '11 at 01:00
  • Note re my earlier comment: DUH. `sqrt(z * conj(z))` **is** the magnitude, which I would have realized if I'd thought about it for ten seconds. The real improvement would be to skip the `sqrt` and compare to 2500. – Tom Zych Oct 20 '11 at 12:42
  • Actually, `cabs(z)` is *more* accurate than `sqrt(z*conj(z))` on many platforms, and should be at least as accurate on the others. – Stephen Canon Oct 28 '11 at 12:34

3 Answers3

3

Try setting cimax to -0.127, I'm also working on this project and it seems to do the trick ;)

geekingreen
  • 982
  • 1
  • 7
  • 8
2

The code looks good. But your starting rectangle doesn't look right!

you are using 

Real ranage [  -.75 ,  -.74 ]
Imag range  [ -.138 ,  -.75 ]

are you sure this is what you intended? It seems like an awfully stretched y-scale to me.

Also, standard mandelbrot algorithms tend to use

magnitude > 2

rather than 50. as an escape check. Though this shouldn't affect the actual shape of the set.

Sanjay Manohar
  • 6,920
  • 3
  • 35
  • 58
0

BTW, there's no point in computing the sqrt of z*conj(z). Simply square the expressions on both sides of the inequality, giving if (z*conj(z) > 2500) and you've boosted the performance.

Jens
  • 69,818
  • 15
  • 125
  • 179