// linalg.h - v2.0 - Single-header public domain linear algebra library
//
// The intent of this library is to provide the bulk of the functionality
// you need to write programs that frequently use small, fixed-size vectors
// and matrices, in domains such as computational geometry or computer
// graphics. It strives for terse, readable source code.
//
// The original author of this software is Sterling Orsten, and its permanent
// home is . If you find this software
// useful, an acknowledgement in your source text and/or product documentation
// is appreciated, but not required.
//
// The author acknowledges significant insights and contributions by:
// Stan Melax
// Dimitri Diakopoulos
// This is free and unencumbered software released into the public domain.
//
// Anyone is free to copy, modify, publish, use, compile, sell, or
// distribute this software, either in source code form or as a compiled
// binary, for any purpose, commercial or non-commercial, and by any
// means.
//
// In jurisdictions that recognize copyright laws, the author or authors
// of this software dedicate any and all copyright interest in the
// software to the public domain. We make this dedication for the benefit
// of the public at large and to the detriment of our heirs and
// successors. We intend this dedication to be an overt act of
// relinquishment in perpetuity of all present and future rights to this
// software under copyright law.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
// IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//
// For more information, please refer to
#pragma once
#ifndef LINALG_H
#define LINALG_H
#include // For various unary math functions, such as std::sqrt
#include // To resolve std::abs ambiguity on clang
#include // For implementing namespace linalg::aliases
#include // For std::array, used in the relational operator overloads
#include // For std::numeric_limits/epsilon
// Visual Studio versions prior to 2015 lack constexpr support
#if defined(_MSC_VER) && _MSC_VER < 1900 && !defined(constexpr)
#define constexpr
#endif
namespace linalg
{
// Small, fixed-length vector type, consisting of exactly M elements of type T, and presumed to be a column-vector unless otherwise noted
template struct vec;
template struct vec
{
T x,y;
constexpr vec() : x(), y() {}
constexpr vec(const T & x_, const T & y_) : x(x_), y(y_) {}
constexpr vec(const std::array & a) : x(a[0]), y(a[1]) {}
constexpr explicit vec(const T & s) : vec(s, s) {}
constexpr explicit vec(const T * p) : vec(p[0], p[1]) {}
template
constexpr explicit vec(const vec & v) : vec(static_cast(v.x), static_cast(v.y)) {}
constexpr const T & operator[] (int i) const { return (&x)[i]; }
T & operator[] (int i) { return (&x)[i]; }
};
template struct vec
{
T x,y,z;
constexpr vec() : x(), y(), z() {}
constexpr vec(const T & x_, const T & y_,
const T & z_) : x(x_), y(y_), z(z_) {}
constexpr vec(const vec & xy,
const T & z_) : vec(xy.x, xy.y, z_) {}
constexpr vec(const std::array & a) : x(a[0]), y(a[1]), z(a[2]) {}
constexpr explicit vec(const T & s) : vec(s, s, s) {}
constexpr explicit vec(const T * p) : vec(p[0], p[1], p[2]) {}
template
constexpr explicit vec(const vec & v) : vec(static_cast(v.x), static_cast(v.y), static_cast(v.z)) {}
constexpr const T & operator[] (int i) const { return (&x)[i]; }
T & operator[] (int i) { return (&x)[i]; }
constexpr const vec & xy() const { return *reinterpret_cast *>(this); }
vec & xy() { return *reinterpret_cast *>(this); }
};
template struct vec
{
T x,y,z,w;
constexpr vec() : x(), y(), z(), w() {}
constexpr vec(const T & x_, const T & y_,
const T & z_, const T & w_) : x(x_), y(y_), z(z_), w(w_) {}
constexpr vec(const vec & xy,
const T & z_, const T & w_) : vec(xy.x, xy.y, z_, w_) {}
constexpr vec(const vec & xyz,
const T & w_) : vec(xyz.x, xyz.y, xyz.z, w_) {}
constexpr vec(const std::array & a) : x(a[0]), y(a[1]), z(a[2]), w(a[3]) {}
constexpr explicit vec(const T & s) : vec(s, s, s, s) {}
constexpr explicit vec(const T * p) : vec(p[0], p[1], p[2], p[3]) {}
template
constexpr explicit vec(const vec & v) : vec(static_cast(v.x), static_cast(v.y), static_cast(v.z), static_cast(v.w)) {}
constexpr const T & operator[] (int i) const { return (&x)[i]; }
T & operator[] (int i) { return (&x)[i]; }
constexpr const vec & xy() const { return *reinterpret_cast *>(this); }
constexpr const vec & xyz() const { return *reinterpret_cast *>(this); }
vec & xy() { return *reinterpret_cast *>(this); }
vec & xyz() { return *reinterpret_cast *>(this); }
};
// Small, fixed-size matrix type, consisting of exactly M rows and N columns of type T, stored in column-major order.
template struct mat;
template struct mat
{
typedef vec V;
V x,y;
constexpr mat() : x(), y() {}
constexpr mat(const V & x_, const V & y_) : x(x_), y(y_) {}
constexpr explicit mat(const T & s) : x(s), y(s) {}
constexpr explicit mat(const T * p) : x(p+M*0), y(p+M*1) {}
template
constexpr explicit mat(const mat & m) : mat(V(m.x), V(m.y)) {}
constexpr vec row(int i) const { return {x[i], y[i]}; }
constexpr const V & operator[] (int j) const { return (&x)[j]; }
V & operator[] (int j) { return (&x)[j]; }
};
template struct mat
{
typedef vec V;
V x,y,z;
constexpr mat() : x(), y(), z() {}
constexpr mat(const V & x_, const V & y_,
const V & z_) : x(x_), y(y_), z(z_) {}
constexpr explicit mat(const T & s) : x(s), y(s), z(s) {}
constexpr explicit mat(const T * p) : x(p+M*0), y(p+M*1), z(p+M*2) {}
template
constexpr explicit mat(const mat & m) : mat(V(m.x), V(m.y), V(m.z)) {}
constexpr vec row(int i) const { return {x[i], y[i], z[i]}; }
constexpr const V & operator[] (int j) const { return (&x)[j]; }
V & operator[] (int j) { return (&x)[j]; }
};
template struct mat
{
typedef vec V;
V x,y,z,w;
constexpr mat() : x(), y(), z(), w() {}
constexpr mat(const V & x_, const V & y_,
const V & z_, const V & w_) : x(x_), y(y_), z(z_), w(w_) {}
constexpr explicit mat(const T & s) : x(s), y(s), z(s), w(s) {}
constexpr explicit mat(const T * p) : x(p+M*0), y(p+M*1), z(p+M*2), w(p+M*3) {}
template
constexpr explicit mat(const mat & m) : mat(V(m.x), V(m.y), V(m.z), V(m.w)) {}
constexpr vec row(int i) const { return {x[i], y[i], z[i], w[i]}; }
constexpr const V & operator[] (int j) const { return (&x)[j]; }
V & operator[] (int j) { return (&x)[j]; }
};
// Type traits for a binary operation involving linear algebra types, used for SFINAE on templated functions and operator overloads
template struct traits {};
template struct traits, vec> { typedef T scalar; typedef vec result; typedef vec bool_result; typedef vec arith_result; typedef std::array compare_as; };
template struct traits, T > { typedef T scalar; typedef vec result; typedef vec bool_result; typedef vec arith_result; };
template struct traits> { typedef T scalar; typedef vec result; typedef vec bool_result; typedef vec arith_result; };
template struct traits, mat> { typedef T scalar; typedef mat result; typedef mat bool_result; typedef mat arith_result; typedef std::array compare_as; };
template struct traits, T > { typedef T scalar; typedef mat result; typedef mat bool_result; typedef mat arith_result; };
template struct traits> { typedef T scalar; typedef mat result; typedef mat bool_result; typedef mat arith_result; };
template using scalar_t = typename traits::scalar; // Underlying scalar type when performing elementwise operations
template using result_t = typename traits::result; // Result of calling a function on linear algebra types
template using bool_result_t = typename traits::bool_result; // Result of a comparison or unary not operation on linear algebra types
template using arith_result_t = typename traits::arith_result; // Result of an arithmetic operation on linear algebra types (accounts for integer promotion)
// Produce a scalar by applying f(T,T) -> T to adjacent pairs of elements from vector/matrix a in left-to-right order (matching the associativity of arithmetic and logical operators)
template constexpr T fold(const vec & a, F f) { return f(a.x,a.y); }
template constexpr T fold(const vec & a, F f) { return f(f(a.x,a.y),a.z); }
template constexpr T fold(const vec & a, F f) { return f(f(f(a.x,a.y),a.z),a.w); }
template constexpr T fold(const mat & a, F f) { return f(fold(a.x,f),fold(a.y,f)); }
template constexpr T fold(const mat & a, F f) { return f(f(fold(a.x,f),fold(a.y,f)),fold(a.z,f)); }
template constexpr T fold(const mat & a, F f) { return f(f(f(fold(a.x,f),fold(a.y,f)),fold(a.z,f)),fold(a.w,f)); }
// Produce a vector/matrix by applying f(T,T) to corresponding pairs of elements from vectors/matrix a and b
template constexpr auto zip(const vec & a, const vec & b, F f) -> vec { return {f(a.x,b.x), f(a.y,b.y)}; }
template constexpr auto zip(const vec & a, const vec & b, F f) -> vec { return {f(a.x,b.x), f(a.y,b.y), f(a.z,b.z)}; }
template constexpr auto zip(const vec & a, const vec & b, F f) -> vec { return {f(a.x,b.x), f(a.y,b.y), f(a.z,b.z), f(a.w,b.w)}; }
template constexpr auto zip(const vec & a, T b, F f) -> vec { return zip(a, vec(b), f); }
template constexpr auto zip( T a, const vec & b, F f) -> vec { return zip(vec(a), b, f); }
template constexpr auto zip(const mat & a, const mat & b, F f) -> mat { return {zip(a.x,b.x,f), zip(a.y,b.y,f)}; }
template constexpr auto zip(const mat & a, const mat & b, F f) -> mat { return {zip(a.x,b.x,f), zip(a.y,b.y,f), zip(a.z,b.z,f)}; }
template constexpr auto zip(const mat & a, const mat & b, F f) -> mat { return {zip(a.x,b.x,f), zip(a.y,b.y,f), zip(a.z,b.z,f), zip(a.w,b.w,f)}; }
template constexpr auto zip(const mat & a, T b, F f) -> mat { return zip(a, mat(b), f); }
template constexpr auto zip( T a, const mat & b, F f) -> mat { return zip(mat(a), b, f); }
// Produce a vector/matrix by applying f(T) to elements from vector/matrix a
template constexpr auto map(const vec & a, F f) -> vec { return zip(a, a, [f](T l, T) { return f(l); }); }
template constexpr auto map(const mat & a, F f) -> mat { return zip(a, a, [f](T l, T) { return f(l); }); }
// Relational operators are defined to compare the elements of two vectors or matrices lexicographically, in column-major order
template::compare_as> constexpr bool operator == (const A & a, const A & b) { return reinterpret_cast(a) == reinterpret_cast(b); }
template::compare_as> constexpr bool operator != (const A & a, const A & b) { return reinterpret_cast(a) != reinterpret_cast(b); }
template::compare_as> constexpr bool operator < (const A & a, const A & b) { return reinterpret_cast(a) < reinterpret_cast(b); }
template::compare_as> constexpr bool operator > (const A & a, const A & b) { return reinterpret_cast(a) > reinterpret_cast(b); }
template::compare_as> constexpr bool operator <= (const A & a, const A & b) { return reinterpret_cast(a) <= reinterpret_cast(b); }
template::compare_as> constexpr bool operator >= (const A & a, const A & b) { return reinterpret_cast(a) >= reinterpret_cast(b); }
// Lambdas are not permitted inside constexpr functions, so we provide explicit function objects instead
namespace op
{
template struct pos { constexpr auto operator() (T r) const -> decltype(+r) { return +r; } };
template struct neg { constexpr auto operator() (T r) const -> decltype(-r) { return -r; } };
template struct add { constexpr auto operator() (T l, T r) const -> decltype(l + r) { return l + r; } };
template struct sub { constexpr auto operator() (T l, T r) const -> decltype(l - r) { return l - r; } };
template struct mul { constexpr auto operator() (T l, T r) const -> decltype(l * r) { return l * r; } };
template struct div { constexpr auto operator() (T l, T r) const -> decltype(l / r) { return l / r; } };
template struct mod { constexpr auto operator() (T l, T r) const -> decltype(l % r) { return l % r; } };
template struct lshift { constexpr auto operator() (T l, T r) const -> decltype(l << r) { return l << r; } };
template struct rshift { constexpr auto operator() (T l, T r) const -> decltype(l >> r) { return l >> r; } };
template struct binary_not { constexpr auto operator() (T r) const -> decltype(+r) { return ~r; } };
template struct binary_or { constexpr auto operator() (T l, T r) const -> decltype(l | r) { return l | r; } };
template struct binary_xor { constexpr auto operator() (T l, T r) const -> decltype(l ^ r) { return l ^ r; } };
template struct binary_and { constexpr auto operator() (T l, T r) const -> decltype(l & r) { return l & r; } };
template struct logical_not { constexpr bool operator() (T r) const { return !r; } };
template struct logical_or { constexpr bool operator() (T l, T r) const { return l || r; } };
template struct logical_and { constexpr bool operator() (T l, T r) const { return l && r; } };
template struct equal { constexpr bool operator() (T l, T r) const { return l == r; } };
template struct nequal { constexpr bool operator() (T l, T r) const { return l != r; } };
template struct less { constexpr bool operator() (T l, T r) const { return l < r; } };
template struct greater { constexpr bool operator() (T l, T r) const { return l > r; } };
template struct lequal { constexpr bool operator() (T l, T r) const { return l <= r; } };
template struct gequal { constexpr bool operator() (T l, T r) const { return l >= r; } };
template struct min { constexpr T operator() (T l, T r) const { return l < r ? l : r; } };
template struct max { constexpr T operator() (T l, T r) const { return l > r ? l : r; } };
}
// Functions for coalescing scalar values
template constexpr scalar_t any (const A & a) { return fold(a, op::logical_or>{}); }
template constexpr scalar_t all (const A & a) { return fold(a, op::logical_and>{}); }
template constexpr scalar_t sum (const A & a) { return fold(a, op::add>{}); }
template constexpr scalar_t product(const A & a) { return fold(a, op::mul>{}); }
template int argmin(const vec & a) { int j=0; for(int i=1; i int argmax(const vec & a) { int j=0; for(int i=1; i a[j]) j = i; return j; }
template T minelem(const vec & a) { return a[argmin(a)]; }
template T maxelem(const vec & a) { return a[argmax(a)]; }
// Overloads for unary operators on vectors are implemented in terms of elementwise application of the operator
template constexpr arith_result_t operator + (const A & a) { return map(a, op::pos>{}); }
template constexpr arith_result_t operator - (const A & a) { return map(a, op::neg>{}); }
template constexpr arith_result_t operator ~ (const A & a) { return map(a, op::binary_not>{}); }
template constexpr bool_result_t operator ! (const A & a) { return map(a, op::logical_not>{}); }
// Mirror the set of unary scalar math functions to apply elementwise to vectors
template result_t abs (const A & a) { return map(a, [](scalar_t l) { return std::abs (l); }); }
template result_t floor(const A & a) { return map(a, [](scalar_t l) { return std::floor(l); }); }
template result_t ceil (const A & a) { return map(a, [](scalar_t l) { return std::ceil (l); }); }
template result_t exp (const A & a) { return map(a, [](scalar_t l) { return std::exp (l); }); }
template result_t log (const A & a) { return map(a, [](scalar_t l) { return std::log (l); }); }
template result_t log10(const A & a) { return map(a, [](scalar_t l) { return std::log10(l); }); }
template result_t