3D fractal

#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>
#define L 400
int cnt[L][L];
#define SCALE 0.7
#define CMAX (9)

void dump()
{
        FILE *fp=fopen("a.pgm","w");
        fprintf(fp, "P5\n%d %d\n255\n",L,L);
        int ix,iy;
        for(ix=0;ix<L;ix++){
                for(iy=0;iy<L;iy++){
                        fputc( (cnt[ix][iy]&255),fp);
                }}
        fclose(fp);
}


#define ISQ3 0.577350269189626
double n4vec[4][3]={
        { ISQ3,  ISQ3, ISQ3},
        {-ISQ3, -ISQ3, ISQ3},
        {-ISQ3, ISQ3, -ISQ3},
        { ISQ3,-ISQ3, -ISQ3}
};
double nvec[6][3]={
        {1,0,0},
        {0,1,0},
        {0,0,1},
        {-1,0,0},
        {0,-1,0},
        {0,0,-1}
};

double n3vec[3][3]={
        {1,0,0},{0,1,0},{0,0,1}};

// simple 3D version of complex square, w.r.t some axis v
// rnew = 2(r.v) r - (r.r)v
void conv(double r[3], int ax)
{
        double *v = nvec[ax];
        double rv = r[0]*v[0]+r[1]*v[1]+r[2]*v[2];
        double rr = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];

        r[0] = 2*rv * r[0] - rr *v[0] ;
        r[1] = 2*rv * r[1] - rr *v[1] ;
        r[2] = 2*rv * r[2] - rr *v[2] ;
}

double rfunc(double c0[3], int c)
{
        int cc;
        double rr[3]={c0[0],c0[1],c0[2]};
        for(cc=0;cc<c;cc++)
        {
                conv(rr, 0);
                conv(rr, 1);            
                conv(rr, 2);
                rr[0]+=c0[0]; rr[1]+=c0[1]; rr[2]+=c0[2];
        }
        return rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2] - 4.0;
}

double bisec_z(int lvl, double p0[3], double dz, double eps, double zmax)
{
        while(p0[2]> -zmax)
        {
                double r2 = rfunc(p0,lvl);
                if (r2<0.0) break;
                p0[2] -= dz;
        }
        if (p0[2]< -zmax)return -999.0;
        p0[2] += dz;
        dz/=2;

        while(dz>eps)
        {
                p0[2] -= dz;
                double r2 = rfunc(p0,lvl);
                if (r2<0.0) p0[2]+=dz;
                dz/=2;
        }
        return p0[2];
}

int  main(int argc, char **argv)
{
        int ix,iy;
        double c0[3],rr[3];
        for(ix=0;ix<L;ix++){
                printf("%d/%d\n",ix,L);
                for(iy=0;iy<L;iy++){

                        c0[0]=SCALE*(2.0*ix/L -1.0);
                        c0[1]=SCALE*(2.0*iy/L -1.0);

                        double rr=c0[0]*c0[0]+c0[1]*c0[1];
                        if(rr>=4.0){ cnt[ix][iy]=0;continue;}

                        double zmax=sqrt(4-c0[0]*c0[0]-c0[1]*c0[1]);
                        double dz=0.1;
                        double zz=zmax+0.01;


                        c0[2]=zmax;

                        int c;

                        int skip=0;
                        for(c=0;c<=CMAX;c++)
                        {
                                double z0=bisec_z(c, c0, dz, 1e-7, zmax);
                                if (z0<-888.0)
                                {
                                        skip=1;
                                        break;
                                }
                        }

                        if(skip) { cnt[ix][iy]=0;continue;}

#define EPS 1e-9
                        double r2=rfunc(c0,CMAX);
                        double cz[3]={c0[0],c0[1],c0[2]+EPS};
                        double r2z=rfunc(cz,CMAX)-r2;
                        double cx[3]={c0[0]+EPS,c0[1],c0[2]};
                        double r2x=rfunc(cx,CMAX)-r2;
                        double cy[3]={c0[0],c0[1]+EPS,c0[2]};
                        double r2y=rfunc(cy,CMAX)-r2;
                        double nn=sqrt(r2z*r2z+r2x*r2x+r2y*r2y);

                        cnt[ix][iy]=(int)((0.2+0.8*fabs( r2z/nn))*255);

                }}
                dump();
}