Z^Z^Z.. Julia set

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

// color stuff
void hue2rgb(int hue, int vmax, int *r, int *g, int *b)
{
    int *max=r;
    int *p1=g;
    int *p2=b;

    int h0= hue+60;
    if(h0>=360)h0-=360;

    if(h0>=240)
    {
        max=b;
        p1=r;
        p2=g;
        h0-=240;
    }
    else if(h0>=120)
    {
        max=g;
        p1=b;
        p2=r;
        h0-=120;
    }

    *max= vmax;
    if(h0<60)
    {
        *p1=0;
        *p2=(int)(vmax*(60-h0)/60);
    }
    else
    {
        *p2=0;
        *p1=(int)(vmax*(h0-60)/60);
    }
}

// complex op. iz=1/z
void cinv(double zr, double zi, double *izr, double *izi)
{
    double r2=zr*zr+zi*zi;
    *izr = zr/r2;
    *izi = -zi/r2;
}

// newZ = exp(a z)
void ite(double ar, double ai, double zr, double zi, double *nzr, double *nzi)
{
    double logr= ar*zr-ai*zi;
    double th=   ar*zi+ai*zr;
    *nzr =exp(logr)*cos(th);
    *nzi =exp(logr)*sin(th);
}

//solve f(z)=exp(az)-z =0
void newton(double ar, double ai, double *zr0, double *zi0)
{

    //initial val
    double zr, zi;

    zr=*zr0;
    zi=*zi0;

    int i;
    for(i=0;i<1000;i++)
    {
        // exp(az)
        double exr,exi;
        ite(ar,ai, zr, zi, &exr, &exi);

        // f(z)=exp(az)-z
        double fr=exr-zr;
        double fi=exi-zi;

        // df/dz= a exp(az)-1
        double dfr = exr*ar - exi*ai -1.0;
        double dfi = exr*ai + exi*ar;

        // t=1/(a exp(az)-1)
        double tr,ti;
        cinv(dfr,dfi, &tr,&ti);

        // f*t
        double dzr= fr*tr-fi*ti;
        double dzi= fr*ti+fi*tr;

        double r2=dzr*dzr+dzi*dzi;
        if(r2<1e-20)break;
        zr -= dzr;
        zi -= dzi;
    }

    *zr0=zr;
    *zi0=zi;
}

char *buf;
int *ctab;
int d2i(double d, double lmax, int div)
{
    int x=(int)( (d+lmax)/(lmax*2)*div);
    if(x<0 || x>=div)return -1;
    return x;
}

void putwh(double x, double y, double lmax, int div)
{
    int ix=d2i(x, lmax, div);
    int iy=d2i(y, lmax, div);
    if(ix==-1 || iy==-1)return;

    int xx,yy;
    for(xx=0;xx<3;xx++){
    for(yy=0;yy<3;yy++){
	int ix2=ix+xx;
	int iy2=iy+yy;
	if(ix2>=div)continue;
	if(iy2>=div)continue;
	char *p = &(buf[(ix2+(div-1-iy2)*div)*3]);
	p[0]=p[1]=p[2]=255;
    }}
}

int  main(int argc, char **argv)
{
    FILE *fp;

    double r2max=1e+20;
    int div=800;
    double lmax=4.0;

    int x,y,c;
    buf=(char *)malloc(div*div*3);

    int cmax=720;
    double *xtab=(double *)malloc(cmax*sizeof(double));
    double *ytab=(double *)malloc(cmax*sizeof(double));

    ctab = (int *)malloc(sizeof(int)*div*div);

    double ar=atof(argv[1]);
    double ai=atof(argv[2]);

    int cmin = -1;
    double xmin,ymin;

    for(x=0;x<div;x++){
        for(y=0;y<div;y++){

            int c;
            int r,g,b;
            char *p=&buf[(x+ (div-1-y)*div)*3];

            double x0= x*lmax*2/div -lmax;
            double y0= y*lmax*2/div -lmax;

            //  printf("%e %e\n",cr,ci);
            for(c=0;c<cmax;c++)
            {
                xtab[c]=x0;
                ytab[c]=y0;

                // check periodic convergence
                int pp;
                int done=0;
#define PMAX 199
                for(pp=1;pp<=PMAX;pp++)
                {
                    if(c-pp<0) break;
                    double dx=x0-xtab[c-pp];
                    double dy=y0-ytab[c-pp];
                    double rr=dx*dx+dy*dy;
                    if(rr<1e-30)
                    {
                        // use 12 colors
                        hue2rgb(((c-pp)*4)%180, 255, &r, &g, &b);
                        done=1;
			printf("#Period %d\n",pp);

			if((cmin==-1) || (cmin>c))
			{
			    cmin=c;
			    xmin=x0;
			    ymin=y0;
			    //printf("conv %f %f  ini %f %f\n",x0,y0, xtab[0],ytab[0]);
			}
                        break;
                    }
                }
                if(done==1)break;

                double nx,ny;
                ite(ar,ai, x0,y0,&nx,&ny);

                //diverge?
                if(nx*nx+ny*ny>r2max)
                {
                    hue2rgb(180+(c*4)%180, 192, &r, &g, &b);
                    break;
                }
                x0=nx;
                y0=ny;
            }

            // neither converge or diverge
            if(c==cmax) r=g=b=128;

            // put
            p[0]=r;
            p[1]=g;
            p[2]=b;
        }}//xy

    double x0=xmin;
    double y0=ymin;
    for(c=0;c<2000+100;c++)
    {
	double nx,ny;
	ite(ar,ai, x0,y0,&nx,&ny);
	x0=nx;
	y0=ny;
	if(c>2000) putwh(x0,y0,lmax,div);
    }
    //printf("#conv %f %f\n",x0,y0);


    // axis
    for(x=0;x<div;x++)
    {
	char *p=&buf[(x+ div/2*div)*3];
	p[0]=p[1]=p[2]=0;
	p=&buf[(div/2+ x*div)*3];
	p[0]=p[1]=p[2]=0;
    }



    fp=fopen("a.ppm","w");
    fprintf(fp,"P6\n%d %d\n255\n",div,div);
    fwrite(buf,1,3*div*div,fp);
    fclose(fp);
    free(xtab);
    free(ytab);

}