Manifold 3.0
Robust geometry
 
Loading...
Searching...
No Matches
linalg.h
1// Copyright 2024 The Manifold Authors.
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7// http://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15// Based on linalg.h - 2.2 - Single-header public domain linear algebra library
16//
17// The intent of this library is to provide the bulk of the functionality
18// you need to write programs that frequently use small, fixed-size vectors
19// and matrices, in domains such as computational geometry or computer
20// graphics. It strives for terse, readable source code.
21//
22// The original author of this software is Sterling Orsten, and its permanent
23// home is <http://github.com/sgorsten/linalg/>. If you find this software
24// useful, an acknowledgement in your source text and/or product documentation
25// is appreciated, but not required.
26//
27// The author acknowledges significant insights and contributions by:
28// Stan Melax <http://github.com/melax/>
29// Dimitri Diakopoulos <http://github.com/ddiakopoulos/>
30//
31// Some features are deprecated. Define LINALG_FORWARD_COMPATIBLE to remove
32// them.
33
34#pragma once
35#ifndef LINALG_H
36#define LINALG_H
37
38#include <array> // For std::array
39#include <cmath> // For various unary math functions, such as std::sqrt
40#include <cstdlib> // To resolve std::abs ambiguity on clang
41#include <functional> // For std::hash declaration
42#include <type_traits> // For std::enable_if, std::is_same, std::declval
43
44// In Visual Studio 2015, `constexpr` applied to a member function implies
45// `const`, which causes ambiguous overload resolution
46#if defined(_MSC_VER) && (_MSC_VER <= 1900)
47#define LINALG_CONSTEXPR14
48#else
49#define LINALG_CONSTEXPR14 constexpr
50#endif
51
52namespace linalg {
53// Small, fixed-length vector type, consisting of exactly M elements of type T,
54// and presumed to be a column-vector unless otherwise noted.
55template <class T, int M>
56struct vec;
57
58// Small, fixed-size matrix type, consisting of exactly M rows and N columns of
59// type T, stored in column-major order.
60template <class T, int M, int N>
61struct mat;
62
63// Specialize converter<T,U> with a function application operator that converts
64// type U to type T to enable implicit conversions
65template <class T, class U>
66struct converter {};
67namespace detail {
68template <class T, class U>
69using conv_t = typename std::enable_if<!std::is_same<T, U>::value,
70 decltype(converter<T, U>{}(
71 std::declval<U>()))>::type;
72
73// Trait for retrieving scalar type of any linear algebra object
74template <class A>
75struct scalar_type {};
76template <class T, int M>
77struct scalar_type<vec<T, M>> {
78 using type = T;
79};
80template <class T, int M, int N>
81struct scalar_type<mat<T, M, N>> {
82 using type = T;
83};
84
85// Type returned by the compare(...) function which supports all six comparison
86// operators against 0
87template <class T>
88struct ord {
89 T a, b;
90};
91template <class T>
92constexpr bool operator==(const ord<T>& o, std::nullptr_t) {
93 return o.a == o.b;
94}
95template <class T>
96constexpr bool operator!=(const ord<T>& o, std::nullptr_t) {
97 return !(o.a == o.b);
98}
99template <class T>
100constexpr bool operator<(const ord<T>& o, std::nullptr_t) {
101 return o.a < o.b;
102}
103template <class T>
104constexpr bool operator>(const ord<T>& o, std::nullptr_t) {
105 return o.b < o.a;
106}
107template <class T>
108constexpr bool operator<=(const ord<T>& o, std::nullptr_t) {
109 return !(o.b < o.a);
110}
111template <class T>
112constexpr bool operator>=(const ord<T>& o, std::nullptr_t) {
113 return !(o.a < o.b);
114}
115
116// Patterns which can be used with the compare(...) function
117template <class A, class B>
118struct any_compare {};
119template <class T>
120struct any_compare<vec<T, 1>, vec<T, 1>> {
121 using type = ord<T>;
122 constexpr ord<T> operator()(const vec<T, 1>& a, const vec<T, 1>& b) const {
123 return ord<T>{a.x, b.x};
124 }
125};
126template <class T>
127struct any_compare<vec<T, 2>, vec<T, 2>> {
128 using type = ord<T>;
129 constexpr ord<T> operator()(const vec<T, 2>& a, const vec<T, 2>& b) const {
130 return !(a.x == b.x) ? ord<T>{a.x, b.x} : ord<T>{a.y, b.y};
131 }
132};
133template <class T>
134struct any_compare<vec<T, 3>, vec<T, 3>> {
135 using type = ord<T>;
136 constexpr ord<T> operator()(const vec<T, 3>& a, const vec<T, 3>& b) const {
137 return !(a.x == b.x) ? ord<T>{a.x, b.x}
138 : !(a.y == b.y) ? ord<T>{a.y, b.y}
139 : ord<T>{a.z, b.z};
140 }
141};
142template <class T>
143struct any_compare<vec<T, 4>, vec<T, 4>> {
144 using type = ord<T>;
145 constexpr ord<T> operator()(const vec<T, 4>& a, const vec<T, 4>& b) const {
146 return !(a.x == b.x) ? ord<T>{a.x, b.x}
147 : !(a.y == b.y) ? ord<T>{a.y, b.y}
148 : !(a.z == b.z) ? ord<T>{a.z, b.z}
149 : ord<T>{a.w, b.w};
150 }
151};
152template <class T, int M>
153struct any_compare<mat<T, M, 1>, mat<T, M, 1>> {
154 using type = ord<T>;
155 constexpr ord<T> operator()(const mat<T, M, 1>& a,
156 const mat<T, M, 1>& b) const {
157 return compare(a.x, b.x);
158 }
159};
160template <class T, int M>
161struct any_compare<mat<T, M, 2>, mat<T, M, 2>> {
162 using type = ord<T>;
163 constexpr ord<T> operator()(const mat<T, M, 2>& a,
164 const mat<T, M, 2>& b) const {
165 return a.x != b.x ? compare(a.x, b.x) : compare(a.y, b.y);
166 }
167};
168template <class T, int M>
169struct any_compare<mat<T, M, 3>, mat<T, M, 3>> {
170 using type = ord<T>;
171 constexpr ord<T> operator()(const mat<T, M, 3>& a,
172 const mat<T, M, 3>& b) const {
173 return a.x != b.x ? compare(a.x, b.x)
174 : a.y != b.y ? compare(a.y, b.y)
175 : compare(a.z, b.z);
176 }
177};
178template <class T, int M>
179struct any_compare<mat<T, M, 4>, mat<T, M, 4>> {
180 using type = ord<T>;
181 constexpr ord<T> operator()(const mat<T, M, 4>& a,
182 const mat<T, M, 4>& b) const {
183 return a.x != b.x ? compare(a.x, b.x)
184 : a.y != b.y ? compare(a.y, b.y)
185 : a.z != b.z ? compare(a.z, b.z)
186 : compare(a.w, b.w);
187 }
188};
189
190// Helper for compile-time index-based access to members of vector and matrix
191// types
192template <int I>
193struct getter;
194template <>
195struct getter<0> {
196 template <class A>
197 constexpr auto operator()(A& a) const -> decltype(a.x) {
198 return a.x;
199 }
200};
201template <>
202struct getter<1> {
203 template <class A>
204 constexpr auto operator()(A& a) const -> decltype(a.y) {
205 return a.y;
206 }
207};
208template <>
209struct getter<2> {
210 template <class A>
211 constexpr auto operator()(A& a) const -> decltype(a.z) {
212 return a.z;
213 }
214};
215template <>
216struct getter<3> {
217 template <class A>
218 constexpr auto operator()(A& a) const -> decltype(a.w) {
219 return a.w;
220 }
221};
222
223// Stand-in for std::integer_sequence/std::make_integer_sequence
224template <int... I>
225struct seq {};
226template <int A, int N>
228template <int A>
229struct make_seq_impl<A, 0> {
230 using type = seq<>;
231};
232template <int A>
233struct make_seq_impl<A, 1> {
234 using type = seq<A + 0>;
235};
236template <int A>
237struct make_seq_impl<A, 2> {
238 using type = seq<A + 0, A + 1>;
239};
240template <int A>
241struct make_seq_impl<A, 3> {
242 using type = seq<A + 0, A + 1, A + 2>;
243};
244template <int A>
245struct make_seq_impl<A, 4> {
247};
248template <int A, int B>
249using make_seq = typename make_seq_impl<A, B - A>::type;
250template <class T, int M, int... I>
251vec<T, sizeof...(I)> constexpr swizzle(const vec<T, M>& v, seq<I...>) {
252 return {getter<I>{}(v)...};
253}
254template <class T, int M, int N, int... I, int... J>
255mat<T, sizeof...(I), sizeof...(J)> constexpr swizzle(const mat<T, M, N>& m,
256 seq<I...> i, seq<J...>) {
257 return {swizzle(getter<J>{}(m), i)...};
258}
259
260// SFINAE helpers to determine result of function application
261template <class F, class... T>
262using ret_t = decltype(std::declval<F>()(std::declval<T>()...));
263
264// SFINAE helper which is defined if all provided types are scalars
265struct empty {};
266template <class... T>
267struct scalars;
268template <>
269struct scalars<> {
270 using type = void;
271};
272template <class T, class... U>
273struct scalars<T, U...> : std::conditional<std::is_arithmetic<T>::value,
274 scalars<U...>, empty>::type {};
275template <class... T>
276using scalars_t = typename scalars<T...>::type;
277
278// Helpers which indicate how apply(F, ...) should be called for various
279// arguments
280template <class F, class Void, class... T>
281struct apply {}; // Patterns which contain only vectors or scalars
282template <class F, int M, class A>
283struct apply<F, scalars_t<>, vec<A, M>> {
284 using type = vec<ret_t<F, A>, M>;
285 enum { size = M, mm = 0 };
286 template <int... I>
287 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a) {
288 return {f(getter<I>{}(a))...};
289 }
290};
291template <class F, int M, class A, class B>
292struct apply<F, scalars_t<>, vec<A, M>, vec<B, M>> {
293 using type = vec<ret_t<F, A, B>, M>;
294 enum { size = M, mm = 0 };
295 template <int... I>
296 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a,
297 const vec<B, M>& b) {
298 return {f(getter<I>{}(a), getter<I>{}(b))...};
299 }
300};
301template <class F, int M, class A, class B>
302struct apply<F, scalars_t<B>, vec<A, M>, B> {
303 using type = vec<ret_t<F, A, B>, M>;
304 enum { size = M, mm = 0 };
305 template <int... I>
306 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a, B b) {
307 return {f(getter<I>{}(a), b)...};
308 }
309};
310template <class F, int M, class A, class B>
311struct apply<F, scalars_t<A>, A, vec<B, M>> {
312 using type = vec<ret_t<F, A, B>, M>;
313 enum { size = M, mm = 0 };
314 template <int... I>
315 static constexpr type impl(seq<I...>, F f, A a, const vec<B, M>& b) {
316 return {f(a, getter<I>{}(b))...};
317 }
318};
319template <class F, int M, class A, class B, class C>
320struct apply<F, scalars_t<>, vec<A, M>, vec<B, M>, vec<C, M>> {
321 using type = vec<ret_t<F, A, B, C>, M>;
322 enum { size = M, mm = 0 };
323 template <int... I>
324 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a,
325 const vec<B, M>& b, const vec<C, M>& c) {
326 return {f(getter<I>{}(a), getter<I>{}(b), getter<I>{}(c))...};
327 }
328};
329template <class F, int M, class A, class B, class C>
330struct apply<F, scalars_t<C>, vec<A, M>, vec<B, M>, C> {
331 using type = vec<ret_t<F, A, B, C>, M>;
332 enum { size = M, mm = 0 };
333 template <int... I>
334 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a,
335 const vec<B, M>& b, C c) {
336 return {f(getter<I>{}(a), getter<I>{}(b), c)...};
337 }
338};
339template <class F, int M, class A, class B, class C>
340struct apply<F, scalars_t<B>, vec<A, M>, B, vec<C, M>> {
341 using type = vec<ret_t<F, A, B, C>, M>;
342 enum { size = M, mm = 0 };
343 template <int... I>
344 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a, B b,
345 const vec<C, M>& c) {
346 return {f(getter<I>{}(a), b, getter<I>{}(c))...};
347 }
348};
349template <class F, int M, class A, class B, class C>
350struct apply<F, scalars_t<B, C>, vec<A, M>, B, C> {
351 using type = vec<ret_t<F, A, B, C>, M>;
352 enum { size = M, mm = 0 };
353 template <int... I>
354 static constexpr type impl(seq<I...>, F f, const vec<A, M>& a, B b, C c) {
355 return {f(getter<I>{}(a), b, c)...};
356 }
357};
358template <class F, int M, class A, class B, class C>
359struct apply<F, scalars_t<A>, A, vec<B, M>, vec<C, M>> {
360 using type = vec<ret_t<F, A, B, C>, M>;
361 enum { size = M, mm = 0 };
362 template <int... I>
363 static constexpr type impl(seq<I...>, F f, A a, const vec<B, M>& b,
364 const vec<C, M>& c) {
365 return {f(a, getter<I>{}(b), getter<I>{}(c))...};
366 }
367};
368template <class F, int M, class A, class B, class C>
369struct apply<F, scalars_t<A, C>, A, vec<B, M>, C> {
370 using type = vec<ret_t<F, A, B, C>, M>;
371 enum { size = M, mm = 0 };
372 template <int... I>
373 static constexpr type impl(seq<I...>, F f, A a, const vec<B, M>& b, C c) {
374 return {f(a, getter<I>{}(b), c)...};
375 }
376};
377template <class F, int M, class A, class B, class C>
378struct apply<F, scalars_t<A, B>, A, B, vec<C, M>> {
379 using type = vec<ret_t<F, A, B, C>, M>;
380 enum { size = M, mm = 0 };
381 template <int... I>
382 static constexpr type impl(seq<I...>, F f, A a, B b, const vec<C, M>& c) {
383 return {f(a, b, getter<I>{}(c))...};
384 }
385};
386template <class F, int M, int N, class A>
387struct apply<F, scalars_t<>, mat<A, M, N>> {
388 using type = mat<ret_t<F, A>, M, N>;
389 enum { size = N, mm = 0 };
390 template <int... J>
391 static constexpr type impl(seq<J...>, F f, const mat<A, M, N>& a) {
392 return {apply<F, void, vec<A, M>>::impl(make_seq<0, M>{}, f,
393 getter<J>{}(a))...};
394 }
395};
396template <class F, int M, int N, class A, class B>
397struct apply<F, scalars_t<>, mat<A, M, N>, mat<B, M, N>> {
398 using type = mat<ret_t<F, A, B>, M, N>;
399 enum { size = N, mm = 1 };
400 template <int... J>
401 static constexpr type impl(seq<J...>, F f, const mat<A, M, N>& a,
402 const mat<B, M, N>& b) {
403 return {apply<F, void, vec<A, M>, vec<B, M>>::impl(
404 make_seq<0, M>{}, f, getter<J>{}(a), getter<J>{}(b))...};
405 }
406};
407template <class F, int M, int N, class A, class B>
408struct apply<F, scalars_t<B>, mat<A, M, N>, B> {
409 using type = mat<ret_t<F, A, B>, M, N>;
410 enum { size = N, mm = 0 };
411 template <int... J>
412 static constexpr type impl(seq<J...>, F f, const mat<A, M, N>& a, B b) {
413 return {apply<F, void, vec<A, M>, B>::impl(make_seq<0, M>{}, f,
414 getter<J>{}(a), b)...};
415 }
416};
417template <class F, int M, int N, class A, class B>
418struct apply<F, scalars_t<A>, A, mat<B, M, N>> {
419 using type = mat<ret_t<F, A, B>, M, N>;
420 enum { size = N, mm = 0 };
421 template <int... J>
422 static constexpr type impl(seq<J...>, F f, A a, const mat<B, M, N>& b) {
423 return {apply<F, void, A, vec<B, M>>::impl(make_seq<0, M>{}, f, a,
424 getter<J>{}(b))...};
425 }
426};
427template <class F, class... A>
428struct apply<F, scalars_t<A...>, A...> {
429 using type = ret_t<F, A...>;
430 enum { size = 0, mm = 0 };
431 static constexpr type impl(seq<>, F f, A... a) { return f(a...); }
432};
433
434// Function objects for selecting between alternatives
435struct min {
436 template <class A, class B>
437 constexpr auto operator()(A a, B b) const ->
438 typename std::remove_reference<decltype(a < b ? a : b)>::type {
439 return a < b ? a : b;
440 }
441};
442struct max {
443 template <class A, class B>
444 constexpr auto operator()(A a, B b) const ->
445 typename std::remove_reference<decltype(a < b ? b : a)>::type {
446 return a < b ? b : a;
447 }
448};
449struct clamp {
450 template <class A, class B, class C>
451 constexpr auto operator()(A a, B b, C c) const ->
452 typename std::remove_reference<decltype(a < b ? b
453 : a < c ? a
454 : c)>::type {
455 return a < b ? b : a < c ? a : c;
456 }
457};
458struct select {
459 template <class A, class B, class C>
460 constexpr auto operator()(A a, B b, C c) const ->
461 typename std::remove_reference<decltype(a ? b : c)>::type {
462 return a ? b : c;
463 }
464};
465struct lerp {
466 template <class A, class B, class C>
467 constexpr auto operator()(A a, B b, C c) const
468 -> decltype(a * (1 - c) + b * c) {
469 return a * (1 - c) + b * c;
470 }
471};
472
473// Function objects for applying operators
474struct op_pos {
475 template <class A>
476 constexpr auto operator()(A a) const -> decltype(+a) {
477 return +a;
478 }
479};
480struct op_neg {
481 template <class A>
482 constexpr auto operator()(A a) const -> decltype(-a) {
483 return -a;
484 }
485};
486struct op_not {
487 template <class A>
488 constexpr auto operator()(A a) const -> decltype(!a) {
489 return !a;
490 }
491};
492struct op_cmp {
493 template <class A>
494 constexpr auto operator()(A a) const -> decltype(~(a)) {
495 return ~a;
496 }
497};
498struct op_mul {
499 template <class A, class B>
500 constexpr auto operator()(A a, B b) const -> decltype(a * b) {
501 return a * b;
502 }
503};
504struct op_div {
505 template <class A, class B>
506 constexpr auto operator()(A a, B b) const -> decltype(a / b) {
507 return a / b;
508 }
509};
510struct op_mod {
511 template <class A, class B>
512 constexpr auto operator()(A a, B b) const -> decltype(a % b) {
513 return a % b;
514 }
515};
516struct op_add {
517 template <class A, class B>
518 constexpr auto operator()(A a, B b) const -> decltype(a + b) {
519 return a + b;
520 }
521};
522struct op_sub {
523 template <class A, class B>
524 constexpr auto operator()(A a, B b) const -> decltype(a - b) {
525 return a - b;
526 }
527};
528struct op_lsh {
529 template <class A, class B>
530 constexpr auto operator()(A a, B b) const -> decltype(a << b) {
531 return a << b;
532 }
533};
534struct op_rsh {
535 template <class A, class B>
536 constexpr auto operator()(A a, B b) const -> decltype(a >> b) {
537 return a >> b;
538 }
539};
540struct op_lt {
541 template <class A, class B>
542 constexpr auto operator()(A a, B b) const -> decltype(a < b) {
543 return a < b;
544 }
545};
546struct op_gt {
547 template <class A, class B>
548 constexpr auto operator()(A a, B b) const -> decltype(a > b) {
549 return a > b;
550 }
551};
552struct op_le {
553 template <class A, class B>
554 constexpr auto operator()(A a, B b) const -> decltype(a <= b) {
555 return a <= b;
556 }
557};
558struct op_ge {
559 template <class A, class B>
560 constexpr auto operator()(A a, B b) const -> decltype(a >= b) {
561 return a >= b;
562 }
563};
564struct op_eq {
565 template <class A, class B>
566 constexpr auto operator()(A a, B b) const -> decltype(a == b) {
567 return a == b;
568 }
569};
570struct op_ne {
571 template <class A, class B>
572 constexpr auto operator()(A a, B b) const -> decltype(a != b) {
573 return a != b;
574 }
575};
576struct op_int {
577 template <class A, class B>
578 constexpr auto operator()(A a, B b) const -> decltype(a & b) {
579 return a & b;
580 }
581};
582struct op_xor {
583 template <class A, class B>
584 constexpr auto operator()(A a, B b) const -> decltype(a ^ b) {
585 return a ^ b;
586 }
587};
588struct op_un {
589 template <class A, class B>
590 constexpr auto operator()(A a, B b) const -> decltype(a | b) {
591 return a | b;
592 }
593};
594struct op_and {
595 template <class A, class B>
596 constexpr auto operator()(A a, B b) const -> decltype(a && b) {
597 return a && b;
598 }
599};
600struct op_or {
601 template <class A, class B>
602 constexpr auto operator()(A a, B b) const -> decltype(a || b) {
603 return a || b;
604 }
605};
606
607// Function objects for applying standard library math functions
609 template <class A>
610 constexpr auto operator()(A a) const -> decltype(std::isfinite(a)) {
611 return std::isfinite(a);
612 }
613};
614struct std_abs {
615 template <class A>
616 constexpr auto operator()(A a) const -> decltype(std::abs(a)) {
617 return std::abs(a);
618 }
619};
620struct std_floor {
621 template <class A>
622 constexpr auto operator()(A a) const -> decltype(std::floor(a)) {
623 return std::floor(a);
624 }
625};
626struct std_ceil {
627 template <class A>
628 constexpr auto operator()(A a) const -> decltype(std::ceil(a)) {
629 return std::ceil(a);
630 }
631};
632struct std_exp {
633 template <class A>
634 constexpr auto operator()(A a) const -> decltype(std::exp(a)) {
635 return std::exp(a);
636 }
637};
638struct std_log {
639 template <class A>
640 constexpr auto operator()(A a) const -> decltype(std::log(a)) {
641 return std::log(a);
642 }
643};
644struct std_log2 {
645 template <class A>
646 constexpr auto operator()(A a) const -> decltype(std::log2(a)) {
647 return std::log2(a);
648 }
649};
650struct std_log10 {
651 template <class A>
652 constexpr auto operator()(A a) const -> decltype(std::log10(a)) {
653 return std::log10(a);
654 }
655};
656struct std_sqrt {
657 template <class A>
658 constexpr auto operator()(A a) const -> decltype(std::sqrt(a)) {
659 return std::sqrt(a);
660 }
661};
662struct std_sin {
663 template <class A>
664 constexpr auto operator()(A a) const -> decltype(std::sin(a)) {
665 return std::sin(a);
666 }
667};
668struct std_cos {
669 template <class A>
670 constexpr auto operator()(A a) const -> decltype(std::cos(a)) {
671 return std::cos(a);
672 }
673};
674struct std_tan {
675 template <class A>
676 constexpr auto operator()(A a) const -> decltype(std::tan(a)) {
677 return std::tan(a);
678 }
679};
680struct std_asin {
681 template <class A>
682 constexpr auto operator()(A a) const -> decltype(std::asin(a)) {
683 return std::asin(a);
684 }
685};
686struct std_acos {
687 template <class A>
688 constexpr auto operator()(A a) const -> decltype(std::acos(a)) {
689 return std::acos(a);
690 }
691};
692struct std_atan {
693 template <class A>
694 constexpr auto operator()(A a) const -> decltype(std::atan(a)) {
695 return std::atan(a);
696 }
697};
698struct std_sinh {
699 template <class A>
700 constexpr auto operator()(A a) const -> decltype(std::sinh(a)) {
701 return std::sinh(a);
702 }
703};
704struct std_cosh {
705 template <class A>
706 constexpr auto operator()(A a) const -> decltype(std::cosh(a)) {
707 return std::cosh(a);
708 }
709};
710struct std_tanh {
711 template <class A>
712 constexpr auto operator()(A a) const -> decltype(std::tanh(a)) {
713 return std::tanh(a);
714 }
715};
716struct std_round {
717 template <class A>
718 constexpr auto operator()(A a) const -> decltype(std::round(a)) {
719 return std::round(a);
720 }
721};
722struct std_fmod {
723 template <class A, class B>
724 constexpr auto operator()(A a, B b) const -> decltype(std::fmod(a, b)) {
725 return std::fmod(a, b);
726 }
727};
728struct std_pow {
729 template <class A, class B>
730 constexpr auto operator()(A a, B b) const -> decltype(std::pow(a, b)) {
731 return std::pow(a, b);
732 }
733};
734struct std_atan2 {
735 template <class A, class B>
736 constexpr auto operator()(A a, B b) const -> decltype(std::atan2(a, b)) {
737 return std::atan2(a, b);
738 }
739};
741 template <class A, class B>
742 constexpr auto operator()(A a, B b) const -> decltype(std::copysign(a, b)) {
743 return std::copysign(a, b);
744 }
745};
746} // namespace detail
747
751
847template <class T>
848struct vec<T, 1> {
849 T x;
850 constexpr vec() : x() {}
851 constexpr vec(const T& x_) : x(x_) {}
852 // NOTE: vec<T,1> does NOT have a constructor from pointer, this can conflict
853 // with initializing its single element from zero
854 template <class U>
855 constexpr explicit vec(const vec<U, 1>& v) : vec(static_cast<T>(v.x)) {}
856 constexpr const T& operator[](int) const { return x; }
857 LINALG_CONSTEXPR14 T& operator[](int) { return x; }
858
859 template <class U, class = detail::conv_t<vec, U>>
860 constexpr vec(const U& u) : vec(converter<vec, U>{}(u)) {}
861 template <class U, class = detail::conv_t<U, vec>>
862 constexpr operator U() const {
863 return converter<U, vec>{}(*this);
864 }
865};
866template <class T>
867struct vec<T, 2> {
868 T x, y;
869 constexpr vec() : x(), y() {}
870 constexpr vec(const T& x_, const T& y_) : x(x_), y(y_) {}
871 constexpr explicit vec(const T& s) : vec(s, s) {}
872 constexpr explicit vec(const T* p) : vec(p[0], p[1]) {}
873 template <class U, int N>
874 constexpr explicit vec(const vec<U, N>& v)
875 : vec(static_cast<T>(v.x), static_cast<T>(v.y)) {
876 static_assert(
877 N >= 2,
878 "You must give extra arguments if your input vector is shorter.");
879 }
880 constexpr const T& operator[](int i) const { return i == 0 ? x : y; }
881 LINALG_CONSTEXPR14 T& operator[](int i) { return i == 0 ? x : y; }
882
883 template <class U, class = detail::conv_t<vec, U>>
884 constexpr vec(const U& u) : vec(converter<vec, U>{}(u)) {}
885 template <class U, class = detail::conv_t<U, vec>>
886 constexpr operator U() const {
887 return converter<U, vec>{}(*this);
888 }
889};
890template <class T>
891struct vec<T, 3> {
892 T x, y, z;
893 constexpr vec() : x(), y(), z() {}
894 constexpr vec(const T& x_, const T& y_, const T& z_) : x(x_), y(y_), z(z_) {}
895 constexpr vec(const vec<T, 2>& xy, const T& z_) : vec(xy.x, xy.y, z_) {}
896 constexpr explicit vec(const T& s) : vec(s, s, s) {}
897 constexpr explicit vec(const T* p) : vec(p[0], p[1], p[2]) {}
898 template <class U, int N>
899 constexpr explicit vec(const vec<U, N>& v)
900 : vec(static_cast<T>(v.x), static_cast<T>(v.y), static_cast<T>(v.z)) {
901 static_assert(
902 N >= 3,
903 "You must give extra arguments if your input vector is shorter.");
904 }
905 constexpr const T& operator[](int i) const {
906 return i == 0 ? x : i == 1 ? y : z;
907 }
908 LINALG_CONSTEXPR14 T& operator[](int i) {
909 return i == 0 ? x : i == 1 ? y : z;
910 }
911 constexpr const vec<T, 2>& xy() const {
912 return *reinterpret_cast<const vec<T, 2>*>(this);
913 }
914 vec<T, 2>& xy() { return *reinterpret_cast<vec<T, 2>*>(this); }
915
916 template <class U, class = detail::conv_t<vec, U>>
917 constexpr vec(const U& u) : vec(converter<vec, U>{}(u)) {}
918 template <class U, class = detail::conv_t<U, vec>>
919 constexpr operator U() const {
920 return converter<U, vec>{}(*this);
921 }
922};
923template <class T>
924struct vec<T, 4> {
925 T x, y, z, w;
926 constexpr vec() : x(), y(), z(), w() {}
927 constexpr vec(const T& x_, const T& y_, const T& z_, const T& w_)
928 : x(x_), y(y_), z(z_), w(w_) {}
929 constexpr vec(const vec<T, 2>& xy, const T& z_, const T& w_)
930 : vec(xy.x, xy.y, z_, w_) {}
931 constexpr vec(const vec<T, 3>& xyz, const T& w_)
932 : vec(xyz.x, xyz.y, xyz.z, w_) {}
933 constexpr explicit vec(const T& s) : vec(s, s, s, s) {}
934 constexpr explicit vec(const T* p) : vec(p[0], p[1], p[2], p[3]) {}
935 template <class U, int N>
936 constexpr explicit vec(const vec<U, N>& v)
937 : vec(static_cast<T>(v.x), static_cast<T>(v.y), static_cast<T>(v.z),
938 static_cast<T>(v.w)) {
939 static_assert(
940 N >= 4,
941 "You must give extra arguments if your input vector is shorter.");
942 }
943 constexpr const T& operator[](int i) const {
944 return i == 0 ? x : i == 1 ? y : i == 2 ? z : w;
945 }
946 LINALG_CONSTEXPR14 T& operator[](int i) {
947 return i == 0 ? x : i == 1 ? y : i == 2 ? z : w;
948 }
949 constexpr const vec<T, 2>& xy() const {
950 return *reinterpret_cast<const vec<T, 2>*>(this);
951 }
952 constexpr const vec<T, 3>& xyz() const {
953 return *reinterpret_cast<const vec<T, 3>*>(this);
954 }
955 vec<T, 2>& xy() { return *reinterpret_cast<vec<T, 2>*>(this); }
956 vec<T, 3>& xyz() { return *reinterpret_cast<vec<T, 3>*>(this); }
957
958 template <class U, class = detail::conv_t<vec, U>>
959 constexpr vec(const U& u) : vec(converter<vec, U>{}(u)) {}
960 template <class U, class = detail::conv_t<U, vec>>
961 constexpr operator U() const {
962 return converter<U, vec>{}(*this);
963 }
964};
965
966
1061template <class T, int M>
1062struct mat<T, M, 1> {
1063 using V = vec<T, M>;
1064 V x;
1065 constexpr mat() : x() {}
1066 constexpr mat(const V& x_) : x(x_) {}
1067 constexpr explicit mat(const T& s) : x(s) {}
1068 constexpr explicit mat(const T* p) : x(p + M * 0) {}
1069 template <class U>
1070 constexpr explicit mat(const mat<U, M, 1>& m) : mat(V(m.x)) {}
1071 constexpr vec<T, 1> row(int i) const { return {x[i]}; }
1072 constexpr const V& operator[](int) const { return x; }
1073 LINALG_CONSTEXPR14 V& operator[](int) { return x; }
1074
1075 template <class U, class = detail::conv_t<mat, U>>
1076 constexpr mat(const U& u) : mat(converter<mat, U>{}(u)) {}
1077 template <class U, class = detail::conv_t<U, mat>>
1078 constexpr operator U() const {
1079 return converter<U, mat>{}(*this);
1080 }
1081};
1082template <class T, int M>
1083struct mat<T, M, 2> {
1084 using V = vec<T, M>;
1085 V x, y;
1086 constexpr mat() : x(), y() {}
1087 constexpr mat(const V& x_, const V& y_) : x(x_), y(y_) {}
1088 constexpr explicit mat(const T& s) : x(s), y(s) {}
1089 constexpr explicit mat(const T* p) : x(p + M * 0), y(p + M * 1) {}
1090 template <class U, int N, int P>
1091 constexpr explicit mat(const mat<U, N, P>& m) : mat(V(m.x), V(m.y)) {
1092 static_assert(P >= 2, "Input matrix dimensions must be at least as big.");
1093 }
1094 constexpr vec<T, 2> row(int i) const { return {x[i], y[i]}; }
1095 constexpr const V& operator[](int j) const { return j == 0 ? x : y; }
1096 LINALG_CONSTEXPR14 V& operator[](int j) { return j == 0 ? x : y; }
1097
1098 template <class U, class = detail::conv_t<mat, U>>
1099 constexpr mat(const U& u) : mat(converter<mat, U>{}(u)) {}
1100 template <class U, class = detail::conv_t<U, mat>>
1101 constexpr operator U() const {
1102 return converter<U, mat>{}(*this);
1103 }
1104};
1105template <class T, int M>
1106struct mat<T, M, 3> {
1107 using V = vec<T, M>;
1108 V x, y, z;
1109 constexpr mat() : x(), y(), z() {}
1110 constexpr mat(const V& x_, const V& y_, const V& z_) : x(x_), y(y_), z(z_) {}
1111 constexpr mat(const mat<T, M, 2>& m_, const V& z_)
1112 : x(m_.x), y(m_.y), z(z_) {}
1113 constexpr explicit mat(const T& s) : x(s), y(s), z(s) {}
1114 constexpr explicit mat(const T* p)
1115 : x(p + M * 0), y(p + M * 1), z(p + M * 2) {}
1116 template <class U, int N, int P>
1117 constexpr explicit mat(const mat<U, N, P>& m) : mat(V(m.x), V(m.y), V(m.z)) {
1118 static_assert(P >= 3, "Input matrix dimensions must be at least as big.");
1119 }
1120 constexpr vec<T, 3> row(int i) const { return {x[i], y[i], z[i]}; }
1121 constexpr const V& operator[](int j) const {
1122 return j == 0 ? x : j == 1 ? y : z;
1123 }
1124 LINALG_CONSTEXPR14 V& operator[](int j) {
1125 return j == 0 ? x : j == 1 ? y : z;
1126 }
1127
1128 template <class U, class = detail::conv_t<mat, U>>
1129 constexpr mat(const U& u) : mat(converter<mat, U>{}(u)) {}
1130 template <class U, class = detail::conv_t<U, mat>>
1131 constexpr operator U() const {
1132 return converter<U, mat>{}(*this);
1133 }
1134};
1135template <class T, int M>
1136struct mat<T, M, 4> {
1137 using V = vec<T, M>;
1138 V x, y, z, w;
1139 constexpr mat() : x(), y(), z(), w() {}
1140 constexpr mat(const V& x_, const V& y_, const V& z_, const V& w_)
1141 : x(x_), y(y_), z(z_), w(w_) {}
1142 constexpr mat(const mat<T, M, 3>& m_, const V& w_)
1143 : x(m_.x), y(m_.y), z(m_.z), w(w_) {}
1144 constexpr explicit mat(const T& s) : x(s), y(s), z(s), w(s) {}
1145 constexpr explicit mat(const T* p)
1146 : x(p + M * 0), y(p + M * 1), z(p + M * 2), w(p + M * 3) {}
1147 template <class U, int N, int P>
1148 constexpr explicit mat(const mat<U, N, P>& m)
1149 : mat(V(m.x), V(m.y), V(m.z), V(m.w)) {
1150 static_assert(P >= 4, "Input matrix dimensions must be at least as big.");
1151 }
1152
1153 constexpr vec<T, 4> row(int i) const { return {x[i], y[i], z[i], w[i]}; }
1154 constexpr const V& operator[](int j) const {
1155 return j == 0 ? x : j == 1 ? y : j == 2 ? z : w;
1156 }
1157 LINALG_CONSTEXPR14 V& operator[](int j) {
1158 return j == 0 ? x : j == 1 ? y : j == 2 ? z : w;
1159 }
1160
1161 template <class U, class = detail::conv_t<mat, U>>
1162 constexpr mat(const U& u) : mat(converter<mat, U>{}(u)) {}
1163 template <class U, class = detail::conv_t<U, mat>>
1164 constexpr operator U() const {
1165 return converter<U, mat>{}(*this);
1166 }
1167};
1168
1169
1176struct identity_t {
1177 constexpr explicit identity_t(int) {}
1178};
1179template <class T>
1180struct converter<mat<T, 1, 1>, identity_t> {
1181 constexpr mat<T, 1, 1> operator()(identity_t) const { return {vec<T, 1>{1}}; }
1182};
1183template <class T>
1184struct converter<mat<T, 2, 2>, identity_t> {
1185 constexpr mat<T, 2, 2> operator()(identity_t) const {
1186 return {{1, 0}, {0, 1}};
1187 }
1188};
1189template <class T>
1190struct converter<mat<T, 2, 3>, identity_t> {
1191 constexpr mat<T, 2, 3> operator()(identity_t) const {
1192 return {{1, 0}, {0, 1}, {0, 0}};
1193 }
1194};
1195template <class T>
1196struct converter<mat<T, 3, 3>, identity_t> {
1197 constexpr mat<T, 3, 3> operator()(identity_t) const {
1198 return {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
1199 }
1200};
1201template <class T>
1202struct converter<mat<T, 3, 4>, identity_t> {
1203 constexpr mat<T, 3, 4> operator()(identity_t) const {
1204 return {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}};
1205 }
1206};
1207template <class T>
1208struct converter<mat<T, 4, 4>, identity_t> {
1209 constexpr mat<T, 4, 4> operator()(identity_t) const {
1210 return {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
1211 }
1212};
1213constexpr identity_t identity{1};
1215
1223template <class F, class A, class B>
1224constexpr A fold(F f, A a, const vec<B, 1>& b) {
1225 return f(a, b.x);
1226}
1227template <class F, class A, class B>
1228constexpr A fold(F f, A a, const vec<B, 2>& b) {
1229 return f(f(a, b.x), b.y);
1230}
1231template <class F, class A, class B>
1232constexpr A fold(F f, A a, const vec<B, 3>& b) {
1233 return f(f(f(a, b.x), b.y), b.z);
1234}
1235template <class F, class A, class B>
1236constexpr A fold(F f, A a, const vec<B, 4>& b) {
1237 return f(f(f(f(a, b.x), b.y), b.z), b.w);
1238}
1239template <class F, class A, class B, int M>
1240constexpr A fold(F f, A a, const mat<B, M, 1>& b) {
1241 return fold(f, a, b.x);
1242}
1243template <class F, class A, class B, int M>
1244constexpr A fold(F f, A a, const mat<B, M, 2>& b) {
1245 return fold(f, fold(f, a, b.x), b.y);
1246}
1247template <class F, class A, class B, int M>
1248constexpr A fold(F f, A a, const mat<B, M, 3>& b) {
1249 return fold(f, fold(f, fold(f, a, b.x), b.y), b.z);
1250}
1251template <class F, class A, class B, int M>
1252constexpr A fold(F f, A a, const mat<B, M, 4>& b) {
1253 return fold(f, fold(f, fold(f, fold(f, a, b.x), b.y), b.z), b.w);
1254}
1256
1263
1264// Type aliases for the result of calling apply(...) with various arguments, can
1265// be used with return type SFINAE to constrain overload sets
1266template <class F, class... A>
1267using apply_t = typename detail::apply<F, void, A...>::type;
1268template <class F, class... A>
1269using mm_apply_t = typename std::enable_if<detail::apply<F, void, A...>::mm,
1270 apply_t<F, A...>>::type;
1271template <class F, class... A>
1272using no_mm_apply_t = typename std::enable_if<!detail::apply<F, void, A...>::mm,
1273 apply_t<F, A...>>::type;
1274template <class A>
1275using scalar_t =
1276 typename detail::scalar_type<A>::type; // Underlying scalar type when
1277 // performing elementwise operations
1278
1279// apply(f,...) applies the provided function in an elementwise fashion to its
1280// arguments, producing an object of the same dimensions
1281template <class F, class... A>
1282constexpr apply_t<F, A...> apply(F func, const A&... args) {
1284 detail::make_seq<0, detail::apply<F, void, A...>::size>{}, func, args...);
1285}
1286
1287// map(a,f) is equivalent to apply(f,a)
1288template <class A, class F>
1289constexpr apply_t<F, A> map(const A& a, F func) {
1290 return apply(func, a);
1291}
1292
1293// zip(a,b,f) is equivalent to apply(f,a,b)
1294template <class A, class B, class F>
1295constexpr apply_t<F, A, B> zip(const A& a, const B& b, F func) {
1296 return apply(func, a, b);
1297}
1299
1306template <class A, class B>
1307constexpr typename detail::any_compare<A, B>::type compare(const A& a,
1308 const B& b) {
1309 return detail::any_compare<A, B>()(a, b);
1310}
1311template <class A, class B>
1312constexpr auto operator==(const A& a, const B& b)
1313 -> decltype(compare(a, b) == 0) {
1314 return compare(a, b) == 0;
1315}
1316template <class A, class B>
1317constexpr auto operator!=(const A& a, const B& b)
1318 -> decltype(compare(a, b) != 0) {
1319 return compare(a, b) != 0;
1320}
1321template <class A, class B>
1322constexpr auto operator<(const A& a, const B& b)
1323 -> decltype(compare(a, b) < 0) {
1324 return compare(a, b) < 0;
1325}
1326template <class A, class B>
1327constexpr auto operator>(const A& a, const B& b)
1328 -> decltype(compare(a, b) > 0) {
1329 return compare(a, b) > 0;
1330}
1331template <class A, class B>
1332constexpr auto operator<=(const A& a, const B& b)
1333 -> decltype(compare(a, b) <= 0) {
1334 return compare(a, b) <= 0;
1335}
1336template <class A, class B>
1337constexpr auto operator>=(const A& a, const B& b)
1338 -> decltype(compare(a, b) >= 0) {
1339 return compare(a, b) >= 0;
1340}
1342
1348template <class A>
1349constexpr bool any(const A& a) {
1350 return fold(detail::op_or{}, false, a);
1351}
1352template <class A>
1353constexpr bool all(const A& a) {
1354 return fold(detail::op_and{}, true, a);
1355}
1356template <class A>
1357constexpr scalar_t<A> sum(const A& a) {
1358 return fold(detail::op_add{}, scalar_t<A>(0), a);
1359}
1360template <class A>
1361constexpr scalar_t<A> product(const A& a) {
1362 return fold(detail::op_mul{}, scalar_t<A>(1), a);
1363}
1364template <class A>
1365constexpr scalar_t<A> minelem(const A& a) {
1366 return fold(detail::min{}, a.x, a);
1367}
1368template <class A>
1369constexpr scalar_t<A> maxelem(const A& a) {
1370 return fold(detail::max{}, a.x, a);
1371}
1372template <class T, int M>
1373int argmin(const vec<T, M>& a) {
1374 int j = 0;
1375 for (int i = 1; i < M; ++i)
1376 if (a[i] < a[j]) j = i;
1377 return j;
1378}
1379template <class T, int M>
1380int argmax(const vec<T, M>& a) {
1381 int j = 0;
1382 for (int i = 1; i < M; ++i)
1383 if (a[i] > a[j]) j = i;
1384 return j;
1385}
1387
1393template <class A>
1394constexpr apply_t<detail::op_pos, A> operator+(const A& a) {
1395 return apply(detail::op_pos{}, a);
1396}
1397template <class A>
1398constexpr apply_t<detail::op_neg, A> operator-(const A& a) {
1399 return apply(detail::op_neg{}, a);
1400}
1401template <class A>
1402constexpr apply_t<detail::op_cmp, A> operator~(const A& a) {
1403 return apply(detail::op_cmp{}, a);
1404}
1405template <class A>
1406constexpr apply_t<detail::op_not, A> operator!(const A& a) {
1407 return apply(detail::op_not{}, a);
1408}
1410
1419template <class A, class B>
1420constexpr apply_t<detail::op_add, A, B> operator+(const A& a, const B& b) {
1421 return apply(detail::op_add{}, a, b);
1422}
1423template <class A, class B>
1424constexpr apply_t<detail::op_sub, A, B> operator-(const A& a, const B& b) {
1425 return apply(detail::op_sub{}, a, b);
1426}
1427template <class A, class B>
1428constexpr apply_t<detail::op_mul, A, B> cmul(const A& a, const B& b) {
1429 return apply(detail::op_mul{}, a, b);
1430}
1431template <class A, class B>
1432constexpr apply_t<detail::op_div, A, B> operator/(const A& a, const B& b) {
1433 return apply(detail::op_div{}, a, b);
1434}
1435template <class A, class B>
1436constexpr apply_t<detail::op_mod, A, B> operator%(const A& a, const B& b) {
1437 return apply(detail::op_mod{}, a, b);
1438}
1439template <class A, class B>
1440constexpr apply_t<detail::op_un, A, B> operator|(const A& a, const B& b) {
1441 return apply(detail::op_un{}, a, b);
1442}
1443template <class A, class B>
1444constexpr apply_t<detail::op_xor, A, B> operator^(const A& a, const B& b) {
1445 return apply(detail::op_xor{}, a, b);
1446}
1447template <class A, class B>
1448constexpr apply_t<detail::op_int, A, B> operator&(const A& a, const B& b) {
1449 return apply(detail::op_int{}, a, b);
1450}
1451template <class A, class B>
1452constexpr apply_t<detail::op_lsh, A, B> operator<<(const A& a, const B& b) {
1453 return apply(detail::op_lsh{}, a, b);
1454}
1455template <class A, class B>
1456constexpr apply_t<detail::op_rsh, A, B> operator>>(const A& a, const B& b) {
1457 return apply(detail::op_rsh{}, a, b);
1458}
1459
1460// Binary `operator *` represents the algebraic matrix product - use cmul(a, b)
1461// for the Hadamard (component-wise) product.
1462template <class A, class B>
1463constexpr auto operator*(const A& a, const B& b) {
1464 return mul(a, b);
1465}
1466
1467// Binary assignment operators a $= b is always defined as though it were
1468// explicitly written a = a $ b
1469template <class A, class B>
1470constexpr auto operator+=(A& a, const B& b) -> decltype(a = a + b) {
1471 return a = a + b;
1472}
1473template <class A, class B>
1474constexpr auto operator-=(A& a, const B& b) -> decltype(a = a - b) {
1475 return a = a - b;
1476}
1477template <class A, class B>
1478constexpr auto operator*=(A& a, const B& b) -> decltype(a = a * b) {
1479 return a = a * b;
1480}
1481template <class A, class B>
1482constexpr auto operator/=(A& a, const B& b) -> decltype(a = a / b) {
1483 return a = a / b;
1484}
1485template <class A, class B>
1486constexpr auto operator%=(A& a, const B& b) -> decltype(a = a % b) {
1487 return a = a % b;
1488}
1489template <class A, class B>
1490constexpr auto operator|=(A& a, const B& b) -> decltype(a = a | b) {
1491 return a = a | b;
1492}
1493template <class A, class B>
1494constexpr auto operator^=(A& a, const B& b) -> decltype(a = a ^ b) {
1495 return a = a ^ b;
1496}
1497template <class A, class B>
1498constexpr auto operator&=(A& a, const B& b) -> decltype(a = a & b) {
1499 return a = a & b;
1500}
1501template <class A, class B>
1502constexpr auto operator<<=(A& a, const B& b) -> decltype(a = a << b) {
1503 return a = a << b;
1504}
1505template <class A, class B>
1506constexpr auto operator>>=(A& a, const B& b) -> decltype(a = a >> b) {
1507 return a = a >> b;
1508}
1510
1520template <int... I, class T, int M>
1521constexpr vec<T, sizeof...(I)> swizzle(const vec<T, M>& a) {
1522 return {detail::getter<I>{}(a)...};
1523}
1524
1528template <int I0, int I1, class T, int M>
1529constexpr vec<T, I1 - I0> subvec(const vec<T, M>& a) {
1530 return detail::swizzle(a, detail::make_seq<I0, I1>{});
1531}
1532
1536template <int I0, int J0, int I1, int J1, class T, int M, int N>
1537constexpr mat<T, I1 - I0, J1 - J0> submat(const mat<T, M, N>& a) {
1538 return detail::swizzle(a, detail::make_seq<I0, I1>{},
1539 detail::make_seq<J0, J1>{});
1540}
1541
1542
1548template <class A>
1549constexpr apply_t<detail::std_isfinite, A> isfinite(const A& a) {
1550 return apply(detail::std_isfinite{}, a);
1551}
1552template <class A>
1553constexpr apply_t<detail::std_abs, A> abs(const A& a) {
1554 return apply(detail::std_abs{}, a);
1555}
1556template <class A>
1557constexpr apply_t<detail::std_floor, A> floor(const A& a) {
1558 return apply(detail::std_floor{}, a);
1559}
1560template <class A>
1561constexpr apply_t<detail::std_ceil, A> ceil(const A& a) {
1562 return apply(detail::std_ceil{}, a);
1563}
1564template <class A>
1565constexpr apply_t<detail::std_exp, A> exp(const A& a) {
1566 return apply(detail::std_exp{}, a);
1567}
1568template <class A>
1569constexpr apply_t<detail::std_log, A> log(const A& a) {
1570 return apply(detail::std_log{}, a);
1571}
1572template <class A>
1573constexpr apply_t<detail::std_log2, A> log2(const A& a) {
1574 return apply(detail::std_log2{}, a);
1575}
1576template <class A>
1577constexpr apply_t<detail::std_log10, A> log10(const A& a) {
1578 return apply(detail::std_log10{}, a);
1579}
1580template <class A>
1581constexpr apply_t<detail::std_sqrt, A> sqrt(const A& a) {
1582 return apply(detail::std_sqrt{}, a);
1583}
1584template <class A>
1585constexpr apply_t<detail::std_sin, A> sin(const A& a) {
1586 return apply(detail::std_sin{}, a);
1587}
1588template <class A>
1589constexpr apply_t<detail::std_cos, A> cos(const A& a) {
1590 return apply(detail::std_cos{}, a);
1591}
1592template <class A>
1593constexpr apply_t<detail::std_tan, A> tan(const A& a) {
1594 return apply(detail::std_tan{}, a);
1595}
1596template <class A>
1597constexpr apply_t<detail::std_asin, A> asin(const A& a) {
1598 return apply(detail::std_asin{}, a);
1599}
1600template <class A>
1601constexpr apply_t<detail::std_acos, A> acos(const A& a) {
1602 return apply(detail::std_acos{}, a);
1603}
1604template <class A>
1605constexpr apply_t<detail::std_atan, A> atan(const A& a) {
1606 return apply(detail::std_atan{}, a);
1607}
1608template <class A>
1609constexpr apply_t<detail::std_sinh, A> sinh(const A& a) {
1610 return apply(detail::std_sinh{}, a);
1611}
1612template <class A>
1613constexpr apply_t<detail::std_cosh, A> cosh(const A& a) {
1614 return apply(detail::std_cosh{}, a);
1615}
1616template <class A>
1617constexpr apply_t<detail::std_tanh, A> tanh(const A& a) {
1618 return apply(detail::std_tanh{}, a);
1619}
1620template <class A>
1621constexpr apply_t<detail::std_round, A> round(const A& a) {
1622 return apply(detail::std_round{}, a);
1623}
1625
1632template <class A, class B>
1633constexpr apply_t<detail::std_fmod, A, B> fmod(const A& a, const B& b) {
1634 return apply(detail::std_fmod{}, a, b);
1635}
1636template <class A, class B>
1637constexpr apply_t<detail::std_pow, A, B> pow(const A& a, const B& b) {
1638 return apply(detail::std_pow{}, a, b);
1639}
1640template <class A, class B>
1641constexpr apply_t<detail::std_atan2, A, B> atan2(const A& a, const B& b) {
1642 return apply(detail::std_atan2{}, a, b);
1643}
1644template <class A, class B>
1645constexpr apply_t<detail::std_copysign, A, B> copysign(const A& a, const B& b) {
1646 return apply(detail::std_copysign{}, a, b);
1647}
1649
1656template <class A, class B>
1657constexpr apply_t<detail::op_eq, A, B> equal(const A& a, const B& b) {
1658 return apply(detail::op_eq{}, a, b);
1659}
1660template <class A, class B>
1661constexpr apply_t<detail::op_ne, A, B> nequal(const A& a, const B& b) {
1662 return apply(detail::op_ne{}, a, b);
1663}
1664template <class A, class B>
1665constexpr apply_t<detail::op_lt, A, B> less(const A& a, const B& b) {
1666 return apply(detail::op_lt{}, a, b);
1667}
1668template <class A, class B>
1669constexpr apply_t<detail::op_gt, A, B> greater(const A& a, const B& b) {
1670 return apply(detail::op_gt{}, a, b);
1671}
1672template <class A, class B>
1673constexpr apply_t<detail::op_le, A, B> lequal(const A& a, const B& b) {
1674 return apply(detail::op_le{}, a, b);
1675}
1676template <class A, class B>
1677constexpr apply_t<detail::op_ge, A, B> gequal(const A& a, const B& b) {
1678 return apply(detail::op_ge{}, a, b);
1679}
1681
1688template <class A, class B>
1689constexpr apply_t<detail::min, A, B> min(const A& a, const B& b) {
1690 return apply(detail::min{}, a, b);
1691}
1692template <class A, class B>
1693constexpr apply_t<detail::max, A, B> max(const A& a, const B& b) {
1694 return apply(detail::max{}, a, b);
1695}
1699template <class X, class L, class H>
1700constexpr apply_t<detail::clamp, X, L, H> clamp(const X& x, const L& l,
1701 const H& h) {
1702 return apply(detail::clamp{}, x, l, h);
1703}
1704
1708template <class P, class A, class B>
1709constexpr apply_t<detail::select, P, A, B> select(const P& p, const A& a,
1710 const B& b) {
1711 return apply(detail::select{}, p, a, b);
1712}
1713
1717template <class A, class B, class T>
1718constexpr apply_t<detail::lerp, A, B, T> lerp(const A& a, const B& b,
1719 const T& t) {
1720 return apply(detail::lerp{}, a, b, t);
1721}
1722
1723
1732template <class T>
1733constexpr T cross(const vec<T, 2>& a, const vec<T, 2>& b) {
1734 return a.x * b.y - a.y * b.x;
1735}
1736
1739template <class T>
1740constexpr vec<T, 2> cross(T a, const vec<T, 2>& b) {
1741 return {-a * b.y, a * b.x};
1742}
1743
1746template <class T>
1747constexpr vec<T, 2> cross(const vec<T, 2>& a, T b) {
1748 return {a.y * b, -a.x * b};
1749}
1750
1754template <class T>
1755constexpr vec<T, 3> cross(const vec<T, 3>& a, const vec<T, 3>& b) {
1756 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};
1757}
1758
1762template <class T, int M>
1763constexpr T dot(const vec<T, M>& a, const vec<T, M>& b) {
1764 return sum(a * b);
1765}
1766
1769template <class T, int M>
1770constexpr T length2(const vec<T, M>& a) {
1771 return dot(a, a);
1772}
1773
1776template <class T, int M>
1777T length(const vec<T, M>& a) {
1778 return std::sqrt(length2(a));
1779}
1780
1785template <class T, int M>
1787 return a / length(a);
1788}
1789
1794template <class T, int M>
1795constexpr T distance2(const vec<T, M>& a, const vec<T, M>& b) {
1796 return length2(b - a);
1797}
1798
1802template <class T, int M>
1803T distance(const vec<T, M>& a, const vec<T, M>& b) {
1804 return length(b - a);
1805}
1806
1809template <class T, int M>
1810T uangle(const vec<T, M>& a, const vec<T, M>& b) {
1811 T d = dot(a, b);
1812 return d > 1 ? 0 : std::acos(d < -1 ? -1 : d);
1813}
1814
1817template <class T, int M>
1818T angle(const vec<T, M>& a, const vec<T, M>& b) {
1819 return uangle(normalize(a), normalize(b));
1820}
1821
1825template <class T>
1826vec<T, 2> rot(T a, const vec<T, 2>& v) {
1827 const T s = std::sin(a), c = std::cos(a);
1828 return {v.x * c - v.y * s, v.x * s + v.y * c};
1829}
1830
1834template <class T>
1835vec<T, 3> rotx(T a, const vec<T, 3>& v) {
1836 const T s = std::sin(a), c = std::cos(a);
1837 return {v.x, v.y * c - v.z * s, v.y * s + v.z * c};
1838}
1839
1843template <class T>
1844vec<T, 3> roty(T a, const vec<T, 3>& v) {
1845 const T s = std::sin(a), c = std::cos(a);
1846 return {v.x * c + v.z * s, v.y, -v.x * s + v.z * c};
1847}
1848
1852template <class T>
1853vec<T, 3> rotz(T a, const vec<T, 3>& v) {
1854 const T s = std::sin(a), c = std::cos(a);
1855 return {v.x * c - v.y * s, v.x * s + v.y * c, v.z};
1856}
1857
1860template <class T, int M>
1861vec<T, M> nlerp(const vec<T, M>& a, const vec<T, M>& b, T t) {
1862 return normalize(lerp(a, b, t));
1863}
1864
1869template <class T, int M>
1870vec<T, M> slerp(const vec<T, M>& a, const vec<T, M>& b, T t) {
1871 T th = uangle(a, b);
1872 return th == 0 ? a
1873 : a * (std::sin(th * (1 - t)) / std::sin(th)) +
1874 b * (std::sin(th * t) / std::sin(th));
1875}
1876
1877
1889template <class T>
1890constexpr vec<T, 4> qconj(const vec<T, 4>& q) {
1891 return {-q.x, -q.y, -q.z, q.w};
1892}
1893
1898template <class T>
1900 return qconj(q) / length2(q);
1901}
1902
1907template <class T>
1909 const auto v = q.xyz();
1910 const auto vv = length(v);
1911 return std::exp(q.w) *
1912 vec<T, 4>{v * (vv > 0 ? std::sin(vv) / vv : 0), std::cos(vv)};
1913}
1914
1919template <class T>
1921 const auto v = q.xyz();
1922 const auto vv = length(v), qq = length(q);
1923 return {v * (vv > 0 ? std::acos(q.w / qq) / vv : 0), std::log(qq)};
1924}
1925
1928template <class T>
1929vec<T, 4> qpow(const vec<T, 4>& q, const T& p) {
1930 const auto v = q.xyz();
1931 const auto vv = length(v), qq = length(q), th = std::acos(q.w / qq);
1932 return std::pow(qq, p) *
1933 vec<T, 4>{v * (vv > 0 ? std::sin(p * th) / vv : 0), std::cos(p * th)};
1934}
1935
1940template <class T>
1941constexpr vec<T, 4> qmul(const vec<T, 4>& a, const vec<T, 4>& b) {
1942 return {a.x * b.w + a.w * b.x + a.y * b.z - a.z * b.y,
1943 a.y * b.w + a.w * b.y + a.z * b.x - a.x * b.z,
1944 a.z * b.w + a.w * b.z + a.x * b.y - a.y * b.x,
1945 a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z};
1946}
1947
1950template <class T, class... R>
1951constexpr vec<T, 4> qmul(const vec<T, 4>& a, R... r) {
1952 return qmul(a, qmul(r...));
1953}
1954
1955
1964template <class T>
1965constexpr vec<T, 3> qxdir(const vec<T, 4>& q) {
1966 return {q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
1967 (q.x * q.y + q.z * q.w) * 2, (q.z * q.x - q.y * q.w) * 2};
1968}
1969
1972template <class T>
1973constexpr vec<T, 3> qydir(const vec<T, 4>& q) {
1974 return {(q.x * q.y - q.z * q.w) * 2,
1975 q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
1976 (q.y * q.z + q.x * q.w) * 2};
1977}
1978
1981template <class T>
1982constexpr vec<T, 3> qzdir(const vec<T, 4>& q) {
1983 return {(q.z * q.x + q.y * q.w) * 2, (q.y * q.z - q.x * q.w) * 2,
1984 q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z};
1985}
1986
1989template <class T>
1990constexpr mat<T, 3, 3> qmat(const vec<T, 4>& q) {
1991 return {qxdir(q), qydir(q), qzdir(q)};
1992}
1993
1996template <class T>
1997constexpr vec<T, 3> qrot(const vec<T, 4>& q, const vec<T, 3>& v) {
1998 return qxdir(q) * v.x + qydir(q) * v.y + qzdir(q) * v.z;
1999}
2000
2004template <class T>
2005T qangle(const vec<T, 4>& q) {
2006 return std::atan2(length(q.xyz()), q.w) * 2;
2007}
2008
2012template <class T>
2014 return normalize(q.xyz());
2015}
2016
2020template <class T>
2021vec<T, 4> qnlerp(const vec<T, 4>& a, const vec<T, 4>& b, T t) {
2022 return nlerp(a, dot(a, b) < 0 ? -b : b, t);
2023}
2024
2027template <class T>
2028vec<T, 4> qslerp(const vec<T, 4>& a, const vec<T, 4>& b, T t) {
2029 return slerp(a, dot(a, b) < 0 ? -b : b, t);
2030}
2031
2035template <class T>
2036vec<T, 4> constexpr rotation_quat(const vec<T, 3>& axis, T angle) {
2037 return {axis * std::sin(angle / 2), std::cos(angle / 2)};
2038}
2039
2043template <class T>
2044vec<T, 4> rotation_quat(const vec<T, 3>& orig, const vec<T, 3>& dest);
2049template <class T>
2052
2058template <class T, int M>
2059constexpr vec<T, M> mul(const vec<T, M>& a, const T& b) {
2060 return cmul(a, b);
2061}
2062template <class T, int M>
2063constexpr vec<T, M> mul(const T& b, const vec<T, M>& a) {
2064 return cmul(b, a);
2065}
2066template <class T, int M, int N>
2067constexpr mat<T, M, N> mul(const mat<T, M, N>& a, const T& b) {
2068 return cmul(a, b);
2069}
2070template <class T, int M, int N>
2071constexpr mat<T, M, N> mul(const T& b, const mat<T, M, N>& a) {
2072 return cmul(b, a);
2073}
2074template <class T, int M>
2075constexpr vec<T, M> mul(const vec<T, M>& a, const vec<T, M>& b) {
2076 return cmul(a, b);
2077}
2078template <class T, int M>
2079constexpr vec<T, M> mul(const mat<T, M, 1>& a, const vec<T, 1>& b) {
2080 return a.x * b.x;
2081}
2082template <class T, int M>
2083constexpr vec<T, M> mul(const mat<T, M, 2>& a, const vec<T, 2>& b) {
2084 return a.x * b.x + a.y * b.y;
2085}
2086template <class T, int M>
2087constexpr vec<T, M> mul(const mat<T, M, 3>& a, const vec<T, 3>& b) {
2088 return a.x * b.x + a.y * b.y + a.z * b.z;
2089}
2090template <class T, int M>
2091constexpr vec<T, M> mul(const mat<T, M, 4>& a, const vec<T, 4>& b) {
2092 return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
2093}
2094template <class T, int M, int N>
2095constexpr mat<T, M, 1> mul(const mat<T, M, N>& a, const mat<T, N, 1>& b) {
2096 return {mul(a, b.x)};
2097}
2098template <class T, int M, int N>
2099constexpr mat<T, M, 2> mul(const mat<T, M, N>& a, const mat<T, N, 2>& b) {
2100 return {mul(a, b.x), mul(a, b.y)};
2101}
2102template <class T, int M, int N>
2103constexpr mat<T, M, 3> mul(const mat<T, M, N>& a, const mat<T, N, 3>& b) {
2104 return {mul(a, b.x), mul(a, b.y), mul(a, b.z)};
2105}
2106template <class T, int M, int N>
2107constexpr mat<T, M, 4> mul(const mat<T, M, N>& a, const mat<T, N, 4>& b) {
2108 return {mul(a, b.x), mul(a, b.y), mul(a, b.z), mul(a, b.w)};
2109}
2110template <class T, int M, int N, int P>
2111constexpr vec<T, M> mul(const mat<T, M, N>& a, const mat<T, N, P>& b,
2112 const vec<T, P>& c) {
2113 return mul(mul(a, b), c);
2114}
2115template <class T, int M, int N, int P, int Q>
2116constexpr mat<T, M, Q> mul(const mat<T, M, N>& a, const mat<T, N, P>& b,
2117 const mat<T, P, Q>& c) {
2118 return mul(mul(a, b), c);
2119}
2120template <class T, int M, int N, int P, int Q>
2121constexpr vec<T, M> mul(const mat<T, M, N>& a, const mat<T, N, P>& b,
2122 const mat<T, P, Q>& c, const vec<T, Q>& d) {
2123 return mul(mul(a, b, c), d);
2124}
2125template <class T, int M, int N, int P, int Q, int R>
2126constexpr mat<T, M, R> mul(const mat<T, M, N>& a, const mat<T, N, P>& b,
2127 const mat<T, P, Q>& c, const mat<T, Q, R>& d) {
2128 return mul(mul(a, b, c), d);
2129}
2130template <class T, int M>
2131constexpr mat<T, M, 1> outerprod(const vec<T, M>& a, const vec<T, 1>& b) {
2132 return {a * b.x};
2133}
2134template <class T, int M>
2135constexpr mat<T, M, 2> outerprod(const vec<T, M>& a, const vec<T, 2>& b) {
2136 return {a * b.x, a * b.y};
2137}
2138template <class T, int M>
2139constexpr mat<T, M, 3> outerprod(const vec<T, M>& a, const vec<T, 3>& b) {
2140 return {a * b.x, a * b.y, a * b.z};
2141}
2142template <class T, int M>
2143constexpr mat<T, M, 4> outerprod(const vec<T, M>& a, const vec<T, 4>& b) {
2144 return {a * b.x, a * b.y, a * b.z, a * b.w};
2145}
2146template <class T>
2147constexpr vec<T, 1> diagonal(const mat<T, 1, 1>& a) {
2148 return {a.x.x};
2149}
2150template <class T>
2151constexpr vec<T, 2> diagonal(const mat<T, 2, 2>& a) {
2152 return {a.x.x, a.y.y};
2153}
2154template <class T>
2155constexpr vec<T, 3> diagonal(const mat<T, 3, 3>& a) {
2156 return {a.x.x, a.y.y, a.z.z};
2157}
2158template <class T>
2159constexpr vec<T, 4> diagonal(const mat<T, 4, 4>& a) {
2160 return {a.x.x, a.y.y, a.z.z, a.w.w};
2161}
2162template <class T, int N>
2163constexpr T trace(const mat<T, N, N>& a) {
2164 return sum(diagonal(a));
2165}
2166template <class T, int M>
2167constexpr mat<T, M, 1> transpose(const mat<T, 1, M>& m) {
2168 return {m.row(0)};
2169}
2170template <class T, int M>
2171constexpr mat<T, M, 2> transpose(const mat<T, 2, M>& m) {
2172 return {m.row(0), m.row(1)};
2173}
2174template <class T, int M>
2175constexpr mat<T, M, 3> transpose(const mat<T, 3, M>& m) {
2176 return {m.row(0), m.row(1), m.row(2)};
2177}
2178template <class T, int M>
2179constexpr mat<T, M, 4> transpose(const mat<T, 4, M>& m) {
2180 return {m.row(0), m.row(1), m.row(2), m.row(3)};
2181}
2182template <class T, int M>
2183constexpr mat<T, 1, M> transpose(const vec<T, M>& m) {
2184 return transpose(mat<T, M, 1>(m));
2185}
2186template <class T>
2187constexpr mat<T, 1, 1> adjugate(const mat<T, 1, 1>&) {
2188 return {vec<T, 1>{1}};
2189}
2190template <class T>
2191constexpr mat<T, 2, 2> adjugate(const mat<T, 2, 2>& a) {
2192 return {{a.y.y, -a.x.y}, {-a.y.x, a.x.x}};
2193}
2194template <class T>
2195constexpr mat<T, 3, 3> adjugate(const mat<T, 3, 3>& a);
2196template <class T>
2197constexpr mat<T, 4, 4> adjugate(const mat<T, 4, 4>& a);
2198template <class T, int N>
2199constexpr mat<T, N, N> comatrix(const mat<T, N, N>& a) {
2200 return transpose(adjugate(a));
2201}
2202template <class T>
2203constexpr T determinant(const mat<T, 1, 1>& a) {
2204 return a.x.x;
2205}
2206template <class T>
2207constexpr T determinant(const mat<T, 2, 2>& a) {
2208 return a.x.x * a.y.y - a.x.y * a.y.x;
2209}
2210template <class T>
2211constexpr T determinant(const mat<T, 3, 3>& a) {
2212 return a.x.x * (a.y.y * a.z.z - a.z.y * a.y.z) +
2213 a.x.y * (a.y.z * a.z.x - a.z.z * a.y.x) +
2214 a.x.z * (a.y.x * a.z.y - a.z.x * a.y.y);
2215}
2216template <class T>
2217constexpr T determinant(const mat<T, 4, 4>& a);
2218template <class T, int N>
2219constexpr mat<T, N, N> inverse(const mat<T, N, N>& a) {
2220 return adjugate(a) / determinant(a);
2221}
2223
2229template <class T, int M>
2230T* begin(vec<T, M>& a) {
2231 return &a.x;
2232}
2233template <class T, int M>
2234const T* begin(const vec<T, M>& a) {
2235 return &a.x;
2236}
2237template <class T, int M>
2238T* end(vec<T, M>& a) {
2239 return begin(a) + M;
2240}
2241template <class T, int M>
2242const T* end(const vec<T, M>& a) {
2243 return begin(a) + M;
2244}
2245template <class T, int M, int N>
2246vec<T, M>* begin(mat<T, M, N>& a) {
2247 return &a.x;
2248}
2249template <class T, int M, int N>
2250const vec<T, M>* begin(const mat<T, M, N>& a) {
2251 return &a.x;
2252}
2253template <class T, int M, int N>
2254vec<T, M>* end(mat<T, M, N>& a) {
2255 return begin(a) + N;
2256}
2257template <class T, int M, int N>
2258const vec<T, M>* end(const mat<T, M, N>& a) {
2259 return begin(a) + N;
2260}
2262
2268enum fwd_axis {
2269 neg_z,
2270 pos_z
2271}; // Should projection matrices be generated assuming forward is {0,0,-1} or
2272 // {0,0,1}
2273enum z_range {
2274 neg_one_to_one,
2275 zero_to_one
2276}; // Should projection matrices map z into the range of [-1,1] or [0,1]?
2277template <class T>
2278mat<T, 4, 4> translation_matrix(const vec<T, 3>& translation) {
2279 return {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {translation, 1}};
2280}
2281template <class T>
2282mat<T, 4, 4> rotation_matrix(const vec<T, 4>& rotation) {
2283 return {{qxdir(rotation), 0},
2284 {qydir(rotation), 0},
2285 {qzdir(rotation), 0},
2286 {0, 0, 0, 1}};
2287}
2288template <class T>
2289mat<T, 4, 4> scaling_matrix(const vec<T, 3>& scaling) {
2290 return {{scaling.x, 0, 0, 0},
2291 {0, scaling.y, 0, 0},
2292 {0, 0, scaling.z, 0},
2293 {0, 0, 0, 1}};
2294}
2295template <class T>
2296mat<T, 4, 4> pose_matrix(const vec<T, 4>& q, const vec<T, 3>& p) {
2297 return {{qxdir(q), 0}, {qydir(q), 0}, {qzdir(q), 0}, {p, 1}};
2298}
2299template <class T>
2300mat<T, 4, 4> lookat_matrix(const vec<T, 3>& eye, const vec<T, 3>& center,
2301 const vec<T, 3>& view_y_dir, fwd_axis fwd = neg_z);
2302template <class T>
2303mat<T, 4, 4> frustum_matrix(T x0, T x1, T y0, T y1, T n, T f,
2304 fwd_axis a = neg_z, z_range z = neg_one_to_one);
2305template <class T>
2306mat<T, 4, 4> perspective_matrix(T fovy, T aspect, T n, T f, fwd_axis a = neg_z,
2307 z_range z = neg_one_to_one) {
2308 T y = n * std::tan(fovy / 2), x = y * aspect;
2309 return frustum_matrix(-x, x, -y, y, n, f, a, z);
2310}
2312
2319template <class T>
2320struct converter<vec<T, 1>, std::array<T, 1>> {
2321 vec<T, 1> operator()(const std::array<T, 1>& a) const { return {a[0]}; }
2322};
2323template <class T>
2324struct converter<vec<T, 2>, std::array<T, 2>> {
2325 vec<T, 2> operator()(const std::array<T, 2>& a) const { return {a[0], a[1]}; }
2326};
2327template <class T>
2328struct converter<vec<T, 3>, std::array<T, 3>> {
2329 vec<T, 3> operator()(const std::array<T, 3>& a) const {
2330 return {a[0], a[1], a[2]};
2331 }
2332};
2333template <class T>
2334struct converter<vec<T, 4>, std::array<T, 4>> {
2335 vec<T, 4> operator()(const std::array<T, 4>& a) const {
2336 return {a[0], a[1], a[2], a[3]};
2337 }
2338};
2339
2340template <class T>
2341struct converter<std::array<T, 1>, vec<T, 1>> {
2342 std::array<T, 1> operator()(const vec<T, 1>& a) const { return {a[0]}; }
2343};
2344template <class T>
2345struct converter<std::array<T, 2>, vec<T, 2>> {
2346 std::array<T, 2> operator()(const vec<T, 2>& a) const { return {a[0], a[1]}; }
2347};
2348template <class T>
2349struct converter<std::array<T, 3>, vec<T, 3>> {
2350 std::array<T, 3> operator()(const vec<T, 3>& a) const {
2351 return {a[0], a[1], a[2]};
2352 }
2353};
2354template <class T>
2355struct converter<std::array<T, 4>, vec<T, 4>> {
2356 std::array<T, 4> operator()(const vec<T, 4>& a) const {
2357 return {a[0], a[1], a[2], a[3]};
2358 }
2359};
2360
2361} // namespace linalg
2362
2363namespace std {
2369template <class T>
2370struct hash<linalg::vec<T, 1>> {
2371 std::size_t operator()(const linalg::vec<T, 1>& v) const {
2372 std::hash<T> h;
2373 return h(v.x);
2374 }
2375};
2376template <class T>
2377struct hash<linalg::vec<T, 2>> {
2378 std::size_t operator()(const linalg::vec<T, 2>& v) const {
2379 std::hash<T> h;
2380 return h(v.x) ^ (h(v.y) << 1);
2381 }
2382};
2383template <class T>
2384struct hash<linalg::vec<T, 3>> {
2385 std::size_t operator()(const linalg::vec<T, 3>& v) const {
2386 std::hash<T> h;
2387 return h(v.x) ^ (h(v.y) << 1) ^ (h(v.z) << 2);
2388 }
2389};
2390template <class T>
2391struct hash<linalg::vec<T, 4>> {
2392 std::size_t operator()(const linalg::vec<T, 4>& v) const {
2393 std::hash<T> h;
2394 return h(v.x) ^ (h(v.y) << 1) ^ (h(v.z) << 2) ^ (h(v.w) << 3);
2395 }
2396};
2397
2398template <class T, int M>
2399struct hash<linalg::mat<T, M, 1>> {
2400 std::size_t operator()(const linalg::mat<T, M, 1>& v) const {
2401 std::hash<linalg::vec<T, M>> h;
2402 return h(v.x);
2403 }
2404};
2405template <class T, int M>
2406struct hash<linalg::mat<T, M, 2>> {
2407 std::size_t operator()(const linalg::mat<T, M, 2>& v) const {
2408 std::hash<linalg::vec<T, M>> h;
2409 return h(v.x) ^ (h(v.y) << M);
2410 }
2411};
2412template <class T, int M>
2413struct hash<linalg::mat<T, M, 3>> {
2414 std::size_t operator()(const linalg::mat<T, M, 3>& v) const {
2415 std::hash<linalg::vec<T, M>> h;
2416 return h(v.x) ^ (h(v.y) << M) ^ (h(v.z) << (M * 2));
2417 }
2418};
2419template <class T, int M>
2420struct hash<linalg::mat<T, M, 4>> {
2421 std::size_t operator()(const linalg::mat<T, M, 4>& v) const {
2422 std::hash<linalg::vec<T, M>> h;
2423 return h(v.x) ^ (h(v.y) << M) ^ (h(v.z) << (M * 2)) ^ (h(v.w) << (M * 3));
2424 }
2425};
2426
2427} // namespace std
2428
2429// Definitions of functions too long to be defined inline
2430template <class T>
2431constexpr linalg::mat<T, 3, 3> linalg::adjugate(const mat<T, 3, 3>& a) {
2432 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,
2433 a.x.y * a.y.z - a.y.y * a.x.z},
2434 {a.y.z * a.z.x - a.z.z * a.y.x, a.z.z * a.x.x - a.x.z * a.z.x,
2435 a.x.z * a.y.x - a.y.z * a.x.x},
2436 {a.y.x * a.z.y - a.z.x * a.y.y, a.z.x * a.x.y - a.x.x * a.z.y,
2437 a.x.x * a.y.y - a.y.x * a.x.y}};
2438}
2439
2440template <class T>
2441constexpr linalg::mat<T, 4, 4> linalg::adjugate(const mat<T, 4, 4>& a) {
2442 return {{a.y.y * a.z.z * a.w.w + a.w.y * a.y.z * a.z.w +
2443 a.z.y * a.w.z * a.y.w - a.y.y * a.w.z * a.z.w -
2444 a.z.y * a.y.z * a.w.w - a.w.y * a.z.z * a.y.w,
2445 a.x.y * a.w.z * a.z.w + a.z.y * a.x.z * a.w.w +
2446 a.w.y * a.z.z * a.x.w - a.w.y * a.x.z * a.z.w -
2447 a.z.y * a.w.z * a.x.w - a.x.y * a.z.z * a.w.w,
2448 a.x.y * a.y.z * a.w.w + a.w.y * a.x.z * a.y.w +
2449 a.y.y * a.w.z * a.x.w - a.x.y * a.w.z * a.y.w -
2450 a.y.y * a.x.z * a.w.w - a.w.y * a.y.z * a.x.w,
2451 a.x.y * a.z.z * a.y.w + a.y.y * a.x.z * a.z.w +
2452 a.z.y * a.y.z * a.x.w - a.x.y * a.y.z * a.z.w -
2453 a.z.y * a.x.z * a.y.w - a.y.y * a.z.z * a.x.w},
2454 {a.y.z * a.w.w * a.z.x + a.z.z * a.y.w * a.w.x +
2455 a.w.z * a.z.w * a.y.x - a.y.z * a.z.w * a.w.x -
2456 a.w.z * a.y.w * a.z.x - a.z.z * a.w.w * a.y.x,
2457 a.x.z * a.z.w * a.w.x + a.w.z * a.x.w * a.z.x +
2458 a.z.z * a.w.w * a.x.x - a.x.z * a.w.w * a.z.x -
2459 a.z.z * a.x.w * a.w.x - a.w.z * a.z.w * a.x.x,
2460 a.x.z * a.w.w * a.y.x + a.y.z * a.x.w * a.w.x +
2461 a.w.z * a.y.w * a.x.x - a.x.z * a.y.w * a.w.x -
2462 a.w.z * a.x.w * a.y.x - a.y.z * a.w.w * a.x.x,
2463 a.x.z * a.y.w * a.z.x + a.z.z * a.x.w * a.y.x +
2464 a.y.z * a.z.w * a.x.x - a.x.z * a.z.w * a.y.x -
2465 a.y.z * a.x.w * a.z.x - a.z.z * a.y.w * a.x.x},
2466 {a.y.w * a.z.x * a.w.y + a.w.w * a.y.x * a.z.y +
2467 a.z.w * a.w.x * a.y.y - a.y.w * a.w.x * a.z.y -
2468 a.z.w * a.y.x * a.w.y - a.w.w * a.z.x * a.y.y,
2469 a.x.w * a.w.x * a.z.y + a.z.w * a.x.x * a.w.y +
2470 a.w.w * a.z.x * a.x.y - a.x.w * a.z.x * a.w.y -
2471 a.w.w * a.x.x * a.z.y - a.z.w * a.w.x * a.x.y,
2472 a.x.w * a.y.x * a.w.y + a.w.w * a.x.x * a.y.y +
2473 a.y.w * a.w.x * a.x.y - a.x.w * a.w.x * a.y.y -
2474 a.y.w * a.x.x * a.w.y - a.w.w * a.y.x * a.x.y,
2475 a.x.w * a.z.x * a.y.y + a.y.w * a.x.x * a.z.y +
2476 a.z.w * a.y.x * a.x.y - a.x.w * a.y.x * a.z.y -
2477 a.z.w * a.x.x * a.y.y - a.y.w * a.z.x * a.x.y},
2478 {a.y.x * a.w.y * a.z.z + a.z.x * a.y.y * a.w.z +
2479 a.w.x * a.z.y * a.y.z - a.y.x * a.z.y * a.w.z -
2480 a.w.x * a.y.y * a.z.z - a.z.x * a.w.y * a.y.z,
2481 a.x.x * a.z.y * a.w.z + a.w.x * a.x.y * a.z.z +
2482 a.z.x * a.w.y * a.x.z - a.x.x * a.w.y * a.z.z -
2483 a.z.x * a.x.y * a.w.z - a.w.x * a.z.y * a.x.z,
2484 a.x.x * a.w.y * a.y.z + a.y.x * a.x.y * a.w.z +
2485 a.w.x * a.y.y * a.x.z - a.x.x * a.y.y * a.w.z -
2486 a.w.x * a.x.y * a.y.z - a.y.x * a.w.y * a.x.z,
2487 a.x.x * a.y.y * a.z.z + a.z.x * a.x.y * a.y.z +
2488 a.y.x * a.z.y * a.x.z - a.x.x * a.z.y * a.y.z -
2489 a.y.x * a.x.y * a.z.z - a.z.x * a.y.y * a.x.z}};
2490}
2491
2492template <class T>
2493constexpr T linalg::determinant(const mat<T, 4, 4>& a) {
2494 return a.x.x * (a.y.y * a.z.z * a.w.w + a.w.y * a.y.z * a.z.w +
2495 a.z.y * a.w.z * a.y.w - a.y.y * a.w.z * a.z.w -
2496 a.z.y * a.y.z * a.w.w - a.w.y * a.z.z * a.y.w) +
2497 a.x.y * (a.y.z * a.w.w * a.z.x + a.z.z * a.y.w * a.w.x +
2498 a.w.z * a.z.w * a.y.x - a.y.z * a.z.w * a.w.x -
2499 a.w.z * a.y.w * a.z.x - a.z.z * a.w.w * a.y.x) +
2500 a.x.z * (a.y.w * a.z.x * a.w.y + a.w.w * a.y.x * a.z.y +
2501 a.z.w * a.w.x * a.y.y - a.y.w * a.w.x * a.z.y -
2502 a.z.w * a.y.x * a.w.y - a.w.w * a.z.x * a.y.y) +
2503 a.x.w * (a.y.x * a.w.y * a.z.z + a.z.x * a.y.y * a.w.z +
2504 a.w.x * a.z.y * a.y.z - a.y.x * a.z.y * a.w.z -
2505 a.w.x * a.y.y * a.z.z - a.z.x * a.w.y * a.y.z);
2506}
2507
2508template <class T>
2509linalg::vec<T, 4> linalg::rotation_quat(const vec<T, 3>& orig,
2510 const vec<T, 3>& dest) {
2511 T cosTheta = dot(orig, dest);
2512 if (cosTheta >= 1 - std::numeric_limits<T>::epsilon()) {
2513 return {0, 0, 0, 1};
2514 }
2515 if (cosTheta < -1 + std::numeric_limits<T>::epsilon()) {
2516 vec<T, 3> axis = cross(vec<T, 3>(0, 0, 1), orig);
2517 if (length2(axis) < std::numeric_limits<T>::epsilon())
2518 axis = cross(vec<T, 3>(1, 0, 0), orig);
2519 return rotation_quat(normalize(axis),
2520 3.14159265358979323846264338327950288);
2521 }
2522 vec<T, 3> axis = cross(orig, dest);
2523 T s = std::sqrt((1 + cosTheta) * 2);
2524 return {axis * (1 / s), s * 0.5};
2525}
2526
2527template <class T>
2528linalg::vec<T, 4> linalg::rotation_quat(const mat<T, 3, 3>& m) {
2529 const vec<T, 4> q{m.x.x - m.y.y - m.z.z, m.y.y - m.x.x - m.z.z,
2530 m.z.z - m.x.x - m.y.y, m.x.x + m.y.y + m.z.z},
2531 s[]{{1, m.x.y + m.y.x, m.z.x + m.x.z, m.y.z - m.z.y},
2532 {m.x.y + m.y.x, 1, m.y.z + m.z.y, m.z.x - m.x.z},
2533 {m.x.z + m.z.x, m.y.z + m.z.y, 1, m.x.y - m.y.x},
2534 {m.y.z - m.z.y, m.z.x - m.x.z, m.x.y - m.y.x, 1}};
2535 return copysign(normalize(sqrt(max(T(0), T(1) + q))), s[argmax(q)]);
2536}
2537
2538template <class T>
2539linalg::mat<T, 4, 4> linalg::lookat_matrix(const vec<T, 3>& eye,
2540 const vec<T, 3>& center,
2541 const vec<T, 3>& view_y_dir,
2542 fwd_axis a) {
2543 const vec<T, 3> f = normalize(center - eye), z = a == pos_z ? f : -f,
2544 x = normalize(cross(view_y_dir, z)), y = cross(z, x);
2545 return inverse(mat<T, 4, 4>{{x, 0}, {y, 0}, {z, 0}, {eye, 1}});
2546}
2547
2548template <class T>
2549linalg::mat<T, 4, 4> linalg::frustum_matrix(T x0, T x1, T y0, T y1, T n, T f,
2550 fwd_axis a, z_range z) {
2551 const T s = a == pos_z ? T(1) : T(-1), o = z == neg_one_to_one ? n : 0;
2552 return {{2 * n / (x1 - x0), 0, 0, 0},
2553 {0, 2 * n / (y1 - y0), 0, 0},
2554 {-s * (x0 + x1) / (x1 - x0), -s * (y0 + y1) / (y1 - y0),
2555 s * (f + o) / (f - n), s},
2556 {0, 0, -(n + o) * f / (f - n), 0}};
2557}
2558#endif
constexpr mat< T, 3, 3 > qmat(const vec< T, 4 > &q)
Create an equivalent mat3 rotation matrix from the input quaternion.
Definition linalg.h:1990
vec< T, 4 > qslerp(const vec< T, 4 > &a, const vec< T, 4 > &b, T t)
Spherical linear interpolation that takes the shortest path.
Definition linalg.h:2028
T qangle(const vec< T, 4 > &q)
Return the angle in radians of the axis-angle representation of the input normalized quaternion.
Definition linalg.h:2005
constexpr vec< T, 3 > qzdir(const vec< T, 4 > &q)
efficient shorthand for qrot(q, {0,0,1})
Definition linalg.h:1982
vec< T, 4 > qnlerp(const vec< T, 4 > &a, const vec< T, 4 > &b, T t)
Linear interpolation that takes the shortest path - this is not geometrically sensible,...
Definition linalg.h:2021
vec< T, 3 > qaxis(const vec< T, 4 > &q)
Return the normalized axis of the axis-angle representation of the input normalized quaternion.
Definition linalg.h:2013
vec< T, 4 > constexpr rotation_quat(const vec< T, 3 > &axis, T angle)
Returns a normalized quaternion representing a rotation by angle in radians about the provided axis.
Definition linalg.h:2036
constexpr vec< T, 3 > qxdir(const vec< T, 4 > &q)
efficient shorthand for qrot(q, {1,0,0})
Definition linalg.h:1965
constexpr vec< T, 3 > qrot(const vec< T, 4 > &q, const vec< T, 3 > &v)
Rotate a vector by a quaternion.
Definition linalg.h:1997
constexpr vec< T, 3 > qydir(const vec< T, 4 > &q)
efficient shorthand for qrot(q, {0,1,0})
Definition linalg.h:1973
vec< T, 4 > qinv(const vec< T, 4 > &q)
inverse or reciprocal of quaternion q (undefined for zero-length quaternions)
Definition linalg.h:1899
vec< T, 4 > qexp(const vec< T, 4 > &q)
exponential of quaternion q
Definition linalg.h:1908
vec< T, 4 > qpow(const vec< T, 4 > &q, const T &p)
quaternion q raised to the exponent p
Definition linalg.h:1929
vec< T, 4 > qlog(const vec< T, 4 > &q)
logarithm of quaternion q
Definition linalg.h:1920
constexpr vec< T, 4 > qmul(const vec< T, 4 > &a, const vec< T, 4 > &b)
Hamilton product of quaternions a and b
Definition linalg.h:1941
constexpr vec< T, 4 > qconj(const vec< T, 4 > &q)
conjugate of quaternion q
Definition linalg.h:1890
constexpr apply_t< detail::clamp, X, L, H > clamp(const X &x, const L &l, const H &h)
Clamps the components of x between l and h, provided l[i] < h[i].
Definition linalg.h:1700
constexpr apply_t< detail::select, P, A, B > select(const P &p, const A &a, const B &b)
Returns the component from a if the corresponding component of p is true and from b otherwise.
Definition linalg.h:1709
constexpr apply_t< detail::lerp, A, B, T > lerp(const A &a, const B &b, const T &t)
Linear interpolation from a to b as t goes from 0 -> 1. Values beyond [a, b] will result if t is outs...
Definition linalg.h:1718
constexpr vec< T, sizeof...(I)> swizzle(const vec< T, M > &a)
Returns a vector containing the specified ordered indices, e.g. linalg::swizzle<1,...
Definition linalg.h:1521
constexpr mat< T, I1 - I0, J1 - J0 > submat(const mat< T, M, N > &a)
Returns a matrix containing the specified row and column range: linalg::submat<rowStart,...
Definition linalg.h:1537
constexpr vec< T, I1 - I0 > subvec(const vec< T, M > &a)
Returns a vector containing the specified index range, e.g. linalg::subvec<1, 4>(vec4(4,...
Definition linalg.h:1529
vec< T, 3 > roty(T a, const vec< T, 3 > &v)
vector v rotated counter-clockwise by the angle a in radians around the Y axis
Definition linalg.h:1844
vec< T, 2 > rot(T a, const vec< T, 2 > &v)
vector v rotated counter-clockwise by the angle a in radians
Definition linalg.h:1826
constexpr T distance2(const vec< T, M > &a, const vec< T, M > &b)
square of the Euclidean distance between points a and b
Definition linalg.h:1795
vec< T, 3 > rotx(T a, const vec< T, 3 > &v)
vector v rotated counter-clockwise by the angle a in radians around the X axis
Definition linalg.h:1835
constexpr T length2(const vec< T, M > &a)
square of the length or magnitude of vector a
Definition linalg.h:1770
constexpr T dot(const vec< T, M > &a, const vec< T, M > &b)
dot or inner product of vectors a and b
Definition linalg.h:1763
vec< T, M > nlerp(const vec< T, M > &a, const vec< T, M > &b, T t)
shorthand for normalize(lerp(a,b,t))
Definition linalg.h:1861
vec< T, M > normalize(const vec< T, M > &a)
unit length vector in the same direction as a (undefined for zero-length vectors)
Definition linalg.h:1786
vec< T, M > slerp(const vec< T, M > &a, const vec< T, M > &b, T t)
spherical linear interpolation between unit vectors a and b (undefined for non-unit vectors) by param...
Definition linalg.h:1870
T angle(const vec< T, M > &a, const vec< T, M > &b)
Return the angle in radians between two non-unit vectors.
Definition linalg.h:1818
T uangle(const vec< T, M > &a, const vec< T, M > &b)
Return the angle in radians between two unit vectors.
Definition linalg.h:1810
constexpr T cross(const vec< T, 2 > &a, const vec< T, 2 > &b)
shorthand for cross({a.x,a.y,0}, {b.x,b.y,0}).z
Definition linalg.h:1733
vec< T, 3 > rotz(T a, const vec< T, 3 > &v)
vector v rotated counter-clockwise by the angle a in radians around the Z axis
Definition linalg.h:1853
T length(const vec< T, M > &a)
length or magnitude of a vector a
Definition linalg.h:1777
T distance(const vec< T, M > &a, const vec< T, M > &b)
Euclidean distance between points a and b
Definition linalg.h:1803
Definition linalg.h:66
Definition linalg.h:118
Definition linalg.h:281
Definition linalg.h:449
Definition linalg.h:265
Definition linalg.h:193
Definition linalg.h:465
Definition linalg.h:227
Definition linalg.h:442
Definition linalg.h:435
Definition linalg.h:516
Definition linalg.h:594
Definition linalg.h:492
Definition linalg.h:504
Definition linalg.h:564
Definition linalg.h:558
Definition linalg.h:546
Definition linalg.h:576
Definition linalg.h:552
Definition linalg.h:528
Definition linalg.h:540
Definition linalg.h:510
Definition linalg.h:498
Definition linalg.h:570
Definition linalg.h:480
Definition linalg.h:486
Definition linalg.h:600
Definition linalg.h:474
Definition linalg.h:534
Definition linalg.h:522
Definition linalg.h:588
Definition linalg.h:582
Definition linalg.h:88
Definition linalg.h:75
Definition linalg.h:267
Definition linalg.h:458
Definition linalg.h:225
Definition linalg.h:614
Definition linalg.h:686
Definition linalg.h:680
Definition linalg.h:734
Definition linalg.h:692
Definition linalg.h:626
Definition linalg.h:740
Definition linalg.h:668
Definition linalg.h:704
Definition linalg.h:632
Definition linalg.h:620
Definition linalg.h:722
Definition linalg.h:608
Definition linalg.h:650
Definition linalg.h:644
Definition linalg.h:638
Definition linalg.h:728
Definition linalg.h:716
Definition linalg.h:662
Definition linalg.h:698
Definition linalg.h:656
Definition linalg.h:674
Definition linalg.h:710
Definition linalg.h:1176
Definition linalg.h:61
Definition linalg.h:56