fruit/libguf
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add guf_nearly_zero/one

Browse Source
This commit is contained in:
jun 2025-03-07 11:09:18 +01:00
parent 7b433cd776
commit 9096a4f9c9
4 changed files with 100 additions and 95 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists.txt
View File

@ -18,10 +18,10 @@ if (NOT DEFINED CMAKE_RUNTIME_OUTPUT_DIRECTORY)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin) set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin)
endif () endif ()
add_executable(libguf_example src/test/example.c src/test/guf_str_impl.c src/test/guf_dict_impl.c) add_executable(libguf_example src/test/example.c src/test/guf_str_impl.c src/test/guf_dict_impl.c src/test/guf_linalg_impl.c)
target_include_directories(libguf_example PRIVATE src src/test) target_include_directories(libguf_example PRIVATE src src/test)
add_executable(libguf_test src/test/test.cpp src/test/guf_init_impl.c src/test/guf_dbuf_impl.c src/test/guf_str_impl.c src/test/guf_dict_impl.c src/test/guf_rand_impl.c src/test/guf_sort_impl.c) add_executable(libguf_test src/test/test.cpp src/test/guf_init_impl.c src/test/guf_dbuf_impl.c src/test/guf_str_impl.c src/test/guf_dict_impl.c src/test/guf_rand_impl.c src/test/guf_sort_impl.c src/test/guf_linalg_impl.c)
target_include_directories(libguf_test PRIVATE src src/test) target_include_directories(libguf_test PRIVATE src src/test)
if (NOT DEFINED MSVC) if (NOT DEFINED MSVC)

181
src/guf_math_linalg.h → src/guf_linalg.h
View File

@ -1,11 +1,11 @@
#if defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(GUF_LINALG_IMPL_STATIC)
#define GUF_MATH_LINALG_KWRDS static #define GUF_LINALG_KWRDS static
#else #else
#define GUF_MATH_LINALG_KWRDS #define GUF_LINALG_KWRDS
#endif #endif
#ifndef GUF_MATH_LINALG_H #ifndef GUF_LINALG_H
#define GUF_MATH_LINALG_H #define GUF_LINALG_H
#include "guf_assert.h" #include "guf_assert.h"
#include "guf_math.h" #include "guf_math.h"
@ -51,44 +51,44 @@
float angle; float angle;
} guf_axis_angle; } guf_axis_angle;
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m); GUF_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m);
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m); GUF_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m);
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis); GUF_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis);
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q); GUF_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q);
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m); GUF_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m);
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m); GUF_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m);
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis); GUF_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis);
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q); GUF_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q);
#ifdef __cplusplus #ifdef __cplusplus
// C++ has no restrict keyword like C99... // C++ has no restrict keyword like C99...
GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed); GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed);
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed); GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed);
#endif #endif
#ifdef __cplusplus #ifdef __cplusplus
// C++ has no restrict keyword like C99... // C++ has no restrict keyword like C99...
GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed); GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed);
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed); GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed);
#endif #endif
#ifdef __cplusplus #ifdef __cplusplus
// C++ has no restrict keyword like C99... // C++ has no restrict keyword like C99...
GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result); GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result);
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4f *a, const guf_mat4x4f *b, guf_mat4x4f * restrict result); GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 * restrict result);
#endif #endif
GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted); GUF_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted);
#ifdef __cplusplus #ifdef __cplusplus
// C++ has no restrict keyword like C99... // C++ has no restrict keyword like C99...
GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result); GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result);
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result); GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result);
#endif #endif
static inline void guf_mat4x4_set_trans(guf_mat4x4 *m, guf_vec3 trans) static inline void guf_mat4x4_set_trans(guf_mat4x4 *m, guf_vec3 trans)
@ -99,12 +99,12 @@
m->data[2][3] = trans.z; m->data[2][3] = trans.z;
} }
GUF_MATH_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3); GUF_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3);
GUF_MATH_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps); GUF_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps);
GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v); GUF_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v);
GUF_MATH_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v); GUF_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v);
static inline guf_vec2 guf_vec2_add(guf_vec2 a, guf_vec2 b) static inline guf_vec2 guf_vec2_add(guf_vec2 a, guf_vec2 b)
{ {
@ -174,13 +174,13 @@
return a.x * b.y - a.y * b.x; // The z-component of a 3d cross-product of two vectors in a plane. return a.x * b.y - a.y * b.x; // The z-component of a 3d cross-product of two vectors in a plane.
} }
GUF_MATH_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v); GUF_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v);
GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v); GUF_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v);
static inline guf_vec2 guf_vec2_normalised(guf_vec2 v) {return *guf_vec2_normalise(&v);}; static inline guf_vec2 guf_vec2_normalised(guf_vec2 v) {return *guf_vec2_normalise(&v);}
static inline guf_vec3 guf_vec3_normalised(guf_vec3 v) {return *guf_vec3_normalise(&v);}; static inline guf_vec3 guf_vec3_normalised(guf_vec3 v) {return *guf_vec3_normalise(&v);}
GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v); GUF_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v);
static inline void guf_quaternion_init_identity(guf_quaternion *q) static inline void guf_quaternion_init_identity(guf_quaternion *q)
{ {
@ -199,7 +199,7 @@
return q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z; return q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z;
} }
GUF_MATH_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q); GUF_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q);
static inline guf_quaternion guf_quaternion_normalised(guf_quaternion q) static inline guf_quaternion guf_quaternion_normalised(guf_quaternion q)
{ {
@ -216,28 +216,29 @@
#ifdef __cplusplus #ifdef __cplusplus
// C++ has no restrict keyword like C99... // C++ has no restrict keyword like C99...
GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result); GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result);
#else #else
GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result); GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result);
#endif #endif
GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q); GUF_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q);
static inline guf_vec3 guf_vec3_rotated_by_quat(guf_vec3 p, const guf_quaternion *q) static inline guf_vec3 guf_vec3_rotated_by_quat(guf_vec3 p, const guf_quaternion *q)
{ {
return *guf_vec3_rotate_by_quat(&p, q); return *guf_vec3_rotate_by_quat(&p, q);
} }
GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler); GUF_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler);
GUF_MATH_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q); GUF_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q);
#endif #endif
#define GUF_LINALG_IMPL_STATIC /* debug */
#if defined(GUF_MATH_LINALG_IMPL) || defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(GUF_LINALG_IMPL) || defined(GUF_LINALG_IMPL_STATIC)
#include <string.h> #include <string.h>
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m) GUF_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m)
{ {
GUF_ASSERT(m); GUF_ASSERT(m);
for (int row = 0; row < 4; ++row) { for (int row = 0; row < 4; ++row) {
@ -247,7 +248,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m)
} }
} }
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m) GUF_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m)
{ {
GUF_ASSERT(m); GUF_ASSERT(m);
for (int row = 0; row < 3; ++row) { for (int row = 0; row < 3; ++row) {
@ -257,7 +258,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m)
} }
} }
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m) GUF_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m)
{ {
GUF_ASSERT(m); GUF_ASSERT(m);
guf_mat4x4_init_zero(m); guf_mat4x4_init_zero(m);
@ -266,7 +267,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m)
} }
} }
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m) GUF_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m)
{ {
GUF_ASSERT(m); GUF_ASSERT(m);
guf_mat3x3_init_zero(m); guf_mat3x3_init_zero(m);
@ -275,7 +276,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m)
} }
} }
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis) GUF_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis)
{ {
GUF_ASSERT(m); GUF_ASSERT(m);
@ -298,7 +299,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_
m->data[3][3] = 1; m->data[3][3] = 1;
} }
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis) GUF_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis)
{ {
GUF_ASSERT(m); GUF_ASSERT(m);
@ -315,10 +316,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_
m->data[2][2] = z_basis.z; m->data[2][2] = z_basis.z;
} }
#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result) GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result)
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 * restrict result) GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 * restrict result)
#endif #endif
{ {
GUF_ASSERT(a && b && result); GUF_ASSERT(a && b && result);
@ -333,10 +334,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4
} }
} }
#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result) GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result)
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result) GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result)
#endif #endif
{ {
GUF_ASSERT(a && b && result); GUF_ASSERT(a && b && result);
@ -351,10 +352,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3
} }
} }
#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed) GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed)
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed) GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed)
#endif #endif
{ {
GUF_ASSERT(m && transposed); GUF_ASSERT(m && transposed);
@ -366,10 +367,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4
} }
} }
#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed) GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed)
#else #else
GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed) GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed)
#endif #endif
{ {
GUF_ASSERT(m && transposed); GUF_ASSERT(m && transposed);
@ -381,7 +382,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3
} }
} }
GUF_MATH_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3) GUF_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3)
{ {
GUF_ASSERT(m && rotmat_3x3) GUF_ASSERT(m && rotmat_3x3)
for (int row = 0; row < 3; ++row) { for (int row = 0; row < 3; ++row) {
@ -391,16 +392,16 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3
} }
} }
GUF_MATH_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps) GUF_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps)
{ {
for (int row = 0; row < 4; ++row) { for (int row = 0; row < 4; ++row) {
for (int col = 0; col < 4; ++col) { for (int col = 0; col < 4; ++col) {
if (row == col) { // Check 1 entries. if (row == col) { // Check 1 entries.
if (fabs(m->data[row][col] - 1) > eps) { if (!guf_nearly_one_f32(m->data[row][col], eps)) {
return false; return false;
} }
} else { // Check 0 entries. } else { // Check 0 entries.
if (fabs(m->data[row][col]) > eps) { if (!guf_nearly_one_f32(m->data[row][col], eps)) {
return false; return false;
} }
} }
@ -431,7 +432,7 @@ static void guf_mat4x4_swap_rows(guf_mat4x4 *m, int row_a_idx, int row_b_idx)
Note: In the special case of square orthonormal matrices, the matrix inverse can be calculated more efficiently Note: In the special case of square orthonormal matrices, the matrix inverse can be calculated more efficiently
and accurately using matrix transposition (e.g. valid 3x3 rotation matrices can be inverted by transposing them). and accurately using matrix transposition (e.g. valid 3x3 rotation matrices can be inverted by transposing them).
*/ */
GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted) GUF_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted)
{ {
GUF_ASSERT(m && inverted); GUF_ASSERT(m && inverted);
@ -441,9 +442,6 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
const float EPS = 1e-6f; const float EPS = 1e-6f;
const float EPS_ASSERT = 0.05f; const float EPS_ASSERT = 0.05f;
#define guf_nearly_zero(x, eps) (fabsf((x)) < (eps))
#define guf_nearly_one(x, eps) (fabsf((x) - 1) < (eps))
for (int row_start = 0; row_start < 3; ++row_start) { // 1.) Try to bring into row echelon form. for (int row_start = 0; row_start < 3; ++row_start) { // 1.) Try to bring into row echelon form.
// Find pivot, i.e. the element with row in range [row_start, 3] which has // Find pivot, i.e. the element with row in range [row_start, 3] which has
// the the largest value in the current_column (current_column = row_start) // the the largest value in the current_column (current_column = row_start)
@ -451,7 +449,7 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
float largest = -INFINITY; float largest = -INFINITY;
for (int row = row_start; row < 4; ++row) { for (int row = row_start; row < 4; ++row) {
const float candidate = fabsf(rref.data[row][row_start]); const float candidate = fabsf(rref.data[row][row_start]);
if (!guf_nearly_zero(candidate, EPS) && candidate > largest) { if (!guf_nearly_zero_f32(candidate, EPS) && candidate > largest) {
largest = candidate; largest = candidate;
first_nonzero_row = row; first_nonzero_row = row;
} }
@ -469,11 +467,11 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
const int pivot_row = row_start; const int pivot_row = row_start;
const int pivot_col = row_start; const int pivot_col = row_start;
const float pivot_val = rref.data[pivot_row][pivot_col]; const float pivot_val = rref.data[pivot_row][pivot_col];
GUF_ASSERT(!guf_nearly_zero(pivot_val, EPS)); GUF_ASSERT(!guf_nearly_zero_f32(pivot_val, EPS));
for (int row = pivot_row + 1; row < 4; ++row) { // Make all values below the pivot element zero. for (int row = pivot_row + 1; row < 4; ++row) { // Make all values below the pivot element zero.
const float val = rref.data[row][pivot_col]; // val: the current value we want to turn into zero. const float val = rref.data[row][pivot_col]; // val: the current value we want to turn into zero.
if (guf_nearly_zero(val, EPS)) { if (guf_nearly_zero_f32(val, EPS)) {
rref.data[row][pivot_col] = 0; rref.data[row][pivot_col] = 0;
continue; continue;
} }
@ -492,21 +490,21 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
} }
for (int i = 3; i >= 0; --i) { // 2.) Try to bring into *reduced* row echelon form. for (int i = 3; i >= 0; --i) { // 2.) Try to bring into *reduced* row echelon form.
if (guf_nearly_zero(rref.data[i][i], EPS)) { if (guf_nearly_zero_f32(rref.data[i][i], EPS)) {
return false; return false;
} }
for (int col = 0; col < 4; ++col) { // Scale row_i to make elem[i][i] 1 for (int col = 0; col < 4; ++col) { // Scale row_i to make elem[i][i] 1
inverted->data[i][col] *= 1.f / rref.data[i][i]; inverted->data[i][col] *= 1.f / rref.data[i][i];
GUF_ASSERT((col != i) ? rref.data[i][col] == 0 : rref.data[i][col] != 0); GUF_ASSERT((col != i) ? rref.data[i][col] == 0 : rref.data[i][col] != 0);
} }
GUF_ASSERT(guf_nearly_one(rref.data[i][i], EPS_ASSERT)); GUF_ASSERT(guf_nearly_one_f32(rref.data[i][i], EPS_ASSERT));
rref.data[i][i] = 1; // *= 1.f / rref.data[i][i] (== 1) rref.data[i][i] = 1; // *= 1.f / rref.data[i][i] (== 1)
for (int row = i - 1; row >= 0; --row) { // Make all elements in column above rref.data[i][i] zero (by adding a scaled row_i to row). for (int row = i - 1; row >= 0; --row) { // Make all elements in column above rref.data[i][i] zero (by adding a scaled row_i to row).
for (int col = 0; col < 4; ++col) { for (int col = 0; col < 4; ++col) {
inverted->data[row][col] += inverted->data[i][col] * -rref.data[row][i]; inverted->data[row][col] += inverted->data[i][col] * -rref.data[row][i];
} }
GUF_ASSERT(guf_nearly_zero(rref.data[row][i], EPS_ASSERT)); GUF_ASSERT(guf_nearly_zero_f32(rref.data[row][i], EPS_ASSERT));
rref.data[row][i] = 0; // += rref.data[i][i] * -rref.data[row][i] (== 0) rref.data[row][i] = 0; // += rref.data[i][i] * -rref.data[row][i] (== 0)
} }
} }
@ -514,11 +512,10 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
GUF_ASSERT(guf_mat4x4_is_identity(&rref, EPS)); GUF_ASSERT(guf_mat4x4_is_identity(&rref, EPS));
(void)EPS_ASSERT; (void)EPS_ASSERT;
return true; return true;
#undef guf_nearly_zero
} }
// Assumes the last row of mat is [0, 0, 0, 1] and v->w implicitly is 1 (affine transformation) // Assumes the last row of mat is [0, 0, 0, 1] and v->w implicitly is 1 (affine transformation)
GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v) GUF_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v)
{ {
GUF_ASSERT(mat && v); GUF_ASSERT(mat && v);
guf_vec3 v_transformed = { guf_vec3 v_transformed = {
@ -529,7 +526,7 @@ GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const
return v_transformed; return v_transformed;
} }
GUF_MATH_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v) GUF_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v)
{ {
GUF_ASSERT(mat && v); GUF_ASSERT(mat && v);
guf_vec4 v_transformed = { guf_vec4 v_transformed = {
@ -541,7 +538,7 @@ GUF_MATH_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const
return v_transformed; return v_transformed;
} }
GUF_MATH_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v) GUF_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v)
{ {
GUF_ASSERT(v); GUF_ASSERT(v);
const float mag = guf_vec2_mag(*v); const float mag = guf_vec2_mag(*v);
@ -554,7 +551,7 @@ GUF_MATH_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v)
return v; return v;
} }
GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v) GUF_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v)
{ {
GUF_ASSERT(v); GUF_ASSERT(v);
const float mag = guf_vec3_mag(*v); const float mag = guf_vec3_mag(*v);
@ -568,10 +565,10 @@ GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v)
return v; return v;
} }
GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v) GUF_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v)
{ {
if (guf_isclose_abstol_f32(v.w, 0, 1e-21f)) { // If w is extremely close 0, we treat it as zero. if (guf_isclose_abstol_f32(v.w, 0, 1e-21f)) { // If w is extremely close 0, we treat it as zero.
guf_vec3 inf_vec = {(float)INFINITY, (float)INFINITY, float(INFINITY)}; guf_vec3 inf_vec = {(float)INFINITY, (float)INFINITY, (float)INFINITY};
return inf_vec; return inf_vec;
} else if (guf_isclose_reltol_f32(v.w, 1, 1e-9f)) { // If w is almost 1, we treat it as exactly 1. } else if (guf_isclose_reltol_f32(v.w, 1, 1e-9f)) { // If w is almost 1, we treat it as exactly 1.
guf_vec3 vec = {.x = v.x, .y = v.y, .z = v.z}; guf_vec3 vec = {.x = v.x, .y = v.y, .z = v.z};
@ -580,14 +577,14 @@ GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v)
guf_vec3 vec = {.x = v.x / v.w, .y = v.y / v.w, .z = v.z / v.w}; guf_vec3 vec = {.x = v.x / v.w, .y = v.y / v.w, .z = v.z / v.w};
return vec; return vec;
} }
}; }
/* /*
(Unit)quaternion operations: (Unit)quaternion operations:
cf. https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html (last-retrieved 2025-03-06) cf. https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html (last-retrieved 2025-03-06)
*/ */
GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_axis_angle(float angle, guf_vec3 unit_axis) GUF_LINALG_KWRDS guf_quaternion guf_quaternion_from_axis_angle(float angle, guf_vec3 unit_axis)
{ {
const float sin_halfangle = sinf(angle / 2.f); const float sin_halfangle = sinf(angle / 2.f);
guf_quaternion q = { guf_quaternion q = {
@ -599,7 +596,7 @@ GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_axis_angle(float angle,
return q; return q;
} }
GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler) GUF_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler)
{ {
const float r = euler.roll / 2; const float r = euler.roll / 2;
const float p = euler.pitch / 2; const float p = euler.pitch / 2;
@ -613,7 +610,7 @@ GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles
return q; return q;
} }
GUF_MATH_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q) GUF_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q)
{ {
guf_euler_angles euler = {0}; guf_euler_angles euler = {0};
euler.pitch = asinf(2 * (q.w * q.y - q.x * q.z)); euler.pitch = asinf(2 * (q.w * q.y - q.x * q.z));
@ -632,7 +629,7 @@ GUF_MATH_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q)
return euler; return euler;
} }
GUF_MATH_LINALG_KWRDS float guf_quaternion_to_axis_angle(const guf_quaternion *q, guf_vec3 *axis_result) GUF_LINALG_KWRDS float guf_quaternion_to_axis_angle(const guf_quaternion *q, guf_vec3 *axis_result)
{ {
GUF_ASSERT(q); GUF_ASSERT(q);
if (fabsf(q->w - 1) <= GUF_QUAT_TO_AXIS_ANGLE_EPS) { // Singularity: q->w is very close to 1 (identity-quaternions produce no rotation). if (fabsf(q->w - 1) <= GUF_QUAT_TO_AXIS_ANGLE_EPS) { // Singularity: q->w is very close to 1 (identity-quaternions produce no rotation).
@ -653,7 +650,7 @@ GUF_MATH_LINALG_KWRDS float guf_quaternion_to_axis_angle(const guf_quaternion *q
return angle; return angle;
} }
GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q) GUF_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q)
{ {
GUF_ASSERT(rotmat); GUF_ASSERT(rotmat);
@ -676,7 +673,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, g
rotmat->data[3][3] = 1; rotmat->data[3][3] = 1;
} }
GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q) GUF_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q)
{ {
GUF_ASSERT(rotmat); GUF_ASSERT(rotmat);
@ -693,7 +690,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, g
rotmat->data[2][2] = 1 - 2 * q.x * q.x - 2 * q.y * q.y; rotmat->data[2][2] = 1 - 2 * q.x * q.x - 2 * q.y * q.y;
} }
GUF_MATH_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q) GUF_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q)
{ {
GUF_ASSERT(q); GUF_ASSERT(q);
const float sq_mag = guf_quaternion_mag_squared(*q); const float sq_mag = guf_quaternion_mag_squared(*q);
@ -712,10 +709,10 @@ GUF_MATH_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q
return q; return q;
} }
#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC) #if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result) GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result)
#else #else
GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result) GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result)
#endif #endif
{ {
GUF_ASSERT(a && b && result); GUF_ASSERT(a && b && result);
@ -726,7 +723,7 @@ GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf
result->z = a->w * b->z - a->x * b->y + a->y * b->x + a->z * b->w; // (r0s3 − r1s2 + r2s1 + r3s0) result->z = a->w * b->z - a->x * b->y + a->y * b->x + a->z * b->w; // (r0s3 − r1s2 + r2s1 + r3s0)
} }
GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q) GUF_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q)
{ {
// Active rotation of the point/vector p by the quaternion q. // Active rotation of the point/vector p by the quaternion q.
GUF_ASSERT(p && q); GUF_ASSERT(p && q);
@ -745,9 +742,9 @@ GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_q
return p; return p;
} }
#undef GUF_MATH_LINALG_IMPL #undef GUF_LINALG_IMPL
#undef GUF_MATH_LINALG_IMPL_STATIC #undef GUF_LINALG_IMPL_STATIC
#endif /* end impl */ #endif /* end impl */
#undef GUF_MATH_LINALG_KWRDS #undef GUF_LINALG_KWRDS

8
src/guf_math.h
View File

@ -72,8 +72,14 @@ static inline uint16_t guf_absdiff_i16(int16_t a, int16_t b) {return a > b ? (u
static inline uint32_t guf_absdiff_i32(int32_t a, int32_t b) {return a > b ? (uint32_t)a - (uint32_t)b : (uint32_t)b - (uint32_t)a;} static inline uint32_t guf_absdiff_i32(int32_t a, int32_t b) {return a > b ? (uint32_t)a - (uint32_t)b : (uint32_t)b - (uint32_t)a;}
static inline uint64_t guf_absdiff_i64(int64_t a, int64_t b) {return a > b ? (uint64_t)a - (uint64_t)b : (uint64_t)b - (uint64_t)a;} static inline uint64_t guf_absdiff_i64(int64_t a, int64_t b) {return a > b ? (uint64_t)a - (uint64_t)b : (uint64_t)b - (uint64_t)a;}
static bool guf_nearly_zero_f32(float x, float eps) {return fabsf(x) <= eps;}
static bool guf_nearly_one_f32(float x, float eps) {return fabsf(x - 1) <= eps;}
static bool guf_nearly_zero_f64(double x, double eps) {return fabs(x) <= eps;}
static bool guf_nearly_one_f64(double x, double eps) {return fabs(x - 1) <= eps;}
/* /*
Floating-point comparison: General floating-point comparison:
cf. https://peps.python.org/pep-0485/ (last-retrieved 2025-03-04) cf. https://peps.python.org/pep-0485/ (last-retrieved 2025-03-04)
https://github.com/PythonCHB/close_pep/blob/master/is_close.py (last-retrieved 2025-03-04) https://github.com/PythonCHB/close_pep/blob/master/is_close.py (last-retrieved 2025-03-04)

2
src/test/guf_linalg_impl.c Normal file
View File

@ -0,0 +1,2 @@
#define GUF_MATH_LINALG_IMPL
#include "guf_linalg.h"