#include "mathlink.h"
#define MAXCOUNT 100
#define MAXVALUE 100000

typedef struct {
  double re;
  double im;
} complex;

complex f(complex z,complex mu) {
  complex ret;
  ret.re=mu.re-z.im*z.im+z.re*z.re;
  ret.im=-mu.im+2*z.im*z.re;
  return ret;
}

void mandelbrotC(double x0, double x1, double y0, double y1, int steps) {
  long dimensions[2]={(long)steps,(long)steps};
  char *heads[2]={"List","List"};
  double x,y;
  int i,j,count;
  complex z,mu;
  int *res = (int*) malloc(steps*steps*sizeof(int));
  double stepx=(x1-x0)/steps;
  double stepy=(y1-y0)/steps;
  for(x=x0,i=0;i<steps;i++,x+=stepx)      
    for(y=y0,j=0;j<steps;j++,y+=stepy) {
      mu.re=x;
      mu.im=y;
      z.re=0.0;
      z.im=0.0;
      for(count=0;count<MAXCOUNT;count++) {
        z=f(z,mu);
        if(fabs(z.im)>MAXVALUE) break;
      }
      res[i+j*steps]=count;
    }
  MLPutIntegerArray(stdlink,res,dimensions,heads,2);
  free(res);
}

int main(int argc, char *argv[]) {
  return MLMain(argc, argv);
}
