#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