// sphere flake bvh raytracer (c) 2005, thierry berger-perrin <tbptbp@gmail.com>
// this code is released under the GNU Public License.
#include <cmath>       // see http://ompf.org/ray/sphereflake/
#include <cstdlib>
#include <iostream>    // compile with ie g++ -O2 -ffast-math sphereflake.cc
#define GIMME_SHADOWS  // usage: ./sphereflake [lvl=6] >pix.ppm

enum { childs = 9, ss= 2, ss_sqr = ss*ss }; /* not really tweakable anymore */
static const double infinity = 1./0, epsilon = 1e-12;

// LLVM LOCAL begin
// Implementations of mathematical functions may vary slightly in the accuracy of
// their results, typically only in the least significant bit.  Round to make
// the results consistent across platforms.
//
// We don't care much about precision, and we do *really*
// want this *NOT* to the the costliest part of any benchmark
#define FACT3 6
#define FACT5 120
#define PI    3.141592656
#define PI2   6.283185307
#define PI_2  1.570796327
#define PI3_2 4.712388984

static int LLVMdiff(double a, double b) {
  double d = a - b;
  return (d > epsilon || d < -epsilon);
}

// Simple iterative multiplication
static double LLVMpow(double d, int n) {
  int i;
  double res = d;
  if (n == 0)
    return 1;
  for (i=1; i<n; i++)
    res *= d;
  return res;
}

// Simplified Taylor series
static double LLVMsin(double d) {
  double sign = 1.0;

  /* move into 2PI area */
  while (d < 0)
    d += PI2;
  while (d > PI2)
    d -= PI2;
  /* move into PI/2 area */
  if (d > PI3_2) {
    d = PI2 - d;
    sign = -1.0;
  } else if (d > PI) {
    d -= PI;
    sign = -1.0;
  } else if (d > PI_2) {
    d = PI - d;
  }
  /* series terms */
  double f3 = LLVMpow(d, 3)/FACT3;
  double f5 = LLVMpow(d, 5)/FACT5;
  d = sign * (d - f3 + f5);
  /* saturate */
  if (d > 1.0)
    d = 1.0;
  if (d < -1.0)
    d = -1.0;
  return d;
}

// Sine displacement
static double LLVMcos(double d) {
  return LLVMsin(d + PI_2);
}

// Simple Newton iteration (for values < 10^10)
static double LLVMsqrt(double x) {
  double xn=0, xk=1;
  /* special case where we'd return nan */
  if (x == infinity)
    return infinity;
  int it=100;
  while (it--) {
    xn = (xk + x/xk)/2;
    if (!LLVMdiff(xn, xk))
      return xn;
    xk = xn;
  }
  /* Return our best approximation */
  return xn;
}

// LLVM LOCAL end

struct v_t{ double x,y,z;v_t(){}
	v_t(const double a,const double b,const double c):x(a),y(b),z(c){}
	v_t operator+(const v_t&v)const{return v_t(x+v.x,y+v.y,z+v.z);}
	v_t operator-(const v_t&v)const{return v_t(x-v.x,y-v.y,z-v.z);}
	v_t operator-()const{return v_t(-x,-y,-z);}
	v_t operator*(const double d)const{return v_t(x*d,y*d,z*d);}
	v_t cross(const v_t&v)const{return v_t(y*v.z-z*v.y,z*v.x-x*v.z,x*v.y-y*v.x);}
	v_t norm()const{return*this*(1./LLVMsqrt(magsqr()));}
	double dot(const v_t&v)const{return x*v.x+y*v.y+z*v.z;}
	double magsqr()const{return dot(*this);}
};

//static const v_t light(v_t(0.5,-.95,1.775).norm()); /*pick one*/
static const v_t light(v_t(-0.5,-.65,.9).norm()); /*fiat lux*/

struct ray_t{
	v_t o,d;
	ray_t(const v_t&v):o(v){}
	 ray_t(const v_t&v,const v_t&w):o(v),d(w){}
};
struct hit_t {
	v_t n;
	double t;
	hit_t():n(v_t(0,0,0)),t(infinity){}
};

struct sphere_t{ 
	v_t o;
	double r;
	sphere_t(){}
	sphere_t(const v_t&v,double d):o(v),r(d){}
	v_t get_normal(const v_t&v)const{return(v-o)*(1./r);}
	double intersect(const ray_t&ray)const{
		const v_t v(o-ray.o); const double b=ray.d.dot(v),disc=b*b-v.magsqr()+r*r;
		if(disc < 0.)
			return infinity; /*branch away from the square root*/
		const double d=LLVMsqrt(disc), t2=b+d, t1=b-d; /*cond. move*/
		if(t2 < 0.)
			return infinity;
		else
			return(t1 > 0.? t1 : t2);
	}
};

struct node_t; 
static node_t *pool=0, *end=0;

struct node_t { /*a bvh in array form+skip for navigation.*/
	sphere_t bound,leaf;
	long diff;/*far from optimal*/
	node_t(){} node_t(const sphere_t&b,const sphere_t&l,const long jump) :bound(b),leaf(l),diff(jump){}
	template<bool shadow> static void intersect(const ray_t &ray,hit_t &hit){
		const node_t*p=pool;
		while(p < end) {
			if(p->bound.intersect(ray)>=hit.t) /*missed bound*/
				p+=p->diff; /*skip subtree*/
			else{
				const double t=p->leaf.intersect(ray);
				if(t < hit.t) { /*if hit, update, then break for shadows*/
					 hit.t=t;
					if(shadow) break;
					hit.n=p->leaf.get_normal(ray.o+ray.d*t);
				}
				++p; /*next!*/
			}
		}
	}
};

static double ray_trace(const node_t*const scene,const ray_t&ray) {
	hit_t hit;
	scene->intersect<false>(ray,hit);// trace primary
	const double diffuse = hit.t==infinity ? 0. : -hit.n.dot(light);
	#ifdef GIMME_SHADOWS
		if (diffuse <= 0.)
			return 0.;
		const ray_t sray(ray.o+(ray.d*hit.t)+(hit.n*epsilon),-light);
		hit_t shit;
		scene->intersect<true>(sray,shit);// trace shadow
		return shit.t==infinity ? diffuse : 0.;
	#else
		return diffuse > 0. ? diffuse : 0.;
	#endif
}

static const double grid[ss_sqr][2]={ /*our rotated grid*/
	{-3/3.,-1/3.},{+1/3.,-3/3.},
	{-1/3.,+3/3.},{+3/3.,+1/3.}
};
static void trace_rgss(const int width,const int height) {
	const double w=width,h=height,rcp=1/double(ss),scale=256./double(ss_sqr);
	ray_t ray(v_t(0,0,-4.5)); /* eye, looking into Z */
	v_t rgss[ss_sqr];
	for(int i=0;i<ss_sqr;++i) /*precomp.*/
		rgss[i]=v_t(grid[i][0]*rcp-w/2.,grid[i][1]*rcp-h/2.,0);
	
	v_t scan(0,w-1,std::max(w,h)); /*scan line*/
	for(int i=height;i;--i) {
		for(int j=width;j;--j) {
			double g=0;
			for(int idx=0;idx < ss_sqr;++idx){ /*AA*/
				ray.d=(scan+rgss[idx]).norm();
				g+=ray_trace(pool,ray); /*trace*/
			}
			std::cout << int(scale*g)<< " ";
			scan.x+=1; /*next pixel*/
		}
		scan.x=0;scan.y-=1; /*next line*/
	}
	std::cout << std::endl; 
}
	
struct basis_t{ /* bogus and compact, exactly what we need */
	v_t up,b1,b2;
	basis_t(const v_t&v){ const v_t n(v.norm());
		if ((n.x*n.x !=1.)&(n.y*n.y !=1.)&(n.z*n.z !=1.)) {/*cough*/
			b1=n;
			if(n.y*n.y>n.x*n.x) {
				if(n.y*n.y>n.z*n.z)
					b1.y=-b1.y;
				else b1.z=-b1.z;
			}
			else if(n.z*n.z > n.x*n.x)
				b1.z=-b1.z;
			else b1.x=-b1.x;
		}
		else
			b1=v_t(n.z,n.x,n.y);/*leaves some cases out,dodge them*/
		
		up=n;
		b2=up.cross(b1);
		b1=up.cross(b2);
	}
};

static node_t *create(node_t*n,const int lvl,int dist,v_t c,v_t d,double r) {
	n = 1 + new (n) node_t(sphere_t(c,2.*r),sphere_t(c,r), lvl > 1 ? dist : 1);
	if (lvl <= 1)
		return n; /*if not at the bottom, recurse a bit more*/
	
	dist=std::max((dist-childs)/childs,1); const basis_t b(d); 
	const double nr=r*1/3.,daL=2.*M_PI/6.,daU=2.*M_PI/3.; double a=0;
	for(int i=0;i<6;++i){ /*lower ring*/
		const v_t ndir((d*-.2+b.b1*LLVMsin(a)+b.b2*LLVMcos(a)).norm()); /*transcendentals?!*/
		n=create(n,lvl-1,dist,c+ndir*(r+nr),ndir,nr);
		a+=daL;
	}
	a-=daL/3.;/*tweak*/
	for(int i=0;i<3;++i){ /*upper ring*/
		const v_t ndir((d*+.6+b.b1*LLVMsin(a)+b.b2*LLVMcos(a)).norm());
		n=create(n,lvl-1,dist,c+ndir*(r+nr),ndir,nr); a+=daU;
	}
	return n;
}

#ifdef SMALL_PROBLEM_SIZE
#define TEST_SIZE 256
#else
#define TEST_SIZE 1024
#endif

int main(int argc,char*argv[]){
	enum{ w = TEST_SIZE, h = w }; /* resolution */
	const int lvl=(argc==2?std::max(atoi(argv[1]),2):6);
	int count=childs, dec=lvl;
	while(--dec > 1) count=(count*childs)+childs;
	++count;
	pool=new node_t[count];  /* raw */
	end=pool+count;
	create(pool,lvl,count,v_t(0,0,0),v_t(+.25,+1,-.5).norm(),1.); /* cooked */
	std::cout << "P2\n" << w << " " << h << "\n256\n";
	trace_rgss(w,h); /* served */
	return 0;
}