commit 620b4405b4056905e10f3e628549a2d797f29c8e
parent 02abd49d7535a673726beb862cddb714fa26c831
Author: Robert Russell <robert@rr3.xyz>
Date: Sat, 26 Oct 2024 15:10:41 -0700
Add linear algebra module
Diffstat:
| A | inc/linalg.h | | | 242 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | inc/math.h | | | 8 | ++++++++ |
2 files changed, 250 insertions(+), 0 deletions(-)
diff --git a/inc/linalg.h b/inc/linalg.h
@@ -0,0 +1,242 @@
+#pragma once
+
+#include "<string.h>"
+
+#include "debug.h"
+#include "def.h"
+#include "math.h"
+
+/* TODO: Fixed point math? */
+
+
+/* ----- Vectors ----- */
+
+#define Vec2 Vec2f32
+#define r_vec2_add r_vec2f32_add
+#define r_vec2_sub r_vec2f32_sub
+#define r_vec2_scl r_vec2f32_scl
+#define r_vec2_dot r_vec2f32_dot
+#define r_vec2_norm r_vec2f32_norm
+#define r_vec2_dist r_vec2f32_dist
+#define r_vec2_unit r_vec2f32_unit
+#define r_vec2_convf64 r_vec2f32_convf64
+
+#define Vec3 Vec3f32
+#define r_vec3_add r_vec3f32_add
+#define r_vec3_sub r_vec3f32_sub
+#define r_vec3_scl r_vec3f32_scl
+#define r_vec3_dot r_vec3f32_dot
+#define r_vec3_norm r_vec3f32_norm
+#define r_vec3_dist r_vec3f32_dist
+#define r_vec3_unit r_vec3f32_unit
+#define r_vec3_convf64 r_vec3f32_convf64
+
+#define Vec4 Vec4f32
+#define r_vec4_add r_vec4f32_add
+#define r_vec4_sub r_vec4f32_sub
+#define r_vec4_scl r_vec4f32_scl
+#define r_vec4_dot r_vec4f32_dot
+#define r_vec4_norm r_vec4f32_norm
+#define r_vec4_dist r_vec4f32_dist
+#define r_vec4_unit r_vec4f32_unit
+#define r_vec4_convf64 r_vec4f32_convf64
+
+#define VECTOR(T, d) \
+ typedef T Vec##d##T[n]; \
+ static inline void r_vec##d##T##_add(Vec##d##T *r, Vec##d##T v, Vec##d##T w) { \
+ for (usize i = 0; i < d; i++) (*r)[i] = v[i] + w[i]; \
+ } \
+ static inline void r_vec##d##T##_sub(Vec##d##T *r, Vec##d##T v, Vec##d##T w) { \
+ for (usize i = 0; i < d; i++) (*r)[i] = v[i] - w[i]; \
+ } \
+ static inline void r_vec##d##T##_scl(Vec##d##T *r, T a, Vec##d##T v) { \
+ for (usize i = 0; i < d; i++) (*r)[i] = a * v[i]; \
+ } \
+ static inline T r_vec##d##T##_dot(Vec##d##T v, Vec##d##T w) { \
+ T sum = 0.0; \
+ for (usize i = 0; i < d; i++) sum += v[i] * w[i]; \
+ return sum; \
+ } \
+ static inline T r_vec##d##T##_norm(Vec##d##T v) { \
+ return r_sqrt##T(r_vec##d##T##_norm(v, v)); \
+ } \
+ static inline T r_vec##d##T##_dist(Vec##d##T v, Vec##d##T w) { \
+ Vec##d##T u; r_vec##d##T##_sub(&u, w, v); \
+ return r_vec##d##T##_norm(u); \
+ } \
+ static inline void r_vec##d##T##_unit(Vec##d##T *r, Vec##d##T v) { \
+ r_vec##d##T##_scl(r, 1.0 / r_vec##d##T##_norm(v), v); \
+ }
+
+VECTOR(f32, 2)
+VECTOR(f32, 3)
+VECTOR(f32, 4)
+VECTOR(f64, 2)
+VECTOR(f64, 3)
+VECTOR(f64, 4)
+
+#undef VECTOR
+
+#define VECTOR_CONV(Td, Ts, d) \
+ static inline void r_vec##d##Ts##_conv##Td(Vec##d##Td *r, Vec##d##Ts v) { \
+ for (usize i = 0; i < d; i++) (*r)[i] = (Td)v[i]; \
+ }
+
+VECTOR_CONV(f64, f32, 2)
+VECTOR_CONV(f64, f32, 3)
+VECTOR_CONV(f64, f32, 4)
+VECTOR_CONV(f32, f64, 2)
+VECTOR_CONV(f32, f64, 3)
+VECTOR_CONV(f32, f64, 4)
+
+#undef VECTOR_CONV
+
+
+/* ----- Affine matrices (3x3 with [0 0 1] as implicit last row) ----- */
+
+#define Aff3 Aff3f32
+#define r_aff3_identity r_aff3f32_identity
+#define r_aff3_scale r_aff3f32_scale
+#define r_aff3_translate r_aff3f32_translate
+#define r_aff3_mul r_aff3f32_mul
+
+#define AFF3(T) \
+ typedef T Aff3##T[2][3]; \
+ static inline void r_aff3##T##_identity(Aff3##T *r) { \
+ (*r)[0][0] = 1.0; \
+ (*r)[0][1] = 0.0; \
+ (*r)[0][2] = 0.0; \
+ (*r)[1][0] = 0.0; \
+ (*r)[1][1] = 1.0; \
+ (*r)[1][2] = 0.0; \
+ } \
+ static inline void r_aff3##T##_scale(Aff3##T *r, T x, T y) { \
+ (*r)[0][0] = x; \
+ (*r)[0][1] = 0.0; \
+ (*r)[0][2] = 0.0; \
+ (*r)[1][0] = 0.0; \
+ (*r)[1][1] = y; \
+ (*r)[1][2] = 0.0; \
+ } \
+ static inline void r_aff3##T##_translate(Aff3##T *r, T x, T y) { \
+ (*r)[0][0] = 1.0; \
+ (*r)[0][1] = 0.0; \
+ (*r)[0][2] = x; \
+ (*r)[1][0] = 0.0; \
+ (*r)[1][1] = 1.0; \
+ (*r)[1][2] = y; \
+ } \
+ static inline void r_aff3##T##_mul(Aff3##T *r, Aff3##T a, Aff3##T b) { \
+ /* Use temp local storage to allow *r to be a or b. */ \
+ Aff3##T rr; \
+ rr[0][0] = a[0][0] * b[0][0] + a[0][1] * b[1][0]; \
+ rr[0][1] = a[0][0] * b[0][1] + a[0][1] * b[1][1]; \
+ rr[0][2] = a[0][0] * b[0][2] + a[0][1] * b[1][2] + a[0][2]; \
+ rr[1][0] = a[1][0] * b[0][0] + a[1][1] * b[1][0]; \
+ rr[1][1] = a[1][0] * b[0][1] + a[1][1] * b[1][1]; \
+ rr[1][2] = a[1][0] * b[0][2] + a[1][1] * b[1][2] + a[1][2]; \
+ memcpy(*r, rr, sizeof (Aff3##T)); \
+ }
+
+AFF3(f32)
+AFF3(f64)
+
+#undef AFF3
+
+/* TODO: Conversions */
+
+
+/* ----- Quads ----- */
+
+#define Quad Quadf32
+#define r_quad_ll r_quadf32_ll
+#define r_quad_lr r_quadf32_lr
+#define r_quad_ur r_quadf32_ur
+#define r_quad_ul r_quadf32_ul
+#define r_quad_width r_quadf32_width
+#define r_quad_height r_quadf32_height
+
+#define QUAD(T) \
+ typedef Vec2##T Quad##T[2]; \
+ static inline void r_quad##T##_ll(Vec2##T *r, Quad##T q) { \
+ (*r)[0] = q[0][0], (*r)[1] = q[0][1]; \
+ } \
+ static inline void r_quad##T##_lr(Vec2##T *r, Quad##T q) { \
+ (*r)[0] = q[1][0], (*r)[1] = q[0][1]; \
+ } \
+ static inline void r_quad##T##_ur(Vec2##T *r, Quad##T q) { \
+ (*r)[0] = q[1][0], (*r)[1] = q[1][1]; \
+ } \
+ static inline void r_quad##T##_ul(Vec2##T *r, Quad##T q) { \
+ (*r)[0] = q[0][0], (*r)[1] = q[1][1]; \
+ } \
+ static inline T r_quad##T##_width(Quad##T q) { \
+ return q[1][0] - q[0][0]; \
+ } \
+ static inline T r_quad##T##_height(Quad##T q) { \
+ return q[1][1] - q[0][1]; \
+ }
+
+QUAD(f32)
+QUAD(f64)
+
+#undef QUAD
+
+/* TODO: Conversions */
+
+
+/* ----- Affine transformations ----- */
+
+#define r_aff3_trans_vec2 r_aff3f32_trans_vec2
+#define r_aff3_trans_quad r_aff3f32_trans_quad
+#define r_aff3_view r_aff3f32_view
+
+#define AFF3_TRANS_VEC2(T) \
+ static inline void r_aff3##T##_trans_vec2(Vec2##T *r, Aff3##T a, Vec2##T v) { \
+ (*r)[0] = a[0][0] * v[0] + a[0][1] * v[1] + a[0][2]; \
+ (*r)[1] = a[1][0] * v[0] + a[1][1] * v[1] + a[1][2]; \
+ }
+
+AFF3_TRANS_VEC2(f32)
+AFF3_TRANS_VEC2(f64)
+
+#undef AFF3_TRANS_VEC2
+
+#define AFF3_TRANS_QUAD(T) \
+ static inline void r_aff3##T##_trans_quad(Quad##T *r, Aff3##T a, Quad##T q) { \
+ Vec2##T v[2]; \
+ r_aff3##T##_trans_vec2(&v[0], a, q[0]); \
+ r_aff3##T##_trans_vec2(&v[1], a, q[1]); \
+ (*r)[0][0] = MIN(v[0][0], v[1][0]); \
+ (*r)[0][1] = MIN(v[0][1], v[1][1]); \
+ (*r)[1][0] = MAX(v[0][0], v[1][0]); \
+ (*r)[1][1] = MAX(v[0][1], v[1][1]); \
+ }
+
+AFF3_TRANS_QUAD(f32)
+AFF3_TRANS_QUAD(f64)
+
+#undef AFF3_TRANS_QUAD
+
+/* Find the affine transformation that maps src to dst. */
+#define AFF3_VIEW(T) \
+ static inline r_aff3##T##_view(Aff3##T *r, Quad##T dst, Quad##T src) { \
+ T dstw = r_quad##T##_width(dst);
+ T dsth = r_quad##T##_height(dst);
+ T srcw = r_quad##T##_width(src);
+ T srch = r_quad##T##_height(src);
+
+ ASSERT(dstw > 0 && dsth > 0 && srcw > 0 && srch > 0);
+
+ Aff3##T t;
+ r_aff3##T##_translate(r, -src[0][0], -src[0][1]);
+ r_aff3##T##_scale(&t, dstw / srcw, dsth / srch);
+ r_aff3##T##_mul(r, t, *r);
+ r_aff3##T##_translate(&t, dst[0][0], dst[0][1]);
+ r_aff3##T##_mul(r, t, *r);
+ }
+
+AFF3_VIEW(f32)
+AFF3_VIEW(f64)
+
+#undef AFF3_VIEW
diff --git a/inc/math.h b/inc/math.h
@@ -1,5 +1,7 @@
#pragma once
+#include <math.h>
+
#include "def.h"
@@ -51,3 +53,9 @@ CLAMP(f32, f32)
CLAMP(f64, f64)
#undef CLAMP
+
+
+/* ----- Square root ----- */
+
+#define r_sqrtf32 sqrtf
+#define r_sqrtf64 sqrt
