Blame geometry.h

d2b2b5
d2b2b5
#ifndef GEOMETRY_H
d2b2b5
#define GEOMETRY_H
d2b2b5
d2b2b5
d2b2b5
#include <cmath></cmath>
d2b2b5
d2b2b5
d2b2b5
typedef double Real;
d2b2b5
b2ca59
const Real precision = 1e-10;
e31ea0
e31ea0
inline int solve(Real* roots, Real k0, Real k1) {
e31ea0
    if (fabs(k1) <= precision) return 0;
e31ea0
    if (roots) roots[0] = -k0/k1;
e31ea0
    return 1;
e31ea0
}
e31ea0
e31ea0
inline int solve(Real* roots, Real k0, Real k1, Real k2) {
b2ca59
    if (fabs(k2) <= precision*precision) return solve(roots, k0, k1);
e31ea0
    Real D = k1*k1 - 4*k2*k0;
b2ca59
    if (fabs(D) <= precision*precision) {
e31ea0
        if (roots) roots[0] = -0.5*k1/k2;
e31ea0
        return 1;
e31ea0
    } else
e31ea0
    if (D > 0) {
e31ea0
        if (roots) {
e31ea0
            Real a = sqrt(D);
e31ea0
            Real b = -0.5/k2;
e31ea0
            roots[0] = (k1 - a)*b;
e31ea0
            roots[1] = (k1 + a)*b;
e31ea0
        }
e31ea0
        return 2;
e31ea0
    }
e31ea0
    return 0;
e31ea0
}
e31ea0
e31ea0
e31ea0
inline Real sign(Real x, Real precision) {
e31ea0
    return x < -precision ? -1
e31ea0
         : x >  precision ?  1
e31ea0
         : 0;
e31ea0
}
e31ea0
d2b2b5
d2b2b5
class Vector2 {
d2b2b5
public:
d2b2b5
    union {
d2b2b5
        struct { Real x, y; };
d2b2b5
        struct { Real c[2]; };
d2b2b5
    };
d2b2b5
    
d2b2b5
    explicit Vector2(Real x = 0, Real y = 0):
d2b2b5
        x(x), y(y) { }
d2b2b5
        
d2b2b5
    Real len_sqr() const { return x*x + y*y; }
d2b2b5
    Real len() const { return sqrt(len_sqr()); }
d2b2b5
};
d2b2b5
d2b2b5
d2b2b5
class Vector3 {
d2b2b5
public:
d2b2b5
    union {
d2b2b5
        struct { Real x, y, z; };
d2b2b5
        struct { Real r, g, b; };
d2b2b5
        struct { Real c[3]; };
d2b2b5
    };
d2b2b5
    
d2b2b5
    explicit Vector3(Real x = 0, Real y = 0, Real z = 0):
d2b2b5
        x(x), y(y), z(z) { }
8ee194
    
8ee194
    Vector3 operator+(const Vector3 &v) const
8ee194
        { return Vector3(x+v.x, y+v.y, z+v.z); }
8ee194
    Vector3 operator-(const Vector3 &v) const
8ee194
        { return Vector3(x-v.x, y-v.y, z-v.z); }
8ee194
8ee194
    Real operator*(const Vector3 &v) const
8ee194
        { return x*v.x + y*v.y + z*v.z; }
8ee194
    Vector3 operator*(Real k) const
8ee194
        { return Vector3(x*k, y*k, z*k); }
e31ea0
    Vector3 operator/(Real k) const
e31ea0
        { return *this*(1/k); }
8ee194
8ee194
    Vector3 cross(const Vector3 &v) const
8ee194
        { return Vector3(y*v.z - z*v.y, z*v.x - x*v.z, x*v.y - y*v.x); }
e31ea0
    Vector3 perp() const
e31ea0
        { return fabs(y) > fabs(x) ? Vector3(0, z, -y) : Vector3(-z, 0, x); }
8ee194
8ee194
    Real len_sqr() const { return x*x + y*y + z*z; }
8ee194
    Real len() const { return sqrt(len_sqr()); }
e31ea0
e31ea0
    Vector3 norm() const
e31ea0
        { Real l = len(); return l > precision ? *this/l : Vector3(); }
d2b2b5
};
d2b2b5
d2b2b5
d2b2b5
class Vector4 {
d2b2b5
public:
d2b2b5
    union {
d2b2b5
        struct { Real x, y, z, w; };
d2b2b5
        struct { Real r, g, b, a; };
d2b2b5
        struct { Real c[4]; };
d2b2b5
    };
d2b2b5
    
d2b2b5
    explicit Vector4(Real x = 0, Real y = 0, Real z = 0, Real w = 0):
d2b2b5
        x(x), y(y), z(z), w(w) { }
d2b2b5
};
d2b2b5
d2b2b5
e31ea0
class Matrix4 {
e31ea0
public:
e31ea0
    union {
e31ea0
        struct {
e31ea0
            Real m00, m01, m02, m03,
e31ea0
                 m10, m11, m12, m13,
e31ea0
                 m20, m21, m22, m23,
e31ea0
                 m30, m31, m32, m33;
e31ea0
        };
e31ea0
        struct { Real m[4][4]; };
e31ea0
        struct { Real a[16]; };
e31ea0
    };
e31ea0
e31ea0
    explicit Matrix4(
e31ea0
        Real m00 = 1, Real m01 = 0, Real m02 = 0, Real m03 = 0,
e31ea0
        Real m10 = 0, Real m11 = 1, Real m12 = 0, Real m13 = 0,
e31ea0
        Real m20 = 0, Real m21 = 0, Real m22 = 1, Real m23 = 0,
e31ea0
        Real m30 = 0, Real m31 = 0, Real m32 = 0, Real m33 = 1
e31ea0
    ):
e31ea0
        m00(m00), m01(m01), m02(m02), m03(m03),
e31ea0
        m10(m10), m11(m11), m12(m12), m13(m13),
e31ea0
        m20(m20), m21(m21), m22(m22), m23(m23),
e31ea0
        m30(m30), m31(m31), m32(m32), m33(m33)
e31ea0
    { }
e31ea0
    
e31ea0
    Vector4 operator*(const Vector4 &v) const {
e31ea0
        return Vector4(
e31ea0
            m00*v.x + m10*v.y + m20*v.z + m30*v.w,
e31ea0
            m01*v.x + m11*v.y + m21*v.z + m31*v.w,
e31ea0
            m02*v.x + m12*v.y + m22*v.z + m32*v.w,
e31ea0
            m03*v.x + m13*v.y + m23*v.z + m33*v.w );
e31ea0
    }
e31ea0
e31ea0
    Vector3 transform(const Vector3 &v, Real w = 1) const {
e31ea0
        return Vector3(
e31ea0
            m00*v.x + m10*v.y + m20*v.z + m30*w,
e31ea0
            m01*v.x + m11*v.y + m21*v.z + m31*w,
e31ea0
            m02*v.x + m12*v.y + m22*v.z + m32*w );
e31ea0
    }
e31ea0
e31ea0
    Matrix4 operator*(const Matrix4 &m) const;
e31ea0
    
e31ea0
    Real det() const;
e31ea0
    Matrix4 inv() const;
e31ea0
    
e31ea0
    static Matrix4 zero();
e31ea0
    static Matrix4 translation(const Vector3 &translate);
e31ea0
    static Matrix4 scaling(const Vector3 &scale);
e31ea0
    static Matrix4 rotation(const Vector3 &axis, Real angle);
1bcd85
    static Matrix4 perspective(Real fovy, Real aspect, Real z_near, Real z_far);
e31ea0
};
e31ea0
d2b2b5
#endif