// #if defined(GUF_MATH_IMPL_STATIC)
//     #define GUF_MATH_KWRDS static
// #else
//     #define GUF_MATH_KWRDS
// #endif

#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 (x << k) | (x >> (64 - k));}
static inline uint32_t guf_rotl_u32(uint32_t x, int k) {return (x << k) | (x >> (32 - k));}

#define GUF_DEFINE_MIN_MAX_CLAMP(int_type, int_type_name)\
    static inline int_type GUF_CAT(guf_min_, int_type_name)(int_type a, int_type b) {return a < b ? a : b;}\
    static inline int_type GUF_CAT(guf_max_, int_type_name)(int_type a, int_type b) {return a > b ? a : b;}\
    static inline int_type GUF_CAT(guf_clamp_, int_type_name)(int_type x, int_type min, int_type max) {if (x < min) {return min;} if (x > max) {return max;} return x;}

GUF_DEFINE_MIN_MAX_CLAMP(char, char)
GUF_DEFINE_MIN_MAX_CLAMP(int, int)
GUF_DEFINE_MIN_MAX_CLAMP(long, long)
GUF_DEFINE_MIN_MAX_CLAMP(long long, long_long)
GUF_DEFINE_MIN_MAX_CLAMP(int8_t, i8)
GUF_DEFINE_MIN_MAX_CLAMP(int16_t, i16)
GUF_DEFINE_MIN_MAX_CLAMP(int32_t, i32)
GUF_DEFINE_MIN_MAX_CLAMP(int64_t, i64)
GUF_DEFINE_MIN_MAX_CLAMP(ptrdiff_t, ptrdiff_t)

GUF_DEFINE_MIN_MAX_CLAMP(unsigned char, uchar)
GUF_DEFINE_MIN_MAX_CLAMP(unsigned, unsigned)
GUF_DEFINE_MIN_MAX_CLAMP(unsigned long, ulong)
GUF_DEFINE_MIN_MAX_CLAMP(unsigned long long, ulong_long)
GUF_DEFINE_MIN_MAX_CLAMP(uint8_t, u8)
GUF_DEFINE_MIN_MAX_CLAMP(uint16_t, u16)
GUF_DEFINE_MIN_MAX_CLAMP(uint32_t, u32)
GUF_DEFINE_MIN_MAX_CLAMP(uint64_t, u64)
GUF_DEFINE_MIN_MAX_CLAMP(size_t, size_t)

GUF_DEFINE_MIN_MAX_CLAMP(float, f32)
GUF_DEFINE_MIN_MAX_CLAMP(double, f64)

#undef GUF_DEFINE_MIN_MAX_CLAMP

static inline int       guf_abs_int(int x)             {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(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(x > INT8_MIN);    return -x;}
static inline int16_t   guf_abs_i16(int16_t x)         {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT16_MIN);   return -x;}
static inline int32_t   guf_abs_i32(int32_t x)         {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT32_MIN);   return -x;}
static inline int64_t   guf_abs_i64(int64_t x)         {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT64_MIN);   return -x;}
static inline ptrdiff_t guf_abs_ptrdiff(ptrdiff_t x)   {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > PTRDIFF_MIN); return -x;}

static inline unsigned char guf_uabs_char(char x) {if (x >= 0) {return x;} else if (x == CHAR_MIN)  {return (unsigned char)CHAR_MAX + 1;} else {return -x;}}
static inline unsigned  guf_uabs_int(int x)       {if (x >= 0) {return x;} else if (x == INT_MIN)   {return (unsigned)INT_MAX   + 1;} else {return -x;}}
static inline uint8_t   guf_uabs_i8(int8_t x)     {if (x >= 0) {return x;} else if (x == INT8_MIN)  {return (uint8_t)INT8_MAX   + 1;} else {return -x;}}
static inline uint16_t  guf_uabs_i16(int16_t x)   {if (x >= 0) {return x;} else if (x == INT16_MIN) {return (uint16_t)INT16_MAX + 1;} else {return -x;}}
static inline uint32_t  guf_uabs_i32(int32_t x)   {if (x >= 0) {return x;} else if (x == INT32_MIN) {return (uint32_t)INT32_MAX + 1;} else {return -x;}}
static inline uint64_t  guf_uabs_i64(int64_t x)   {if (x >= 0) {return x;} else if (x == INT64_MIN) {return (uint64_t)INT64_MAX + 1;} else {return -x;}}

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 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_add_is_overflow_u32(uint32_t a, uint32_t b)
{
    return (a + b) < a;
}
static inline bool guf_sub_is_overflow_u32(uint32_t a, uint32_t b)
{
    return (a - b) > a;
}
static inline bool guf_mul_is_overflow_u32(uint32_t a, uint32_t b) 
{
    const uint32_t c = a * b; 
    return a != 0 && ((c / a) != b);
}

// cf. https://stackoverflow.com/questions/199333/how-do-i-detect-unsigned-integer-overflow (last-retrieved 2025-03-17)
static inline bool guf_add_is_overflow_ptrdiff(ptrdiff_t a, ptrdiff_t b)
{
    return (b > 0 && a > PTRDIFF_MAX - b) || // a + b overflow
           (b < 0 && a < PTRDIFF_MIN - b);   // a + b underflow
}
static inline bool guf_sub_is_overflow_ptrdiff(ptrdiff_t a, ptrdiff_t b)
{
    return (b < 0 && a > PTRDIFF_MAX + b) || // a - b overflow
           (b > 0 && a < PTRDIFF_MIN + b);   // a - b underflow
}

static inline bool guf_add_is_overflow_i32(int32_t a, int32_t b)
{
    return (b > 0 && a > INT32_MAX - b) || // a + b overflow
           (b < 0 && a < INT32_MIN - b);   // a + b underflow
}
static inline bool guf_sub_is_overflow_i32(int32_t a, int32_t b)
{
    return (b < 0 && a > INT32_MAX + b) || // a - b overflow
           (b > 0 && a < INT32_MIN + b);   // a - b underflow
}

static inline bool guf_add_is_overflow_i64(int64_t a, int64_t b)
{
    return (b > 0 && a > INT64_MAX - b) || // a + b overflow
           (b < 0 && a < INT64_MIN - b);   // a + b underflow
}
static inline bool guf_sub_is_overflow_i64(int64_t a, int64_t b)
{
    return (b < 0 && a > INT64_MAX + b) || // a - b overflow
           (b > 0 && a < INT64_MIN + b);   // a - b underflow
}

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