#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<fstream>
#include<vector>
#include<iostream>
#include<cassert>
#ifdef __linux__
#else
#define M_PI 3.141592653589793
#define INFINITY 1e8
#endif
template<typename T>
class Vec3
{
public:
    T x, y, z;
    Vec3() :x(T(0)), y(T(0)), z(T(0)){}
    Vec3(T xx) :x(xx), y(xx), z(xx){}
    Vec3(T xx, T yy, T zz) :x(xx), y(yy), z(zz){}
    Vec3&normalize()
    {
        T nor2 = length2();
        if (nor2>0){
            T invNor = 1 / sqrt(nor2);
            x *= invNor, y *= invNor, z *= invNor;
        }
        return*this;
    }
    Vec3<T>operator*(const T&f)const{ return Vec3<T>(x*f, y*f, z*f); }
    Vec3<T>operator*(const Vec3<T>&v)const{ return Vec3<T>(x*v.x, y*v.y, z*v.z); }
    T dot(const Vec3<T>&v)const{ return x*v.x + y*v.y + z*v.z; }
    Vec3<T>operator-(const Vec3<T>&v)const{ return Vec3<T>(x - v.x, y - v.y, z - v.z); }
    Vec3<T>operator+(const Vec3<T>&v)const{ return Vec3<T>(x + v.x, y + v.y, z + v.z); }
    Vec3<T>&operator+=(const Vec3<T>&v){ x += v.x, y += v.y, z += v.z; return*this; }
    Vec3<T>&operator*=(const Vec3<T>&v){ x *= v.x, y *= v.y, z *= v.z; return*this; }
    Vec3<T>operator-()const{ return Vec3<T>(-x, -y, -z); }
    T length2()const{ return x*x + y*y + z*z; }
    T length()const{ return sqrt(length2()); }
    friend std::ostream&operator<<(std::ostream&os, const Vec3<T>&v)
    {
        os << "[" << v.x << " " << v.y << " " << v.z << "]";
        return os;
    }
};
template<typename T>
class Sphere
{
public:
    Vec3<T>center;
    T radius, radius2;
    Vec3<T>surfaceColor, emissionColor;
    T transparency, reflection;
    Sphere(const Vec3<T>&c, const T&r, const Vec3<T>&sc,
        const T&refl = 0, const T&transp = 0, const Vec3<T>&ec = 0) :
        center(c), radius(r), radius2(r*r), surfaceColor(sc), emissionColor(ec),
        transparency(transp), reflection(refl)
    {}
    bool intersect(const Vec3<T>&rayorig, const Vec3<T>&raydir, T*t0 = NULL, T*t1 = NULL)const
    {
        Vec3<T>l = center - rayorig;
        T tca = l.dot(raydir);
        if (tca<0)return false;
        T d2 = l.dot(l) - tca*tca;
        if (d2>radius2)return false;
        T thc = sqrt(radius2 - d2);
        if (t0 != NULL&&t1 != NULL){
            *t0 = tca - thc;
            *t1 = tca + thc;
        }
        return true;
    }
};
#define MAX_RAY_DEPTH 5
template<typename T>
T mix(const T&a, const T&b, const T&mix)
{
    return b*mix + a*(T(1) - mix);
}
float min(float i, float j){ return(i<j ? i : j); }
float max(float i, float j){ return(i>j ? i : j); }
template<typename T>
Vec3<T>trace(const Vec3<T>&rayorig, const Vec3<T>&raydir,
    const std::vector<Sphere<T>*>&spheres, const int&depth)
{
    T tnear = INFINITY;
    const Sphere<T>*sphere = NULL;
    for (unsigned i = 0; i<spheres.size(); ++i){
        T t0 = INFINITY, t1 = INFINITY;
        if (spheres[i]->intersect(rayorig, raydir, &t0, &t1)){
            if (t0<0)t0 = t1;
            if (t0<tnear){
                tnear = t0;
                sphere = spheres[i];
            }
        }
    }
    if (!sphere)return Vec3<T>(2);
    Vec3<T>surfaceColor = 0;
    Vec3<T>phit = rayorig + raydir*tnear;
    Vec3<T>nhit = phit - sphere->center;
    nhit.normalize();
    T bias = 1e-4;
    bool inside = false;
    if (raydir.dot(nhit)>0)nhit = -nhit, inside = true;
    if ((sphere->transparency>0 || sphere->reflection>0) && depth<MAX_RAY_DEPTH){
        T facingratio = -raydir.dot(nhit);
        T fresneleffect = mix<T>(pow(1 - facingratio, 3), 1, 0.1);
        Vec3<T>refldir = raydir - nhit * 2 * raydir.dot(nhit);
        refldir.normalize();
        Vec3<T>reflection = trace(phit + nhit*bias, refldir, spheres, depth + 1);
        Vec3<T>refraction = 0;
        if (sphere->transparency){
            T ior = 1.1, eta = (inside) ? ior : 1 / ior;
            T cosi = -nhit.dot(raydir);
            T k = 1 - eta*eta*(1 - cosi*cosi);
            Vec3<T>refrdir = raydir*eta + nhit*(eta*cosi - sqrt(k));
            refrdir.normalize();
            refraction = trace(phit - nhit*bias, refrdir, spheres, depth + 1);
        }
        surfaceColor = (reflection*fresneleffect +
            refraction*(1 - fresneleffect)*sphere->transparency)*sphere->surfaceColor;
    }
    else{
        for (unsigned i = 0; i<spheres.size(); ++i){
            if (spheres[i]->emissionColor.x>0){
                Vec3<T>transmission = 1;
                Vec3<T>lightDirection = spheres[i]->center - phit;
                lightDirection.normalize();
                for (unsigned j = 0; j<spheres.size(); ++j){
                    if (i != j){
                        T t0, t1;
                        if (spheres[j]->intersect(phit + nhit*bias, lightDirection, &t0, &t1)){
                            transmission = 0;
                            break;
                        }
                    }
                }
                surfaceColor += sphere->surfaceColor*transmission*
                    max(T(0), nhit.dot(lightDirection))*spheres[i]->emissionColor;
            }
        }
    }
    return surfaceColor + sphere->emissionColor;
}
template<typename T>
void render(const std::vector<Sphere<T>*>&spheres)
{
    unsigned width = 64, height = 64;
    Vec3<T>*image = new Vec3<T>[width*height], *pixel = image;
    T invWidth = 1 / T(width), invHeight = 1 / T(height);
    T fov = 90, aspectratio = width / T(height);
    T angle = tan(M_PI*0.5*fov / T(180));
    for (unsigned y = 0; y<height; ++y){
        for (unsigned x = 0; x<width; ++x, ++pixel){
            T xx = (2 * ((x + 0.5)*invWidth) - 1)*angle*aspectratio;
            T yy = (1 - 2 * ((y + 0.5)*invHeight))*angle;
            Vec3<T>raydir(xx, yy, -1);
            raydir.normalize();
            *pixel = trace(Vec3<T>(0), raydir, spheres, 0);
        }
    }
    std::ofstream ofs("./untitled.ppm", std::ios::out | std::ios::binary);
    ofs << "P6\n" << width << " " << height << "\n255\n";
    for (unsigned i = 0; i<width*height; ++i){
        ofs << (unsigned char)(min(T(1), image[i].x) * 255) <<
            (unsigned char)(min(T(1), image[i].y) * 255) <<
            (unsigned char)(min(T(1), image[i].z) * 255);
    }
    ofs.close();
    delete[]image;
}
int main(int argc, char**argv)
{
    std::vector<Sphere<float>*>spheres;
    spheres.push_back(new Sphere<float>(Vec3<float>(0, -10004, -20), 10000, Vec3<float>(0.2), 0, 0.0));
    spheres.push_back(new Sphere<float>(Vec3<float>(0, 0, -20), 4, Vec3<float>(1.00, 0.32, 0.36), 0, 0.0));
    spheres.push_back(new Sphere<float>(Vec3<float>(5, -1, -15), 2, Vec3<float>(0.90, 0.76, 0.46), 0, 0.0));
    spheres.push_back(new Sphere<float>(Vec3<float>(5, 0, -25), 3, Vec3<float>(0.65, 0.77, 0.97), 0, 0.0));
    spheres.push_back(new Sphere<float>(Vec3<float>(-5.5, 0, -15), 3, Vec3<float>(0.90, 0.90, 0.90), 0, 0.0));
    spheres.push_back(new Sphere<float>(Vec3<float>(0, 20, -30), 3, Vec3<float>(0), 0, 0, Vec3<float>(3)));
    render<float>(spheres);
    while (!spheres.empty()){
        Sphere<float>*sph = spheres.back();
        spheres.pop_back();
        delete sph;
    }
    return 0;
}