59#if defined(_MSC_VER) && (_MSC_VER <= 1900)
60#define LINALG_CONSTEXPR14
62#define LINALG_CONSTEXPR14 constexpr
68template <
class T,
int M>
73template <
class T,
int M,
int N>
78template <
class T,
class U>
81template <
class T,
class U>
82using conv_t =
typename std::enable_if<!std::is_same<T, U>::value,
84 std::declval<U>()))>::type;
89template <
class T,
int M>
93template <
class T,
int M,
int N>
105constexpr bool operator==(
const ord<T> &o, std::nullptr_t) {
109constexpr bool operator!=(
const ord<T> &o, std::nullptr_t) {
110 return !(o.a == o.b);
113constexpr bool operator<(
const ord<T> &o, std::nullptr_t) {
117constexpr bool operator>(
const ord<T> &o, std::nullptr_t) {
121constexpr bool operator<=(
const ord<T> &o, std::nullptr_t) {
125constexpr bool operator>=(
const ord<T> &o, std::nullptr_t) {
130template <
class A,
class B>
143 return !(a.x == b.x) ?
ord<T>{a.x, b.x} :
ord<T>{a.y, b.y};
150 return !(a.x == b.x) ?
ord<T>{a.x, b.x}
151 : !(a.y == b.y) ?
ord<T>{a.y, b.y}
159 return !(a.x == b.x) ?
ord<T>{a.x, b.x}
160 : !(a.y == b.y) ?
ord<T>{a.y, b.y}
161 : !(a.z == b.z) ?
ord<T>{a.z, b.z}
165template <
class T,
int M>
170 return compare(a.x, b.x);
173template <
class T,
int M>
178 return a.x != b.x ? compare(a.x, b.x) : compare(a.y, b.y);
181template <
class T,
int M>
186 return a.x != b.x ? compare(a.x, b.x)
187 : a.y != b.y ? compare(a.y, b.y)
191template <
class T,
int M>
196 return a.x != b.x ? compare(a.x, b.x)
197 : a.y != b.y ? compare(a.y, b.y)
198 : a.z != b.z ? compare(a.z, b.z)
210 constexpr auto operator()(A &a)
const ->
decltype(a.x) {
217 constexpr auto operator()(A &a)
const ->
decltype(a.y) {
224 constexpr auto operator()(A &a)
const ->
decltype(a.z) {
231 constexpr auto operator()(A &a)
const ->
decltype(a.w) {
239template <
int A,
int N>
261template <
int A,
int B>
263template <
class T,
int M,
int... I>
267template <
class T,
int M,
int N,
int... I,
int... J>
268mat<T,
sizeof...(I),
sizeof...(J)>
constexpr swizzle(
const mat<T, M, N> &m,
269 seq<I...> i, seq<J...> j) {
270 return {swizzle(getter<J>{}(m), i)...};
274template <
class F,
class... T>
275using ret_t =
decltype(std::declval<F>()(std::declval<T>()...));
285template <
class T,
class... U>
286struct scalars<T, U...> : std::conditional<std::is_arithmetic<T>::value,
287 scalars<U...>, empty>::type {};
289using scalars_t =
typename scalars<T...>::type;
293template <
class F,
class Void,
class... T>
295template <
class F,
int M,
class A>
298 enum { size = M, mm = 0 };
304template <
class F,
int M,
class A,
class B>
307 enum { size = M, mm = 0 };
314template <
class F,
int M,
class A,
class B>
317 enum { size = M, mm = 0 };
323template <
class F,
int M,
class A,
class B>
326 enum { size = M, mm = 0 };
332template <
class F,
int M,
class A,
class B,
class C>
335 enum { size = M, mm = 0 };
342template <
class F,
int M,
class A,
class B,
class C>
345 enum { size = M, mm = 0 };
352template <
class F,
int M,
class A,
class B,
class C>
355 enum { size = M, mm = 0 };
362template <
class F,
int M,
class A,
class B,
class C>
365 enum { size = M, mm = 0 };
371template <
class F,
int M,
class A,
class B,
class C>
374 enum { size = M, mm = 0 };
381template <
class F,
int M,
class A,
class B,
class C>
384 enum { size = M, mm = 0 };
390template <
class F,
int M,
class A,
class B,
class C>
393 enum { size = M, mm = 0 };
399template <
class F,
int M,
int N,
class A>
402 enum { size = N, mm = 0 };
409template <
class F,
int M,
int N,
class A,
class B>
412 enum { size = N, mm = 1 };
420template <
class F,
int M,
int N,
class A,
class B>
423 enum { size = N, mm = 0 };
430template <
class F,
int M,
int N,
class A,
class B>
433 enum { size = N, mm = 0 };
440template <
class F,
class... A>
441struct apply<F, scalars_t<A...>, A...> {
442 using type = ret_t<F, A...>;
443 enum { size = 0, mm = 0 };
444 static constexpr type impl(
seq<>, F f, A... a) {
return f(a...); }
449 template <
class A,
class B>
450 constexpr auto operator()(A a, B b)
const ->
451 typename std::remove_reference<
decltype(a < b ? a : b)>::type {
452 return a < b ? a : b;
456 template <
class A,
class B>
457 constexpr auto operator()(A a, B b)
const ->
458 typename std::remove_reference<
decltype(a < b ? b : a)>::type {
459 return a < b ? b : a;
463 template <
class A,
class B,
class C>
464 constexpr auto operator()(A a, B b, C c)
const ->
465 typename std::remove_reference<
decltype(a < b ? b
468 return a < b ? b : a < c ? a : c;
472 template <
class A,
class B,
class C>
473 constexpr auto operator()(A a, B b, C c)
const ->
474 typename std::remove_reference<
decltype(a ? b : c)>::type {
479 template <
class A,
class B,
class C>
480 constexpr auto operator()(A a, B b,
481 C c)
const ->
decltype(a * (1 - c) + b * c) {
482 return a * (1 - c) + b * c;
489 constexpr auto operator()(A a)
const ->
decltype(+a) {
495 constexpr auto operator()(A a)
const ->
decltype(-a) {
501 constexpr auto operator()(A a)
const ->
decltype(!a) {
507 constexpr auto operator()(A a)
const ->
decltype(~(a)) {
512 template <
class A,
class B>
513 constexpr auto operator()(A a, B b)
const ->
decltype(a * b) {
518 template <
class A,
class B>
519 constexpr auto operator()(A a, B b)
const ->
decltype(a / b) {
524 template <
class A,
class B>
525 constexpr auto operator()(A a, B b)
const ->
decltype(a % b) {
530 template <
class A,
class B>
531 constexpr auto operator()(A a, B b)
const ->
decltype(a + b) {
536 template <
class A,
class B>
537 constexpr auto operator()(A a, B b)
const ->
decltype(a - b) {
542 template <
class A,
class B>
543 constexpr auto operator()(A a, B b)
const ->
decltype(a << b) {
548 template <
class A,
class B>
549 constexpr auto operator()(A a, B b)
const ->
decltype(a >> b) {
554 template <
class A,
class B>
555 constexpr auto operator()(A a, B b)
const ->
decltype(a < b) {
560 template <
class A,
class B>
561 constexpr auto operator()(A a, B b)
const ->
decltype(a > b) {
566 template <
class A,
class B>
567 constexpr auto operator()(A a, B b)
const ->
decltype(a <= b) {
572 template <
class A,
class B>
573 constexpr auto operator()(A a, B b)
const ->
decltype(a >= b) {
578 template <
class A,
class B>
579 constexpr auto operator()(A a, B b)
const ->
decltype(a == b) {
584 template <
class A,
class B>
585 constexpr auto operator()(A a, B b)
const ->
decltype(a != b) {
590 template <
class A,
class B>
591 constexpr auto operator()(A a, B b)
const ->
decltype(a & b) {
596 template <
class A,
class B>
597 constexpr auto operator()(A a, B b)
const ->
decltype(a ^ b) {
602 template <
class A,
class B>
603 constexpr auto operator()(A a, B b)
const ->
decltype(a | b) {
608 template <
class A,
class B>
609 constexpr auto operator()(A a, B b)
const ->
decltype(a && b) {
614 template <
class A,
class B>
615 constexpr auto operator()(A a, B b)
const ->
decltype(a || b) {
623 constexpr auto operator()(A a)
const ->
decltype(std::isfinite(a)) {
624 return std::isfinite(a);
629 constexpr auto operator()(A a)
const ->
decltype(std::abs(a)) {
635 constexpr auto operator()(A a)
const ->
decltype(std::floor(a)) {
636 return std::floor(a);
641 constexpr auto operator()(A a)
const ->
decltype(std::ceil(a)) {
647 constexpr auto operator()(A a)
const ->
decltype(std::exp(a)) {
653 constexpr auto operator()(A a)
const ->
decltype(std::log(a)) {
659 constexpr auto operator()(A a)
const ->
decltype(std::log2(a)) {
665 constexpr auto operator()(A a)
const ->
decltype(std::log10(a)) {
666 return std::log10(a);
671 constexpr auto operator()(A a)
const ->
decltype(std::sqrt(a)) {
677 constexpr auto operator()(A a)
const ->
decltype(std::sin(a)) {
683 constexpr auto operator()(A a)
const ->
decltype(std::cos(a)) {
689 constexpr auto operator()(A a)
const ->
decltype(std::tan(a)) {
695 constexpr auto operator()(A a)
const ->
decltype(std::asin(a)) {
701 constexpr auto operator()(A a)
const ->
decltype(std::acos(a)) {
707 constexpr auto operator()(A a)
const ->
decltype(std::atan(a)) {
713 constexpr auto operator()(A a)
const ->
decltype(std::sinh(a)) {
719 constexpr auto operator()(A a)
const ->
decltype(std::cosh(a)) {
725 constexpr auto operator()(A a)
const ->
decltype(std::tanh(a)) {
731 constexpr auto operator()(A a)
const ->
decltype(std::round(a)) {
732 return std::round(a);
736 template <
class A,
class B>
737 constexpr auto operator()(A a, B b)
const ->
decltype(std::fmod(a, b)) {
738 return std::fmod(a, b);
742 template <
class A,
class B>
743 constexpr auto operator()(A a, B b)
const ->
decltype(std::pow(a, b)) {
744 return std::pow(a, b);
748 template <
class A,
class B>
749 constexpr auto operator()(A a, B b)
const ->
decltype(std::atan2(a, b)) {
750 return std::atan2(a, b);
754 template <
class A,
class B>
755 constexpr auto operator()(A a, B b)
const ->
decltype(std::copysign(a, b)) {
756 return std::copysign(a, b);
766 constexpr vec() : x() {}
767 constexpr vec(
const T &x_) : x(x_) {}
771 constexpr explicit vec(
const vec<U, 1> &v) :
vec(
static_cast<T
>(v.x)) {}
772 constexpr const T &operator[](
int i)
const {
return x; }
773 LINALG_CONSTEXPR14 T &operator[](
int i) {
return x; }
775 template <
class U,
class = detail::conv_t<vec, U>>
777 template <
class U,
class = detail::conv_t<U, vec>>
778 constexpr operator U()
const {
785 constexpr vec() : x(), y() {}
786 constexpr vec(
const T &x_,
const T &y_) : x(x_), y(y_) {}
787 constexpr explicit vec(
const T &s) :
vec(s, s) {}
788 constexpr explicit vec(
const T *p) :
vec(p[0], p[1]) {}
789 template <
class U,
int N>
791 :
vec(
static_cast<T
>(v.x),
static_cast<T
>(v.y)) {
794 "You must give extra arguments if your input vector is shorter.");
796 constexpr const T &operator[](
int i)
const {
return i == 0 ? x : y; }
797 LINALG_CONSTEXPR14 T &operator[](
int i) {
return i == 0 ? x : y; }
799 template <
class U,
class = detail::conv_t<vec, U>>
801 template <
class U,
class = detail::conv_t<U, vec>>
802 constexpr operator U()
const {
809 constexpr vec() : x(), y(), z() {}
810 constexpr vec(
const T &x_,
const T &y_,
const T &z_) : x(x_), y(y_), z(z_) {}
811 constexpr vec(
const vec<T, 2> &xy,
const T &z_) :
vec(xy.x, xy.y, z_) {}
812 constexpr explicit vec(
const T &s) :
vec(s, s, s) {}
813 constexpr explicit vec(
const T *p) :
vec(p[0], p[1], p[2]) {}
814 template <
class U,
int N>
816 :
vec(
static_cast<T
>(v.x),
static_cast<T
>(v.y),
static_cast<T
>(v.z)) {
819 "You must give extra arguments if your input vector is shorter.");
821 constexpr const T &operator[](
int i)
const {
822 return i == 0 ? x : i == 1 ? y : z;
824 LINALG_CONSTEXPR14 T &operator[](
int i) {
825 return i == 0 ? x : i == 1 ? y : z;
828 return *
reinterpret_cast<const vec<T, 2> *
>(
this);
832 template <
class U,
class = detail::conv_t<vec, U>>
834 template <
class U,
class = detail::conv_t<U, vec>>
835 constexpr operator U()
const {
842 constexpr vec() : x(), y(), z(), w() {}
843 constexpr vec(
const T &x_,
const T &y_,
const T &z_,
const T &w_)
844 : x(x_), y(y_), z(z_), w(w_) {}
845 constexpr vec(
const vec<T, 2> &xy,
const T &z_,
const T &w_)
846 :
vec(xy.x, xy.y, z_, w_) {}
848 :
vec(xyz.x, xyz.y, xyz.z, w_) {}
849 constexpr explicit vec(
const T &s) :
vec(s, s, s, s) {}
850 constexpr explicit vec(
const T *p) :
vec(p[0], p[1], p[2], p[3]) {}
851 template <
class U,
int N>
853 :
vec(
static_cast<T
>(v.x),
static_cast<T
>(v.y),
static_cast<T
>(v.z),
854 static_cast<T
>(v.w)) {
857 "You must give extra arguments if your input vector is shorter.");
859 constexpr const T &operator[](
int i)
const {
860 return i == 0 ? x : i == 1 ? y : i == 2 ? z : w;
862 LINALG_CONSTEXPR14 T &operator[](
int i) {
863 return i == 0 ? x : i == 1 ? y : i == 2 ? z : w;
866 return *
reinterpret_cast<const vec<T, 2> *
>(
this);
869 return *
reinterpret_cast<const vec<T, 3> *
>(
this);
874 template <
class U,
class = detail::conv_t<vec, U>>
876 template <
class U,
class = detail::conv_t<U, vec>>
877 constexpr operator U()
const {
884template <
class T,
int M>
888 constexpr mat() : x() {}
889 constexpr mat(
const V &x_) : x(x_) {}
890 constexpr explicit mat(
const T &s) : x(s) {}
891 constexpr explicit mat(
const T *p) : x(p + M * 0) {}
894 constexpr vec<T, 1> row(
int i)
const {
return {x[i]}; }
895 constexpr const V &operator[](
int j)
const {
return x; }
896 LINALG_CONSTEXPR14
V &operator[](
int j) {
return x; }
898 template <
class U,
class = detail::conv_t<mat, U>>
900 template <
class U,
class = detail::conv_t<U, mat>>
901 constexpr operator U()
const {
905template <
class T,
int M>
909 constexpr mat() : x(), y() {}
910 constexpr mat(
const V &x_,
const V &y_) : x(x_), y(y_) {}
911 constexpr explicit mat(
const T &s) : x(s), y(s) {}
912 constexpr explicit mat(
const T *p) : x(p + M * 0), y(p + M * 1) {}
913 template <
class U,
int N,
int P>
915 static_assert(P >= 2,
"Input matrix dimensions must be at least as big.");
917 constexpr vec<T, 2> row(
int i)
const {
return {x[i], y[i]}; }
918 constexpr const V &operator[](
int j)
const {
return j == 0 ? x : y; }
919 LINALG_CONSTEXPR14
V &operator[](
int j) {
return j == 0 ? x : y; }
921 template <
class U,
class = detail::conv_t<mat, U>>
923 template <
class U,
class = detail::conv_t<U, mat>>
924 constexpr operator U()
const {
928template <
class T,
int M>
932 constexpr mat() : x(), y(), z() {}
933 constexpr mat(
const V &x_,
const V &y_,
const V &z_) : x(x_), y(y_), z(z_) {}
935 : x(m_.x), y(m_.y), z(z_) {}
936 constexpr explicit mat(
const T &s) : x(s), y(s), z(s) {}
937 constexpr explicit mat(
const T *p)
938 : x(p + M * 0), y(p + M * 1), z(p + M * 2) {}
939 template <
class U,
int N,
int P>
941 static_assert(P >= 3,
"Input matrix dimensions must be at least as big.");
943 constexpr vec<T, 3> row(
int i)
const {
return {x[i], y[i], z[i]}; }
944 constexpr const V &operator[](
int j)
const {
945 return j == 0 ? x : j == 1 ? y : z;
947 LINALG_CONSTEXPR14
V &operator[](
int j) {
948 return j == 0 ? x : j == 1 ? y : z;
951 template <
class U,
class = detail::conv_t<mat, U>>
953 template <
class U,
class = detail::conv_t<U, mat>>
954 constexpr operator U()
const {
958template <
class T,
int M>
962 constexpr mat() : x(), y(), z(), w() {}
963 constexpr mat(
const V &x_,
const V &y_,
const V &z_,
const V &w_)
964 : x(x_), y(y_), z(z_), w(w_) {}
966 : x(m_.x), y(m_.y), z(m_.z), w(w_) {}
967 constexpr explicit mat(
const T &s) : x(s), y(s), z(s), w(s) {}
968 constexpr explicit mat(
const T *p)
969 : x(p + M * 0), y(p + M * 1), z(p + M * 2), w(p + M * 3) {}
970 template <
class U,
int N,
int P>
972 :
mat(
V(m.x),
V(m.y),
V(m.z),
V(m.w)) {
973 static_assert(P >= 4,
"Input matrix dimensions must be at least as big.");
976 constexpr vec<T, 4> row(
int i)
const {
return {x[i], y[i], z[i], w[i]}; }
977 constexpr const V &operator[](
int j)
const {
978 return j == 0 ? x : j == 1 ? y : j == 2 ? z : w;
980 LINALG_CONSTEXPR14
V &operator[](
int j) {
981 return j == 0 ? x : j == 1 ? y : j == 2 ? z : w;
984 template <
class U,
class = detail::conv_t<mat, U>>
986 template <
class U,
class = detail::conv_t<U, mat>>
987 constexpr operator U()
const {
1004 return {{1, 0}, {0, 1}};
1010 return {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
1016 return {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
1024template <
class F,
class A,
class B>
1025constexpr A fold(F f, A a,
const vec<B, 1> &b) {
1028template <
class F,
class A,
class B>
1029constexpr A fold(F f, A a,
const vec<B, 2> &b) {
1030 return f(f(a, b.x), b.y);
1032template <
class F,
class A,
class B>
1033constexpr A fold(F f, A a,
const vec<B, 3> &b) {
1034 return f(f(f(a, b.x), b.y), b.z);
1036template <
class F,
class A,
class B>
1037constexpr A fold(F f, A a,
const vec<B, 4> &b) {
1038 return f(f(f(f(a, b.x), b.y), b.z), b.w);
1040template <
class F,
class A,
class B,
int M>
1041constexpr A fold(F f, A a,
const mat<B, M, 1> &b) {
1042 return fold(f, a, b.x);
1044template <
class F,
class A,
class B,
int M>
1045constexpr A fold(F f, A a,
const mat<B, M, 2> &b) {
1046 return fold(f, fold(f, a, b.x), b.y);
1048template <
class F,
class A,
class B,
int M>
1049constexpr A fold(F f, A a,
const mat<B, M, 3> &b) {
1050 return fold(f, fold(f, fold(f, a, b.x), b.y), b.z);
1052template <
class F,
class A,
class B,
int M>
1053constexpr A fold(F f, A a,
const mat<B, M, 4> &b) {
1054 return fold(f, fold(f, fold(f, fold(f, a, b.x), b.y), b.z), b.w);
1059template <
class F,
class... A>
1060using apply_t =
typename detail::apply<F, void, A...>::type;
1061template <
class F,
class... A>
1062using mm_apply_t =
typename std::enable_if<detail::apply<F, void, A...>::mm,
1063 apply_t<F, A...>>::type;
1064template <
class F,
class... A>
1065using no_mm_apply_t =
typename std::enable_if<!detail::apply<F, void, A...>::mm,
1066 apply_t<F, A...>>::type;
1069 typename detail::scalar_type<A>::type;
1074template <
class F,
class... A>
1075constexpr apply_t<F, A...> apply(F func,
const A &...args) {
1076 return detail::apply<F, void, A...>::impl(
1077 detail::make_seq<0, detail::apply<F, void, A...>::size>{}, func, args...);
1081template <
class A,
class F>
1082constexpr apply_t<F, A> map(
const A &a, F func) {
1083 return apply(func, a);
1087template <
class A,
class B,
class F>
1088constexpr apply_t<F, A, B> zip(
const A &a,
const B &b, F func) {
1089 return apply(func, a, b);
1094template <
class A,
class B>
1095constexpr typename detail::any_compare<A, B>::type compare(
const A &a,
1097 return detail::any_compare<A, B>()(a, b);
1099template <
class A,
class B>
1100constexpr auto operator==(
const A &a,
1101 const B &b) ->
decltype(compare(a, b) == 0) {
1102 return compare(a, b) == 0;
1104template <
class A,
class B>
1105constexpr auto operator!=(
const A &a,
1106 const B &b) ->
decltype(compare(a, b) != 0) {
1107 return compare(a, b) != 0;
1109template <
class A,
class B>
1110constexpr auto operator<(
const A &a,
1111 const B &b) ->
decltype(compare(a, b) < 0) {
1112 return compare(a, b) < 0;
1114template <
class A,
class B>
1115constexpr auto operator>(
const A &a,
1116 const B &b) ->
decltype(compare(a, b) > 0) {
1117 return compare(a, b) > 0;
1119template <
class A,
class B>
1120constexpr auto operator<=(
const A &a,
1121 const B &b) ->
decltype(compare(a, b) <= 0) {
1122 return compare(a, b) <= 0;
1124template <
class A,
class B>
1125constexpr auto operator>=(
const A &a,
1126 const B &b) ->
decltype(compare(a, b) >= 0) {
1127 return compare(a, b) >= 0;
1132constexpr bool any(
const A &a) {
1133 return fold(detail::op_or{},
false, a);
1136constexpr bool all(
const A &a) {
1137 return fold(detail::op_and{},
true, a);
1140constexpr scalar_t<A> sum(
const A &a) {
1141 return fold(detail::op_add{}, scalar_t<A>(0), a);
1144constexpr scalar_t<A> product(
const A &a) {
1145 return fold(detail::op_mul{}, scalar_t<A>(1), a);
1148constexpr scalar_t<A> minelem(
const A &a) {
1149 return fold(detail::min{}, a.x, a);
1152constexpr scalar_t<A> maxelem(
const A &a) {
1153 return fold(detail::max{}, a.x, a);
1155template <
class T,
int M>
1156int argmin(
const vec<T, M> &a) {
1158 for (
int i = 1; i < M; ++i)
1159 if (a[i] < a[j]) j = i;
1162template <
class T,
int M>
1163int argmax(
const vec<T, M> &a) {
1165 for (
int i = 1; i < M; ++i)
1166 if (a[i] > a[j]) j = i;
1172constexpr apply_t<detail::op_pos, A> operator+(
const A &a) {
1173 return apply(detail::op_pos{}, a);
1176constexpr apply_t<detail::op_neg, A> operator-(
const A &a) {
1177 return apply(detail::op_neg{}, a);
1180constexpr apply_t<detail::op_cmp, A> operator~(
const A &a) {
1181 return apply(detail::op_cmp{}, a);
1184constexpr apply_t<detail::op_not, A> operator!(
const A &a) {
1185 return apply(detail::op_not{}, a);
1190template <
class A,
class B>
1191constexpr apply_t<detail::op_add, A, B> operator+(
const A &a,
const B &b) {
1192 return apply(detail::op_add{}, a, b);
1194template <
class A,
class B>
1195constexpr apply_t<detail::op_sub, A, B> operator-(
const A &a,
const B &b) {
1196 return apply(detail::op_sub{}, a, b);
1198template <
class A,
class B>
1199constexpr apply_t<detail::op_mul, A, B> cmul(
const A &a,
const B &b) {
1200 return apply(detail::op_mul{}, a, b);
1202template <
class A,
class B>
1203constexpr apply_t<detail::op_div, A, B> operator/(
const A &a,
const B &b) {
1204 return apply(detail::op_div{}, a, b);
1206template <
class A,
class B>
1207constexpr apply_t<detail::op_mod, A, B> operator%(
const A &a,
const B &b) {
1208 return apply(detail::op_mod{}, a, b);
1210template <
class A,
class B>
1211constexpr apply_t<detail::op_un, A, B> operator|(
const A &a,
const B &b) {
1212 return apply(detail::op_un{}, a, b);
1214template <
class A,
class B>
1215constexpr apply_t<detail::op_xor, A, B> operator^(
const A &a,
const B &b) {
1216 return apply(detail::op_xor{}, a, b);
1218template <
class A,
class B>
1219constexpr apply_t<detail::op_int, A, B> operator&(
const A &a,
const B &b) {
1220 return apply(detail::op_int{}, a, b);
1222template <
class A,
class B>
1223constexpr apply_t<detail::op_lsh, A, B> operator<<(
const A &a,
const B &b) {
1224 return apply(detail::op_lsh{}, a, b);
1226template <
class A,
class B>
1227constexpr apply_t<detail::op_rsh, A, B> operator>>(
const A &a,
const B &b) {
1228 return apply(detail::op_rsh{}, a, b);
1233template <
class A,
class B>
1234constexpr auto operator*(
const A &a,
const B &b) {
1240template <
class A,
class B>
1241constexpr auto operator+=(A &a,
const B &b) ->
decltype(a = a + b) {
1244template <
class A,
class B>
1245constexpr auto operator-=(A &a,
const B &b) ->
decltype(a = a - b) {
1248template <
class A,
class B>
1249constexpr auto operator*=(A &a,
const B &b) ->
decltype(a = a * b) {
1252template <
class A,
class B>
1253constexpr auto operator/=(A &a,
const B &b) ->
decltype(a = a / b) {
1256template <
class A,
class B>
1257constexpr auto operator%=(A &a,
const B &b) ->
decltype(a = a % b) {
1260template <
class A,
class B>
1261constexpr auto operator|=(A &a,
const B &b) ->
decltype(a = a | b) {
1264template <
class A,
class B>
1265constexpr auto operator^=(A &a,
const B &b) ->
decltype(a = a ^ b) {
1268template <
class A,
class B>
1269constexpr auto operator&=(A &a,
const B &b) ->
decltype(a = a & b) {
1272template <
class A,
class B>
1273constexpr auto operator<<=(A &a,
const B &b) ->
decltype(a = a << b) {
1276template <
class A,
class B>
1277constexpr auto operator>>=(A &a,
const B &b) ->
decltype(a = a >> b) {
1282template <
int... I,
class T,
int M>
1283constexpr vec<T,
sizeof...(I)> swizzle(
const vec<T, M> &a) {
1284 return {detail::getter<I>{}(a)...};
1286template <
int I0,
int I1,
class T,
int M>
1287constexpr vec<T, I1 - I0> subvec(
const vec<T, M> &a) {
1288 return detail::swizzle(a, detail::make_seq<I0, I1>{});
1290template <
int I0,
int J0,
int I1,
int J1,
class T,
int M,
int N>
1291constexpr mat<T, I1 - I0, J1 - J0> submat(
const mat<T, M, N> &a) {
1292 return detail::swizzle(a, detail::make_seq<I0, I1>{},
1293 detail::make_seq<J0, J1>{});
1298constexpr apply_t<detail::std_isfinite, A> isfinite(
const A &a) {
1299 return apply(detail::std_isfinite{}, a);
1302constexpr apply_t<detail::std_abs, A> abs(
const A &a) {
1303 return apply(detail::std_abs{}, a);
1306constexpr apply_t<detail::std_floor, A> floor(
const A &a) {
1307 return apply(detail::std_floor{}, a);
1310constexpr apply_t<detail::std_ceil, A> ceil(
const A &a) {
1311 return apply(detail::std_ceil{}, a);
1314constexpr apply_t<detail::std_exp, A> exp(
const A &a) {
1315 return apply(detail::std_exp{}, a);
1318constexpr apply_t<detail::std_log, A> log(
const A &a) {
1319 return apply(detail::std_log{}, a);
1322constexpr apply_t<detail::std_log2, A> log2(
const A &a) {
1323 return apply(detail::std_log2{}, a);
1326constexpr apply_t<detail::std_log10, A> log10(
const A &a) {
1327 return apply(detail::std_log10{}, a);
1330constexpr apply_t<detail::std_sqrt, A> sqrt(
const A &a) {
1331 return apply(detail::std_sqrt{}, a);
1334constexpr apply_t<detail::std_sin, A> sin(
const A &a) {
1335 return apply(detail::std_sin{}, a);
1338constexpr apply_t<detail::std_cos, A> cos(
const A &a) {
1339 return apply(detail::std_cos{}, a);
1342constexpr apply_t<detail::std_tan, A> tan(
const A &a) {
1343 return apply(detail::std_tan{}, a);
1346constexpr apply_t<detail::std_asin, A> asin(
const A &a) {
1347 return apply(detail::std_asin{}, a);
1350constexpr apply_t<detail::std_acos, A> acos(
const A &a) {
1351 return apply(detail::std_acos{}, a);
1354constexpr apply_t<detail::std_atan, A> atan(
const A &a) {
1355 return apply(detail::std_atan{}, a);
1358constexpr apply_t<detail::std_sinh, A> sinh(
const A &a) {
1359 return apply(detail::std_sinh{}, a);
1362constexpr apply_t<detail::std_cosh, A> cosh(
const A &a) {
1363 return apply(detail::std_cosh{}, a);
1366constexpr apply_t<detail::std_tanh, A> tanh(
const A &a) {
1367 return apply(detail::std_tanh{}, a);
1370constexpr apply_t<detail::std_round, A> round(
const A &a) {
1371 return apply(detail::std_round{}, a);
1374template <
class A,
class B>
1375constexpr apply_t<detail::std_fmod, A, B> fmod(
const A &a,
const B &b) {
1376 return apply(detail::std_fmod{}, a, b);
1378template <
class A,
class B>
1379constexpr apply_t<detail::std_pow, A, B> pow(
const A &a,
const B &b) {
1380 return apply(detail::std_pow{}, a, b);
1382template <
class A,
class B>
1383constexpr apply_t<detail::std_atan2, A, B> atan2(
const A &a,
const B &b) {
1384 return apply(detail::std_atan2{}, a, b);
1386template <
class A,
class B>
1387constexpr apply_t<detail::std_copysign, A, B> copysign(
const A &a,
const B &b) {
1388 return apply(detail::std_copysign{}, a, b);
1392template <
class A,
class B>
1393constexpr apply_t<detail::op_eq, A, B> equal(
const A &a,
const B &b) {
1394 return apply(detail::op_eq{}, a, b);
1396template <
class A,
class B>
1397constexpr apply_t<detail::op_ne, A, B> nequal(
const A &a,
const B &b) {
1398 return apply(detail::op_ne{}, a, b);
1400template <
class A,
class B>
1401constexpr apply_t<detail::op_lt, A, B> less(
const A &a,
const B &b) {
1402 return apply(detail::op_lt{}, a, b);
1404template <
class A,
class B>
1405constexpr apply_t<detail::op_gt, A, B> greater(
const A &a,
const B &b) {
1406 return apply(detail::op_gt{}, a, b);
1408template <
class A,
class B>
1409constexpr apply_t<detail::op_le, A, B> lequal(
const A &a,
const B &b) {
1410 return apply(detail::op_le{}, a, b);
1412template <
class A,
class B>
1413constexpr apply_t<detail::op_ge, A, B> gequal(
const A &a,
const B &b) {
1414 return apply(detail::op_ge{}, a, b);
1418template <
class A,
class B>
1419constexpr apply_t<detail::min, A, B> min(
const A &a,
const B &b) {
1420 return apply(detail::min{}, a, b);
1422template <
class A,
class B>
1423constexpr apply_t<detail::max, A, B> max(
const A &a,
const B &b) {
1424 return apply(detail::max{}, a, b);
1426template <
class X,
class L,
class H>
1427constexpr apply_t<detail::clamp, X, L, H> clamp(
const X &x,
const L &l,
1429 return apply(detail::clamp{}, x, l, h);
1431template <
class P,
class A,
class B>
1432constexpr apply_t<detail::select, P, A, B> select(
const P &p,
const A &a,
1434 return apply(detail::select{}, p, a, b);
1436template <
class A,
class B,
class T>
1437constexpr apply_t<detail::lerp, A, B, T> lerp(
const A &a,
const B &b,
1439 return apply(detail::lerp{}, a, b, t);
1444constexpr T cross(
const vec<T, 2> &a,
const vec<T, 2> &b) {
1445 return a.x * b.y - a.y * b.x;
1448constexpr vec<T, 2> cross(T a,
const vec<T, 2> &b) {
1449 return {-a * b.y, a * b.x};
1452constexpr vec<T, 2> cross(
const vec<T, 2> &a, T b) {
1453 return {a.y * b, -a.x * b};
1456constexpr vec<T, 3> cross(
const vec<T, 3> &a,
const vec<T, 3> &b) {
1457 return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
1459template <
class T,
int M>
1460constexpr T dot(
const vec<T, M> &a,
const vec<T, M> &b) {
1463template <
class T,
int M>
1464constexpr T length2(
const vec<T, M> &a) {
1467template <
class T,
int M>
1468T length(
const vec<T, M> &a) {
1469 return std::sqrt(length2(a));
1471template <
class T,
int M>
1472vec<T, M> normalize(
const vec<T, M> &a) {
1473 return a / length(a);
1475template <
class T,
int M>
1476constexpr T distance2(
const vec<T, M> &a,
const vec<T, M> &b) {
1477 return length2(b - a);
1479template <
class T,
int M>
1480T distance(
const vec<T, M> &a,
const vec<T, M> &b) {
1481 return length(b - a);
1483template <
class T,
int M>
1484T uangle(
const vec<T, M> &a,
const vec<T, M> &b) {
1486 return d > 1 ? 0 : std::acos(d < -1 ? -1 : d);
1488template <
class T,
int M>
1489T angle(
const vec<T, M> &a,
const vec<T, M> &b) {
1490 return uangle(normalize(a), normalize(b));
1493vec<T, 2> rot(T a,
const vec<T, 2> &v) {
1494 const T s = std::sin(a), c = std::cos(a);
1495 return {v.x * c - v.y * s, v.x * s + v.y * c};
1498vec<T, 3> rotx(T a,
const vec<T, 3> &v) {
1499 const T s = std::sin(a), c = std::cos(a);
1500 return {v.x, v.y * c - v.z * s, v.y * s + v.z * c};
1503vec<T, 3> roty(T a,
const vec<T, 3> &v) {
1504 const T s = std::sin(a), c = std::cos(a);
1505 return {v.x * c + v.z * s, v.y, -v.x * s + v.z * c};
1508vec<T, 3> rotz(T a,
const vec<T, 3> &v) {
1509 const T s = std::sin(a), c = std::cos(a);
1510 return {v.x * c - v.y * s, v.x * s + v.y * c, v.z};
1512template <
class T,
int M>
1513vec<T, M> nlerp(
const vec<T, M> &a,
const vec<T, M> &b, T t) {
1514 return normalize(lerp(a, b, t));
1516template <
class T,
int M>
1517vec<T, M> slerp(
const vec<T, M> &a,
const vec<T, M> &b, T t) {
1518 T th = uangle(a, b);
1520 : a * (std::sin(th * (1 - t)) / std::sin(th)) +
1521 b * (std::sin(th * t) / std::sin(th));
1527constexpr vec<T, 4> qconj(
const vec<T, 4> &q) {
1528 return {-q.x, -q.y, -q.z, q.w};
1531vec<T, 4> qinv(
const vec<T, 4> &q) {
1532 return qconj(q) / length2(q);
1535vec<T, 4> qexp(
const vec<T, 4> &q) {
1536 const auto v = q.xyz();
1537 const auto vv = length(v);
1538 return std::exp(q.w) *
1539 vec<T, 4>{v * (vv > 0 ? std::sin(vv) / vv : 0), std::cos(vv)};
1542vec<T, 4> qlog(
const vec<T, 4> &q) {
1543 const auto v = q.xyz();
1544 const auto vv = length(v), qq = length(q);
1545 return {v * (vv > 0 ? std::acos(q.w / qq) / vv : 0), std::log(qq)};
1548vec<T, 4> qpow(
const vec<T, 4> &q,
const T &p) {
1549 const auto v = q.xyz();
1550 const auto vv = length(v), qq = length(q), th = std::acos(q.w / qq);
1551 return std::pow(qq, p) *
1552 vec<T, 4>{v * (vv > 0 ? std::sin(p * th) / vv : 0), std::cos(p * th)};
1555constexpr vec<T, 4> qmul(
const vec<T, 4> &a,
const vec<T, 4> &b) {
1556 return {a.x * b.w + a.w * b.x + a.y * b.z - a.z * b.y,
1557 a.y * b.w + a.w * b.y + a.z * b.x - a.x * b.z,
1558 a.z * b.w + a.w * b.z + a.x * b.y - a.y * b.x,
1559 a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z};
1561template <
class T,
class... R>
1562constexpr vec<T, 4> qmul(
const vec<T, 4> &a, R... r) {
1563 return qmul(a, qmul(r...));
1569constexpr vec<T, 3> qxdir(
const vec<T, 4> &q) {
1570 return {q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
1571 (q.x * q.y + q.z * q.w) * 2, (q.z * q.x - q.y * q.w) * 2};
1574constexpr vec<T, 3> qydir(
const vec<T, 4> &q) {
1575 return {(q.x * q.y - q.z * q.w) * 2,
1576 q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
1577 (q.y * q.z + q.x * q.w) * 2};
1580constexpr vec<T, 3> qzdir(
const vec<T, 4> &q) {
1581 return {(q.z * q.x + q.y * q.w) * 2, (q.y * q.z - q.x * q.w) * 2,
1582 q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z};
1585constexpr mat<T, 3, 3> qmat(
const vec<T, 4> &q) {
1586 return {qxdir(q), qydir(q), qzdir(q)};
1589constexpr vec<T, 3> qrot(
const vec<T, 4> &q,
const vec<T, 3> &v) {
1590 return qxdir(q) * v.x + qydir(q) * v.y + qzdir(q) * v.z;
1593T qangle(
const vec<T, 4> &q) {
1594 return std::atan2(length(q.xyz()), q.w) * 2;
1597vec<T, 3> qaxis(
const vec<T, 4> &q) {
1598 return normalize(q.xyz());
1601vec<T, 4> qnlerp(
const vec<T, 4> &a,
const vec<T, 4> &b, T t) {
1602 return nlerp(a, dot(a, b) < 0 ? -b : b, t);
1605vec<T, 4> qslerp(
const vec<T, 4> &a,
const vec<T, 4> &b, T t) {
1606 return slerp(a, dot(a, b) < 0 ? -b : b, t);
1610template <
class T,
int M>
1611constexpr vec<T, M> mul(
const vec<T, M> &a,
const T &b) {
1614template <
class T,
int M>
1615constexpr vec<T, M> mul(
const T &b,
const vec<T, M> &a) {
1618template <
class T,
int M,
int N>
1619constexpr mat<T, M, N> mul(
const mat<T, M, N> &a,
const T &b) {
1622template <
class T,
int M,
int N>
1623constexpr mat<T, M, N> mul(
const T &b,
const mat<T, M, N> &a) {
1626template <
class T,
int M>
1627constexpr vec<T, M> mul(
const vec<T, M> &a,
const vec<T, M> &b) {
1630template <
class T,
int M>
1631constexpr vec<T, M> mul(
const mat<T, M, 1> &a,
const vec<T, 1> &b) {
1634template <
class T,
int M>
1635constexpr vec<T, M> mul(
const mat<T, M, 2> &a,
const vec<T, 2> &b) {
1636 return a.x * b.x + a.y * b.y;
1638template <
class T,
int M>
1639constexpr vec<T, M> mul(
const mat<T, M, 3> &a,
const vec<T, 3> &b) {
1640 return a.x * b.x + a.y * b.y + a.z * b.z;
1642template <
class T,
int M>
1643constexpr vec<T, M> mul(
const mat<T, M, 4> &a,
const vec<T, 4> &b) {
1644 return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
1646template <
class T,
int M,
int N>
1647constexpr mat<T, M, 1> mul(
const mat<T, M, N> &a,
const mat<T, N, 1> &b) {
1648 return {mul(a, b.x)};
1650template <
class T,
int M,
int N>
1651constexpr mat<T, M, 2> mul(
const mat<T, M, N> &a,
const mat<T, N, 2> &b) {
1652 return {mul(a, b.x), mul(a, b.y)};
1654template <
class T,
int M,
int N>
1655constexpr mat<T, M, 3> mul(
const mat<T, M, N> &a,
const mat<T, N, 3> &b) {
1656 return {mul(a, b.x), mul(a, b.y), mul(a, b.z)};
1658template <
class T,
int M,
int N>
1659constexpr mat<T, M, 4> mul(
const mat<T, M, N> &a,
const mat<T, N, 4> &b) {
1660 return {mul(a, b.x), mul(a, b.y), mul(a, b.z), mul(a, b.w)};
1662template <
class T,
int M,
int N,
int P>
1663constexpr vec<T, M> mul(
const mat<T, M, N> &a,
const mat<T, N, P> &b,
1664 const vec<T, P> &c) {
1665 return mul(mul(a, b), c);
1667template <
class T,
int M,
int N,
int P,
int Q>
1668constexpr mat<T, M, Q> mul(
const mat<T, M, N> &a,
const mat<T, N, P> &b,
1669 const mat<T, P, Q> &c) {
1670 return mul(mul(a, b), c);
1672template <
class T,
int M,
int N,
int P,
int Q>
1673constexpr vec<T, M> mul(
const mat<T, M, N> &a,
const mat<T, N, P> &b,
1674 const mat<T, P, Q> &c,
const vec<T, Q> &d) {
1675 return mul(mul(a, b, c), d);
1677template <
class T,
int M,
int N,
int P,
int Q,
int R>
1678constexpr mat<T, M, R> mul(
const mat<T, M, N> &a,
const mat<T, N, P> &b,
1679 const mat<T, P, Q> &c,
const mat<T, Q, R> &d) {
1680 return mul(mul(a, b, c), d);
1682template <
class T,
int M>
1683constexpr mat<T, M, 1> outerprod(
const vec<T, M> &a,
const vec<T, 1> &b) {
1686template <
class T,
int M>
1687constexpr mat<T, M, 2> outerprod(
const vec<T, M> &a,
const vec<T, 2> &b) {
1688 return {a * b.x, a * b.y};
1690template <
class T,
int M>
1691constexpr mat<T, M, 3> outerprod(
const vec<T, M> &a,
const vec<T, 3> &b) {
1692 return {a * b.x, a * b.y, a * b.z};
1694template <
class T,
int M>
1695constexpr mat<T, M, 4> outerprod(
const vec<T, M> &a,
const vec<T, 4> &b) {
1696 return {a * b.x, a * b.y, a * b.z, a * b.w};
1699constexpr vec<T, 1> diagonal(
const mat<T, 1, 1> &a) {
1703constexpr vec<T, 2> diagonal(
const mat<T, 2, 2> &a) {
1704 return {a.x.x, a.y.y};
1707constexpr vec<T, 3> diagonal(
const mat<T, 3, 3> &a) {
1708 return {a.x.x, a.y.y, a.z.z};
1711constexpr vec<T, 4> diagonal(
const mat<T, 4, 4> &a) {
1712 return {a.x.x, a.y.y, a.z.z, a.w.w};
1714template <
class T,
int N>
1715constexpr T trace(
const mat<T, N, N> &a) {
1716 return sum(diagonal(a));
1718template <
class T,
int M>
1719constexpr mat<T, M, 1> transpose(
const mat<T, 1, M> &m) {
1722template <
class T,
int M>
1723constexpr mat<T, M, 2> transpose(
const mat<T, 2, M> &m) {
1724 return {m.row(0), m.row(1)};
1726template <
class T,
int M>
1727constexpr mat<T, M, 3> transpose(
const mat<T, 3, M> &m) {
1728 return {m.row(0), m.row(1), m.row(2)};
1730template <
class T,
int M>
1731constexpr mat<T, M, 4> transpose(
const mat<T, 4, M> &m) {
1732 return {m.row(0), m.row(1), m.row(2), m.row(3)};
1734template <
class T,
int M>
1735constexpr mat<T, 1, M> transpose(
const vec<T, M> &m) {
1736 return transpose(mat<T, M, 1>(m));
1739constexpr mat<T, 1, 1> adjugate(
const mat<T, 1, 1> &a) {
1740 return {vec<T, 1>{1}};
1743constexpr mat<T, 2, 2> adjugate(
const mat<T, 2, 2> &a) {
1744 return {{a.y.y, -a.x.y}, {-a.y.x, a.x.x}};
1747constexpr mat<T, 3, 3> adjugate(
const mat<T, 3, 3> &a);
1749constexpr mat<T, 4, 4> adjugate(
const mat<T, 4, 4> &a);
1750template <
class T,
int N>
1751constexpr mat<T, N, N> comatrix(
const mat<T, N, N> &a) {
1752 return transpose(adjugate(a));
1755constexpr T determinant(
const mat<T, 1, 1> &a) {
1759constexpr T determinant(
const mat<T, 2, 2> &a) {
1760 return a.x.x * a.y.y - a.x.y * a.y.x;
1763constexpr T determinant(
const mat<T, 3, 3> &a) {
1764 return a.x.x * (a.y.y * a.z.z - a.z.y * a.y.z) +
1765 a.x.y * (a.y.z * a.z.x - a.z.z * a.y.x) +
1766 a.x.z * (a.y.x * a.z.y - a.z.x * a.y.y);
1769constexpr T determinant(
const mat<T, 4, 4> &a);
1770template <
class T,
int N>
1771constexpr mat<T, N, N> inverse(
const mat<T, N, N> &a) {
1772 return adjugate(a) / determinant(a);
1776template <
class T,
int M>
1777T *begin(vec<T, M> &a) {
1780template <
class T,
int M>
1781const T *begin(
const vec<T, M> &a) {
1784template <
class T,
int M>
1785T *end(vec<T, M> &a) {
1786 return begin(a) + M;
1788template <
class T,
int M>
1789const T *end(
const vec<T, M> &a) {
1790 return begin(a) + M;
1792template <
class T,
int M,
int N>
1793vec<T, M> *begin(mat<T, M, N> &a) {
1796template <
class T,
int M,
int N>
1797const vec<T, M> *begin(
const mat<T, M, N> &a) {
1800template <
class T,
int M,
int N>
1801vec<T, M> *end(mat<T, M, N> &a) {
1802 return begin(a) + N;
1804template <
class T,
int M,
int N>
1805const vec<T, M> *end(
const mat<T, M, N> &a) {
1806 return begin(a) + N;
1821vec<T, 4>
constexpr rotation_quat(
const vec<T, 3> &axis, T angle) {
1822 return {axis * std::sin(angle / 2), std::cos(angle / 2)};
1825vec<T, 4>
constexpr rotation_quat(
const vec<T, 3> &orig,
1826 const vec<T, 3> &dest) {
1827 T cosTheta = dot(orig, dest);
1828 if (cosTheta >= 1 - std::numeric_limits<T>::epsilon()) {
1829 return {0, 0, 0, 1};
1831 if (cosTheta < -1 + std::numeric_limits<T>::epsilon()) {
1832 vec<T, 3> axis = cross(vec<T, 3>(0, 0, 1), orig);
1833 if (length2(axis) < std::numeric_limits<T>::epsilon())
1834 axis = cross(vec<T, 3>(1, 0, 0), orig);
1835 return rotation_quat(normalize(axis),
1836 3.14159265358979323846264338327950288);
1838 vec<T, 3> axis = cross(orig, dest);
1839 T s = std::sqrt((1 + cosTheta) * 2);
1840 return {axis * (1 / s), s * 0.5};
1843vec<T, 4> rotation_quat(
const mat<T, 3, 3> &m);
1845mat<T, 4, 4> translation_matrix(
const vec<T, 3> &translation) {
1846 return {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {translation, 1}};
1849mat<T, 4, 4> rotation_matrix(
const vec<T, 4> &rotation) {
1850 return {{qxdir(rotation), 0},
1851 {qydir(rotation), 0},
1852 {qzdir(rotation), 0},
1856mat<T, 4, 4> scaling_matrix(
const vec<T, 3> &scaling) {
1857 return {{scaling.x, 0, 0, 0},
1858 {0, scaling.y, 0, 0},
1859 {0, 0, scaling.z, 0},
1863mat<T, 4, 4> pose_matrix(
const vec<T, 4> &q,
const vec<T, 3> &p) {
1864 return {{qxdir(q), 0}, {qydir(q), 0}, {qzdir(q), 0}, {p, 1}};
1867mat<T, 4, 4> lookat_matrix(
const vec<T, 3> &eye,
const vec<T, 3> ¢er,
1868 const vec<T, 3> &view_y_dir, fwd_axis fwd = neg_z);
1870mat<T, 4, 4> frustum_matrix(T x0, T x1, T y0, T y1, T n, T f,
1871 fwd_axis a = neg_z, z_range z = neg_one_to_one);
1873mat<T, 4, 4> perspective_matrix(T fovy, T aspect, T n, T f, fwd_axis a = neg_z,
1874 z_range z = neg_one_to_one) {
1875 T y = n * std::tan(fovy / 2), x = y * aspect;
1876 return frustum_matrix(-x, x, -y, y, n, f, a, z);
1882 vec<T, 1> operator()(
const std::array<T, 1> &a)
const {
return {a[0]}; }
1886 vec<T, 2> operator()(
const std::array<T, 2> &a)
const {
return {a[0], a[1]}; }
1890 vec<T, 3> operator()(
const std::array<T, 3> &a)
const {
1891 return {a[0], a[1], a[2]};
1896 vec<T, 4> operator()(
const std::array<T, 4> &a)
const {
1897 return {a[0], a[1], a[2], a[3]};
1903 std::array<T, 1> operator()(
const vec<T, 1> &a)
const {
return {a[0]}; }
1907 std::array<T, 2> operator()(
const vec<T, 2> &a)
const {
return {a[0], a[1]}; }
1911 std::array<T, 3> operator()(
const vec<T, 3> &a)
const {
1912 return {a[0], a[1], a[2]};
1917 std::array<T, 4> operator()(
const vec<T, 4> &a)
const {
1918 return {a[0], a[1], a[2], a[3]};
1926struct hash<linalg::vec<T, 1>> {
1933struct hash<linalg::vec<T, 2>> {
1936 return h(v.x) ^ (h(v.y) << 1);
1940struct hash<linalg::vec<T, 3>> {
1943 return h(v.x) ^ (h(v.y) << 1) ^ (h(v.z) << 2);
1947struct hash<linalg::vec<T, 4>> {
1950 return h(v.x) ^ (h(v.y) << 1) ^ (h(v.z) << 2) ^ (h(v.w) << 3);
1954template <
class T,
int M>
1955struct hash<linalg::mat<T, M, 1>> {
1957 std::hash<linalg::vec<T, M>> h;
1961template <
class T,
int M>
1962struct hash<linalg::mat<T, M, 2>> {
1964 std::hash<linalg::vec<T, M>> h;
1965 return h(v.x) ^ (h(v.y) << M);
1968template <
class T,
int M>
1969struct hash<linalg::mat<T, M, 3>> {
1971 std::hash<linalg::vec<T, M>> h;
1972 return h(v.x) ^ (h(v.y) << M) ^ (h(v.z) << (M * 2));
1975template <
class T,
int M>
1976struct hash<linalg::mat<T, M, 4>> {
1978 std::hash<linalg::vec<T, M>> h;
1979 return h(v.x) ^ (h(v.y) << M) ^ (h(v.z) << (M * 2)) ^ (h(v.w) << (M * 3));
1987 return {{a.y.y * a.z.z - a.z.y * a.y.z, a.z.y * a.x.z - a.x.y * a.z.z,
1988 a.x.y * a.y.z - a.y.y * a.x.z},
1989 {a.y.z * a.z.x - a.z.z * a.y.x, a.z.z * a.x.x - a.x.z * a.z.x,
1990 a.x.z * a.y.x - a.y.z * a.x.x},
1991 {a.y.x * a.z.y - a.z.x * a.y.y, a.z.x * a.x.y - a.x.x * a.z.y,
1992 a.x.x * a.y.y - a.y.x * a.x.y}};
1997 return {{a.y.y * a.z.z * a.w.w + a.w.y * a.y.z * a.z.w +
1998 a.z.y * a.w.z * a.y.w - a.y.y * a.w.z * a.z.w -
1999 a.z.y * a.y.z * a.w.w - a.w.y * a.z.z * a.y.w,
2000 a.x.y * a.w.z * a.z.w + a.z.y * a.x.z * a.w.w +
2001 a.w.y * a.z.z * a.x.w - a.w.y * a.x.z * a.z.w -
2002 a.z.y * a.w.z * a.x.w - a.x.y * a.z.z * a.w.w,
2003 a.x.y * a.y.z * a.w.w + a.w.y * a.x.z * a.y.w +
2004 a.y.y * a.w.z * a.x.w - a.x.y * a.w.z * a.y.w -
2005 a.y.y * a.x.z * a.w.w - a.w.y * a.y.z * a.x.w,
2006 a.x.y * a.z.z * a.y.w + a.y.y * a.x.z * a.z.w +
2007 a.z.y * a.y.z * a.x.w - a.x.y * a.y.z * a.z.w -
2008 a.z.y * a.x.z * a.y.w - a.y.y * a.z.z * a.x.w},
2009 {a.y.z * a.w.w * a.z.x + a.z.z * a.y.w * a.w.x +
2010 a.w.z * a.z.w * a.y.x - a.y.z * a.z.w * a.w.x -
2011 a.w.z * a.y.w * a.z.x - a.z.z * a.w.w * a.y.x,
2012 a.x.z * a.z.w * a.w.x + a.w.z * a.x.w * a.z.x +
2013 a.z.z * a.w.w * a.x.x - a.x.z * a.w.w * a.z.x -
2014 a.z.z * a.x.w * a.w.x - a.w.z * a.z.w * a.x.x,
2015 a.x.z * a.w.w * a.y.x + a.y.z * a.x.w * a.w.x +
2016 a.w.z * a.y.w * a.x.x - a.x.z * a.y.w * a.w.x -
2017 a.w.z * a.x.w * a.y.x - a.y.z * a.w.w * a.x.x,
2018 a.x.z * a.y.w * a.z.x + a.z.z * a.x.w * a.y.x +
2019 a.y.z * a.z.w * a.x.x - a.x.z * a.z.w * a.y.x -
2020 a.y.z * a.x.w * a.z.x - a.z.z * a.y.w * a.x.x},
2021 {a.y.w * a.z.x * a.w.y + a.w.w * a.y.x * a.z.y +
2022 a.z.w * a.w.x * a.y.y - a.y.w * a.w.x * a.z.y -
2023 a.z.w * a.y.x * a.w.y - a.w.w * a.z.x * a.y.y,
2024 a.x.w * a.w.x * a.z.y + a.z.w * a.x.x * a.w.y +
2025 a.w.w * a.z.x * a.x.y - a.x.w * a.z.x * a.w.y -
2026 a.w.w * a.x.x * a.z.y - a.z.w * a.w.x * a.x.y,
2027 a.x.w * a.y.x * a.w.y + a.w.w * a.x.x * a.y.y +
2028 a.y.w * a.w.x * a.x.y - a.x.w * a.w.x * a.y.y -
2029 a.y.w * a.x.x * a.w.y - a.w.w * a.y.x * a.x.y,
2030 a.x.w * a.z.x * a.y.y + a.y.w * a.x.x * a.z.y +
2031 a.z.w * a.y.x * a.x.y - a.x.w * a.y.x * a.z.y -
2032 a.z.w * a.x.x * a.y.y - a.y.w * a.z.x * a.x.y},
2033 {a.y.x * a.w.y * a.z.z + a.z.x * a.y.y * a.w.z +
2034 a.w.x * a.z.y * a.y.z - a.y.x * a.z.y * a.w.z -
2035 a.w.x * a.y.y * a.z.z - a.z.x * a.w.y * a.y.z,
2036 a.x.x * a.z.y * a.w.z + a.w.x * a.x.y * a.z.z +
2037 a.z.x * a.w.y * a.x.z - a.x.x * a.w.y * a.z.z -
2038 a.z.x * a.x.y * a.w.z - a.w.x * a.z.y * a.x.z,
2039 a.x.x * a.w.y * a.y.z + a.y.x * a.x.y * a.w.z +
2040 a.w.x * a.y.y * a.x.z - a.x.x * a.y.y * a.w.z -
2041 a.w.x * a.x.y * a.y.z - a.y.x * a.w.y * a.x.z,
2042 a.x.x * a.y.y * a.z.z + a.z.x * a.x.y * a.y.z +
2043 a.y.x * a.z.y * a.x.z - a.x.x * a.z.y * a.y.z -
2044 a.y.x * a.x.y * a.z.z - a.z.x * a.y.y * a.x.z}};
2048constexpr T linalg::determinant(
const mat<T, 4, 4> &a) {
2049 return a.x.x * (a.y.y * a.z.z * a.w.w + a.w.y * a.y.z * a.z.w +
2050 a.z.y * a.w.z * a.y.w - a.y.y * a.w.z * a.z.w -
2051 a.z.y * a.y.z * a.w.w - a.w.y * a.z.z * a.y.w) +
2052 a.x.y * (a.y.z * a.w.w * a.z.x + a.z.z * a.y.w * a.w.x +
2053 a.w.z * a.z.w * a.y.x - a.y.z * a.z.w * a.w.x -
2054 a.w.z * a.y.w * a.z.x - a.z.z * a.w.w * a.y.x) +
2055 a.x.z * (a.y.w * a.z.x * a.w.y + a.w.w * a.y.x * a.z.y +
2056 a.z.w * a.w.x * a.y.y - a.y.w * a.w.x * a.z.y -
2057 a.z.w * a.y.x * a.w.y - a.w.w * a.z.x * a.y.y) +
2058 a.x.w * (a.y.x * a.w.y * a.z.z + a.z.x * a.y.y * a.w.z +
2059 a.w.x * a.z.y * a.y.z - a.y.x * a.z.y * a.w.z -
2060 a.w.x * a.y.y * a.z.z - a.z.x * a.w.y * a.y.z);
2065 const vec<T, 4> q{m.x.x - m.y.y - m.z.z, m.y.y - m.x.x - m.z.z,
2066 m.z.z - m.x.x - m.y.y, m.x.x + m.y.y + m.z.z},
2067 s[]{{1, m.x.y + m.y.x, m.z.x + m.x.z, m.y.z - m.z.y},
2068 {m.x.y + m.y.x, 1, m.y.z + m.z.y, m.z.x - m.x.z},
2069 {m.x.z + m.z.x, m.y.z + m.z.y, 1, m.x.y - m.y.x},
2070 {m.y.z - m.z.y, m.z.x - m.x.z, m.x.y - m.y.x, 1}};
2071 return copysign(normalize(sqrt(max(T(0), T(1) + q))), s[argmax(q)]);
2076 const vec<T, 3> ¢er,
2077 const vec<T, 3> &view_y_dir,
2079 const vec<T, 3> f = normalize(center - eye), z = a == pos_z ? f : -f,
2080 x = normalize(cross(view_y_dir, z)), y = cross(z, x);
2081 return inverse(mat<T, 4, 4>{{x, 0}, {y, 0}, {z, 0}, {eye, 1}});
2086 fwd_axis a, z_range z) {
2087 const T s = a == pos_z ? T(1) : T(-1), o = z == neg_one_to_one ? n : 0;
2088 return {{2 * n / (x1 - x0), 0, 0, 0},
2089 {0, 2 * n / (y1 - y0), 0, 0},
2090 {-s * (x0 + x1) / (x1 - x0), -s * (y0 + y1) / (y1 - y0),
2091 s * (f + o) / (f - n), s},
2092 {0, 0, -(n + o) * f / (f - n), 0}};