libguf/src/guf_math.h

287 lines
16 KiB
C

/*
is parametrized: no
*/
#ifndef GUF_MATH_H
#define GUF_MATH_H
#include "guf_common.h"
#include "guf_assert.h"
#include <math.h>
#include <float.h>
#define GUF_PI 3.14159265358979323846264338327950288
#define GUF_PI_F32 3.14159265358979323846264338327950288f
// Biggest float/double less than one, cf. https://stackoverflow.com/a/33832753 (last-retrieved 2025-03-05)
#define GUF_MAX_F32_LT_ONE (1.f - FLT_EPSILON/FLT_RADIX)
#define GUF_MAX_F64_LT_ONE (1.0 - DBL_EPSILON/FLT_RADIX)
// Rotate left.
static inline uint64_t guf_rotl_u64(uint64_t x, int k) {return (1u*x << k) | (1u*x >> (64 - k));}
static inline uint32_t guf_rotl_u32(uint32_t x, int k) {return (1u*x << k) | (1u*x >> (32 - k));}
// Signed min, max, clamp functions (generated with libguf/tools/min_max_clamp-gen.py)
static inline char guf_min_char(char a, char b) { return a < b ? a : b; }
static inline char guf_max_char(char a, char b) { return a > b ? a : b; }
static inline char guf_clamp_char(char x, char min, char max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline int guf_min_int(int a, int b) { return a < b ? a : b; }
static inline int guf_max_int(int a, int b) { return a > b ? a : b; }
static inline int guf_clamp_int(int x, int min, int max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline long guf_min_long(long a, long b) { return a < b ? a : b; }
static inline long guf_max_long(long a, long b) { return a > b ? a : b; }
static inline long guf_clamp_long(long x, long min, long max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline long long guf_min_long_long(long long a, long long b) { return a < b ? a : b; }
static inline long long guf_max_long_long(long long a, long long b) { return a > b ? a : b; }
static inline long long guf_clamp_long_long(long long x, long long min, long long max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline int8_t guf_min_i8(int8_t a, int8_t b) { return a < b ? a : b; }
static inline int8_t guf_max_i8(int8_t a, int8_t b) { return a > b ? a : b; }
static inline int8_t guf_clamp_i8(int8_t x, int8_t min, int8_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline int16_t guf_min_i16(int16_t a, int16_t b) { return a < b ? a : b; }
static inline int16_t guf_max_i16(int16_t a, int16_t b) { return a > b ? a : b; }
static inline int16_t guf_clamp_i16(int16_t x, int16_t min, int16_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline int32_t guf_min_i32(int32_t a, int32_t b) { return a < b ? a : b; }
static inline int32_t guf_max_i32(int32_t a, int32_t b) { return a > b ? a : b; }
static inline int32_t guf_clamp_i32(int32_t x, int32_t min, int32_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline int64_t guf_min_i64(int64_t a, int64_t b) { return a < b ? a : b; }
static inline int64_t guf_max_i64(int64_t a, int64_t b) { return a > b ? a : b; }
static inline int64_t guf_clamp_i64(int64_t x, int64_t min, int64_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline ptrdiff_t guf_min_ptrdiff_t(ptrdiff_t a, ptrdiff_t b) { return a < b ? a : b; }
static inline ptrdiff_t guf_max_ptrdiff_t(ptrdiff_t a, ptrdiff_t b) { return a > b ? a : b; }
static inline ptrdiff_t guf_clamp_ptrdiff_t(ptrdiff_t x, ptrdiff_t min, ptrdiff_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline float guf_min_f32(float a, float b) { return a < b ? a : b; }
static inline float guf_max_f32(float a, float b) { return a > b ? a : b; }
static inline float guf_clamp_f32(float x, float min, float max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline double guf_min_f64(double a, double b) { return a < b ? a : b; }
static inline double guf_max_f64(double a, double b) { return a > b ? a : b; }
static inline double guf_clamp_f64(double x, double min, double max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
// Unsigned min, max, clamp functions (generated with libguf/tools/min_max_clamp-gen.py)
static inline unsigned char guf_min_uchar(unsigned char a, unsigned char b) { return a < b ? a : b; }
static inline unsigned char guf_max_uchar(unsigned char a, unsigned char b) { return a > b ? a : b; }
static inline unsigned char guf_clamp_uchar(unsigned char x, unsigned char min, unsigned char max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline unsigned guf_min_unsigned(unsigned a, unsigned b) { return a < b ? a : b; }
static inline unsigned guf_max_unsigned(unsigned a, unsigned b) { return a > b ? a : b; }
static inline unsigned guf_clamp_unsigned(unsigned x, unsigned min, unsigned max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline unsigned long guf_min_ulong(unsigned long a, unsigned long b) { return a < b ? a : b; }
static inline unsigned long guf_max_ulong(unsigned long a, unsigned long b) { return a > b ? a : b; }
static inline unsigned long guf_clamp_ulong(unsigned long x, unsigned long min, unsigned long max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline unsigned long long guf_min_ulong_long(unsigned long long a, unsigned long long b) { return a < b ? a : b; }
static inline unsigned long long guf_max_ulong_long(unsigned long long a, unsigned long long b) { return a > b ? a : b; }
static inline unsigned long long guf_clamp_ulong_long(unsigned long long x, unsigned long long min, unsigned long long max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline uint8_t guf_min_u8(uint8_t a, uint8_t b) { return a < b ? a : b; }
static inline uint8_t guf_max_u8(uint8_t a, uint8_t b) { return a > b ? a : b; }
static inline uint8_t guf_clamp_u8(uint8_t x, uint8_t min, uint8_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline uint16_t guf_min_u16(uint16_t a, uint16_t b) { return a < b ? a : b; }
static inline uint16_t guf_max_u16(uint16_t a, uint16_t b) { return a > b ? a : b; }
static inline uint16_t guf_clamp_u16(uint16_t x, uint16_t min, uint16_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline uint32_t guf_min_u32(uint32_t a, uint32_t b) { return a < b ? a : b; }
static inline uint32_t guf_max_u32(uint32_t a, uint32_t b) { return a > b ? a : b; }
static inline uint32_t guf_clamp_u32(uint32_t x, uint32_t min, uint32_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline uint64_t guf_min_u64(uint64_t a, uint64_t b) { return a < b ? a : b; }
static inline uint64_t guf_max_u64(uint64_t a, uint64_t b) { return a > b ? a : b; }
static inline uint64_t guf_clamp_u64(uint64_t x, uint64_t min, uint64_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
static inline size_t guf_min_size_t(size_t a, size_t b) { return a < b ? a : b; }
static inline size_t guf_max_size_t(size_t a, size_t b) { return a > b ? a : b; }
static inline size_t guf_clamp_size_t(size_t x, size_t min, size_t max) { if (x < min) {return min;} if (x > max) {return max;} return x; }
// abs functions with signed result (can fail/panic for abs(INT_TYPE_MIN) for platforms with two's complement signed ints; C2X for example guarantees two's complement)
static inline int guf_abs_int(int x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(-INT_MAX == INT_MIN || x > INT_MIN); return -x;} // I would not drink that...
static inline int8_t guf_abs_i8 (int8_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(-INT8_MAX == INT8_MIN || x > INT8_MIN); return -x;}
static inline int16_t guf_abs_i16(int16_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(-INT16_MAX == INT16_MIN || x > INT16_MIN); return -x;}
static inline int32_t guf_abs_i32(int32_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(-INT32_MAX == INT32_MIN || x > INT32_MIN); return -x;}
static inline int64_t guf_abs_i64(int64_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(-INT64_MAX == INT64_MIN || x > INT64_MIN); return -x;}
static inline ptrdiff_t guf_abs_ptrdiff(ptrdiff_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(-PTRDIFF_MAX == PTRDIFF_MIN || x > PTRDIFF_MIN); return -x;}
// abs functions with unsigned result functions (cannot fail)
static inline unsigned guf_uabs_int(int x) {if (x >= 0) {return x;} else if (x == INT_MIN && -INT_MAX != INT_MIN) {return (unsigned)INT_MAX + 1u;} else {return -x;}}
static inline uint8_t guf_uabs_i8(int8_t x) {if (x >= 0) {return x;} else if (x == INT8_MIN && -INT8_MAX != INT8_MIN) {return (uint8_t)INT8_MAX + 1u;} else {return -x;}}
static inline uint16_t guf_uabs_i16(int16_t x) {if (x >= 0) {return x;} else if (x == INT16_MIN && -INT16_MAX != INT16_MIN) {return (uint16_t)INT16_MAX + 1u;} else {return -x;}}
static inline uint32_t guf_uabs_i32(int32_t x) {if (x >= 0) {return x;} else if (x == INT32_MIN && -INT32_MAX != INT32_MIN) {return (uint32_t)INT32_MAX + 1u;} else {return -x;}}
static inline uint64_t guf_uabs_i64(int64_t x) {if (x >= 0) {return x;} else if (x == INT64_MIN && -INT64_MAX != INT64_MIN) {return (uint64_t)INT64_MAX + 1u;} else {return -x;}}
static inline size_t guf_uabs_ptrdiff_t(ptrdiff_t x) {if (x >= 0) {return x;} else if (x == PTRDIFF_MIN && -PTRDIFF_MAX != PTRDIFF_MIN) {return (size_t)PTRDIFF_MAX + 1u;} else {return -x;}}
// absdiff functions with unsigned result (cannot fail)
static inline unsigned char guf_absdiff_char(char a, char b) {return a > b ? (unsigned char)a - (unsigned char)b : (unsigned char)b - (unsigned char)a;}
static inline unsigned guf_absdiff_int(int a, int b) {return a > b ? (unsigned)a - (unsigned)b : (unsigned)b - (unsigned)a;}
static inline uint8_t guf_absdiff_i8(int8_t a, int8_t b) {return a > b ? (uint8_t)a - (uint8_t)b : (uint8_t)b - (uint8_t)a;}
static inline uint16_t guf_absdiff_i16(int16_t a, int16_t b) {return a > b ? (uint16_t)a - (uint16_t)b : (uint16_t)b - (uint16_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 size_t guf_absdiff_ptrdiff_t(ptrdiff_t a, ptrdiff_t b) {return a > b ? (size_t)a - (size_t)b : (size_t)b - (size_t)a;}
static inline bool guf_add_is_overflow_size_t(size_t a, size_t b)
{
return (a + b) < a;
}
static inline bool guf_sub_is_overflow_size_t(size_t a, size_t b)
{
return (a - b) > a;
}
static inline bool guf_mul_is_overflow_size_t(size_t a, size_t b)
{
const size_t c = a * b;
return a != 0 && ((c / a) != b);
}
static inline bool guf_size_calc_safe(ptrdiff_t count, ptrdiff_t sizeof_elem, ptrdiff_t *result)
{
if (count < 0 || sizeof_elem <= 0) {
return false;
}
if (guf_mul_is_overflow_size_t((size_t)count, (size_t)sizeof_elem)) {
return false;
}
const size_t size = (size_t)count * (size_t)sizeof_elem;
const bool is_safe = size <= PTRDIFF_MAX;
if (result) {
*result = is_safe ? (ptrdiff_t)size : -1;
}
return is_safe;
}
// cf. https://graphics.stanford.edu/~seander/bithacks.html#DetermineIfPowerOf2 (last-retrieved 2025-03-19)
static inline bool guf_is_pow2_u8(uint8_t x) { return x && !(x & (x - 1)); }
static inline bool guf_is_pow2_u16(uint16_t x) { return x && !(x & (x - 1)); }
static inline bool guf_is_pow2_u32(uint32_t x) { return x && !(x & (x - 1)); }
static inline bool guf_is_pow2_u64(uint64_t x) { return x && !(x & (x - 1)); }
static inline bool guf_is_pow2_size_t(size_t x) { return x && !(x & (x - 1)); }
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;}
/*
General floating-point comparison:
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)
Based on is_close.py (Copyright: Christopher H. Barker, License: Apache License 2.0 http://opensource.org/licenses/apache2.0.php)
rel_tol: "The relative tolerance -- the amount of error allowed, relative to the magnitude of the input values."
(Example: To set a tolerance of 5%, pass rel_tol=0.05)
(Example: rel_tol=1e-9 assures that the two values are the same within about 9 decimal digits)
abs_tol: "The minimum absolute tolerance level -- useful for comparisons to zero." (or close to zero, I think).
*/
static bool guf_isclose_f32(float a, float b, float rel_tol, float abs_tol)
{
if (a == b) {
return true;
}
if (isinf(a) || isinf(b)) { // "Two infinities of opposite sign, or one infinity and one finite number."
return false;
}
rel_tol = (rel_tol < 0) ? 0 : ((rel_tol > 1) ? 1 : rel_tol); // [0, 1]
abs_tol = (abs_tol < 0) ? 0 : abs_tol; // [0, inf]
// "The relative tolerance is scaled by the larger of the two values."
const float diff = fabsf(b - a);
return ((diff <= fabsf(rel_tol * b)) || (diff <= fabsf(rel_tol * a))) || (diff <= abs_tol);
}
static inline bool guf_isclose_reltol_f32(float a, float b, float rel_tol)
{
return guf_isclose_f32(a, b, rel_tol, 0);
}
static inline bool guf_isclose_abstol_f32(float a, float b, float abs_tol)
{
return guf_isclose_f32(a, b, 0, abs_tol);
}
static bool guf_isclose_f64(double a, double b, double rel_tol, double abs_tol)
{
if (a == b) {
return true;
}
if (isinf(a) || isinf(b)) { // "Two infinities of opposite sign, or one infinity and one finite number."
return false;
}
rel_tol = (rel_tol < 0) ? 0 : ((rel_tol > 1) ? 1 : rel_tol); // [0, 1]
abs_tol = (abs_tol < 0) ? 0 : abs_tol; // [0, inf]
// "The relative tolerance is scaled by the larger of the two values."
const double diff = fabs(b - a);
return ((diff <= fabs(rel_tol * b)) || (diff <= fabs(rel_tol * a))) || (diff <= abs_tol);
}
static inline bool guf_isclose_reltol_f64(double a, double b, double rel_tol)
{
return guf_isclose_f64(a, b, rel_tol, 0);
}
static inline bool guf_isclose_abstol_f64(double a, double b, double abs_tol)
{
return guf_isclose_f64(a, b, 0, abs_tol);
}
// An alternative lerp would be a + alpha * (b - a) (advantage: would be weakly monotonic, disadvantage: would not guarantee a for alpha = 0 and b for alpha = 1)
static inline float guf_lerp_f32(float a, float b, float alpha) {return (1 - alpha) * a + alpha * b;}
static inline double guf_lerp_f64(double a, double b, double alpha) {return (1 - alpha) * a + alpha * b;}
// smoothstep interpolation, cf. https://en.wikipedia.org/wiki/Smoothstep (last-retrieved 2025-02-18)
static inline float guf_smoothstep_f32(float edge0, float edge1, float x)
{
if (edge0 == edge1) { // Prevent division by zero.
return 1;
}
x = guf_clamp_f32((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
return x * x * (3.f - 2.f * x);
}
static inline float guf_smootherstep_f32(float edge0, float edge1, float x)
{
if (edge0 == edge1) { // Prevent division by zero.
return 1;
}
x = guf_clamp_f32((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
return x * x * x * (x * (6.f * x - 15.f) + 10.f);
}
static inline double guf_smoothstep_f64(double edge0, double edge1, double x)
{
if (edge0 == edge1) { // Prevent division by zero.
return 1;
}
x = guf_clamp_f64((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
return x * x * (3.0 - 2.0 * x);
}
static inline double guf_smootherstep_f64(double edge0, double edge1, double x)
{
if (edge0 == edge1) { // Prevent division by zero.
return 1;
}
x = guf_clamp_f64((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
return x * x * x * (x * (6.0 * x - 15.0) + 10.0);
}
#endif