//From:
//  http://codegolf.stackexchange.com/questions/35569/tweetable-mathematical-art

#include <iostream>
#include <cmath>


#define DIM 256
#define DM1 (DIM-1)

#define _sq(x) ((x)*(x))                           // square
#define _cb(x) abs((x)*(x)*(x))                    // absolute value of cube
#define _cr(x) (unsigned short)(pow((x),1.0/3.0))  // cube root

#define L(X,Y,Z) sqrt(X*X+Y*Y+Z*Z)
#define N(X,Y,Z) {float l=L(X,Y,Z);X/=l;Y/=l;Z/=l;}
float mandelbulb_de(float x,float y,float z, int n) {
    if (L(x,y,z)>7.0) return L(x,y,z)-6.0;

    float x2=x,y2=y,z2=z;
    float dr = 1.0;
    float r = 0.0; //initialization only to keep compiler happy

    #define POWER 8.0

    for (int i=0;i<n;++i) {
        r = L(x2,y2,z2);
        if (r>6.0) break;

        //convert to polar coordinates
        float theta = acos(z2/r);
        float phi = atan2(y2,x2);
        dr = pow(r,POWER-1.0)*POWER*dr + 1.0;

        //scale and rotate the point
        float zr = pow(r,POWER);
        theta *= POWER;
        phi *= POWER;

        //convert back to cartesian coordinates
        float s_theta = sin(theta);
        x2 = zr * s_theta*cos(phi) + x;
        y2 = zr * s_theta*sin(phi) + y;
        z2 = zr * cos(theta)       + z;
    }

    return 0.5*log(r)*r/dr;
}
unsigned char RD(int i,int j) {
    float x=0.0f, y=4.0f, z=0.0f;

    float dx = (float)(i)/DIM - 0.5f;
    float dz = (float)(j)/DIM - 0.5f;
    float dy = -1.0f;
    N(dx,dy,dz)

    float total = 0.0f;
    while (total<10.0f) {
        float d = mandelbulb_de(x,y,z,500);
        total += d;
        x += d * dx;
        y += d * dy;
        z += d * dz;
        if (d<=0.01f) {
            //return 0xFFu*total/5.0;
            float c = mandelbulb_de(x,y,z,500);
            float nx = mandelbulb_de(x+0.001,y,z,500) - c;
            float ny = mandelbulb_de(x,y+0.001,z,500) - c;
            float nz = mandelbulb_de(x,y,z+0.001,500) - c;
            N(nx,ny,nz)
            return 0xFFu*fmax(0,-nz);
        }
    }
    return 0x00u;
}
unsigned char GR(int i,int j) {
    return RD(i,j);
}
unsigned char BL(int i,int j) {
    return RD(i,j);
}

FILE* fp;
void pixel_write(int i, int j){
    static unsigned char color[3];
    color[0] = RD(i,j)&DM1;
    color[1] = GR(i,j)&DM1;
    color[2] = BL(i,j)&DM1;
    fwrite(color, sizeof(unsigned char),3, fp);
}
int main(int argc, char* argv[]) {
    fp = fopen("MathPic.ppm","wb");
    fprintf(fp, "P6\n%d %d\n255\n", DIM,DIM);
    for (int j=0;j<DIM;++j) {
        for (int i=0;i<DIM;++i) {
            pixel_write(i,j);
        }
    }
    fclose(fp);
    return 0;
}