--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+#include <stdlib.h>
+
+static CGM_INLINE float cgm_deg_to_rad(float deg)
+{
+ return M_PI * deg / 180.0f;
+}
+
+static CGM_INLINE float cgm_rad_to_deg(float rad)
+{
+ return 180.0f * rad / M_PI;
+}
+
+static CGM_INLINE float cgm_smoothstep(float a, float b, float x)
+{
+ if(x < a) return 0.0f;
+ if(x >= b) return 1.0f;
+
+ x = (x - a) / (b - a);
+ return x * x * (3.0f - 2.0f * x);
+}
+
+static CGM_INLINE float cgm_lerp(float a, float b, float t)
+{
+ return a + (b - a) * t;
+}
+
+static CGM_INLINE float cgm_logerp(float a, float b, float t)
+{
+ if(a == 0.0f) return 0.0f;
+ return a * pow(b / a, t);
+}
+
+static CGM_INLINE float cgm_bezier(float a, float b, float c, float d, float t)
+{
+ float omt, omt3, t3, f;
+ t3 = t * t * t;
+ omt = 1.0f - t;
+ omt3 = omt * omt * omt;
+ f = 3.0f * t * omt;
+
+ return (a * omt3) + (b * f * omt) + (c * f * t) + (d * t3);
+}
+
+static CGM_INLINE float cgm_bspline(float a, float b, float c, float d, float t)
+{
+ static const float mat[] = {
+ -1, 3, -3, 1,
+ 3, -6, 0, 4,
+ -3, 3, 3, 1,
+ 1, 0, 0, 0
+ };
+ float tsq = t * t;
+ cgm_vec4 tmp = cgm_wmul_m4v4(mat, cgm_wcons(a, b, c, d));
+ return cgm_wdot(cgm_wscale(tmp, 1.0f / 6.0f), cgm_wcons(tsq * t, tsq, t, 1.0f));
+}
+
+static CGM_INLINE float cgm_spline(float a, float b, float c, float d, float t)
+{
+ static const float mat[] = {
+ -1, 2, -1, 0,
+ 3, -5, 0, 2,
+ -3, 4, 1, 0,
+ 1, -1, 0, 0
+ };
+ float tsq = t * t;
+ cgm_vec4 tmp = cgm_wmul_m4v4(mat, cgm_wcons(a, b, c, d));
+ return cgm_wdot(cgm_wscale(tmp, 1.0f / 6.0f), cgm_wcons(tsq * t, tsq, t, 1.0f));
+}
+
+static CGM_INLINE cgm_vec3 cgm_discrand(float rad)
+{
+ cgm_vec3 res;
+ float theta = 2.0f * M_PI * (float)rand() / RAND_MAX;
+ float r = sqrt((float)rand() / RAND_MAX) * rad;
+ res.x = cos(theta) * r;
+ res.y = sin(theta) * r;
+ res.z = 0.0f;
+ return res;
+}
+
+static CGM_INLINE cgm_vec3 cgm_sphrand(float rad)
+{
+ cgm_vec3 res;
+ float u, v, theta, phi;
+
+ u = (float)rand() / RAND_MAX;
+ v = (float)rand() / RAND_MAX;
+
+ theta = 2.0f * M_PI * u;
+ phi = acos(2.0f * v - 1.0f);
+
+ res.x = cos(theta) * sin(phi) * rad;
+ res.y = sin(theta) * sin(phi) * rad;
+ res.z = cos(phi) * rad;
+ return res;
+}
+
+static CGM_INLINE cgm_vec3 cgm_unproject(cgm_vec3 norm_scrpos, const float *inv_viewproj)
+{
+ cgm_vec3 res;
+ cgm_vec4 pos;
+
+ pos.x = 2.0f * norm_scrpos.x - 1.0f;
+ pos.y = 2.0f * norm_scrpos.y - 1.0f;
+ pos.z = 2.0f * norm_scrpos.z - 1.0f;
+ pos.w = 1.0f;
+
+ cgm_wmul_m4v4(inv_viewproj, pos);
+
+ res.x = pos.x / pos.w;
+ res.y = pos.y / pos.w;
+ res.z = pos.z / pos.w;
+ return res;
+}
+
+static CGM_INLINE void cgm_glu_unproject(float winx, float winy, float winz,
+ const float *view, const float *proj, const int *vp,
+ float *objx, float *objy, float *objz)
+{
+ cgm_vec3 npos, res;
+ float inv_pv[16];
+
+ cgm_mcopy(inv_pv, view);
+ cgm_mmul(inv_pv, proj);
+ cgm_minverse(inv_pv);
+
+ npos.x = (winx - vp[0]) / vp[2];
+ npos.y = (winy - vp[1]) / vp[4];
+ npos.z = winz;
+
+ res = cgm_unproject(npos, inv_pv);
+
+ *objx = res.x;
+ *objy = res.y;
+ *objz = res.z;
+}
+
+static CGM_INLINE cgm_ray cgm_pick_ray(float nx, float ny, const float *viewmat,
+ const float *projmat)
+{
+ cgm_ray ray;
+ cgm_vec3 farpt;
+ float inv_pv[16];
+
+ cgm_mcopy(inv_pv, viewmat);
+ cgm_mmul(inv_pv, projmat);
+ cgm_minverse(inv_pv);
+
+ ray.origin = cgm_unproject(cgm_vcons(nx, ny, 0.0f), inv_pv);
+ farpt = cgm_unproject(cgm_vcons(nx, ny, 1.0f), inv_pv);
+
+ ray.dir.x = farpt.x - ray.origin.x;
+ ray.dir.y = farpt.y - ray.origin.y;
+ ray.dir.z = farpt.z - ray.origin.z;
+ return ray;
+}
+
+static CGM_INLINE cgm_vec3 cgm_raypos(cgm_ray ray, float t)
+{
+ cgm_vec3 res;
+ res.x = ray.origin.x + ray.dir.x * t;
+ res.y = ray.origin.y + ray.dir.y * t;
+ res.z = ray.origin.z + ray.dir.z * t;
+ return res;
+}
+
+static CGM_INLINE cgm_vec3 cgm_bary(cgm_vec3 a, cgm_vec3 b, cgm_vec3 c, cgm_vec3 pt)
+{
+ float d00, d01, d11, d20, d21, denom;
+ cgm_vec3 res, v0, v1, v2;
+
+ v0 = cgm_vsub(b, a);
+ v1 = cgm_vsub(c, a);
+ v2 = cgm_vsub(pt, a);
+
+ d00 = cgm_vdot(v0, v0);
+ d01 = cgm_vdot(v0, v1);
+ d11 = cgm_vdot(v1, v1);
+ d20 = cgm_vdot(v2, v0);
+ d21 = cgm_vdot(v2, v1);
+ denom = d00 * d11 - d01 * d01;
+
+ res.y = (d11 * d20 - d01 * d21) / denom;
+ res.z = (d00 * d21 - d01 * d20) / denom;
+ res.x = 1.0f - res.y - res.z;
+ return res;
+}
+
+static CGM_INLINE cgm_vec3 cgm_uvec_to_sph(cgm_vec3 v)
+{
+ cgm_vec3 res;
+ res.x = atan2(v.z, v.x);
+ res.y = acos(v.y);
+ res.z = 1.0f;
+ return res;
+}
+
+static CGM_INLINE cgm_vec3 cgm_sph_to_uvec(float theta, float phi)
+{
+ cgm_vec3 res;
+ res.x = sin(theta) * cos(phi);
+ res.y = sin(phi);
+ res.z = cos(theta) * cos(phi);
+ return res;
+}