aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CMakeLists.txt23
-rw-r--r--include/math/camera.h42
-rw-r--r--include/math/defs.h51
-rw-r--r--include/math/float.h9
-rw-r--r--include/math/fwd.h9
-rw-r--r--include/math/mat4.h488
-rw-r--r--include/math/spatial3.h194
-rw-r--r--include/math/vec.h16
-rw-r--r--include/math/vec2.h10
-rw-r--r--include/math/vec3.h143
-rw-r--r--include/math/vec4.h15
-rw-r--r--test/mat4_test.c111
-rw-r--r--test/test.h185
13 files changed, 1296 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..b9102d0
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,23 @@
1cmake_minimum_required(VERSION 3.16)
2
3project(math)
4
5# Library.
6
7add_library(math INTERFACE)
8
9target_include_directories(math INTERFACE
10 include)
11
12target_link_libraries(math INTERFACE
13 -lm)
14
15# Test.
16
17add_executable(math_test
18 test/mat4_test.c)
19
20target_link_libraries(math_test
21 math)
22
23target_compile_options(math_test PRIVATE -DUNIT_TEST -DNDEBUG -Wall -Wextra)
diff --git a/include/math/camera.h b/include/math/camera.h
new file mode 100644
index 0000000..1621927
--- /dev/null
+++ b/include/math/camera.h
@@ -0,0 +1,42 @@
1#pragma once
2
3#include "mat4.h"
4#include "spatial3.h"
5
6typedef struct Camera {
7 Spatial3 spatial;
8 mat4 projection;
9} Camera;
10
11/// Create an orthographic camera.
12///
13/// The camera is positioned at the origin with canonical right/up/forward
14/// vectors.
15///
16/// \param left The coordinate for the left vertical clipping plane.
17/// \param right The coordinate for the right vertical clipping plane.
18/// \param bottom The coordinate for the bottom horizontal clipping plane.
19/// \param top The coordinate for the top horizontal clipping plane.
20/// \param near The distance to the near clipping plane.
21/// \param far The distance to the far clipping plane.
22static inline Camera camera_orthographic(R left, R right, R bottom, R top,
23 R near, R far) {
24 return (Camera){.spatial = spatial3_make(),
25 .projection =
26 mat4_ortho(left, right, bottom, top, near, far)};
27}
28
29/// Create a perspective camera.
30///
31/// The camera is positioned at the origin with canonical right/up/forward
32/// vectors.
33///
34/// \param fovy The vertical field of view angle in degrees.
35/// \param aspect The aspect ratio that determines the field of view in the
36/// x-direction.
37/// \param near Distance to the near clipping plane.
38/// \param far Distance to the far clipping plane.
39static inline Camera camera_perspective(R fovy, R aspect, R near, R far) {
40 return (Camera){.spatial = spatial3_make(),
41 .projection = mat4_perspective(fovy, aspect, near, far)};
42}
diff --git a/include/math/defs.h b/include/math/defs.h
new file mode 100644
index 0000000..acfb037
--- /dev/null
+++ b/include/math/defs.h
@@ -0,0 +1,51 @@
1#pragma once
2
3//
4// Configuration macros:
5//
6// MATH_USE_FLOAT
7// - Use floats instead of doubles for scalar values.
8
9#include <float.h>
10#include <math.h>
11
12#ifdef MATH_USE_DOUBLE
13typedef double R;
14#define R_MIN DBL_MIN
15#define R_MAX DBL_MAX
16#else // floats
17typedef float R;
18#define R_MIN FLT_MIN
19#define R_MAX FLT_MAX
20#endif
21
22#define PI 3.14159265359
23#define INV_PI 0.31830988618
24
25/// Radians per degree.
26#define TO_RAD (PI / 180.0)
27
28/// Degrees per radian.
29#define TO_DEG (180.0 / PI)
30
31#ifdef MATH_USE_DOUBLE
32static inline double min(R a, R b) { return fmin(a, b); }
33static inline double max(R a, R b) { return fmax(a, b); }
34#else // floats
35static inline float min(R a, R b) { return fminf(a, b); }
36static inline float max(R a, R b) { return fmaxf(a, b); }
37#endif
38
39static inline R rabs(R x) { return x >= 0.0 ? x : -x; }
40static inline R clamp(R x, R low, R high) { return max(low, min(high, x)); }
41static inline R sq(R x) { return x * x; }
42static inline R sign(R x) {
43 if (x < 0) {
44 return -1;
45 } else if (x > 0) {
46 return 1;
47 } else {
48 return 0;
49 }
50}
51static inline R lerp(R a, R b, R t) { return a + (b - a) * t; }
diff --git a/include/math/float.h b/include/math/float.h
new file mode 100644
index 0000000..9d289b9
--- /dev/null
+++ b/include/math/float.h
@@ -0,0 +1,9 @@
1#pragma once
2
3#include <math.h>
4#include <stdbool.h>
5
6/// Compare two floats for equality.
7static inline bool float_eq(float a, float b, float eps) {
8 return fabsf(a - b) <= eps;
9}
diff --git a/include/math/fwd.h b/include/math/fwd.h
new file mode 100644
index 0000000..6c72ded
--- /dev/null
+++ b/include/math/fwd.h
@@ -0,0 +1,9 @@
1/// Forward declarations for all math objects.
2#pragma once
3
4typedef struct Camera Camera;
5typedef struct mat4 mat4;
6typedef struct spatial3 spatial3;
7typedef struct vec2 vec2;
8typedef struct vec3 vec3;
9typedef struct vec4 vec4;
diff --git a/include/math/mat4.h b/include/math/mat4.h
new file mode 100644
index 0000000..e6c707a
--- /dev/null
+++ b/include/math/mat4.h
@@ -0,0 +1,488 @@
1#pragma once
2
3#include "defs.h"
4#include "float.h"
5#include "vec3.h"
6#include "vec4.h"
7
8#include <stdbool.h>
9
10/// A 4x4 column-major matrix.
11typedef struct mat4 {
12 R val[4][4];
13} mat4;
14
15/// Construct a matrix from 16 values.
16static inline mat4 mat4_make(R m00, R m10, R m20, R m30, R m01, R m11, R m21,
17 R m31, R m02, R m12, R m22, R m32, R m03, R m13,
18 R m23, R m33) {
19 mat4 m;
20 m.val[0][0] = m00;
21 m.val[0][1] = m01;
22 m.val[0][2] = m02;
23 m.val[0][3] = m03;
24
25 m.val[1][0] = m10;
26 m.val[1][1] = m11;
27 m.val[1][2] = m12;
28 m.val[1][3] = m13;
29
30 m.val[2][0] = m20;
31 m.val[2][1] = m21;
32 m.val[2][2] = m22;
33 m.val[2][3] = m23;
34
35 m.val[3][0] = m30;
36 m.val[3][1] = m31;
37 m.val[3][2] = m32;
38 m.val[3][3] = m33;
39 return m;
40}
41
42/// Construct a matrix from a column-major matrix array.
43static inline mat4 mat4_from_array(const R M[16]) {
44 mat4 m;
45 m.val[0][0] = M[0];
46 m.val[0][1] = M[1];
47 m.val[0][2] = M[2];
48 m.val[0][3] = M[3];
49
50 m.val[1][0] = M[4];
51 m.val[1][1] = M[5];
52 m.val[1][2] = M[6];
53 m.val[1][3] = M[7];
54
55 m.val[2][0] = M[8];
56 m.val[2][1] = M[9];
57 m.val[2][2] = M[10];
58 m.val[2][3] = M[11];
59
60 m.val[3][0] = M[12];
61 m.val[3][1] = M[13];
62 m.val[3][2] = M[14];
63 m.val[3][3] = M[15];
64 return m;
65}
66
67/// Construct the identity matrix.
68static inline mat4 mat4_id() {
69 return mat4_make(1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
70 0.0, 0.0, 0.0, 1.0);
71}
72
73/// Construct a matrix from 4 column vectors.
74static inline mat4 mat4_from_vec4(vec4 v0, vec4 v1, vec4 v2, vec4 v3) {
75 return mat4_make(v0.x, v0.y, v0.z, v0.w, v1.x, v1.y, v1.z, v1.w, v2.x, v2.y,
76 v2.z, v2.w, v3.x, v3.y, v3.z, v3.w);
77}
78
79/// Construct a transformation matrix from 4 vectors.
80static inline mat4 mat4_from_vec3(vec3 right, vec3 up, vec3 forward,
81 vec3 position) {
82 return mat4_make(right.x, right.y, right.z, 0.0, up.x, up.y, up.z, 0.0,
83 forward.x, forward.y, forward.z, 0.0, position.x, position.y,
84 position.z, 1.0);
85}
86
87/// Return the value at the specified position.
88static inline R mat4_at(mat4 m, int row, int col) { return m.val[col][row]; }
89
90/// Return the matrix's first column.
91static inline vec3 mat4_v0(mat4 m) { return *((vec3*)m.val[0]); }
92
93/// Return the matrix's second column.
94static inline vec3 mat4_v1(mat4 m) { return *((vec3*)m.val[1]); }
95
96/// Return the matrix's third column.
97static inline vec3 mat4_v2(mat4 m) { return *((vec3*)m.val[2]); }
98
99/// Return the matrix's fourth column.
100static inline vec3 mat4_v3(mat4 m) { return *((vec3*)m.val[3]); }
101
102/// Set the matrix's first column.
103static inline void mat4_set_v0(mat4* m, vec3 v) { *((vec3*)m->val[0]) = v; }
104
105/// Set the matrix's second column.
106static inline void mat4_set_v1(mat4* m, vec3 v) { *((vec3*)m->val[1]) = v; }
107
108/// Set the matrix's third column.
109static inline void mat4_set_v2(mat4* m, vec3 v) { *((vec3*)m->val[2]) = v; }
110
111/// Set the matrix's fourth column.
112static inline void mat4_set_v3(mat4* m, vec3 v) { *((vec3*)m->val[3]) = v; }
113
114/// Set the 3x3 portion of the first matrix equal to the 3x3 portion of the
115/// second matrix.
116static inline void mat4_set_3x3(mat4* m, mat4 n) {
117 m->val[0][0] = n.val[0][0];
118 m->val[0][1] = n.val[0][1];
119 m->val[0][2] = n.val[0][2];
120
121 m->val[1][0] = n.val[1][0];
122 m->val[1][1] = n.val[1][1];
123 m->val[1][2] = n.val[1][2];
124
125 m->val[2][0] = n.val[2][0];
126 m->val[2][1] = n.val[2][1];
127 m->val[2][2] = n.val[2][2];
128}
129
130/// Multiply two matrices.
131/// A * B = AB.
132static inline mat4 mat4_mul(mat4 A, mat4 B) {
133 R m00 = mat4_at(A, 0, 0) * mat4_at(B, 0, 0) +
134 mat4_at(A, 0, 1) * mat4_at(B, 1, 0) +
135 mat4_at(A, 0, 2) * mat4_at(B, 2, 0) +
136 mat4_at(A, 0, 3) * mat4_at(B, 3, 0);
137 R m01 = mat4_at(A, 0, 0) * mat4_at(B, 0, 1) +
138 mat4_at(A, 0, 1) * mat4_at(B, 1, 1) +
139 mat4_at(A, 0, 2) * mat4_at(B, 2, 1) +
140 mat4_at(A, 0, 3) * mat4_at(B, 3, 1);
141 R m02 = mat4_at(A, 0, 0) * mat4_at(B, 0, 2) +
142 mat4_at(A, 0, 1) * mat4_at(B, 1, 2) +
143 mat4_at(A, 0, 2) * mat4_at(B, 2, 2) +
144 mat4_at(A, 0, 3) * mat4_at(B, 3, 2);
145 R m03 = mat4_at(A, 0, 0) * mat4_at(B, 0, 3) +
146 mat4_at(A, 0, 1) * mat4_at(B, 1, 3) +
147 mat4_at(A, 0, 2) * mat4_at(B, 2, 3) +
148 mat4_at(A, 0, 3) * mat4_at(B, 3, 3);
149
150 R m10 = mat4_at(A, 1, 0) * mat4_at(B, 0, 0) +
151 mat4_at(A, 1, 1) * mat4_at(B, 1, 0) +
152 mat4_at(A, 1, 2) * mat4_at(B, 2, 0) +
153 mat4_at(A, 1, 3) * mat4_at(B, 3, 0);
154 R m11 = mat4_at(A, 1, 0) * mat4_at(B, 0, 1) +
155 mat4_at(A, 1, 1) * mat4_at(B, 1, 1) +
156 mat4_at(A, 1, 2) * mat4_at(B, 2, 1) +
157 mat4_at(A, 1, 3) * mat4_at(B, 3, 1);
158 R m12 = mat4_at(A, 1, 0) * mat4_at(B, 0, 2) +
159 mat4_at(A, 1, 1) * mat4_at(B, 1, 2) +
160 mat4_at(A, 1, 2) * mat4_at(B, 2, 2) +
161 mat4_at(A, 1, 3) * mat4_at(B, 3, 2);
162 R m13 = mat4_at(A, 1, 0) * mat4_at(B, 0, 3) +
163 mat4_at(A, 1, 1) * mat4_at(B, 1, 3) +
164 mat4_at(A, 1, 2) * mat4_at(B, 2, 3) +
165 mat4_at(A, 1, 3) * mat4_at(B, 3, 3);
166
167 R m20 = mat4_at(A, 2, 0) * mat4_at(B, 0, 0) +
168 mat4_at(A, 2, 1) * mat4_at(B, 1, 0) +
169 mat4_at(A, 2, 2) * mat4_at(B, 2, 0) +
170 mat4_at(A, 2, 3) * mat4_at(B, 3, 0);
171 R m21 = mat4_at(A, 2, 0) * mat4_at(B, 0, 1) +
172 mat4_at(A, 2, 1) * mat4_at(B, 1, 1) +
173 mat4_at(A, 2, 2) * mat4_at(B, 2, 1) +
174 mat4_at(A, 2, 3) * mat4_at(B, 3, 1);
175 R m22 = mat4_at(A, 2, 0) * mat4_at(B, 0, 2) +
176 mat4_at(A, 2, 1) * mat4_at(B, 1, 2) +
177 mat4_at(A, 2, 2) * mat4_at(B, 2, 2) +
178 mat4_at(A, 2, 3) * mat4_at(B, 3, 2);
179 R m23 = mat4_at(A, 2, 0) * mat4_at(B, 0, 3) +
180 mat4_at(A, 2, 1) * mat4_at(B, 1, 3) +
181 mat4_at(A, 2, 2) * mat4_at(B, 2, 3) +
182 mat4_at(A, 2, 3) * mat4_at(B, 3, 3);
183
184 R m30 = mat4_at(A, 3, 0) * mat4_at(B, 0, 0) +
185 mat4_at(A, 3, 1) * mat4_at(B, 1, 0) +
186 mat4_at(A, 3, 2) * mat4_at(B, 2, 0) +
187 mat4_at(A, 3, 3) * mat4_at(B, 3, 0);
188 R m31 = mat4_at(A, 3, 0) * mat4_at(B, 0, 1) +
189 mat4_at(A, 3, 1) * mat4_at(B, 1, 1) +
190 mat4_at(A, 3, 2) * mat4_at(B, 2, 1) +
191 mat4_at(A, 3, 3) * mat4_at(B, 3, 1);
192 R m32 = mat4_at(A, 3, 0) * mat4_at(B, 0, 2) +
193 mat4_at(A, 3, 1) * mat4_at(B, 1, 2) +
194 mat4_at(A, 3, 2) * mat4_at(B, 2, 2) +
195 mat4_at(A, 3, 3) * mat4_at(B, 3, 2);
196 R m33 = mat4_at(A, 3, 0) * mat4_at(B, 0, 3) +
197 mat4_at(A, 3, 1) * mat4_at(B, 1, 3) +
198 mat4_at(A, 3, 2) * mat4_at(B, 2, 3) +
199 mat4_at(A, 3, 3) * mat4_at(B, 3, 3);
200
201 return mat4_make(m00, m01, m02, m03, m10, m11, m12, m13, m20, m21, m22, m23,
202 m30, m31, m32, m33);
203}
204
205/// Return the translation component of the matrix.
206static inline mat4 mat4_translation(mat4 m) {
207 return mat4_make(1.0, 0.0, 0.0, mat4_at(m, 0, 3), 0.0, 1.0, 0.0,
208 mat4_at(m, 1, 3), 0.0, 0.0, 1.0, mat4_at(m, 2, 3), 0.0, 0.0,
209 0.0, 1.0);
210}
211
212/// Return the rotation component of the matrix.
213static inline mat4 mat4_rotation(mat4 m) {
214 return mat4_make(mat4_at(m, 0, 0), mat4_at(m, 0, 1), mat4_at(m, 0, 2), 0.0,
215 mat4_at(m, 1, 0), mat4_at(m, 1, 1), mat4_at(m, 1, 2), 0.0,
216 mat4_at(m, 2, 0), mat4_at(m, 2, 1), mat4_at(m, 2, 2), 0.0,
217 0.0, 0.0, 0.0, 1.0);
218}
219
220/// Create an X-axis rotation matrix.
221static inline mat4 mat4_rotx(R angle) {
222 const R s = sin(angle);
223 const R c = cos(angle);
224 return mat4_make(1, 0, 0, 0, 0, c, -s, 0, 0, s, c, 0, 0, 0, 0, 1);
225}
226
227/// Create a Y-axis rotation matrix.
228static inline mat4 mat4_roty(R angle) {
229 const R s = sin(angle);
230 const R c = cos(angle);
231 return mat4_make(c, 0, s, 0, 0, 1, 0, 0, -s, 0, c, 0, 0, 0, 0, 1);
232}
233
234/// Create a Z-axis rotation matrix.
235static inline mat4 mat4_rotz(R angle) {
236 const R s = sin(angle);
237 const R c = cos(angle);
238 return mat4_make(c, -s, 0, 0, s, c, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
239}
240
241/// Create a rotation matrix.
242static inline mat4 mat4_rot(vec3 axis, R angle) {
243 const R s = sin(angle);
244 const R c = cos(angle);
245 const R x = axis.x;
246 const R y = axis.y;
247 const R z = axis.z;
248 const R xy = x * y;
249 const R xz = x * z;
250 const R yz = y * z;
251 const R sx = s * x;
252 const R sy = s * y;
253 const R sz = s * z;
254 const R omc = 1.0 - c;
255 return mat4_make(c + omc * x * x, omc * xy - sz, omc * xz + sy, 0,
256 omc * xy + sz, c + omc * y * y, omc * yz - sx, 0,
257 omc * xz - sy, omc * yz + sx, c + omc * z * z, 0, 0, 0, 0,
258 1);
259}
260
261/// Create a scaling matrix.
262static inline mat4 mat4_scale(vec3 s) {
263 return mat4_make(s.x, 0, 0, 0, 0, s.y, 0, 0, 0, 0, s.z, 0, 0, 0, 0, 1);
264}
265
266/// Create a translation matrix.
267static inline mat4 mat4_translate(vec3 v) {
268 return mat4_make(1, 0, 0, v.x, 0, 1, 0, v.y, 0, 0, 1, v.z, 0, 0, 0, 1);
269}
270
271/// The X-axis reflection matrix.
272static inline mat4 mat4_reflectx() {
273 return mat4_make(-1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
274}
275
276/// The Y-axis reflection matrix.
277static inline mat4 mat4_reflecty() {
278 return mat4_make(1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
279}
280
281/// The Z-axis reflection matrix.
282static inline mat4 mat4_reflectz() {
283 return mat4_make(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1);
284}
285
286/// Create a transformation matrix from the given forward vector.
287static inline mat4 mat4_from_forward(vec3 forward) {
288 const vec3 f = vec3_normalize(forward);
289 const vec3 r = vec3_normalize(vec3_cross(f, up3()));
290 const vec3 u = vec3_normalize(vec3_cross(r, f));
291 return mat4_make(r.x, u.x, -f.x, 0.0, r.y, u.y, -f.y, 0.0, r.z, u.z, -f.z,
292 0.0, 0.0, 0.0, 0.0, 1.0);
293}
294
295/// Create a transformation matrix.
296static inline mat4 mat4_lookat(vec3 position, vec3 target, vec3 up) {
297 const vec3 fwd = vec3_normalize(vec3_sub(target, position));
298 const vec3 right = vec3_normalize(vec3_cross(fwd, up));
299 up = vec3_normalize(vec3_cross(right, fwd));
300 return mat4_from_vec3(right, up, fwd, position);
301}
302
303/// Create an orthographic projection matrix.
304/// \param left The coordinate for the left vertical clipping plane.
305/// \param right The coordinate for the right vertical clipping plane.
306/// \param bottom The coordinate for the bottom horizontal clipping plane.
307/// \param top The coordinate for the top horizontal clipping plane.
308/// \param near The distance to the near clipping plane.
309/// \param far The distance to the far clipping plane.
310static inline mat4 mat4_ortho(R left, R right, R bottom, R top, R near, R far) {
311 const R tx = -(right + left) / (right - left);
312 const R ty = -(top + bottom) / (top - bottom);
313 const R tz = -(far + near) / (far - near);
314 return mat4_make(2 / (right - left), 0, 0, tx, 0, 2 / (top - bottom), 0, ty,
315 0, 0, -2 / (far - near), tz, 0, 0, 0, 1);
316}
317
318/// Create a perspective projection matrix.
319/// \param fovy The vertical field of view angle in degrees.
320/// \param aspect The aspect ratio that determines the field of view in the
321/// x-direction.
322/// \param near Distance to the near clipping plane.
323/// \param far Distance to the far clipping plane.
324static inline mat4 mat4_perspective(R fovy, R aspect, R near, R far) {
325 R f = tan(fovy / 2.0);
326 assert(f > 0.0);
327 f = 1.0 / f;
328 const R a = near - far;
329 return mat4_make(f / aspect, 0, 0, 0, 0, f, 0, 0, 0, 0, (far + near) / a,
330 (2 * far * near / a), 0, 0, -1, 0);
331}
332
333/// Create the inverse of a perspective projection matrix.
334/// \param fovy The vertical field of view angle in degrees.
335/// \param aspect The aspect ratio that determines the field of view in the
336/// x-direction.
337/// \param near Distance to the near clipping plane.
338/// \param far Distance to the far clipping plane.
339static inline mat4 mat4_perspective_inverse(R fovy, R aspect, R near, R far) {
340 R f = tan(fovy / 2.0);
341 assert(f > 0.0);
342 f = 1.0 / f;
343 const R a = far * near;
344 const R P32 = 0.5 * (near - far) / a;
345 const R P33 = 0.5 * (far + near) / a;
346 return mat4_make(aspect / f, 0, 0, 0, 0, 1.0f / f, 0, 0, 0, 0, 0, -1, 0, 0,
347 P32, P33);
348}
349
350/// Return the matrix's determinant.
351static inline R mat4_det(mat4 m) {
352 const R* M = (const R*)(m.val);
353 const R inv0 = M[5] * M[10] * M[15] - M[5] * M[11] * M[14] -
354 M[9] * M[6] * M[15] + M[9] * M[7] * M[14] +
355 M[13] * M[6] * M[11] - M[13] * M[7] * M[10];
356 const R inv1 = -M[4] * M[10] * M[15] + M[4] * M[11] * M[14] +
357 M[8] * M[6] * M[15] - M[8] * M[7] * M[14] -
358 M[12] * M[6] * M[11] + M[12] * M[7] * M[10];
359 const R inv2 = M[4] * M[9] * M[15] - M[4] * M[11] * M[13] -
360 M[8] * M[5] * M[15] + M[8] * M[7] * M[13] +
361 M[12] * M[5] * M[11] - M[12] * M[7] * M[9];
362 const R inv3 = -M[4] * M[9] * M[14] + M[4] * M[10] * M[13] +
363 M[8] * M[5] * M[14] - M[8] * M[6] * M[13] -
364 M[12] * M[5] * M[10] + M[12] * M[6] * M[9];
365 return M[0] * inv0 + M[1] * inv1 + M[2] * inv2 + M[3] * inv3;
366}
367
368/// Invert the matrix.
369static inline mat4 mat4_inverse(mat4 m) {
370 const R* M = (const R*)(m.val);
371
372 R inv[16];
373 inv[0] = M[5] * M[10] * M[15] - M[5] * M[11] * M[14] - M[9] * M[6] * M[15] +
374 M[9] * M[7] * M[14] + M[13] * M[6] * M[11] - M[13] * M[7] * M[10];
375 inv[4] = -M[4] * M[10] * M[15] + M[4] * M[11] * M[14] + M[8] * M[6] * M[15] -
376 M[8] * M[7] * M[14] - M[12] * M[6] * M[11] + M[12] * M[7] * M[10];
377 inv[8] = M[4] * M[9] * M[15] - M[4] * M[11] * M[13] - M[8] * M[5] * M[15] +
378 M[8] * M[7] * M[13] + M[12] * M[5] * M[11] - M[12] * M[7] * M[9];
379 inv[12] = -M[4] * M[9] * M[14] + M[4] * M[10] * M[13] + M[8] * M[5] * M[14] -
380 M[8] * M[6] * M[13] - M[12] * M[5] * M[10] + M[12] * M[6] * M[9];
381 inv[1] = -M[1] * M[10] * M[15] + M[1] * M[11] * M[14] + M[9] * M[2] * M[15] -
382 M[9] * M[3] * M[14] - M[13] * M[2] * M[11] + M[13] * M[3] * M[10];
383 inv[5] = M[0] * M[10] * M[15] - M[0] * M[11] * M[14] - M[8] * M[2] * M[15] +
384 M[8] * M[3] * M[14] + M[12] * M[2] * M[11] - M[12] * M[3] * M[10];
385 inv[9] = -M[0] * M[9] * M[15] + M[0] * M[11] * M[13] + M[8] * M[1] * M[15] -
386 M[8] * M[3] * M[13] - M[12] * M[1] * M[11] + M[12] * M[3] * M[9];
387 inv[13] = M[0] * M[9] * M[14] - M[0] * M[10] * M[13] - M[8] * M[1] * M[14] +
388 M[8] * M[2] * M[13] + M[12] * M[1] * M[10] - M[12] * M[2] * M[9];
389 inv[2] = M[1] * M[6] * M[15] - M[1] * M[7] * M[14] - M[5] * M[2] * M[15] +
390 M[5] * M[3] * M[14] + M[13] * M[2] * M[7] - M[13] * M[3] * M[6];
391 inv[6] = -M[0] * M[6] * M[15] + M[0] * M[7] * M[14] + M[4] * M[2] * M[15] -
392 M[4] * M[3] * M[14] - M[12] * M[2] * M[7] + M[12] * M[3] * M[6];
393 inv[10] = M[0] * M[5] * M[15] - M[0] * M[7] * M[13] - M[4] * M[1] * M[15] +
394 M[4] * M[3] * M[13] + M[12] * M[1] * M[7] - M[12] * M[3] * M[5];
395 inv[14] = -M[0] * M[5] * M[14] + M[0] * M[6] * M[13] + M[4] * M[1] * M[14] -
396 M[4] * M[2] * M[13] - M[12] * M[1] * M[6] + M[12] * M[2] * M[5];
397 inv[3] = -M[1] * M[6] * M[11] + M[1] * M[7] * M[10] + M[5] * M[2] * M[11] -
398 M[5] * M[3] * M[10] - M[9] * M[2] * M[7] + M[9] * M[3] * M[6];
399 inv[7] = M[0] * M[6] * M[11] - M[0] * M[7] * M[10] - M[4] * M[2] * M[11] +
400 M[4] * M[3] * M[10] + M[8] * M[2] * M[7] - M[8] * M[3] * M[6];
401 inv[11] = -M[0] * M[5] * M[11] + M[0] * M[7] * M[9] + M[4] * M[1] * M[11] -
402 M[4] * M[3] * M[9] - M[8] * M[1] * M[7] + M[8] * M[3] * M[5];
403 inv[15] = M[0] * M[5] * M[10] - M[0] * M[6] * M[9] - M[4] * M[1] * M[10] +
404 M[4] * M[2] * M[9] + M[8] * M[1] * M[6] - M[8] * M[2] * M[5];
405
406 R det = M[0] * inv[0] + M[1] * inv[4] + M[2] * inv[8] + M[3] * inv[12];
407 assert(det != 0.0);
408 det = 1.0 / det;
409 return mat4_make(inv[0] * det, inv[4] * det, inv[8] * det, inv[12] * det,
410 inv[1] * det, inv[5] * det, inv[9] * det, inv[13] * det,
411 inv[2] * det, inv[6] * det, inv[10] * det, inv[14] * det,
412 inv[3] * det, inv[7] * det, inv[11] * det, inv[15] * det);
413}
414
415/// Invert the transformation matrix.
416/// This is much faster than the more general inverse() function, but assumes
417/// that the matrix is of the form TR, where T is a translation and R a
418/// rotation.
419static inline mat4 mat4_inverse_transform(mat4 m) {
420 const vec3 r = mat4_v0(m);
421 const vec3 u = mat4_v1(m);
422 const vec3 f = mat4_v2(m);
423 const vec3 t = mat4_v3(m);
424 return mat4_make(r.x, r.y, r.z, -vec3_dot(r, t), u.x, u.y, u.z,
425 -vec3_dot(u, t), f.x, f.y, f.z, -vec3_dot(f, t), 0.0, 0.0,
426 0.0, 1.0);
427}
428
429/// Transpose the matrix.
430static inline mat4 mat4_transpose(mat4 m) {
431 return mat4_make(
432 mat4_at(m, 0, 0), mat4_at(m, 1, 0), mat4_at(m, 2, 0), mat4_at(m, 3, 0),
433 mat4_at(m, 0, 1), mat4_at(m, 1, 1), mat4_at(m, 2, 1), mat4_at(m, 3, 1),
434 mat4_at(m, 0, 2), mat4_at(m, 1, 2), mat4_at(m, 2, 2), mat4_at(m, 3, 2),
435 mat4_at(m, 0, 3), mat4_at(m, 1, 3), mat4_at(m, 2, 3), mat4_at(m, 3, 3));
436}
437
438/// Transform the vector with the matrix.
439static inline vec3 mat4_mul_vec3(mat4 m, vec3 v, R w) {
440 vec3 u;
441 u.x = mat4_at(m, 0, 0) * v.x + mat4_at(m, 0, 1) * v.y +
442 mat4_at(m, 0, 2) * v.z + mat4_at(m, 0, 3) * w;
443 u.y = mat4_at(m, 1, 0) * v.x + mat4_at(m, 1, 1) * v.y +
444 mat4_at(m, 1, 2) * v.z + mat4_at(m, 1, 3) * w;
445 u.z = mat4_at(m, 2, 0) * v.x + mat4_at(m, 2, 1) * v.y +
446 mat4_at(m, 2, 2) * v.z + mat4_at(m, 2, 3) * w;
447 return u;
448}
449
450/// Return the vector multiplied by the matrix.
451static inline vec4 mat4_mul_vec4(mat4 m, vec4 v) {
452 vec4 u;
453 u.x = mat4_at(m, 0, 0) * v.x + mat4_at(m, 0, 1) * v.y +
454 mat4_at(m, 0, 2) * v.z + mat4_at(m, 0, 3) * v.w;
455 u.y = mat4_at(m, 1, 0) * v.x + mat4_at(m, 1, 1) * v.y +
456 mat4_at(m, 1, 2) * v.z + mat4_at(m, 1, 3) * v.w;
457 u.z = mat4_at(m, 2, 0) * v.x + mat4_at(m, 2, 1) * v.y +
458 mat4_at(m, 2, 2) * v.z + mat4_at(m, 2, 3) * v.w;
459 u.w = mat4_at(m, 3, 0) * v.x + mat4_at(m, 3, 1) * v.y +
460 mat4_at(m, 3, 2) * v.z + mat4_at(m, 3, 3) * v.w;
461 return u;
462}
463
464/// Compare two matrices for equality.
465/// Returns true if the difference between each ij-value across matrices is
466/// within |eps|, false if there is at least one ij-value difference that is
467/// greater than eps.
468static inline bool mat4_eq(mat4 m, mat4 w, float eps) {
469 return (float_eq(mat4_at(m, 0, 0), mat4_at(w, 0, 0), eps) &&
470 float_eq(mat4_at(m, 0, 1), mat4_at(w, 0, 1), eps) &&
471 float_eq(mat4_at(m, 0, 2), mat4_at(w, 0, 2), eps) &&
472 float_eq(mat4_at(m, 0, 3), mat4_at(w, 0, 3), eps) &&
473
474 float_eq(mat4_at(m, 1, 0), mat4_at(w, 1, 0), eps) &&
475 float_eq(mat4_at(m, 1, 1), mat4_at(w, 1, 1), eps) &&
476 float_eq(mat4_at(m, 1, 2), mat4_at(w, 1, 2), eps) &&
477 float_eq(mat4_at(m, 1, 3), mat4_at(w, 1, 3), eps) &&
478
479 float_eq(mat4_at(m, 2, 0), mat4_at(w, 2, 0), eps) &&
480 float_eq(mat4_at(m, 2, 1), mat4_at(w, 2, 1), eps) &&
481 float_eq(mat4_at(m, 2, 2), mat4_at(w, 2, 2), eps) &&
482 float_eq(mat4_at(m, 2, 3), mat4_at(w, 2, 3), eps) &&
483
484 float_eq(mat4_at(m, 3, 0), mat4_at(w, 3, 0), eps) &&
485 float_eq(mat4_at(m, 3, 1), mat4_at(w, 3, 1), eps) &&
486 float_eq(mat4_at(m, 3, 2), mat4_at(w, 3, 2), eps) &&
487 float_eq(mat4_at(m, 3, 3), mat4_at(w, 3, 3), eps));
488}
diff --git a/include/math/spatial3.h b/include/math/spatial3.h
new file mode 100644
index 0000000..f8caf5d
--- /dev/null
+++ b/include/math/spatial3.h
@@ -0,0 +1,194 @@
1#pragma once
2
3#include "mat4.h"
4#include "vec3.h"
5
6/// An object in 3D space.
7typedef struct Spatial3 {
8 vec3 p; // Position.
9 vec3 r; // Right vector.
10 vec3 u; // Up vector.
11 vec3 f; // Forward vector.
12} Spatial3;
13
14/// Construct a spatial with position 0 and canonical direction vectors.
15static inline Spatial3 spatial3_make() {
16 return (Spatial3){.p = zero3(), .r = right3(), .u = up3(), .f = forward3()};
17}
18
19/// Return the spatial's transformation matrix (from local to world
20/// coordinates).
21static inline mat4 spatial3_transform(const Spatial3* spatial) {
22 const vec3 p = spatial->p;
23 const vec3 r = spatial->r;
24 const vec3 u = spatial->u;
25 const vec3 f = spatial->f;
26 return mat4_make(r.x, u.x, -f.x, p.x, r.y, u.y, -f.y, p.y, r.z, u.z, -f.z,
27 p.z, 0.0, 0.0, 0.0, 1.0f);
28}
29
30/// Return the spatial's inverse transformation matrix (from world to local
31/// coordinates).
32static inline mat4 spatial3_inverse_transform(const Spatial3* spatial) {
33 return mat4_inverse_transform(spatial3_transform(spatial));
34}
35
36/// Move the spatial by the given vector.
37static inline void spatial3_move(Spatial3* spatial, vec3 v) {
38 spatial->p = vec3_add(spatial->p, v);
39}
40
41/// Move the spatial along its forward vector.
42static inline void spatial3_move_forwards(Spatial3* spatial, R speed) {
43 spatial->p = vec3_add(spatial->p, vec3_scale(spatial->f, speed));
44}
45
46/// Move the spatial along its backwards vector.
47static inline void spatial3_move_backwards(Spatial3* spatial, R speed) {
48 spatial->p = vec3_add(spatial->p, vec3_scale(spatial->f, -speed));
49}
50
51/// Move the spatial along its left vector.
52static inline void spatial3_move_left(Spatial3* spatial, R speed) {
53 spatial->p = vec3_add(spatial->p, vec3_scale(spatial->r, -speed));
54}
55
56/// Move the spatial along its right vector.
57static inline void spatial3_move_right(Spatial3* spatial, R speed) {
58 spatial->p = vec3_add(spatial->p, vec3_scale(spatial->r, speed));
59}
60
61/// Move the spatial along the global up vector.
62static inline void spatial3_move_up(Spatial3* spatial, R speed) {
63 spatial->p = vec3_add(spatial->p, vec3_scale(up3(), speed));
64}
65
66/// Move the spatial along the global down vector.
67static inline void spatial3_move_down(Spatial3* spatial, R speed) {
68 spatial->p = vec3_add(spatial->p, vec3_scale(up3(), -speed));
69}
70
71/// Rotate the spatial about the given axis by the given angle.
72static inline void spatial3_rotate(Spatial3* spatial, vec3 axis, R angle) {
73 mat4 transf = spatial3_transform(spatial);
74 const vec3 local_axis =
75 vec3_normalize(mat4_mul_vec3(mat4_inverse_transform(transf), axis, 0.0));
76 transf = mat4_mul(transf, mat4_rot(local_axis, angle));
77 spatial->r = vec3_normalize(mat4_v0(transf));
78 spatial->u = vec3_normalize(mat4_v1(transf));
79 spatial->f = vec3_normalize(vec3_neg(mat4_v2(transf)));
80}
81
82/// Rotate the spatial about its local y axis.
83static inline void spatial3_yaw(Spatial3* spatial, const R angle) {
84 const R sa = sin(angle);
85 const R ca = cos(angle);
86 spatial->f = vec3_normalize(
87 vec3_sub(vec3_scale(spatial->f, ca), vec3_scale(spatial->r, sa)));
88 spatial->r = vec3_normalize(vec3_cross(spatial->f, spatial->u));
89}
90
91/// Rotate the spatial about its local x axis.
92static inline void spatial3_pitch(Spatial3* spatial, const R angle) {
93 const R sa = sin(angle);
94 const R ca = cos(angle);
95 spatial->f = vec3_normalize(
96 vec3_add(vec3_scale(spatial->f, ca), vec3_scale(spatial->u, sa)));
97 spatial->u = vec3_normalize(vec3_cross(spatial->r, spatial->f));
98}
99
100/// Rotate the spatial about its local z axis.
101static inline void spatial3_roll(Spatial3* spatial, const R angle) {
102 const R sa = sin(angle);
103 const R ca = cos(angle);
104 spatial->u = vec3_normalize(
105 vec3_sub(vec3_scale(spatial->u, ca), vec3_scale(spatial->r, sa)));
106 spatial->r = vec3_normalize(vec3_cross(spatial->f, spatial->u));
107}
108
109/// Set the spatial's transformation matrix.
110static inline void spatial3_set_transform(Spatial3* spatial, mat4 transform) {
111 spatial->r = mat4_v0(transform);
112 spatial->u = mat4_v1(transform);
113 spatial->f = mat4_v2(transform);
114 spatial->p = mat4_v3(transform);
115}
116
117static inline void spatial3_set_forward(Spatial3* spatial, vec3 forward) {
118 spatial->f = vec3_normalize(forward);
119 // Use aux vector to define right vector orthogonal to forward.
120 if (vec3_eq(vec3_abs(spatial->f), up3())) {
121 spatial->r = vec3_normalize(vec3_cross(spatial->f, forward3()));
122 } else {
123 spatial->r = vec3_normalize(vec3_cross(spatial->f, up3()));
124 }
125 spatial->u = vec3_normalize(vec3_cross(spatial->r, spatial->f));
126}
127
128/// Make the spatial look at the given target.
129static inline void spatial3_lookat(Spatial3* spatial, vec3 target) {
130 spatial3_set_forward(spatial, vec3_sub(target, spatial->p));
131}
132
133/// Make the spatial look at the given target.
134static inline void spatial3_lookat_spatial(Spatial3* spatial,
135 const Spatial3* target) {
136 spatial3_set_forward(spatial, vec3_sub(target->p, spatial->p));
137}
138
139/// Make the spatial orbit around the given target.
140/// \param target Target position.
141/// \param radius Radial distance.
142/// \param azimuth Azimuthal (horizontal) angle.
143/// \param zenith Polar (vertical) angle.
144static inline void spatial3_orbit(Spatial3* spatial, vec3 target, R radius,
145 R azimuth, R zenith) {
146 const R sx = sin(azimuth);
147 const R sy = sin(zenith);
148 const R cx = cos(azimuth);
149 const R cy = cos(zenith);
150 spatial->p = (vec3){target.x + radius * cy * sx, target.y + radius * sy,
151 target.z + radius * cx * cy};
152}
153
154/// Make the spatial orbit around the given target.
155/// \param target Target spatial.
156/// \param radius Radial distance.
157/// \param azimuth Azimuthal (horizontal) angle.
158/// \param zenith Polar (vertical) angle.
159static inline void spatial3_orbit_spatial(Spatial3* spatial,
160 const Spatial3* target, R radius,
161 R azimuth, R zenith) {
162 spatial3_orbit(spatial, target->p, radius, azimuth, zenith);
163}
164
165// The functions below are provided so that client code can work with a Spatial3
166// with no assumptions on the data type's definition.
167
168/// Return the spatial's position.
169static inline vec3 spatial3_position(const Spatial3* spatial) {
170 return spatial->p;
171}
172
173/// Return the spatial's right vector.
174static inline vec3 spatial3_right(const Spatial3* spatial) {
175 return spatial->r;
176}
177
178/// Return the spatial's up vector.
179static inline vec3 spatial3_up(const Spatial3* spatial) { return spatial->u; }
180
181/// Return the spatial's forward vector.
182static inline vec3 spatial3_forward(const Spatial3* spatial) {
183 return spatial->f;
184}
185
186static inline void spatial3_set_position(Spatial3* spatial, vec3 p) {
187 spatial->p = p;
188}
189
190static inline void spatial3_setx(Spatial3* spatial, R x) { spatial->p.x = x; }
191
192static inline void spatial3_sety(Spatial3* spatial, R y) { spatial->p.y = y; }
193
194static inline void spatial3_setz(Spatial3* spatial, R z) { spatial->p.z = z; }
diff --git a/include/math/vec.h b/include/math/vec.h
new file mode 100644
index 0000000..7b870e1
--- /dev/null
+++ b/include/math/vec.h
@@ -0,0 +1,16 @@
1// Common functions for vectors.
2//
3// This header file contains functions that operate with different kinds of
4// vectors that would result in circular dependencies if defined in any of the
5// vec2/3/4 header files.
6#pragma once
7
8#include "vec2.h"
9#include "vec3.h"
10#include "vec4.h"
11
12/// Construct a 3D vector from a 2D one with z = 0.
13static inline vec3 vec3_from_vec2(vec2 v) { return vec3{v.x, v.y, 0, 0}; }
14
15/// Project a 4D vector onto w=0.
16static inline vec3 vec3_from_vec4(vec4 v) { return vec3{v.x, v.y, v.z, 0}; }
diff --git a/include/math/vec2.h b/include/math/vec2.h
new file mode 100644
index 0000000..0a3f692
--- /dev/null
+++ b/include/math/vec2.h
@@ -0,0 +1,10 @@
1#pragma once
2
3#include "defs.h"
4
5typedef struct vec2 {
6 R x, y;
7} vec2;
8
9/// Construct a vector from 2 coordinates.
10static inline vec2 vec2_make(R x, R y) { return (vec2){x, y}; }
diff --git a/include/math/vec3.h b/include/math/vec3.h
new file mode 100644
index 0000000..3c3b053
--- /dev/null
+++ b/include/math/vec3.h
@@ -0,0 +1,143 @@
1#pragma once
2
3#include "defs.h"
4
5#include <assert.h>
6#include <stdbool.h>
7
8/// A 3D vector.
9typedef struct vec3 {
10 R x, y, z;
11} vec3;
12
13/// Construct a vector from 3 coordinates.
14static inline vec3 vec3_make(R x, R y, R z) { return (vec3){x, y, z}; }
15
16/// Construct a vector from an array.
17static inline vec3 vec3_from_array(const R xyz[3]) {
18 return (vec3){xyz[0], xyz[1], xyz[2]};
19}
20
21/// Construct a vector from a single scalar value.
22/// x = y = z = val.
23static inline vec3 vec3_from_scalar(R val) { return (vec3){val, val, val}; }
24
25/// Normalize the vector.
26static inline vec3 vec3_normalize(vec3 v) {
27 R n = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
28 assert(n > 0);
29 return (vec3){v.x / n, v.y / n, v.z / n};
30}
31
32/// Return the vector's ith coordinate.
33static inline R vec3_ith(vec3 v, int i) {
34 assert(i >= 0 && i < 3);
35 return ((const R*)&v)[i];
36}
37
38/// Negate the given vector.
39static inline vec3 vec3_neg(vec3 v) { return (vec3){-v.x, -v.y, -v.z}; }
40
41/// Add two vectors.
42static inline vec3 vec3_add(vec3 a, vec3 b) {
43 return (vec3){a.x + b.x, a.y + b.y, a.z + b.z};
44}
45
46/// Subtract two vectors.
47static inline vec3 vec3_sub(vec3 a, vec3 b) {
48 return (vec3){a.x - b.x, a.y - b.y, a.z - b.z};
49}
50
51/// Modulate two vectors (component-wise multiplication).
52static inline vec3 vec3_mul(vec3 a, vec3 b) {
53 return (vec3){a.x * b.x, a.y * b.y, a.z * b.z};
54}
55
56/// Divide two vectors component-wise.
57static inline vec3 vec3_div(vec3 a, vec3 b) {
58 return (vec3){a.x / b.x, a.y / b.y, a.z / b.z};
59}
60
61/// Scale a vector by a scalar value.
62static inline vec3 vec3_scale(vec3 v, R s) {
63 return (vec3){v.x * s, v.y * s, v.z * s};
64}
65
66/// Compare two vectors for equality.
67static inline bool vec3_eq(vec3 a, vec3 b) {
68 return a.x == b.x && a.y == b.y && a.z == b.z;
69}
70
71/// Return the absolute value of the vector.
72static inline vec3 vec3_abs(vec3 v) {
73 return (vec3){rabs(v.x), rabs(v.y), rabs(v.z)};
74}
75
76/// Compare two vectors for inequality.
77static inline bool vec3_ne(vec3 a, vec3 b) { return !(vec3_eq(a, b)); }
78
79/// Return the vector's squared magnitude.
80static inline R vec3_norm2(vec3 v) { return v.x * v.x + v.y * v.y + v.z * v.z; }
81
82/// Return the vector's magnitude.
83static inline R vec3_norm(vec3 v) { return sqrt(vec3_norm2(v)); }
84
85/// Return the squared distance between two points.
86static inline R vec3_dist2(vec3 a, vec3 b) {
87 const vec3 v = vec3_sub(b, a);
88 return vec3_norm2(v);
89}
90
91/// Return the distance between two points.
92static inline R vec3_dist(vec3 a, vec3 b) { return sqrt(vec3_dist2(a, b)); }
93
94/// Return the given vector divided by its magnitude.
95static inline vec3 normalize(vec3 v) {
96 const R n = vec3_norm(v);
97 assert(n > 0);
98 return (vec3){v.x / n, v.y / n, v.z / n};
99}
100
101/// Return the dot product of two vectors.
102static inline R vec3_dot(vec3 a, vec3 b) {
103 return a.x * b.x + a.y * b.y + a.z * b.z;
104}
105
106/// Return the cross product of two vectors.
107static inline vec3 vec3_cross(vec3 a, vec3 b) {
108 return (vec3){a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z,
109 a.x * b.y - a.y * b.x};
110}
111
112/// Reflect the vector about the normal.
113static inline vec3 vec3_reflect(vec3 v, vec3 n) {
114 // r = v - 2 * dot(v, n) * n
115 return vec3_sub(v, vec3_scale(n, 2 * vec3_dot(v, n)));
116}
117
118/// Refract the vector about the normal.
119static inline vec3 refract(vec3 v, vec3 n, R e) {
120 // k = 1 - e^2(1 - dot(n,v) * dot(n,v))
121 const R k = 1.0 - e * e * (1.0 - vec3_dot(n, v) * vec3_dot(n, v));
122 assert(k >= 0);
123 // r = e*v - (e * dot(n,v) + sqrt(k)) * n
124 return vec3_sub(vec3_scale(v, e),
125 vec3_scale(n, e * vec3_dot(n, v) * sqrt(k)));
126}
127
128/// Elevate the vector to a power.
129static inline vec3 vec3_pow(vec3 v, R p) {
130 return (vec3){pow(v.x, p), pow(v.y, p), pow(v.z, p)};
131}
132
133/// The (1, 0, 0) vector.
134static inline vec3 right3() { return (vec3){1.0, 0.0, 0.0}; }
135
136/// The (0, 1, 0) vector.
137static inline vec3 up3() { return (const vec3){0.0, 1.0, 0.0}; }
138
139/// The (0, 0, -1) vector.
140static inline vec3 forward3() { return (const vec3){0.0, 0.0, -1.0}; }
141
142/// The (0, 0, 0) vector.
143static inline vec3 zero3() { return (const vec3){0.0, 0.0, 0.0}; }
diff --git a/include/math/vec4.h b/include/math/vec4.h
new file mode 100644
index 0000000..4ab843b
--- /dev/null
+++ b/include/math/vec4.h
@@ -0,0 +1,15 @@
1#pragma once
2
3#include "defs.h"
4
5typedef struct vec4 {
6 R x, y, w, z;
7} vec4;
8
9/// Construct a vector from 4 coordinates.
10static inline vec4 vec4_make(R x, R y, R z, R w) { return (vec4){x, y, z, w}; }
11
12/// Construct a vector from an array.
13static inline vec4 vec4_from_array(const R xyzw[4]) {
14 return (vec4){xyzw[0], xyzw[1], xyzw[2], xyzw[3]};
15}
diff --git a/test/mat4_test.c b/test/mat4_test.c
new file mode 100644
index 0000000..d4f61bc
--- /dev/null
+++ b/test/mat4_test.c
@@ -0,0 +1,111 @@
1#include <math/mat4.h>
2
3#include <math/float.h>
4
5#include "test.h"
6
7#include <stdio.h>
8
9static const float eps = 1e-9;
10
11static inline void print_mat4(mat4 m) {
12 printf("\n");
13 for (int i = 0; i < 4; ++i) {
14 for (int j = 0; j < 4; ++j) {
15 printf("%6.3f ", mat4_at(m, i, j));
16 }
17 printf("\n");
18 }
19}
20
21TEST_CASE(lookat_right) {
22 const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
23 /*target=*/vec3_make(1, 0, 0),
24 /*up=*/vec3_make(0, 1, 0));
25
26 const mat4 expected = mat4_from_array((const float*)(float[4][4]){
27 {0, 0, 1, 0},
28 {0, 1, 0, 0},
29 {1, 0, 0, 0},
30 {0, 0, 0, 1},
31 });
32
33 TEST_TRUE(mat4_eq(m, expected, eps));
34}
35
36TEST_CASE(lookat_left) {
37 const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
38 /*target=*/vec3_make(-1, 0, 0),
39 /*up=*/vec3_make(0, 1, 0));
40
41 const mat4 expected = mat4_from_array((const float*)(float[4][4]){
42 {0, 0, -1, 0},
43 {0, 1, 0, 0},
44 {-1, 0, 0, 0},
45 {0, 0, 0, 1},
46 });
47
48 TEST_TRUE(mat4_eq(m, expected, eps));
49}
50
51TEST_CASE(lookat_up) {
52 const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
53 /*target=*/vec3_make(0, 1, 0),
54 /*up=*/vec3_make(0, 0, 1));
55
56 const mat4 expected = mat4_from_array((const float*)(float[4][4]){
57 {1, 0, 0, 0},
58 {0, 0, 1, 0},
59 {0, 1, 0, 0},
60 {0, 0, 0, 1},
61 });
62
63 TEST_TRUE(mat4_eq(m, expected, eps));
64}
65
66TEST_CASE(lookat_down) {
67 const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
68 /*target=*/vec3_make(0, -1, 0),
69 /*up=*/vec3_make(0, 0, -1));
70
71 const mat4 expected = mat4_from_array((const float*)(float[4][4]){
72 {1, 0, 0, 0},
73 {0, 0, -1, 0},
74 {0, -1, 0, 0},
75 {0, 0, 0, 1},
76 });
77
78 TEST_TRUE(mat4_eq(m, expected, eps));
79}
80
81TEST_CASE(lookat_forward) {
82 const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
83 /*target=*/vec3_make(0, 0, -1),
84 /*up=*/vec3_make(0, 1, 0));
85
86 const mat4 expected = mat4_from_array((const float*)(float[4][4]){
87 {1, 0, 0, 0},
88 {0, 1, 0, 0},
89 {0, 0, -1, 0},
90 {0, 0, 0, 1},
91 });
92
93 TEST_TRUE(mat4_eq(m, expected, eps));
94}
95
96TEST_CASE(lookat_back) {
97 const mat4 m = mat4_lookat(/*position=*/vec3_make(0, 0, 0),
98 /*target=*/vec3_make(0, 0, 1),
99 /*up=*/vec3_make(0, 1, 0));
100
101 const mat4 expected = mat4_from_array((const float*)(float[4][4]){
102 {-1, 0, 0, 0},
103 {0, 1, 0, 0},
104 {0, 0, 1, 0},
105 {0, 0, 0, 1},
106 });
107
108 TEST_TRUE(mat4_eq(m, expected, eps));
109}
110
111int main() { return 0; }
diff --git a/test/test.h b/test/test.h
new file mode 100644
index 0000000..fd8dc22
--- /dev/null
+++ b/test/test.h
@@ -0,0 +1,185 @@
1// SPDX-License-Identifier: MIT
2#pragma once
3
4#ifdef UNIT_TEST
5
6#include <stdbool.h>
7#include <stdio.h>
8#include <stdlib.h>
9#include <string.h>
10
11#if defined(__DragonFly__) || defined(__FreeBSD__) || defined(__FreeBSD_kernel__) || \
12 defined(__NetBSD__) || defined(__OpenBSD__)
13#define USE_SYSCTL_FOR_ARGS 1
14// clang-format off
15#include <sys/types.h>
16#include <sys/sysctl.h>
17// clang-format on
18#include <unistd.h> // getpid
19#endif
20
21struct test_file_metadata;
22
23struct test_failure {
24 bool present;
25 const char *message;
26 const char *file;
27 int line;
28};
29
30struct test_case_metadata {
31 void (*fn)(struct test_case_metadata *, struct test_file_metadata *);
32 struct test_failure failure;
33 const char *name;
34 struct test_case_metadata *next;
35};
36
37struct test_file_metadata {
38 bool registered;
39 const char *name;
40 struct test_file_metadata *next;
41 struct test_case_metadata *tests;
42};
43
44struct test_file_metadata __attribute__((weak)) * test_file_head;
45
46#define SET_FAILURE(_message) \
47 metadata->failure = (struct test_failure) { \
48 .message = _message, .file = __FILE__, .line = __LINE__, .present = true, \
49 }
50
51#define TEST_EQUAL(a, b) \
52 do { \
53 if ((a) != (b)) { \
54 SET_FAILURE(#a " != " #b); \
55 return; \
56 } \
57 } while (0)
58
59#define TEST_TRUE(a) \
60 do { \
61 if (!(a)) { \
62 SET_FAILURE(#a " is not true"); \
63 return; \
64 } \
65 } while (0)
66
67#define TEST_STREQUAL(a, b) \
68 do { \
69 if (strcmp(a, b) != 0) { \
70 SET_FAILURE(#a " != " #b); \
71 return; \
72 } \
73 } while (0)
74
75#define TEST_CASE(_name) \
76 static void __test_h_##_name(struct test_case_metadata *, \
77 struct test_file_metadata *); \
78 static struct test_file_metadata __test_h_file; \
79 static struct test_case_metadata __test_h_meta_##_name = { \
80 .name = #_name, \
81 .fn = __test_h_##_name, \
82 }; \
83 static void __attribute__((constructor(101))) __test_h_##_name##_register(void) { \
84 __test_h_meta_##_name.next = __test_h_file.tests; \
85 __test_h_file.tests = &__test_h_meta_##_name; \
86 if (!__test_h_file.registered) { \
87 __test_h_file.name = __FILE__; \
88 __test_h_file.next = test_file_head; \
89 test_file_head = &__test_h_file; \
90 __test_h_file.registered = true; \
91 } \
92 } \
93 static void __test_h_##_name( \
94 struct test_case_metadata *metadata __attribute__((unused)), \
95 struct test_file_metadata *file_metadata __attribute__((unused)))
96
97extern void __attribute__((weak)) (*test_h_unittest_setup)(void);
98/// Run defined tests, return true if all tests succeeds
99/// @param[out] tests_run if not NULL, set to whether tests were run
100static inline void __attribute__((constructor(102))) run_tests(void) {
101 bool should_run = false;
102#ifdef USE_SYSCTL_FOR_ARGS
103 int mib[] = {
104 CTL_KERN,
105#if defined(__NetBSD__) || defined(__OpenBSD__)
106 KERN_PROC_ARGS,
107 getpid(),
108 KERN_PROC_ARGV,
109#else
110 KERN_PROC,
111 KERN_PROC_ARGS,
112 getpid(),
113#endif
114 };
115 char *arg = NULL;
116 size_t arglen;
117 sysctl(mib, sizeof(mib) / sizeof(mib[0]), NULL, &arglen, NULL, 0);
118 arg = malloc(arglen);
119 sysctl(mib, sizeof(mib) / sizeof(mib[0]), arg, &arglen, NULL, 0);
120#else
121 FILE *cmdlinef = fopen("/proc/self/cmdline", "r");
122 char *arg = NULL;
123 int arglen;
124 fscanf(cmdlinef, "%ms%n", &arg, &arglen);
125 fclose(cmdlinef);
126#endif
127 for (char *pos = arg; pos < arg + arglen; pos += strlen(pos) + 1) {
128 if (strcmp(pos, "--unittest") == 0) {
129 should_run = true;
130 break;
131 }
132 }
133 free(arg);
134
135 if (!should_run) {
136 return;
137 }
138
139 if (&test_h_unittest_setup) {
140 test_h_unittest_setup();
141 }
142
143 struct test_file_metadata *i = test_file_head;
144 int failed = 0, success = 0;
145 while (i) {
146 fprintf(stderr, "Running tests from %s:\n", i->name);
147 struct test_case_metadata *j = i->tests;
148 while (j) {
149 fprintf(stderr, "\t%s ... ", j->name);
150 j->failure.present = false;
151 j->fn(j, i);
152 if (j->failure.present) {
153 fprintf(stderr, "failed (%s at %s:%d)\n", j->failure.message,
154 j->failure.file, j->failure.line);
155 failed++;
156 } else {
157 fprintf(stderr, "passed\n");
158 success++;
159 }
160 j = j->next;
161 }
162 fprintf(stderr, "\n");
163 i = i->next;
164 }
165 int total = failed + success;
166 fprintf(stderr, "Test results: passed %d/%d, failed %d/%d\n", success, total,
167 failed, total);
168 exit(failed == 0 ? EXIT_SUCCESS : EXIT_FAILURE);
169}
170
171#else
172
173#include <stdbool.h>
174
175#define TEST_CASE(name) static void __attribute__((unused)) __test_h_##name(void)
176
177#define TEST_EQUAL(a, b) \
178 (void)(a); \
179 (void)(b)
180#define TEST_TRUE(a) (void)(a)
181#define TEST_STREQUAL(a, b) \
182 (void)(a); \
183 (void)(b)
184
185#endif