ImathMatrix.h

//
// SPDX-License-Identifier: BSD-3-Clause
// Copyright Contributors to the OpenEXR Project.
//
 
//
// 2x2, 3x3, and 4x4 transformation matrix templates
//
 
#ifndef INCLUDED_ATOMSMATHMATRIX_H
#define INCLUDED_ATOMSMATHMATRIX_H
 
#include <AtomsMath/ImathFun.h>
#include <AtomsMath/ImathNamespace.h>
#include <AtomsMath/ImathPlatform.h>
#include <AtomsMath/ImathShear.h>
#include <AtomsMath/ImathVec.h>
 
#include <cstring>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string.h>
 
#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
// suppress exception specification warnings
#    pragma warning(disable : 4290)
#endif
 
ATOMSMATH_INTERNAL_NAMESPACE_HEADER_ENTER
 
enum Uninitialized
{
    UNINITIALIZED
};
 
 
template <class T> class Matrix22
{
  public:
 
    
    T x[2][2];
 
    
    ATOMSMATH_HOSTDEVICE T* operator[] (int i) noexcept;
 
    ATOMSMATH_HOSTDEVICE const T* operator[] (int i) const noexcept;
 
 
    //ATOMSMATH_HOSTDEVICE Matrix22 (Uninitialized) noexcept {}
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22 (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22 (const T a[2][2]) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22 (T a, T b, T c, T d) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22 (const Matrix22& v) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 explicit Matrix22 (const Matrix22<S>& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator= (const Matrix22& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator= (T a) noexcept;
 
    ~Matrix22() noexcept = default;
 
 
 
    ATOMSMATH_HOSTDEVICE T* getValue() noexcept;
 
    ATOMSMATH_HOSTDEVICE const T* getValue() const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void getValue (Matrix22<S>& v) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22& setValue (const Matrix22<S>& v) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22& setTheMatrix (const Matrix22<S>& v) noexcept;
 
 
    
    ATOMSMATH_HOSTDEVICE constexpr bool operator== (const Matrix22& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr bool operator!= (const Matrix22& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 bool equalWithAbsError (const Matrix22<T>& v, T e) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 bool equalWithRelError (const Matrix22<T>& v, T e) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator+= (const Matrix22& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator+= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix22 operator+ (const Matrix22& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator-= (const Matrix22& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator-= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix22 operator- (const Matrix22& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix22 operator-() const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& negate() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator*= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix22 operator* (T a) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator/= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix22 operator/ (T a) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& operator*= (const Matrix22& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22 operator* (const Matrix22& v) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const noexcept;
 
 
 
    ATOMSMATH_HOSTDEVICE void makeIdentity() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& transpose() noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix22 transposed() const noexcept;
 
    ATOMSMATH_CONSTEXPR14 const Matrix22& invert (bool singExc);
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& invert() noexcept;
 
    ATOMSMATH_CONSTEXPR14 Matrix22<T> inverse (bool singExc) const;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix22<T> inverse() const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr T determinant() const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE const Matrix22& setRotation (S r) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& rotate (S r) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& setScale (T s) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& setScale (const Vec2<S>& s) noexcept;
 
    // Scale the matrix by s
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix22& scale (const Vec2<S>& s) noexcept;
 
    
    
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeLowest() noexcept { return std::numeric_limits<T>::lowest(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeMax() noexcept { return std::numeric_limits<T>::max(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeSmallest() noexcept { return std::numeric_limits<T>::min(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeEpsilon() noexcept { return std::numeric_limits<T>::epsilon(); }
 
    
    ATOMSMATH_HOSTDEVICE constexpr static unsigned int dimensions() noexcept { return 2; }
 
    typedef T BaseType;
 
    typedef Vec2<T> BaseVecType;
 
  private:
    template <typename R, typename S> struct isSameType
    {
        enum
        {
            value = 0
        };
    };
 
    template <typename R> struct isSameType<R, R>
    {
        enum
        {
            value = 1
        };
    };
};
 
 
template <class T> class Matrix33
{
  public:
 
    
    T x[3][3];
 
    
    ATOMSMATH_HOSTDEVICE T* operator[] (int i) noexcept;
 
    ATOMSMATH_HOSTDEVICE const T* operator[] (int i) const noexcept;
 
 
    //ATOMSMATH_HOSTDEVICE Matrix33 (Uninitialized) noexcept {}
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33 (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33 (const T a[3][3]) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33 (const Matrix33& v) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 explicit Matrix33 (const Matrix33<S>& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator= (const Matrix33& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator= (T a) noexcept;
 
    ~Matrix33() noexcept = default;
 
    
 
    ATOMSMATH_HOSTDEVICE T* getValue() noexcept;
 
    ATOMSMATH_HOSTDEVICE const T* getValue() const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void getValue (Matrix33<S>& v) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33& setValue (const Matrix33<S>& v) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33& setTheMatrix (const Matrix33<S>& v) noexcept;
 
    
    
    ATOMSMATH_HOSTDEVICE constexpr bool operator== (const Matrix33& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr bool operator!= (const Matrix33& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 bool equalWithAbsError (const Matrix33<T>& v, T e) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 bool equalWithRelError (const Matrix33<T>& v, T e) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator+= (const Matrix33& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator+= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix33 operator+ (const Matrix33& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator-= (const Matrix33& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator-= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix33 operator- (const Matrix33& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix33 operator-() const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& negate() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator*= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix33 operator* (T a) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator/= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix33 operator/ (T a) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& operator*= (const Matrix33& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33 operator* (const Matrix33& v) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void multVecMatrix (const Vec2<S>& src, Vec2<S>& dst) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const noexcept;
 
 
 
    ATOMSMATH_HOSTDEVICE void makeIdentity() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& transpose() noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix33 transposed() const noexcept;
 
    ATOMSMATH_CONSTEXPR14 const Matrix33& invert (bool singExc);
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& invert() noexcept;
 
    ATOMSMATH_CONSTEXPR14 Matrix33<T> inverse (bool singExc) const;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix33<T> inverse() const noexcept;
 
    const Matrix33& gjInvert (bool singExc);
    
    ATOMSMATH_HOSTDEVICE const Matrix33& gjInvert() noexcept;
 
    Matrix33<T> gjInverse (bool singExc) const;
 
    ATOMSMATH_HOSTDEVICE Matrix33<T> gjInverse() const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 T minorOf (const int r, const int c) const noexcept;
 
    ATOMSMATH_HOSTDEVICE
    constexpr T fastMinor (const int r0, const int r1, const int c0, const int c1) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr T determinant() const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE const Matrix33& setRotation (S r) noexcept;
 
    // Rotate the given matrix by r (in radians)
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& rotate (S r) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& setScale (T s) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& setScale (const Vec2<S>& s) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& scale (const Vec2<S>& s) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& setTranslation (const Vec2<S>& t) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Vec2<T> translation() const noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& translate (const Vec2<S>& t) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& setShear (const S& h) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& setShear (const Vec2<S>& h) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& shear (const S& xy) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix33& shear (const Vec2<S>& h) noexcept;
 
    
    
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeLowest() noexcept { return std::numeric_limits<T>::lowest(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeMax() noexcept { return std::numeric_limits<T>::max(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeSmallest() noexcept { return std::numeric_limits<T>::min(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeEpsilon() noexcept { return std::numeric_limits<T>::epsilon(); }
 
    
    ATOMSMATH_HOSTDEVICE constexpr static unsigned int dimensions() noexcept { return 3; }
 
    typedef T BaseType;
 
    typedef Vec3<T> BaseVecType;
 
  private:
    template <typename R, typename S> struct isSameType
    {
        enum
        {
            value = 0
        };
    };
 
    template <typename R> struct isSameType<R, R>
    {
        enum
        {
            value = 1
        };
    };
};
 
 
template <class T> class Matrix44
{
  public:
 
    
    T x[4][4];
 
    
    ATOMSMATH_HOSTDEVICE T* operator[] (int i) noexcept;
 
    ATOMSMATH_HOSTDEVICE const T* operator[] (int i) const noexcept;
 
 
    //ATOMSMATH_HOSTDEVICE constexpr Matrix44 (Uninitialized) noexcept {}
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44 (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44 (const T a[4][4]) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14
    Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, T i, T j, T k, T l, T m, T n, T o, T p) noexcept;
 
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44 (Matrix33<T> r, Vec3<T> t) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44 (const Matrix44& v) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 explicit Matrix44 (const Matrix44<S>& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator= (const Matrix44& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator= (T a) noexcept;
 
    ~Matrix44() noexcept = default;
 
    
 
    ATOMSMATH_HOSTDEVICE T* getValue() noexcept;
 
    ATOMSMATH_HOSTDEVICE const T* getValue() const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void getValue (Matrix44<S>& v) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44& setValue (const Matrix44<S>& v) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44& setTheMatrix (const Matrix44<S>& v) noexcept;
 
 
    
    ATOMSMATH_HOSTDEVICE constexpr bool operator== (const Matrix44& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr bool operator!= (const Matrix44& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 bool equalWithAbsError (const Matrix44<T>& v, T e) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 bool equalWithRelError (const Matrix44<T>& v, T e) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator+= (const Matrix44& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator+= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix44 operator+ (const Matrix44& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator-= (const Matrix44& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator-= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix44 operator- (const Matrix44& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix44 operator-() const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& negate() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator*= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix44 operator* (T a) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator/= (T a) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix44 operator/ (T a) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& operator*= (const Matrix44& v) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44 operator* (const Matrix44& v) const noexcept;
 
    ATOMSMATH_HOSTDEVICE
    static void multiply (const Matrix44& a,     // assumes that
                          const Matrix44& b,     // &a != &c and
                          Matrix44& c) noexcept; // &b != &c.
 
    template <class S> ATOMSMATH_HOSTDEVICE void multVecMatrix (const Vec3<S>& src, Vec3<S>& dst) const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE void multDirMatrix (const Vec3<S>& src, Vec3<S>& dst) const noexcept;
 
 
 
    ATOMSMATH_HOSTDEVICE void makeIdentity() noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& transpose() noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr Matrix44 transposed() const noexcept;
 
    ATOMSMATH_CONSTEXPR14 const Matrix44& invert (bool singExc);
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& invert() noexcept;
 
    ATOMSMATH_CONSTEXPR14 Matrix44<T> inverse (bool singExc) const;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 Matrix44<T> inverse() const noexcept;
 
    ATOMSMATH_CONSTEXPR14 const Matrix44& gjInvert (bool singExc);
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& gjInvert() noexcept;
 
    Matrix44<T> gjInverse (bool singExc) const;
 
    ATOMSMATH_HOSTDEVICE Matrix44<T> gjInverse() const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 T minorOf (const int r, const int c) const noexcept;
 
    ATOMSMATH_HOSTDEVICE
    constexpr T fastMinor (const int r0,
                           const int r1,
                           const int r2,
                           const int c0,
                           const int c1,
                           const int c2) const noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 T determinant() const noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE const Matrix44& setEulerAngles (const Vec3<S>& r) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& setAxisAngle (const Vec3<S>& ax, S ang) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE const Matrix44& rotate (const Vec3<S>& r) noexcept;
 
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& setScale (T s) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& setScale (const Vec3<S>& s) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& scale (const Vec3<S>& s) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& setTranslation (const Vec3<S>& t) noexcept;
 
    ATOMSMATH_HOSTDEVICE constexpr const Vec3<T> translation() const noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& translate (const Vec3<S>& t) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& setShear (const Vec3<S>& h) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& setShear (const Shear6<S>& h) noexcept;
 
    template <class S> ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& shear (const Vec3<S>& h) noexcept;
 
    template <class S>
    ATOMSMATH_HOSTDEVICE ATOMSMATH_CONSTEXPR14 const Matrix44& shear (const Shear6<S>& h) noexcept;
 
    
    
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeLowest() noexcept { return std::numeric_limits<T>::lowest(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeMax() noexcept { return std::numeric_limits<T>::max(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeSmallest() noexcept { return std::numeric_limits<T>::min(); }
 
    ATOMSMATH_HOSTDEVICE constexpr static T baseTypeEpsilon() noexcept { return std::numeric_limits<T>::epsilon(); }
 
    
    ATOMSMATH_HOSTDEVICE constexpr static unsigned int dimensions() noexcept { return 4; }
 
    typedef T BaseType;
 
    typedef Vec4<T> BaseVecType;
 
  private:
    template <typename R, typename S> struct isSameType
    {
        enum
        {
            value = 0
        };
    };
 
    template <typename R> struct isSameType<R, R>
    {
        enum
        {
            value = 1
        };
    };
};
 
template <class T> std::ostream& operator<< (std::ostream& s, const Matrix22<T>& m);
 
template <class T> std::ostream& operator<< (std::ostream& s, const Matrix33<T>& m);
 
template <class T> std::ostream& operator<< (std::ostream& s, const Matrix44<T>& m);
 
//---------------------------------------------
// Vector-times-matrix multiplication operators
//---------------------------------------------
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline const Vec2<S>& operator*= (Vec2<S>& v, const Matrix22<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline Vec2<S> operator* (const Vec2<S>& v, const Matrix22<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline const Vec2<S>& operator*= (Vec2<S>& v, const Matrix33<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline Vec2<S> operator* (const Vec2<S>& v, const Matrix33<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline const Vec3<S>& operator*= (Vec3<S>& v, const Matrix33<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline Vec3<S> operator* (const Vec3<S>& v, const Matrix33<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline const Vec3<S>& operator*= (Vec3<S>& v, const Matrix44<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline Vec3<S> operator* (const Vec3<S>& v, const Matrix44<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline const Vec4<S>& operator*= (Vec4<S>& v, const Matrix44<T>& m) noexcept;
 
template <class S, class T>
ATOMSMATH_HOSTDEVICE inline Vec4<S> operator* (const Vec4<S>& v, const Matrix44<T>& m) noexcept;
 
//-------------------------
// Typedefs for convenience
//-------------------------
 
typedef Matrix22<float> M22f;
 
typedef Matrix22<double> M22d;
 
typedef Matrix33<float> M33f;
 
typedef Matrix33<double> M33d;
 
typedef Matrix44<float> M44f;
 
typedef Matrix44<double> M44d;
 
//---------------------------
// Implementation of Matrix22
//---------------------------
 
template <class T>
inline T*
Matrix22<T>::operator[] (int i) noexcept
{
    return x[i];
}
 
template <class T>
inline const T*
Matrix22<T>::operator[] (int i) const noexcept
{
    return x[i];
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22() noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[1][0] = 0;
    x[1][1] = 1;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (T a) noexcept
{
    x[0][0] = a;
    x[0][1] = a;
    x[1][0] = a;
    x[1][1] = a;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (const T a[2][2]) noexcept
{
    // Function calls and aliasing issues can inhibit vectorization versus
    // straight assignment of data members, so instead of this:
    //     memcpy (x, a, sizeof (x));
    // we do this:
    x[0][0] = a[0][0];
    x[0][1] = a[0][1];
    x[1][0] = a[1][0];
    x[1][1] = a[1][1];
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (T a, T b, T c, T d) noexcept
{
    x[0][0] = a;
    x[0][1] = b;
    x[1][0] = c;
    x[1][1] = d;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (const Matrix22& v) noexcept
{
    // Function calls and aliasing issues can inhibit vectorization versus
    // straight assignment of data members, so we don't do this:
    //     memcpy (x, v.x, sizeof (x));
    // we do this:
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (const Matrix22<S>& v) noexcept
{
    x[0][0] = T (v.x[0][0]);
    x[0][1] = T (v.x[0][1]);
    x[1][0] = T (v.x[1][0]);
    x[1][1] = T (v.x[1][1]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator= (const Matrix22& v) noexcept
{
    // Function calls and aliasing issues can inhibit vectorization versus
    // straight assignment of data members, so we don't do this:
    //     memcpy (x, v.x, sizeof (x));
    // we do this:
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator= (T a) noexcept
{
    x[0][0] = a;
    x[0][1] = a;
    x[1][0] = a;
    x[1][1] = a;
    return *this;
}
 
template <class T>
inline T*
Matrix22<T>::getValue() noexcept
{
    return (T*) &x[0][0];
}
 
template <class T>
inline const T*
Matrix22<T>::getValue() const noexcept
{
    return (const T*) &x[0][0];
}
 
template <class T>
template <class S>
inline void
Matrix22<T>::getValue (Matrix22<S>& v) const noexcept
{
    v.x[0][0] = x[0][0];
    v.x[0][1] = x[0][1];
    v.x[1][0] = x[1][0];
    v.x[1][1] = x[1][1];
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix22<T>&
Matrix22<T>::setValue (const Matrix22<S>& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix22<T>&
Matrix22<T>::setTheMatrix (const Matrix22<S>& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    return *this;
}
 
template <class T>
inline void
Matrix22<T>::makeIdentity() noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[1][0] = 0;
    x[1][1] = 1;
}
 
template <class T>
constexpr inline bool
Matrix22<T>::operator== (const Matrix22& v) const noexcept
{
    return x[0][0] == v.x[0][0] && x[0][1] == v.x[0][1] && x[1][0] == v.x[1][0] &&
           x[1][1] == v.x[1][1];
}
 
template <class T>
constexpr inline bool
Matrix22<T>::operator!= (const Matrix22& v) const noexcept
{
    return x[0][0] != v.x[0][0] || x[0][1] != v.x[0][1] || x[1][0] != v.x[1][0] ||
           x[1][1] != v.x[1][1];
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline bool
Matrix22<T>::equalWithAbsError (const Matrix22<T>& m, T e) const noexcept
{
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            if (!ATOMSMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this)[i][j], m[i][j], e))
                return false;
 
    return true;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline bool
Matrix22<T>::equalWithRelError (const Matrix22<T>& m, T e) const noexcept
{
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            if (!ATOMSMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this)[i][j], m[i][j], e))
                return false;
 
    return true;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator+= (const Matrix22<T>& v) noexcept
{
    x[0][0] += v.x[0][0];
    x[0][1] += v.x[0][1];
    x[1][0] += v.x[1][0];
    x[1][1] += v.x[1][1];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator+= (T a) noexcept
{
    x[0][0] += a;
    x[0][1] += a;
    x[1][0] += a;
    x[1][1] += a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix22<T>
Matrix22<T>::operator+ (const Matrix22<T>& v) const noexcept
{
    return Matrix22 (x[0][0] + v.x[0][0],
                     x[0][1] + v.x[0][1],
                     x[1][0] + v.x[1][0],
                     x[1][1] + v.x[1][1]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator-= (const Matrix22<T>& v) noexcept
{
    x[0][0] -= v.x[0][0];
    x[0][1] -= v.x[0][1];
    x[1][0] -= v.x[1][0];
    x[1][1] -= v.x[1][1];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator-= (T a) noexcept
{
    x[0][0] -= a;
    x[0][1] -= a;
    x[1][0] -= a;
    x[1][1] -= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix22<T>
Matrix22<T>::operator- (const Matrix22<T>& v) const noexcept
{
    return Matrix22 (x[0][0] - v.x[0][0],
                     x[0][1] - v.x[0][1],
                     x[1][0] - v.x[1][0],
                     x[1][1] - v.x[1][1]);
}
 
template <class T>
constexpr inline Matrix22<T>
Matrix22<T>::operator-() const noexcept
{
    return Matrix22 (-x[0][0], -x[0][1], -x[1][0], -x[1][1]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::negate() noexcept
{
    x[0][0] = -x[0][0];
    x[0][1] = -x[0][1];
    x[1][0] = -x[1][0];
    x[1][1] = -x[1][1];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator*= (T a) noexcept
{
    x[0][0] *= a;
    x[0][1] *= a;
    x[1][0] *= a;
    x[1][1] *= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix22<T>
Matrix22<T>::operator* (T a) const noexcept
{
    return Matrix22 (x[0][0] * a, x[0][1] * a, x[1][0] * a, x[1][1] * a);
}
 
template <class T>
inline Matrix22<T>
operator* (T a, const Matrix22<T>& v) noexcept
{
    return v * a;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator*= (const Matrix22<T>& v) noexcept
{
    Matrix22 tmp (T (0));
 
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            for (int k = 0; k < 2; k++)
                tmp.x[i][j] += x[i][k] * v.x[k][j];
 
    *this = tmp;
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix22<T>
Matrix22<T>::operator* (const Matrix22<T>& v) const noexcept
{
    Matrix22 tmp (T (0));
 
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            for (int k = 0; k < 2; k++)
                tmp.x[i][j] += x[i][k] * v.x[k][j];
 
    return tmp;
}
 
template <class T>
template <class S>
inline void
Matrix22<T>::multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const noexcept
{
    S a, b;
 
    a = src.x * x[0][0] + src.y * x[1][0];
    b = src.x * x[0][1] + src.y * x[1][1];
 
    dst.x = a;
    dst.y = b;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::operator/= (T a) noexcept
{
    x[0][0] /= a;
    x[0][1] /= a;
    x[1][0] /= a;
    x[1][1] /= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix22<T>
Matrix22<T>::operator/ (T a) const noexcept
{
    return Matrix22 (x[0][0] / a, x[0][1] / a, x[1][0] / a, x[1][1] / a);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::transpose() noexcept
{
    Matrix22 tmp (x[0][0], x[1][0], x[0][1], x[1][1]);
    *this = tmp;
    return *this;
}
 
template <class T>
constexpr inline Matrix22<T>
Matrix22<T>::transposed() const noexcept
{
    return Matrix22 (x[0][0], x[1][0], x[0][1], x[1][1]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::invert (bool singExc)
{
    *this = inverse (singExc);
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::invert() noexcept
{
    *this = inverse();
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix22<T>
Matrix22<T>::inverse (bool singExc) const
{
    Matrix22 s (x[1][1], -x[0][1], -x[1][0], x[0][0]);
 
    T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
 
    if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
    {
        for (int i = 0; i < 2; ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                s[i][j] /= r;
            }
        }
    }
    else
    {
        T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
        for (int i = 0; i < 2; ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                {
                    s[i][j] /= r;
                }
                else
                {
                    #ifdef ATOMS_USE_EXCEPTIONS
                    if (singExc)
                        throw std::invalid_argument ("Cannot invert "
                                                     "singular matrix.");
                    #endif
                    return Matrix22();
                }
            }
        }
    }
    return s;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix22<T>
Matrix22<T>::inverse() const noexcept
{
    Matrix22 s (x[1][1], -x[0][1], -x[1][0], x[0][0]);
 
    T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
 
    if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
    {
        for (int i = 0; i < 2; ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                s[i][j] /= r;
            }
        }
    }
    else
    {
        T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
        for (int i = 0; i < 2; ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                {
                    s[i][j] /= r;
                }
                else
                {
                    return Matrix22();
                }
            }
        }
    }
    return s;
}
 
template <class T>
constexpr inline T
Matrix22<T>::determinant() const noexcept
{
    return x[0][0] * x[1][1] - x[1][0] * x[0][1];
}
 
template <class T>
template <class S>
inline const Matrix22<T>&
Matrix22<T>::setRotation (S r) noexcept
{
    S cos_r, sin_r;
 
    cos_r = cos ((T) r);
    sin_r = sin ((T) r);
 
    x[0][0] = cos_r;
    x[0][1] = sin_r;
 
    x[1][0] = -sin_r;
    x[1][1] = cos_r;
 
    return *this;
}
 
template <class T>
template <class S>
 
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::rotate (S r) noexcept
{
    *this *= Matrix22<T>().setRotation (r);
    return *this;
}
 
template <class T>
 
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::setScale (T s) noexcept
{
    //
    // Set the matrix to:
    //  | s 0 |
    //  | 0 s |
    //
 
    x[0][0] = s;
    x[0][1] = static_cast<T> (0);
    x[1][0] = static_cast<T> (0);
    x[1][1] = s;
 
    return *this;
}
 
template <class T>
template <class S>
 
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::setScale (const Vec2<S>& s) noexcept
{
    //
    // Set the matrix to:
    //  | s.x  0  |
    //  |  0  s.y |
    //
 
    x[0][0] = s.x;
    x[0][1] = static_cast<T> (0);
    x[1][0] = static_cast<T> (0);
    x[1][1] = s.y;
 
    return *this;
}
 
template <class T>
template <class S>
 
ATOMSMATH_CONSTEXPR14 inline const Matrix22<T>&
Matrix22<T>::scale (const Vec2<S>& s) noexcept
{
    x[0][0] *= s.x;
    x[0][1] *= s.x;
 
    x[1][0] *= s.y;
    x[1][1] *= s.y;
 
    return *this;
}
 
//---------------------------
// Implementation of Matrix33
//---------------------------
 
template <class T>
inline T*
Matrix33<T>::operator[] (int i) noexcept
{
    return x[i];
}
 
template <class T>
inline const T*
Matrix33<T>::operator[] (int i) const noexcept
{
    return x[i];
}
 
template <class T>
inline ATOMSMATH_CONSTEXPR14
Matrix33<T>::Matrix33() noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (T a) noexcept
{
    x[0][0] = a;
    x[0][1] = a;
    x[0][2] = a;
    x[1][0] = a;
    x[1][1] = a;
    x[1][2] = a;
    x[2][0] = a;
    x[2][1] = a;
    x[2][2] = a;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (const T a[3][3]) noexcept
{
    // Function calls and aliasing issues can inhibit vectorization versus
    // straight assignment of data members, so instead of this:
    //     memcpy (x, a, sizeof (x));
    // we do this:
    x[0][0] = a[0][0];
    x[0][1] = a[0][1];
    x[0][2] = a[0][2];
    x[1][0] = a[1][0];
    x[1][1] = a[1][1];
    x[1][2] = a[1][2];
    x[2][0] = a[2][0];
    x[2][1] = a[2][1];
    x[2][2] = a[2][2];
}
 
template <class T>
 
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i) noexcept
{
    x[0][0] = a;
    x[0][1] = b;
    x[0][2] = c;
    x[1][0] = d;
    x[1][1] = e;
    x[1][2] = f;
    x[2][0] = g;
    x[2][1] = h;
    x[2][2] = i;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (const Matrix33& v) noexcept
{
    // Function calls and aliasing issues can inhibit vectorization versus
    // straight assignment of data members, so instead of this:
    //     memcpy (x, v.x, sizeof (x));
    // we do this:
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
}
 
template <class T>
template <class S>
 
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (const Matrix33<S>& v) noexcept
{
    x[0][0] = T (v.x[0][0]);
    x[0][1] = T (v.x[0][1]);
    x[0][2] = T (v.x[0][2]);
    x[1][0] = T (v.x[1][0]);
    x[1][1] = T (v.x[1][1]);
    x[1][2] = T (v.x[1][2]);
    x[2][0] = T (v.x[2][0]);
    x[2][1] = T (v.x[2][1]);
    x[2][2] = T (v.x[2][2]);
}
 
template <class T>
 
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator= (const Matrix33& v) noexcept
{
    // Function calls and aliasing issues can inhibit vectorization versus
    // straight assignment of data members, so instead of this:
    //     memcpy (x, v.x, sizeof (x));
    // we do this:
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    return *this;
}
 
template <class T>
 
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator= (T a) noexcept
{
    x[0][0] = a;
    x[0][1] = a;
    x[0][2] = a;
    x[1][0] = a;
    x[1][1] = a;
    x[1][2] = a;
    x[2][0] = a;
    x[2][1] = a;
    x[2][2] = a;
    return *this;
}
 
template <class T>
inline T*
Matrix33<T>::getValue() noexcept
{
    return (T*) &x[0][0];
}
 
template <class T>
inline const T*
Matrix33<T>::getValue() const noexcept
{
    return (const T*) &x[0][0];
}
 
template <class T>
template <class S>
inline void
Matrix33<T>::getValue (Matrix33<S>& v) const noexcept
{
    v.x[0][0] = x[0][0];
    v.x[0][1] = x[0][1];
    v.x[0][2] = x[0][2];
    v.x[1][0] = x[1][0];
    v.x[1][1] = x[1][1];
    v.x[1][2] = x[1][2];
    v.x[2][0] = x[2][0];
    v.x[2][1] = x[2][1];
    v.x[2][2] = x[2][2];
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>&
Matrix33<T>::setValue (const Matrix33<S>& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>&
Matrix33<T>::setTheMatrix (const Matrix33<S>& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    return *this;
}
 
template <class T>
inline void
Matrix33<T>::makeIdentity() noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
}
 
template <class T>
constexpr inline bool
Matrix33<T>::operator== (const Matrix33& v) const noexcept
{
    return x[0][0] == v.x[0][0] && x[0][1] == v.x[0][1] && x[0][2] == v.x[0][2] &&
           x[1][0] == v.x[1][0] && x[1][1] == v.x[1][1] && x[1][2] == v.x[1][2] &&
           x[2][0] == v.x[2][0] && x[2][1] == v.x[2][1] && x[2][2] == v.x[2][2];
}
 
template <class T>
constexpr inline bool
Matrix33<T>::operator!= (const Matrix33& v) const noexcept
{
    return x[0][0] != v.x[0][0] || x[0][1] != v.x[0][1] || x[0][2] != v.x[0][2] ||
           x[1][0] != v.x[1][0] || x[1][1] != v.x[1][1] || x[1][2] != v.x[1][2] ||
           x[2][0] != v.x[2][0] || x[2][1] != v.x[2][1] || x[2][2] != v.x[2][2];
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline bool
Matrix33<T>::equalWithAbsError (const Matrix33<T>& m, T e) const noexcept
{
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            if (!ATOMSMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this)[i][j], m[i][j], e))
                return false;
 
    return true;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline bool
Matrix33<T>::equalWithRelError (const Matrix33<T>& m, T e) const noexcept
{
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            if (!ATOMSMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this)[i][j], m[i][j], e))
                return false;
 
    return true;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator+= (const Matrix33<T>& v) noexcept
{
    x[0][0] += v.x[0][0];
    x[0][1] += v.x[0][1];
    x[0][2] += v.x[0][2];
    x[1][0] += v.x[1][0];
    x[1][1] += v.x[1][1];
    x[1][2] += v.x[1][2];
    x[2][0] += v.x[2][0];
    x[2][1] += v.x[2][1];
    x[2][2] += v.x[2][2];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator+= (T a) noexcept
{
    x[0][0] += a;
    x[0][1] += a;
    x[0][2] += a;
    x[1][0] += a;
    x[1][1] += a;
    x[1][2] += a;
    x[2][0] += a;
    x[2][1] += a;
    x[2][2] += a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix33<T>
Matrix33<T>::operator+ (const Matrix33<T>& v) const noexcept
{
    return Matrix33 (x[0][0] + v.x[0][0],
                     x[0][1] + v.x[0][1],
                     x[0][2] + v.x[0][2],
                     x[1][0] + v.x[1][0],
                     x[1][1] + v.x[1][1],
                     x[1][2] + v.x[1][2],
                     x[2][0] + v.x[2][0],
                     x[2][1] + v.x[2][1],
                     x[2][2] + v.x[2][2]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator-= (const Matrix33<T>& v) noexcept
{
    x[0][0] -= v.x[0][0];
    x[0][1] -= v.x[0][1];
    x[0][2] -= v.x[0][2];
    x[1][0] -= v.x[1][0];
    x[1][1] -= v.x[1][1];
    x[1][2] -= v.x[1][2];
    x[2][0] -= v.x[2][0];
    x[2][1] -= v.x[2][1];
    x[2][2] -= v.x[2][2];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator-= (T a) noexcept
{
    x[0][0] -= a;
    x[0][1] -= a;
    x[0][2] -= a;
    x[1][0] -= a;
    x[1][1] -= a;
    x[1][2] -= a;
    x[2][0] -= a;
    x[2][1] -= a;
    x[2][2] -= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix33<T>
Matrix33<T>::operator- (const Matrix33<T>& v) const noexcept
{
    return Matrix33 (x[0][0] - v.x[0][0],
                     x[0][1] - v.x[0][1],
                     x[0][2] - v.x[0][2],
                     x[1][0] - v.x[1][0],
                     x[1][1] - v.x[1][1],
                     x[1][2] - v.x[1][2],
                     x[2][0] - v.x[2][0],
                     x[2][1] - v.x[2][1],
                     x[2][2] - v.x[2][2]);
}
 
template <class T>
constexpr inline Matrix33<T>
Matrix33<T>::operator-() const noexcept
{
    return Matrix33 (-x[0][0],
                     -x[0][1],
                     -x[0][2],
                     -x[1][0],
                     -x[1][1],
                     -x[1][2],
                     -x[2][0],
                     -x[2][1],
                     -x[2][2]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::negate() noexcept
{
    x[0][0] = -x[0][0];
    x[0][1] = -x[0][1];
    x[0][2] = -x[0][2];
    x[1][0] = -x[1][0];
    x[1][1] = -x[1][1];
    x[1][2] = -x[1][2];
    x[2][0] = -x[2][0];
    x[2][1] = -x[2][1];
    x[2][2] = -x[2][2];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator*= (T a) noexcept
{
    x[0][0] *= a;
    x[0][1] *= a;
    x[0][2] *= a;
    x[1][0] *= a;
    x[1][1] *= a;
    x[1][2] *= a;
    x[2][0] *= a;
    x[2][1] *= a;
    x[2][2] *= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix33<T>
Matrix33<T>::operator* (T a) const noexcept
{
    return Matrix33 (x[0][0] * a,
                     x[0][1] * a,
                     x[0][2] * a,
                     x[1][0] * a,
                     x[1][1] * a,
                     x[1][2] * a,
                     x[2][0] * a,
                     x[2][1] * a,
                     x[2][2] * a);
}
 
template <class T>
inline Matrix33<T> constexpr
operator* (T a, const Matrix33<T>& v) noexcept
{
    return v * a;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator*= (const Matrix33<T>& v) noexcept
{
    // Avoid initializing with 0 values before immediately overwriting them,
    // and unroll all loops for the best autovectorization.
    Matrix33 tmp (ATOMSMATH_NAMESPACE::UNINITIALIZED);
 
    tmp.x[0][0] = x[0][0] * v.x[0][0] + x[0][1] * v.x[1][0] + x[0][2] * v.x[2][0];
    tmp.x[0][1] = x[0][0] * v.x[0][1] + x[0][1] * v.x[1][1] + x[0][2] * v.x[2][1];
    tmp.x[0][2] = x[0][0] * v.x[0][2] + x[0][1] * v.x[1][2] + x[0][2] * v.x[2][2];
 
    tmp.x[1][0] = x[1][0] * v.x[0][0] + x[1][1] * v.x[1][0] + x[1][2] * v.x[2][0];
    tmp.x[1][1] = x[1][0] * v.x[0][1] + x[1][1] * v.x[1][1] + x[1][2] * v.x[2][1];
    tmp.x[1][2] = x[1][0] * v.x[0][2] + x[1][1] * v.x[1][2] + x[1][2] * v.x[2][2];
 
    tmp.x[2][0] = x[2][0] * v.x[0][0] + x[2][1] * v.x[1][0] + x[2][2] * v.x[2][0];
    tmp.x[2][1] = x[2][0] * v.x[0][1] + x[2][1] * v.x[1][1] + x[2][2] * v.x[2][1];
    tmp.x[2][2] = x[2][0] * v.x[0][2] + x[2][1] * v.x[1][2] + x[2][2] * v.x[2][2];
 
    *this = tmp;
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>
Matrix33<T>::operator* (const Matrix33<T>& v) const noexcept
{
    // Avoid initializing with 0 values before immediately overwriting them,
    // and unroll all loops for the best autovectorization.
    Matrix33 tmp (ATOMSMATH_NAMESPACE::UNINITIALIZED);
 
    tmp.x[0][0] = x[0][0] * v.x[0][0] + x[0][1] * v.x[1][0] + x[0][2] * v.x[2][0];
    tmp.x[0][1] = x[0][0] * v.x[0][1] + x[0][1] * v.x[1][1] + x[0][2] * v.x[2][1];
    tmp.x[0][2] = x[0][0] * v.x[0][2] + x[0][1] * v.x[1][2] + x[0][2] * v.x[2][2];
 
    tmp.x[1][0] = x[1][0] * v.x[0][0] + x[1][1] * v.x[1][0] + x[1][2] * v.x[2][0];
    tmp.x[1][1] = x[1][0] * v.x[0][1] + x[1][1] * v.x[1][1] + x[1][2] * v.x[2][1];
    tmp.x[1][2] = x[1][0] * v.x[0][2] + x[1][1] * v.x[1][2] + x[1][2] * v.x[2][2];
 
    tmp.x[2][0] = x[2][0] * v.x[0][0] + x[2][1] * v.x[1][0] + x[2][2] * v.x[2][0];
    tmp.x[2][1] = x[2][0] * v.x[0][1] + x[2][1] * v.x[1][1] + x[2][2] * v.x[2][1];
    tmp.x[2][2] = x[2][0] * v.x[0][2] + x[2][1] * v.x[1][2] + x[2][2] * v.x[2][2];
 
    return tmp;
}
 
template <class T>
template <class S>
inline void
Matrix33<T>::multVecMatrix (const Vec2<S>& src, Vec2<S>& dst) const noexcept
{
    S a, b, w;
 
    a = src.x * x[0][0] + src.y * x[1][0] + x[2][0];
    b = src.x * x[0][1] + src.y * x[1][1] + x[2][1];
    w = src.x * x[0][2] + src.y * x[1][2] + x[2][2];
 
    dst.x = a / w;
    dst.y = b / w;
}
 
template <class T>
template <class S>
inline void
Matrix33<T>::multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const noexcept
{
    S a, b;
 
    a = src.x * x[0][0] + src.y * x[1][0];
    b = src.x * x[0][1] + src.y * x[1][1];
 
    dst.x = a;
    dst.y = b;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::operator/= (T a) noexcept
{
    x[0][0] /= a;
    x[0][1] /= a;
    x[0][2] /= a;
    x[1][0] /= a;
    x[1][1] /= a;
    x[1][2] /= a;
    x[2][0] /= a;
    x[2][1] /= a;
    x[2][2] /= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix33<T>
Matrix33<T>::operator/ (T a) const noexcept
{
    return Matrix33 (x[0][0] / a,
                     x[0][1] / a,
                     x[0][2] / a,
                     x[1][0] / a,
                     x[1][1] / a,
                     x[1][2] / a,
                     x[2][0] / a,
                     x[2][1] / a,
                     x[2][2] / a);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::transpose() noexcept
{
    Matrix33 tmp (x[0][0], x[1][0], x[2][0], x[0][1], x[1][1], x[2][1], x[0][2], x[1][2], x[2][2]);
    *this = tmp;
    return *this;
}
 
template <class T>
constexpr inline Matrix33<T>
Matrix33<T>::transposed() const noexcept
{
    return Matrix33 (x[0][0],
                     x[1][0],
                     x[2][0],
                     x[0][1],
                     x[1][1],
                     x[2][1],
                     x[0][2],
                     x[1][2],
                     x[2][2]);
}
 
template <class T>
const inline Matrix33<T>&
Matrix33<T>::gjInvert (bool singExc)
{
    *this = gjInverse (singExc);
    return *this;
}
 
template <class T>
const inline Matrix33<T>&
Matrix33<T>::gjInvert() noexcept
{
    *this = gjInverse();
    return *this;
}
 
template <class T>
inline Matrix33<T>
Matrix33<T>::gjInverse (bool singExc) const
{
    int i, j, k;
    Matrix33 s;
    Matrix33 t (*this);
 
    // Forward elimination
 
    for (i = 0; i < 2; i++)
    {
        int pivot = i;
 
        T pivotsize = t[i][i];
 
        if (pivotsize < 0)
            pivotsize = -pivotsize;
 
        for (j = i + 1; j < 3; j++)
        {
            T tmp = t[j][i];
 
            if (tmp < 0)
                tmp = -tmp;
 
            if (tmp > pivotsize)
            {
                pivot     = j;
                pivotsize = tmp;
            }
        }
 
        if (pivotsize == 0)
        {
            #ifdef ATOMS_USE_EXCEPTIONS
            if (singExc)
                throw std::invalid_argument ("Cannot invert singular matrix.");
            #endif
            return Matrix33();
        }
 
        if (pivot != i)
        {
            for (j = 0; j < 3; j++)
            {
                T tmp;
 
                tmp         = t[i][j];
                t[i][j]     = t[pivot][j];
                t[pivot][j] = tmp;
 
                tmp         = s[i][j];
                s[i][j]     = s[pivot][j];
                s[pivot][j] = tmp;
            }
        }
 
        for (j = i + 1; j < 3; j++)
        {
            T f = t[j][i] / t[i][i];
 
            for (k = 0; k < 3; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    // Backward substitution
 
    for (i = 2; i >= 0; --i)
    {
        T f;
 
        if ((f = t[i][i]) == 0)
        {
            #ifdef ATOMS_USE_EXCEPTIONS
            if (singExc)
                throw std::invalid_argument ("Cannot invert singular matrix.");
            #endif
            return Matrix33();
        }
 
        for (j = 0; j < 3; j++)
        {
            t[i][j] /= f;
            s[i][j] /= f;
        }
 
        for (j = 0; j < i; j++)
        {
            f = t[j][i];
 
            for (k = 0; k < 3; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    return s;
}
 
template <class T>
inline Matrix33<T>
Matrix33<T>::gjInverse() const noexcept
{
    int i, j, k;
    Matrix33 s;
    Matrix33 t (*this);
 
    // Forward elimination
 
    for (i = 0; i < 2; i++)
    {
        int pivot = i;
 
        T pivotsize = t[i][i];
 
        if (pivotsize < 0)
            pivotsize = -pivotsize;
 
        for (j = i + 1; j < 3; j++)
        {
            T tmp = t[j][i];
 
            if (tmp < 0)
                tmp = -tmp;
 
            if (tmp > pivotsize)
            {
                pivot     = j;
                pivotsize = tmp;
            }
        }
 
        if (pivotsize == 0)
        {
            return Matrix33();
        }
 
        if (pivot != i)
        {
            for (j = 0; j < 3; j++)
            {
                T tmp;
 
                tmp         = t[i][j];
                t[i][j]     = t[pivot][j];
                t[pivot][j] = tmp;
 
                tmp         = s[i][j];
                s[i][j]     = s[pivot][j];
                s[pivot][j] = tmp;
            }
        }
 
        for (j = i + 1; j < 3; j++)
        {
            T f = t[j][i] / t[i][i];
 
            for (k = 0; k < 3; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    // Backward substitution
 
    for (i = 2; i >= 0; --i)
    {
        T f;
 
        if ((f = t[i][i]) == 0)
        {
            return Matrix33();
        }
 
        for (j = 0; j < 3; j++)
        {
            t[i][j] /= f;
            s[i][j] /= f;
        }
 
        for (j = 0; j < i; j++)
        {
            f = t[j][i];
 
            for (k = 0; k < 3; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    return s;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::invert (bool singExc)
{
    *this = inverse (singExc);
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::invert() noexcept
{
    *this = inverse();
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>
Matrix33<T>::inverse (bool singExc) const
{
    if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
    {
        Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
                    x[2][1] * x[0][2] - x[0][1] * x[2][2],
                    x[0][1] * x[1][2] - x[1][1] * x[0][2],
 
                    x[2][0] * x[1][2] - x[1][0] * x[2][2],
                    x[0][0] * x[2][2] - x[2][0] * x[0][2],
                    x[1][0] * x[0][2] - x[0][0] * x[1][2],
 
                    x[1][0] * x[2][1] - x[2][0] * x[1][1],
                    x[2][0] * x[0][1] - x[0][0] * x[2][1],
                    x[0][0] * x[1][1] - x[1][0] * x[0][1]);
 
        T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
 
        if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
        {
            for (int i = 0; i < 3; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    s[i][j] /= r;
                }
            }
        }
        else
        {
            T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
            for (int i = 0; i < 3; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                    {
                        s[i][j] /= r;
                    }
                    else
                    {
                        #ifdef ATOMS_USE_EXCEPTIONS
                        if (singExc)
                            throw std::invalid_argument ("Cannot invert "
                                                         "singular matrix.");
                        #endif
                        return Matrix33();
                    }
                }
            }
        }
 
        return s;
    }
    else
    {
        Matrix33 s (x[1][1],
                    -x[0][1],
                    0,
 
                    -x[1][0],
                    x[0][0],
                    0,
 
                    0,
                    0,
                    1);
 
        T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
 
        if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
        {
            for (int i = 0; i < 2; ++i)
            {
                for (int j = 0; j < 2; ++j)
                {
                    s[i][j] /= r;
                }
            }
        }
        else
        {
            T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
            for (int i = 0; i < 2; ++i)
            {
                for (int j = 0; j < 2; ++j)
                {
                    if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                    {
                        s[i][j] /= r;
                    }
                    else
                    {
                        #ifdef ATOMS_USE_EXCEPTIONS
                        if (singExc)
                            throw std::invalid_argument ("Cannot invert "
                                                         "singular matrix.");
                        #endif
                        return Matrix33();
                    }
                }
            }
        }
 
        s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0];
        s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1];
 
        return s;
    }
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix33<T>
Matrix33<T>::inverse() const noexcept
{
    if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
    {
        Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
                    x[2][1] * x[0][2] - x[0][1] * x[2][2],
                    x[0][1] * x[1][2] - x[1][1] * x[0][2],
 
                    x[2][0] * x[1][2] - x[1][0] * x[2][2],
                    x[0][0] * x[2][2] - x[2][0] * x[0][2],
                    x[1][0] * x[0][2] - x[0][0] * x[1][2],
 
                    x[1][0] * x[2][1] - x[2][0] * x[1][1],
                    x[2][0] * x[0][1] - x[0][0] * x[2][1],
                    x[0][0] * x[1][1] - x[1][0] * x[0][1]);
 
        T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
 
        if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
        {
            for (int i = 0; i < 3; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    s[i][j] /= r;
                }
            }
        }
        else
        {
            T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
            for (int i = 0; i < 3; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                    {
                        s[i][j] /= r;
                    }
                    else
                    {
                        return Matrix33();
                    }
                }
            }
        }
 
        return s;
    }
    else
    {
        Matrix33 s (x[1][1],
                    -x[0][1],
                    0,
 
                    -x[1][0],
                    x[0][0],
                    0,
 
                    0,
                    0,
                    1);
 
        T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
 
        if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
        {
            for (int i = 0; i < 2; ++i)
            {
                for (int j = 0; j < 2; ++j)
                {
                    s[i][j] /= r;
                }
            }
        }
        else
        {
            T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
            for (int i = 0; i < 2; ++i)
            {
                for (int j = 0; j < 2; ++j)
                {
                    if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                    {
                        s[i][j] /= r;
                    }
                    else
                    {
                        return Matrix33();
                    }
                }
            }
        }
 
        s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0];
        s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1];
 
        return s;
    }
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline T
Matrix33<T>::minorOf (const int r, const int c) const noexcept
{
    int r0 = 0 + (r < 1 ? 1 : 0);
    int r1 = 1 + (r < 2 ? 1 : 0);
    int c0 = 0 + (c < 1 ? 1 : 0);
    int c1 = 1 + (c < 2 ? 1 : 0);
 
    return x[r0][c0] * x[r1][c1] - x[r1][c0] * x[r0][c1];
}
 
template <class T>
constexpr inline T
Matrix33<T>::fastMinor (const int r0, const int r1, const int c0, const int c1) const noexcept
{
    return x[r0][c0] * x[r1][c1] - x[r0][c1] * x[r1][c0];
}
 
template <class T>
constexpr inline T
Matrix33<T>::determinant() const noexcept
{
    return x[0][0] * (x[1][1] * x[2][2] - x[1][2] * x[2][1]) +
           x[0][1] * (x[1][2] * x[2][0] - x[1][0] * x[2][2]) +
           x[0][2] * (x[1][0] * x[2][1] - x[1][1] * x[2][0]);
}
 
template <class T>
template <class S>
inline const Matrix33<T>&
Matrix33<T>::setRotation (S r) noexcept
{
    S cos_r, sin_r;
 
    cos_r = cos ((T) r);
    sin_r = sin ((T) r);
 
    x[0][0] = cos_r;
    x[0][1] = sin_r;
    x[0][2] = 0;
 
    x[1][0] = -sin_r;
    x[1][1] = cos_r;
    x[1][2] = 0;
 
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::rotate (S r) noexcept
{
    *this *= Matrix33<T>().setRotation (r);
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::setScale (T s) noexcept
{
    //
    // Set the matrix to a 2D homogeneous transform scale:
    //  | s 0 0 |
    //  | 0 s 0 |
    //  | 0 0 1 |
    //
 
    x[0][0] = s;
    x[0][1] = 0;
    x[0][2] = 0;
    x[1][0] = 0;
    x[1][1] = s;
    x[1][2] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::setScale (const Vec2<S>& s) noexcept
{
    //
    // Set the matrix to a 2D homogeneous transform scale:
    //  | s.x  0   0 |
    //  |  0  s.y  0 |
    //  |  0   0   1 |
    //
 
    x[0][0] = s.x;
    x[0][1] = 0;
    x[0][2] = 0;
    x[1][0] = 0;
    x[1][1] = s.y;
    x[1][2] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::scale (const Vec2<S>& s) noexcept
{
    x[0][0] *= s.x;
    x[0][1] *= s.x;
    x[0][2] *= s.x;
 
    x[1][0] *= s.y;
    x[1][1] *= s.y;
    x[1][2] *= s.y;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::setTranslation (const Vec2<S>& t) noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
 
    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
 
    x[2][0] = t.x;
    x[2][1] = t.y;
    x[2][2] = 1;
 
    return *this;
}
 
template <class T>
constexpr inline Vec2<T>
Matrix33<T>::translation() const noexcept
{
    return Vec2<T> (x[2][0], x[2][1]);
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::translate (const Vec2<S>& t) noexcept
{
    x[2][0] += t.x * x[0][0] + t.y * x[1][0];
    x[2][1] += t.x * x[0][1] + t.y * x[1][1];
    x[2][2] += t.x * x[0][2] + t.y * x[1][2];
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::setShear (const S& xy) noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
 
    x[1][0] = xy;
    x[1][1] = 1;
    x[1][2] = 0;
 
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::setShear (const Vec2<S>& h) noexcept
{
    x[0][0] = 1;
    x[0][1] = h.y;
    x[0][2] = 0;
 
    x[1][0] = h.x;
    x[1][1] = 1;
    x[1][2] = 0;
 
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::shear (const S& xy) noexcept
{
    //
    // In this case, we don't need a temp. copy of the matrix
    // because we never use a value on the RHS after we've
    // changed it on the LHS.
    //
 
    x[1][0] += xy * x[0][0];
    x[1][1] += xy * x[0][1];
    x[1][2] += xy * x[0][2];
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix33<T>&
Matrix33<T>::shear (const Vec2<S>& h) noexcept
{
    Matrix33<T> P (*this);
 
    x[0][0] = P[0][0] + h.y * P[1][0];
    x[0][1] = P[0][1] + h.y * P[1][1];
    x[0][2] = P[0][2] + h.y * P[1][2];
 
    x[1][0] = P[1][0] + h.x * P[0][0];
    x[1][1] = P[1][1] + h.x * P[0][1];
    x[1][2] = P[1][2] + h.x * P[0][2];
 
    return *this;
}
 
//---------------------------
// Implementation of Matrix44
//---------------------------
 
template <class T>
inline T*
Matrix44<T>::operator[] (int i) noexcept
{
    return x[i];
}
 
template <class T>
inline const T*
Matrix44<T>::operator[] (int i) const noexcept
{
    return x[i];
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44() noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;
    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
    x[1][3] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
    x[2][3] = 0;
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (T a) noexcept
{
    x[0][0] = a;
    x[0][1] = a;
    x[0][2] = a;
    x[0][3] = a;
    x[1][0] = a;
    x[1][1] = a;
    x[1][2] = a;
    x[1][3] = a;
    x[2][0] = a;
    x[2][1] = a;
    x[2][2] = a;
    x[2][3] = a;
    x[3][0] = a;
    x[3][1] = a;
    x[3][2] = a;
    x[3][3] = a;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (const T a[4][4]) noexcept
{
    x[0][0] = a[0][0];
    x[0][1] = a[0][1];
    x[0][2] = a[0][2];
    x[0][3] = a[0][3];
    x[1][0] = a[1][0];
    x[1][1] = a[1][1];
    x[1][2] = a[1][2];
    x[1][3] = a[1][3];
    x[2][0] = a[2][0];
    x[2][1] = a[2][1];
    x[2][2] = a[2][2];
    x[2][3] = a[2][3];
    x[3][0] = a[3][0];
    x[3][1] = a[3][1];
    x[3][2] = a[3][2];
    x[3][3] = a[3][3];
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix44<
    T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, T i, T j, T k, T l, T m, T n, T o, T p) noexcept
{
    x[0][0] = a;
    x[0][1] = b;
    x[0][2] = c;
    x[0][3] = d;
    x[1][0] = e;
    x[1][1] = f;
    x[1][2] = g;
    x[1][3] = h;
    x[2][0] = i;
    x[2][1] = j;
    x[2][2] = k;
    x[2][3] = l;
    x[3][0] = m;
    x[3][1] = n;
    x[3][2] = o;
    x[3][3] = p;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t) noexcept
{
    x[0][0] = r[0][0];
    x[0][1] = r[0][1];
    x[0][2] = r[0][2];
    x[0][3] = 0;
    x[1][0] = r[1][0];
    x[1][1] = r[1][1];
    x[1][2] = r[1][2];
    x[1][3] = 0;
    x[2][0] = r[2][0];
    x[2][1] = r[2][1];
    x[2][2] = r[2][2];
    x[2][3] = 0;
    x[3][0] = t.x;
    x[3][1] = t.y;
    x[3][2] = t.z;
    x[3][3] = 1;
}
 
template <class T> ATOMSMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (const Matrix44& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[0][3] = v.x[0][3];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[1][3] = v.x[1][3];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    x[2][3] = v.x[2][3];
    x[3][0] = v.x[3][0];
    x[3][1] = v.x[3][1];
    x[3][2] = v.x[3][2];
    x[3][3] = v.x[3][3];
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (const Matrix44<S>& v) noexcept
{
    x[0][0] = T (v.x[0][0]);
    x[0][1] = T (v.x[0][1]);
    x[0][2] = T (v.x[0][2]);
    x[0][3] = T (v.x[0][3]);
    x[1][0] = T (v.x[1][0]);
    x[1][1] = T (v.x[1][1]);
    x[1][2] = T (v.x[1][2]);
    x[1][3] = T (v.x[1][3]);
    x[2][0] = T (v.x[2][0]);
    x[2][1] = T (v.x[2][1]);
    x[2][2] = T (v.x[2][2]);
    x[2][3] = T (v.x[2][3]);
    x[3][0] = T (v.x[3][0]);
    x[3][1] = T (v.x[3][1]);
    x[3][2] = T (v.x[3][2]);
    x[3][3] = T (v.x[3][3]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator= (const Matrix44& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[0][3] = v.x[0][3];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[1][3] = v.x[1][3];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    x[2][3] = v.x[2][3];
    x[3][0] = v.x[3][0];
    x[3][1] = v.x[3][1];
    x[3][2] = v.x[3][2];
    x[3][3] = v.x[3][3];
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator= (T a) noexcept
{
    x[0][0] = a;
    x[0][1] = a;
    x[0][2] = a;
    x[0][3] = a;
    x[1][0] = a;
    x[1][1] = a;
    x[1][2] = a;
    x[1][3] = a;
    x[2][0] = a;
    x[2][1] = a;
    x[2][2] = a;
    x[2][3] = a;
    x[3][0] = a;
    x[3][1] = a;
    x[3][2] = a;
    x[3][3] = a;
    return *this;
}
 
template <class T>
inline T*
Matrix44<T>::getValue() noexcept
{
    return (T*) &x[0][0];
}
 
template <class T>
inline const T*
Matrix44<T>::getValue() const noexcept
{
    return (const T*) &x[0][0];
}
 
template <class T>
template <class S>
inline void
Matrix44<T>::getValue (Matrix44<S>& v) const noexcept
{
    v.x[0][0] = x[0][0];
    v.x[0][1] = x[0][1];
    v.x[0][2] = x[0][2];
    v.x[0][3] = x[0][3];
    v.x[1][0] = x[1][0];
    v.x[1][1] = x[1][1];
    v.x[1][2] = x[1][2];
    v.x[1][3] = x[1][3];
    v.x[2][0] = x[2][0];
    v.x[2][1] = x[2][1];
    v.x[2][2] = x[2][2];
    v.x[2][3] = x[2][3];
    v.x[3][0] = x[3][0];
    v.x[3][1] = x[3][1];
    v.x[3][2] = x[3][2];
    v.x[3][3] = x[3][3];
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix44<T>&
Matrix44<T>::setValue (const Matrix44<S>& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[0][3] = v.x[0][3];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[1][3] = v.x[1][3];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    x[2][3] = v.x[2][3];
    x[3][0] = v.x[3][0];
    x[3][1] = v.x[3][1];
    x[3][2] = v.x[3][2];
    x[3][3] = v.x[3][3];
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline Matrix44<T>&
Matrix44<T>::setTheMatrix (const Matrix44<S>& v) noexcept
{
    x[0][0] = v.x[0][0];
    x[0][1] = v.x[0][1];
    x[0][2] = v.x[0][2];
    x[0][3] = v.x[0][3];
    x[1][0] = v.x[1][0];
    x[1][1] = v.x[1][1];
    x[1][2] = v.x[1][2];
    x[1][3] = v.x[1][3];
    x[2][0] = v.x[2][0];
    x[2][1] = v.x[2][1];
    x[2][2] = v.x[2][2];
    x[2][3] = v.x[2][3];
    x[3][0] = v.x[3][0];
    x[3][1] = v.x[3][1];
    x[3][2] = v.x[3][2];
    x[3][3] = v.x[3][3];
    return *this;
}
 
template <class T>
inline void
Matrix44<T>::makeIdentity() noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;
    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
    x[1][3] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
    x[2][3] = 0;
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
}
 
template <class T>
constexpr inline bool
Matrix44<T>::operator== (const Matrix44& v) const noexcept
{
    return x[0][0] == v.x[0][0] && x[0][1] == v.x[0][1] && x[0][2] == v.x[0][2] &&
           x[0][3] == v.x[0][3] && x[1][0] == v.x[1][0] && x[1][1] == v.x[1][1] &&
           x[1][2] == v.x[1][2] && x[1][3] == v.x[1][3] && x[2][0] == v.x[2][0] &&
           x[2][1] == v.x[2][1] && x[2][2] == v.x[2][2] && x[2][3] == v.x[2][3] &&
           x[3][0] == v.x[3][0] && x[3][1] == v.x[3][1] && x[3][2] == v.x[3][2] &&
           x[3][3] == v.x[3][3];
}
 
template <class T>
constexpr inline bool
Matrix44<T>::operator!= (const Matrix44& v) const noexcept
{
    return x[0][0] != v.x[0][0] || x[0][1] != v.x[0][1] || x[0][2] != v.x[0][2] ||
           x[0][3] != v.x[0][3] || x[1][0] != v.x[1][0] || x[1][1] != v.x[1][1] ||
           x[1][2] != v.x[1][2] || x[1][3] != v.x[1][3] || x[2][0] != v.x[2][0] ||
           x[2][1] != v.x[2][1] || x[2][2] != v.x[2][2] || x[2][3] != v.x[2][3] ||
           x[3][0] != v.x[3][0] || x[3][1] != v.x[3][1] || x[3][2] != v.x[3][2] ||
           x[3][3] != v.x[3][3];
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline bool
Matrix44<T>::equalWithAbsError (const Matrix44<T>& m, T e) const noexcept
{
    for (int i = 0; i < 4; i++)
        for (int j = 0; j < 4; j++)
            if (!ATOMSMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this)[i][j], m[i][j], e))
                return false;
 
    return true;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline bool
Matrix44<T>::equalWithRelError (const Matrix44<T>& m, T e) const noexcept
{
    for (int i = 0; i < 4; i++)
        for (int j = 0; j < 4; j++)
            if (!ATOMSMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this)[i][j], m[i][j], e))
                return false;
 
    return true;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator+= (const Matrix44<T>& v) noexcept
{
    x[0][0] += v.x[0][0];
    x[0][1] += v.x[0][1];
    x[0][2] += v.x[0][2];
    x[0][3] += v.x[0][3];
    x[1][0] += v.x[1][0];
    x[1][1] += v.x[1][1];
    x[1][2] += v.x[1][2];
    x[1][3] += v.x[1][3];
    x[2][0] += v.x[2][0];
    x[2][1] += v.x[2][1];
    x[2][2] += v.x[2][2];
    x[2][3] += v.x[2][3];
    x[3][0] += v.x[3][0];
    x[3][1] += v.x[3][1];
    x[3][2] += v.x[3][2];
    x[3][3] += v.x[3][3];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator+= (T a) noexcept
{
    x[0][0] += a;
    x[0][1] += a;
    x[0][2] += a;
    x[0][3] += a;
    x[1][0] += a;
    x[1][1] += a;
    x[1][2] += a;
    x[1][3] += a;
    x[2][0] += a;
    x[2][1] += a;
    x[2][2] += a;
    x[2][3] += a;
    x[3][0] += a;
    x[3][1] += a;
    x[3][2] += a;
    x[3][3] += a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix44<T>
Matrix44<T>::operator+ (const Matrix44<T>& v) const noexcept
{
    return Matrix44 (x[0][0] + v.x[0][0],
                     x[0][1] + v.x[0][1],
                     x[0][2] + v.x[0][2],
                     x[0][3] + v.x[0][3],
                     x[1][0] + v.x[1][0],
                     x[1][1] + v.x[1][1],
                     x[1][2] + v.x[1][2],
                     x[1][3] + v.x[1][3],
                     x[2][0] + v.x[2][0],
                     x[2][1] + v.x[2][1],
                     x[2][2] + v.x[2][2],
                     x[2][3] + v.x[2][3],
                     x[3][0] + v.x[3][0],
                     x[3][1] + v.x[3][1],
                     x[3][2] + v.x[3][2],
                     x[3][3] + v.x[3][3]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator-= (const Matrix44<T>& v) noexcept
{
    x[0][0] -= v.x[0][0];
    x[0][1] -= v.x[0][1];
    x[0][2] -= v.x[0][2];
    x[0][3] -= v.x[0][3];
    x[1][0] -= v.x[1][0];
    x[1][1] -= v.x[1][1];
    x[1][2] -= v.x[1][2];
    x[1][3] -= v.x[1][3];
    x[2][0] -= v.x[2][0];
    x[2][1] -= v.x[2][1];
    x[2][2] -= v.x[2][2];
    x[2][3] -= v.x[2][3];
    x[3][0] -= v.x[3][0];
    x[3][1] -= v.x[3][1];
    x[3][2] -= v.x[3][2];
    x[3][3] -= v.x[3][3];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator-= (T a) noexcept
{
    x[0][0] -= a;
    x[0][1] -= a;
    x[0][2] -= a;
    x[0][3] -= a;
    x[1][0] -= a;
    x[1][1] -= a;
    x[1][2] -= a;
    x[1][3] -= a;
    x[2][0] -= a;
    x[2][1] -= a;
    x[2][2] -= a;
    x[2][3] -= a;
    x[3][0] -= a;
    x[3][1] -= a;
    x[3][2] -= a;
    x[3][3] -= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix44<T>
Matrix44<T>::operator- (const Matrix44<T>& v) const noexcept
{
    return Matrix44 (x[0][0] - v.x[0][0],
                     x[0][1] - v.x[0][1],
                     x[0][2] - v.x[0][2],
                     x[0][3] - v.x[0][3],
                     x[1][0] - v.x[1][0],
                     x[1][1] - v.x[1][1],
                     x[1][2] - v.x[1][2],
                     x[1][3] - v.x[1][3],
                     x[2][0] - v.x[2][0],
                     x[2][1] - v.x[2][1],
                     x[2][2] - v.x[2][2],
                     x[2][3] - v.x[2][3],
                     x[3][0] - v.x[3][0],
                     x[3][1] - v.x[3][1],
                     x[3][2] - v.x[3][2],
                     x[3][3] - v.x[3][3]);
}
 
template <class T>
constexpr inline Matrix44<T>
Matrix44<T>::operator-() const noexcept
{
    return Matrix44 (-x[0][0],
                     -x[0][1],
                     -x[0][2],
                     -x[0][3],
                     -x[1][0],
                     -x[1][1],
                     -x[1][2],
                     -x[1][3],
                     -x[2][0],
                     -x[2][1],
                     -x[2][2],
                     -x[2][3],
                     -x[3][0],
                     -x[3][1],
                     -x[3][2],
                     -x[3][3]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::negate() noexcept
{
    x[0][0] = -x[0][0];
    x[0][1] = -x[0][1];
    x[0][2] = -x[0][2];
    x[0][3] = -x[0][3];
    x[1][0] = -x[1][0];
    x[1][1] = -x[1][1];
    x[1][2] = -x[1][2];
    x[1][3] = -x[1][3];
    x[2][0] = -x[2][0];
    x[2][1] = -x[2][1];
    x[2][2] = -x[2][2];
    x[2][3] = -x[2][3];
    x[3][0] = -x[3][0];
    x[3][1] = -x[3][1];
    x[3][2] = -x[3][2];
    x[3][3] = -x[3][3];
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator*= (T a) noexcept
{
    x[0][0] *= a;
    x[0][1] *= a;
    x[0][2] *= a;
    x[0][3] *= a;
    x[1][0] *= a;
    x[1][1] *= a;
    x[1][2] *= a;
    x[1][3] *= a;
    x[2][0] *= a;
    x[2][1] *= a;
    x[2][2] *= a;
    x[2][3] *= a;
    x[3][0] *= a;
    x[3][1] *= a;
    x[3][2] *= a;
    x[3][3] *= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix44<T>
Matrix44<T>::operator* (T a) const noexcept
{
    return Matrix44 (x[0][0] * a,
                     x[0][1] * a,
                     x[0][2] * a,
                     x[0][3] * a,
                     x[1][0] * a,
                     x[1][1] * a,
                     x[1][2] * a,
                     x[1][3] * a,
                     x[2][0] * a,
                     x[2][1] * a,
                     x[2][2] * a,
                     x[2][3] * a,
                     x[3][0] * a,
                     x[3][1] * a,
                     x[3][2] * a,
                     x[3][3] * a);
}
 
template <class T>
inline Matrix44<T>
operator* (T a, const Matrix44<T>& v) noexcept
{
    return v * a;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator*= (const Matrix44<T>& v) noexcept
{
    Matrix44 tmp (T (0));
 
    multiply (*this, v, tmp);
    *this = tmp;
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix44<T>
Matrix44<T>::operator* (const Matrix44<T>& v) const noexcept
{
    Matrix44 tmp (T (0));
 
    multiply (*this, v, tmp);
    return tmp;
}
 
template <class T>
inline void
Matrix44<T>::multiply (const Matrix44<T>& a, const Matrix44<T>& b, Matrix44<T>& c) noexcept
{
    const T* ATOMSMATH_RESTRICT ap = &a.x[0][0];
    const T* ATOMSMATH_RESTRICT bp = &b.x[0][0];
    T* ATOMSMATH_RESTRICT cp       = &c.x[0][0];
 
    T a0, a1, a2, a3;
 
    a0 = ap[0];
    a1 = ap[1];
    a2 = ap[2];
    a3 = ap[3];
 
    cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
    cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
    cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
    cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
 
    a0 = ap[4];
    a1 = ap[5];
    a2 = ap[6];
    a3 = ap[7];
 
    cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
    cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
    cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
    cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
 
    a0 = ap[8];
    a1 = ap[9];
    a2 = ap[10];
    a3 = ap[11];
 
    cp[8]  = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
    cp[9]  = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
    cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
    cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
 
    a0 = ap[12];
    a1 = ap[13];
    a2 = ap[14];
    a3 = ap[15];
 
    cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
    cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
    cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
    cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
}
 
template <class T>
template <class S>
inline void
Matrix44<T>::multVecMatrix (const Vec3<S>& src, Vec3<S>& dst) const noexcept
{
    S a, b, c, w;
 
    a = src.x * x[0][0] + src.y * x[1][0] + src.z * x[2][0] + x[3][0];
    b = src.x * x[0][1] + src.y * x[1][1] + src.z * x[2][1] + x[3][1];
    c = src.x * x[0][2] + src.y * x[1][2] + src.z * x[2][2] + x[3][2];
    w = src.x * x[0][3] + src.y * x[1][3] + src.z * x[2][3] + x[3][3];
 
    dst.x = a / w;
    dst.y = b / w;
    dst.z = c / w;
}
 
template <class T>
template <class S>
inline void
Matrix44<T>::multDirMatrix (const Vec3<S>& src, Vec3<S>& dst) const noexcept
{
    S a, b, c;
 
    a = src.x * x[0][0] + src.y * x[1][0] + src.z * x[2][0];
    b = src.x * x[0][1] + src.y * x[1][1] + src.z * x[2][1];
    c = src.x * x[0][2] + src.y * x[1][2] + src.z * x[2][2];
 
    dst.x = a;
    dst.y = b;
    dst.z = c;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::operator/= (T a) noexcept
{
    x[0][0] /= a;
    x[0][1] /= a;
    x[0][2] /= a;
    x[0][3] /= a;
    x[1][0] /= a;
    x[1][1] /= a;
    x[1][2] /= a;
    x[1][3] /= a;
    x[2][0] /= a;
    x[2][1] /= a;
    x[2][2] /= a;
    x[2][3] /= a;
    x[3][0] /= a;
    x[3][1] /= a;
    x[3][2] /= a;
    x[3][3] /= a;
 
    return *this;
}
 
template <class T>
constexpr inline Matrix44<T>
Matrix44<T>::operator/ (T a) const noexcept
{
    return Matrix44 (x[0][0] / a,
                     x[0][1] / a,
                     x[0][2] / a,
                     x[0][3] / a,
                     x[1][0] / a,
                     x[1][1] / a,
                     x[1][2] / a,
                     x[1][3] / a,
                     x[2][0] / a,
                     x[2][1] / a,
                     x[2][2] / a,
                     x[2][3] / a,
                     x[3][0] / a,
                     x[3][1] / a,
                     x[3][2] / a,
                     x[3][3] / a);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::transpose() noexcept
{
    Matrix44 tmp (x[0][0],
                  x[1][0],
                  x[2][0],
                  x[3][0],
                  x[0][1],
                  x[1][1],
                  x[2][1],
                  x[3][1],
                  x[0][2],
                  x[1][2],
                  x[2][2],
                  x[3][2],
                  x[0][3],
                  x[1][3],
                  x[2][3],
                  x[3][3]);
    *this = tmp;
    return *this;
}
 
template <class T>
constexpr inline Matrix44<T>
Matrix44<T>::transposed() const noexcept
{
    return Matrix44 (x[0][0],
                     x[1][0],
                     x[2][0],
                     x[3][0],
                     x[0][1],
                     x[1][1],
                     x[2][1],
                     x[3][1],
                     x[0][2],
                     x[1][2],
                     x[2][2],
                     x[3][2],
                     x[0][3],
                     x[1][3],
                     x[2][3],
                     x[3][3]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::gjInvert (bool singExc)
{
    *this = gjInverse (singExc);
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::gjInvert() noexcept
{
    *this = gjInverse();
    return *this;
}
 
template <class T>
inline Matrix44<T>
Matrix44<T>::gjInverse (bool singExc) const
{
    int i, j, k;
    Matrix44 s;
    Matrix44 t (*this);
 
    // Forward elimination
 
    for (i = 0; i < 3; i++)
    {
        int pivot = i;
 
        T pivotsize = t[i][i];
 
        if (pivotsize < 0)
            pivotsize = -pivotsize;
 
        for (j = i + 1; j < 4; j++)
        {
            T tmp = t[j][i];
 
            if (tmp < 0)
                tmp = -tmp;
 
            if (tmp > pivotsize)
            {
                pivot     = j;
                pivotsize = tmp;
            }
        }
 
        if (pivotsize == 0)
        {
            #ifdef ATOMS_USE_EXCEPTIONS
            if (singExc)
                throw std::invalid_argument ("Cannot invert singular matrix.");
            #endif
            return Matrix44();
        }
 
        if (pivot != i)
        {
            for (j = 0; j < 4; j++)
            {
                T tmp;
 
                tmp         = t[i][j];
                t[i][j]     = t[pivot][j];
                t[pivot][j] = tmp;
 
                tmp         = s[i][j];
                s[i][j]     = s[pivot][j];
                s[pivot][j] = tmp;
            }
        }
 
        for (j = i + 1; j < 4; j++)
        {
            T f = t[j][i] / t[i][i];
 
            for (k = 0; k < 4; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    // Backward substitution
 
    for (i = 3; i >= 0; --i)
    {
        T f;
 
        if ((f = t[i][i]) == 0)
        {
            #ifdef ATOMS_USE_EXCEPTIONS
            if (singExc)
                throw std::invalid_argument ("Cannot invert singular matrix.");
            #endif
            return Matrix44();
        }
 
        for (j = 0; j < 4; j++)
        {
            t[i][j] /= f;
            s[i][j] /= f;
        }
 
        for (j = 0; j < i; j++)
        {
            f = t[j][i];
 
            for (k = 0; k < 4; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    return s;
}
 
template <class T>
inline Matrix44<T>
Matrix44<T>::gjInverse() const noexcept
{
    int i, j, k;
    Matrix44 s;
    Matrix44 t (*this);
 
    // Forward elimination
 
    for (i = 0; i < 3; i++)
    {
        int pivot = i;
 
        T pivotsize = t[i][i];
 
        if (pivotsize < 0)
            pivotsize = -pivotsize;
 
        for (j = i + 1; j < 4; j++)
        {
            T tmp = t[j][i];
 
            if (tmp < 0)
                tmp = -tmp;
 
            if (tmp > pivotsize)
            {
                pivot     = j;
                pivotsize = tmp;
            }
        }
 
        if (pivotsize == 0)
        {
            return Matrix44();
        }
 
        if (pivot != i)
        {
            for (j = 0; j < 4; j++)
            {
                T tmp;
 
                tmp         = t[i][j];
                t[i][j]     = t[pivot][j];
                t[pivot][j] = tmp;
 
                tmp         = s[i][j];
                s[i][j]     = s[pivot][j];
                s[pivot][j] = tmp;
            }
        }
 
        for (j = i + 1; j < 4; j++)
        {
            T f = t[j][i] / t[i][i];
 
            for (k = 0; k < 4; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    // Backward substitution
 
    for (i = 3; i >= 0; --i)
    {
        T f;
 
        if ((f = t[i][i]) == 0)
        {
            return Matrix44();
        }
 
        for (j = 0; j < 4; j++)
        {
            t[i][j] /= f;
            s[i][j] /= f;
        }
 
        for (j = 0; j < i; j++)
        {
            f = t[j][i];
 
            for (k = 0; k < 4; k++)
            {
                t[j][k] -= f * t[i][k];
                s[j][k] -= f * s[i][k];
            }
        }
    }
 
    return s;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::invert (bool singExc)
{
    *this = inverse (singExc);
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::invert() noexcept
{
    *this = inverse();
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix44<T>
Matrix44<T>::inverse (bool singExc) const
{
    if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
        return gjInverse (singExc);
 
    Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
                x[2][1] * x[0][2] - x[0][1] * x[2][2],
                x[0][1] * x[1][2] - x[1][1] * x[0][2],
                0,
 
                x[2][0] * x[1][2] - x[1][0] * x[2][2],
                x[0][0] * x[2][2] - x[2][0] * x[0][2],
                x[1][0] * x[0][2] - x[0][0] * x[1][2],
                0,
 
                x[1][0] * x[2][1] - x[2][0] * x[1][1],
                x[2][0] * x[0][1] - x[0][0] * x[2][1],
                x[0][0] * x[1][1] - x[1][0] * x[0][1],
                0,
 
                0,
                0,
                0,
                1);
 
    T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
 
    if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
    {
        for (int i = 0; i < 3; ++i)
        {
            for (int j = 0; j < 3; ++j)
            {
                s[i][j] /= r;
            }
        }
    }
    else
    {
        T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
        for (int i = 0; i < 3; ++i)
        {
            for (int j = 0; j < 3; ++j)
            {
                if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                {
                    s[i][j] /= r;
                }
                else
                {
                    #ifdef ATOMS_USE_EXCEPTIONS
                    if (singExc)
                        throw std::invalid_argument ("Cannot invert singular matrix.");
                    #endif
                    return Matrix44();
                }
            }
        }
    }
 
    s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
    s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
    s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];
 
    return s;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline Matrix44<T>
Matrix44<T>::inverse() const noexcept
{
    if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
        return gjInverse();
 
    Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
                x[2][1] * x[0][2] - x[0][1] * x[2][2],
                x[0][1] * x[1][2] - x[1][1] * x[0][2],
                0,
 
                x[2][0] * x[1][2] - x[1][0] * x[2][2],
                x[0][0] * x[2][2] - x[2][0] * x[0][2],
                x[1][0] * x[0][2] - x[0][0] * x[1][2],
                0,
 
                x[1][0] * x[2][1] - x[2][0] * x[1][1],
                x[2][0] * x[0][1] - x[0][0] * x[2][1],
                x[0][0] * x[1][1] - x[1][0] * x[0][1],
                0,
 
                0,
                0,
                0,
                1);
 
    T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
 
    if (ATOMSMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
    {
        for (int i = 0; i < 3; ++i)
        {
            for (int j = 0; j < 3; ++j)
            {
                s[i][j] /= r;
            }
        }
    }
    else
    {
        T mr = ATOMSMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
 
        for (int i = 0; i < 3; ++i)
        {
            for (int j = 0; j < 3; ++j)
            {
                if (mr > ATOMSMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
                {
                    s[i][j] /= r;
                }
                else
                {
                    return Matrix44();
                }
            }
        }
    }
 
    s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
    s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
    s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];
 
    return s;
}
 
template <class T>
constexpr inline T
Matrix44<T>::fastMinor (const int r0,
                        const int r1,
                        const int r2,
                        const int c0,
                        const int c1,
                        const int c2) const noexcept
{
    return x[r0][c0] * (x[r1][c1] * x[r2][c2] - x[r1][c2] * x[r2][c1]) +
           x[r0][c1] * (x[r1][c2] * x[r2][c0] - x[r1][c0] * x[r2][c2]) +
           x[r0][c2] * (x[r1][c0] * x[r2][c1] - x[r1][c1] * x[r2][c0]);
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline T
Matrix44<T>::minorOf (const int r, const int c) const noexcept
{
    int r0 = 0 + (r < 1 ? 1 : 0);
    int r1 = 1 + (r < 2 ? 1 : 0);
    int r2 = 2 + (r < 3 ? 1 : 0);
    int c0 = 0 + (c < 1 ? 1 : 0);
    int c1 = 1 + (c < 2 ? 1 : 0);
    int c2 = 2 + (c < 3 ? 1 : 0);
 
    Matrix33<T> working (x[r0][c0],
                         x[r1][c0],
                         x[r2][c0],
                         x[r0][c1],
                         x[r1][c1],
                         x[r2][c1],
                         x[r0][c2],
                         x[r1][c2],
                         x[r2][c2]);
 
    return working.determinant();
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline T
Matrix44<T>::determinant() const noexcept
{
    T sum = (T) 0;
 
    if (x[0][3] != 0.)
        sum -= x[0][3] * fastMinor (1, 2, 3, 0, 1, 2);
    if (x[1][3] != 0.)
        sum += x[1][3] * fastMinor (0, 2, 3, 0, 1, 2);
    if (x[2][3] != 0.)
        sum -= x[2][3] * fastMinor (0, 1, 3, 0, 1, 2);
    if (x[3][3] != 0.)
        sum += x[3][3] * fastMinor (0, 1, 2, 0, 1, 2);
 
    return sum;
}
 
template <class T>
template <class S>
inline const Matrix44<T>&
Matrix44<T>::setEulerAngles (const Vec3<S>& r) noexcept
{
    S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
 
    cos_rz = cos ((T) r.z);
    cos_ry = cos ((T) r.y);
    cos_rx = cos ((T) r.x);
 
    sin_rz = sin ((T) r.z);
    sin_ry = sin ((T) r.y);
    sin_rx = sin ((T) r.x);
 
    x[0][0] = cos_rz * cos_ry;
    x[0][1] = sin_rz * cos_ry;
    x[0][2] = -sin_ry;
    x[0][3] = 0;
 
    x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
    x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
    x[1][2] = cos_ry * sin_rx;
    x[1][3] = 0;
 
    x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
    x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
    x[2][2] = cos_ry * cos_rx;
    x[2][3] = 0;
 
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle) noexcept
{
    Vec3<S> unit (axis.normalized());
    S sine   = std::sin (angle);
    S cosine = std::cos (angle);
 
    x[0][0] = unit.x * unit.x * (1 - cosine) + cosine;
    x[0][1] = unit.x * unit.y * (1 - cosine) + unit.z * sine;
    x[0][2] = unit.x * unit.z * (1 - cosine) - unit.y * sine;
    x[0][3] = 0;
 
    x[1][0] = unit.x * unit.y * (1 - cosine) - unit.z * sine;
    x[1][1] = unit.y * unit.y * (1 - cosine) + cosine;
    x[1][2] = unit.y * unit.z * (1 - cosine) + unit.x * sine;
    x[1][3] = 0;
 
    x[2][0] = unit.x * unit.z * (1 - cosine) + unit.y * sine;
    x[2][1] = unit.y * unit.z * (1 - cosine) - unit.x * sine;
    x[2][2] = unit.z * unit.z * (1 - cosine) + cosine;
    x[2][3] = 0;
 
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
inline const Matrix44<T>&
Matrix44<T>::rotate (const Vec3<S>& r) noexcept
{
    S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
    S m00, m01, m02;
    S m10, m11, m12;
    S m20, m21, m22;
 
    cos_rz = cos ((S) r.z);
    cos_ry = cos ((S) r.y);
    cos_rx = cos ((S) r.x);
 
    sin_rz = sin ((S) r.z);
    sin_ry = sin ((S) r.y);
    sin_rx = sin ((S) r.x);
 
    m00 = cos_rz * cos_ry;
    m01 = sin_rz * cos_ry;
    m02 = -sin_ry;
    m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
    m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
    m12 = cos_ry * sin_rx;
    m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
    m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
    m22 = cos_ry * cos_rx;
 
    Matrix44<T> P (*this);
 
    x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02;
    x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02;
    x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02;
    x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02;
 
    x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12;
    x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12;
    x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12;
    x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12;
 
    x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22;
    x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22;
    x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22;
    x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22;
 
    return *this;
}
 
template <class T>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::setScale (T s) noexcept
{
    //
    // Set the matrix to a 3D homogeneous transform scale:
    //  | s 0 0 0 |
    //  | 0 s 0 0 |
    //  | 0 0 s 0 |
    //  | 0 0 0 1 |
    //
 
    x[0][0] = s;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;
    x[1][0] = 0;
    x[1][1] = s;
    x[1][2] = 0;
    x[1][3] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = s;
    x[2][3] = 0;
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::setScale (const Vec3<S>& s) noexcept
{
    //
    // Set the matrix to a 3D homogeneous transform scale:
    //  | s.x  0   0   0 |
    //  |  0  s.y  0   0 |
    //  |  0   0  s.z  0 |
    //  |  0   0   0   1 |
    //
 
    x[0][0] = s.x;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;
    x[1][0] = 0;
    x[1][1] = s.y;
    x[1][2] = 0;
    x[1][3] = 0;
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = s.z;
    x[2][3] = 0;
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::scale (const Vec3<S>& s) noexcept
{
    x[0][0] *= s.x;
    x[0][1] *= s.x;
    x[0][2] *= s.x;
    x[0][3] *= s.x;
 
    x[1][0] *= s.y;
    x[1][1] *= s.y;
    x[1][2] *= s.y;
    x[1][3] *= s.y;
 
    x[2][0] *= s.z;
    x[2][1] *= s.z;
    x[2][2] *= s.z;
    x[2][3] *= s.z;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::setTranslation (const Vec3<S>& t) noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;
 
    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
    x[1][3] = 0;
 
    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
    x[2][3] = 0;
 
    x[3][0] = t.x;
    x[3][1] = t.y;
    x[3][2] = t.z;
    x[3][3] = 1;
 
    return *this;
}
 
template <class T>
constexpr inline const Vec3<T>
Matrix44<T>::translation() const noexcept
{
    return Vec3<T> (x[3][0], x[3][1], x[3][2]);
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::translate (const Vec3<S>& t) noexcept
{
    x[3][0] += t.x * x[0][0] + t.y * x[1][0] + t.z * x[2][0];
    x[3][1] += t.x * x[0][1] + t.y * x[1][1] + t.z * x[2][1];
    x[3][2] += t.x * x[0][2] + t.y * x[1][2] + t.z * x[2][2];
    x[3][3] += t.x * x[0][3] + t.y * x[1][3] + t.z * x[2][3];
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::setShear (const Vec3<S>& h) noexcept
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;
 
    x[1][0] = h.x;
    x[1][1] = 1;
    x[1][2] = 0;
    x[1][3] = 0;
 
    x[2][0] = h.y;
    x[2][1] = h.z;
    x[2][2] = 1;
    x[2][3] = 0;
 
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::setShear (const Shear6<S>& h) noexcept
{
    x[0][0] = 1;
    x[0][1] = h.yx;
    x[0][2] = h.zx;
    x[0][3] = 0;
 
    x[1][0] = h.xy;
    x[1][1] = 1;
    x[1][2] = h.zy;
    x[1][3] = 0;
 
    x[2][0] = h.xz;
    x[2][1] = h.yz;
    x[2][2] = 1;
    x[2][3] = 0;
 
    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::shear (const Vec3<S>& h) noexcept
{
    //
    // In this case, we don't need a temp. copy of the matrix
    // because we never use a value on the RHS after we've
    // changed it on the LHS.
    //
 
    for (int i = 0; i < 4; i++)
    {
        x[2][i] += h.y * x[0][i] + h.z * x[1][i];
        x[1][i] += h.x * x[0][i];
    }
 
    return *this;
}
 
template <class T>
template <class S>
ATOMSMATH_CONSTEXPR14 inline const Matrix44<T>&
Matrix44<T>::shear (const Shear6<S>& h) noexcept
{
    Matrix44<T> P (*this);
 
    for (int i = 0; i < 4; i++)
    {
        x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i];
        x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i];
        x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i];
    }
 
    return *this;
}
 
//--------------------------------
// Implementation of stream output
//--------------------------------
 
template <class T>
std::ostream&
operator<< (std::ostream& s, const Matrix22<T>& m)
{
    std::ios_base::fmtflags oldFlags = s.flags();
    int width;
 
    if (s.flags() & std::ios_base::fixed)
    {
        s.setf (std::ios_base::showpoint);
        width = static_cast<int> (s.precision()) + 5;
    }
    else
    {
        s.setf (std::ios_base::scientific);
        s.setf (std::ios_base::showpoint);
        width = static_cast<int> (s.precision()) + 8;
    }
 
    s << "(" << std::setw (width) << m[0][0] << ", " << std::setw (width) << m[0][1] << ",\n"
      <<
 
        " " << std::setw (width) << m[1][0] << ", " << std::setw (width) << m[1][1] << ")\n";
 
    s.flags (oldFlags);
    return s;
}
 
template <class T>
std::ostream&
operator<< (std::ostream& s, const Matrix33<T>& m)
{
    std::ios_base::fmtflags oldFlags = s.flags();
    int width;
 
    if (s.flags() & std::ios_base::fixed)
    {
        s.setf (std::ios_base::showpoint);
        width = static_cast<int> (s.precision()) + 5;
    }
    else
    {
        s.setf (std::ios_base::scientific);
        s.setf (std::ios_base::showpoint);
        width = static_cast<int> (s.precision()) + 8;
    }
 
    s << "(" << std::setw (width) << m[0][0] << ", " << std::setw (width) << m[0][1] << ", "
      << std::setw (width) << m[0][2] << ",\n"
      <<
 
        " " << std::setw (width) << m[1][0] << ", " << std::setw (width) << m[1][1] << ", "
      << std::setw (width) << m[1][2] << ",\n"
      <<
 
        " " << std::setw (width) << m[2][0] << ", " << std::setw (width) << m[2][1] << ", "
      << std::setw (width) << m[2][2] << ")\n";
 
    s.flags (oldFlags);
    return s;
}
 
template <class T>
std::ostream&
operator<< (std::ostream& s, const Matrix44<T>& m)
{
    std::ios_base::fmtflags oldFlags = s.flags();
    int width;
 
    if (s.flags() & std::ios_base::fixed)
    {
        s.setf (std::ios_base::showpoint);
        width = static_cast<int> (s.precision()) + 5;
    }
    else
    {
        s.setf (std::ios_base::scientific);
        s.setf (std::ios_base::showpoint);
        width = static_cast<int> (s.precision()) + 8;
    }
 
    s << "(" << std::setw (width) << m[0][0] << ", " << std::setw (width) << m[0][1] << ", "
      << std::setw (width) << m[0][2] << ", " << std::setw (width) << m[0][3] << ",\n"
      <<
 
        " " << std::setw (width) << m[1][0] << ", " << std::setw (width) << m[1][1] << ", "
      << std::setw (width) << m[1][2] << ", " << std::setw (width) << m[1][3] << ",\n"
      <<
 
        " " << std::setw (width) << m[2][0] << ", " << std::setw (width) << m[2][1] << ", "
      << std::setw (width) << m[2][2] << ", " << std::setw (width) << m[2][3] << ",\n"
      <<
 
        " " << std::setw (width) << m[3][0] << ", " << std::setw (width) << m[3][1] << ", "
      << std::setw (width) << m[3][2] << ", " << std::setw (width) << m[3][3] << ")\n";
 
    s.flags (oldFlags);
    return s;
}
 
//---------------------------------------------------------------
// Implementation of vector-times-matrix multiplication operators
//---------------------------------------------------------------
 
template <class S, class T>
inline const Vec2<S>&
operator*= (Vec2<S>& v, const Matrix22<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1]);
 
    v.x = x;
    v.y = y;
 
    return v;
}
 
template <class S, class T>
inline Vec2<S>
operator* (const Vec2<S>& v, const Matrix22<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1]);
 
    return Vec2<S> (x, y);
}
 
template <class S, class T>
inline const Vec2<S>&
operator*= (Vec2<S>& v, const Matrix33<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
    S w = S (v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
 
    v.x = x / w;
    v.y = y / w;
 
    return v;
}
 
template <class S, class T>
inline Vec2<S>
operator* (const Vec2<S>& v, const Matrix33<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
    S w = S (v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
 
    return Vec2<S> (x / w, y / w);
}
 
template <class S, class T>
inline const Vec3<S>&
operator*= (Vec3<S>& v, const Matrix33<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
    S z = S (v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
 
    v.x = x;
    v.y = y;
    v.z = z;
 
    return v;
}
 
template <class S, class T>
inline Vec3<S>
operator* (const Vec3<S>& v, const Matrix33<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
    S z = S (v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
 
    return Vec3<S> (x, y, z);
}
 
template <class S, class T>
inline const Vec3<S>&
operator*= (Vec3<S>& v, const Matrix44<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
    S z = S (v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
    S w = S (v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
 
    v.x = x / w;
    v.y = y / w;
    v.z = z / w;
 
    return v;
}
 
template <class S, class T>
inline Vec3<S>
operator* (const Vec3<S>& v, const Matrix44<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
    S z = S (v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
    S w = S (v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
 
    return Vec3<S> (x / w, y / w, z / w);
}
 
template <class S, class T>
inline const Vec4<S>&
operator*= (Vec4<S>& v, const Matrix44<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]);
    S z = S (v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]);
    S w = S (v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]);
 
    v.x = x;
    v.y = y;
    v.z = z;
    v.w = w;
 
    return v;
}
 
template <class S, class T>
inline Vec4<S>
operator* (const Vec4<S>& v, const Matrix44<T>& m) noexcept
{
    S x = S (v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]);
    S y = S (v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]);
    S z = S (v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]);
    S w = S (v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]);
 
    return Vec4<S> (x, y, z, w);
}
 
ATOMSMATH_INTERNAL_NAMESPACE_HEADER_EXIT
 
#endif // INCLUDED_ATOMSMATHMATRIX_H