#ifndef GEOMETRY_H
#define GEOMETRY_H


#include <cmath>


typedef double Real;

const Real precision = 1e-10;
const Real pi = 3.14159265358979323846;
const Real sqrt2 = 0.7071067811865475244;

inline int solve(Real* roots, Real k0, Real k1) {
    if (fabs(k1) <= precision) return 0;
    if (roots) roots[0] = -k0/k1;
    return 1;
}

inline int solve(Real* roots, Real k0, Real k1, Real k2) {
    if (fabs(k2) <= precision*precision) return solve(roots, k0, k1);
    Real D = k1*k1 - 4*k2*k0;
    if (fabs(D) <= precision*precision) {
        if (roots) roots[0] = -0.5*k1/k2;
        return 1;
    } else
    if (D > 0) {
        if (roots) {
            Real a = sqrt(D);
            Real b = -0.5/k2;
            roots[0] = (k1 - a)*b;
            roots[1] = (k1 + a)*b;
        }
        return 2;
    }
    return 0;
}


inline Real sign(Real x, Real precision) {
    return x < -precision ? -1
         : x >  precision ?  1
         : 0;
}


Real real_random();

inline Real real_random2()
    { return real_random()*2 - 1; }
inline Real real_random_normal()
    { return sqrt(-2*log(real_random())) * sin(2*pi*real_random()); }


class Vector2 {
public:
    union {
        struct { Real x, y; };
        struct { Real c[2]; };
    };
    
    explicit Vector2(Real x = 0, Real y = 0):
        x(x), y(y) { }
        
    Real len_sqr() const { return x*x + y*y; }
    Real len() const { return sqrt(len_sqr()); }
};


class Vector3 {
public:
    union {
        struct { Real x, y, z; };
        struct { Real r, g, b; };
        struct { Real c[3]; };
    };
    
    explicit Vector3(Real x = 0, Real y = 0, Real z = 0):
        x(x), y(y), z(z) { }
    
    Vector3 operator+(const Vector3 &v) const
        { return Vector3(x+v.x, y+v.y, z+v.z); }
    Vector3 operator-(const Vector3 &v) const
        { return Vector3(x-v.x, y-v.y, z-v.z); }

    Real operator*(const Vector3 &v) const
        { return x*v.x + y*v.y + z*v.z; }
    Vector3 operator*(Real k) const
        { return Vector3(x*k, y*k, z*k); }
    Vector3 operator/(Real k) const
        { return *this*(1/k); }

    Vector3& operator+=(const Vector3 &v)
        { return *this = *this + v; }
    Vector3& operator-=(const Vector3 &v)
        { return *this = *this - v; }

    Vector3 operator*=(Real k)
        { return *this = *this*k; }
    Vector3 operator/=(Real k)
        { return *this = *this/k; }
        
    Vector3 cross(const Vector3 &v) const
        { return Vector3(y*v.z - z*v.y, z*v.x - x*v.z, x*v.y - y*v.x); }
    Vector3 perp() const
        { return fabs(y) > fabs(x) ? Vector3(0, z, -y) : Vector3(-z, 0, x); }

    Real len_sqr() const { return x*x + y*y + z*z; }
    Real len() const { return sqrt(len_sqr()); }

    Vector3 norm() const
        { Real l = len(); return l > precision ? *this/l : Vector3(); }
};


class Vector4 {
public:
    union {
        struct { Real x, y, z, w; };
        struct { Real r, g, b, a; };
        struct { Real c[4]; };
    };
    
    explicit Vector4(Real x = 0, Real y = 0, Real z = 0, Real w = 0):
        x(x), y(y), z(z), w(w) { }
};


class Matrix4 {
public:
    union {
        struct {
            Real m00, m01, m02, m03,
                 m10, m11, m12, m13,
                 m20, m21, m22, m23,
                 m30, m31, m32, m33;
        };
        struct { Real m[4][4]; };
        struct { Real a[16]; };
    };

    explicit Matrix4(
        Real m00 = 1, Real m01 = 0, Real m02 = 0, Real m03 = 0,
        Real m10 = 0, Real m11 = 1, Real m12 = 0, Real m13 = 0,
        Real m20 = 0, Real m21 = 0, Real m22 = 1, Real m23 = 0,
        Real m30 = 0, Real m31 = 0, Real m32 = 0, Real m33 = 1
    ):
        m00(m00), m01(m01), m02(m02), m03(m03),
        m10(m10), m11(m11), m12(m12), m13(m13),
        m20(m20), m21(m21), m22(m22), m23(m23),
        m30(m30), m31(m31), m32(m32), m33(m33)
    { }
    
    Vector4 operator*(const Vector4 &v) const {
        return Vector4(
            m00*v.x + m10*v.y + m20*v.z + m30*v.w,
            m01*v.x + m11*v.y + m21*v.z + m31*v.w,
            m02*v.x + m12*v.y + m22*v.z + m32*v.w,
            m03*v.x + m13*v.y + m23*v.z + m33*v.w );
    }

    Vector3 transform(const Vector3 &v, Real w = 1) const {
        return Vector3(
            m00*v.x + m10*v.y + m20*v.z + m30*w,
            m01*v.x + m11*v.y + m21*v.z + m31*w,
            m02*v.x + m12*v.y + m22*v.z + m32*w );
    }

    Matrix4 operator*(const Matrix4 &m) const;
    
    Real det() const;
    Matrix4 inv() const;
    
    static Matrix4 zero();
    static Matrix4 translation(const Vector3 &translate);
    static Matrix4 scaling(const Vector3 &scale);
    static Matrix4 rotation(const Vector3 &axis, Real angle);
    static Matrix4 perspective(Real fovy, Real aspect, Real z_near, Real z_far);
};

#endif