tests: dont test isinf == true on fast math

mat3x2: fix multiplication functions

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.8.2)
 project(cglm
-  VERSION 0.9.3
+  VERSION 0.9.4
   HOMEPAGE_URL https://github.com/recp/cglm
   DESCRIPTION "OpenGL Mathematics (glm) for C"
   LANGUAGES C

CREDITS

@@ -82,3 +82,17 @@ http://github.com/microsoft/DirectXMath
 17. Pick Matrix
 
 glu project -> project.c
+
+18. Ray sphere intersection
+
+RAY TRACING GEMS
+HIGH-QUALITY AND REAL-TIME RENDERING WITH DXR AND OTHER APIS
+
+CHAPTER 7
+Precision Improvements for Ray/Sphere Intersection
+Eric Haines (1), Johannes Günther (2), and Tomas Akenine-Möller (1)
+  (1) NVIDIA
+  (2) Intel
+
+Wyman, C., and Haines, E. Getting Started with RTX Ray Tracing.
+https://github.com/NVIDIAGameWorks/GettingStartedWithRTXRayTracing

Makefile.am

@@ -216,7 +216,8 @@ cglm_struct_HEADERS = include/cglm/struct/mat4.h \
                       include/cglm/struct/project.h \
                       include/cglm/struct/sphere.h \
                       include/cglm/struct/color.h \
-                      include/cglm/struct/curve.h
+                      include/cglm/struct/curve.h \
+                      include/cglm/struct/ray.h
 
 cglm_struct_clipspacedir=$(includedir)/cglm/struct/clipspace
 cglm_struct_clipspace_HEADERS = include/cglm/struct/clipspace/persp_lh_no.h \

cglm.podspec

@@ -2,7 +2,7 @@ Pod::Spec.new do |s|
 
   # Description
   s.name         = "cglm"
-  s.version      = "0.9.2"
+  s.version      = "0.9.3"
   s.summary      = "📽 Highly Optimized Graphics Math (glm) for C"
   s.description  = <<-DESC
 cglm is math library for graphics programming for C. See the documentation or README for all features.

configure.ac

@@ -7,7 +7,7 @@
 #*****************************************************************************
 
 AC_PREREQ([2.69])
-AC_INIT([cglm], [0.9.3], [info@recp.me])
+AC_INIT([cglm], [0.9.4], [info@recp.me])
 AM_INIT_AUTOMAKE([-Wall foreign subdir-objects serial-tests])
 
 # Don't use the default cflags (-O2 -g), we set ours manually in Makefile.am.

docs/source/affine-pre.rst

@@ -110,7 +110,24 @@ Functions documentation
 
 .. c:function:: void  glm_rotate(mat4 m, float angle, vec3 axis)
 
-    rotate existing transform matrix around Z axis by angle and axis
+    rotate existing transform matrix around given axis by angle at ORIGIN (0,0,0)
+
+    **❗️IMPORTANT ❗️**
+    
+    If you need to rotate object around itself e.g. center of object or at
+    some point [of object] then `glm_rotate_at()` would be better choice to do so.
+    
+    Even if object's model transform is identiy, rotation may not be around
+    center of object if object does not lay out at ORIGIN perfectly.
+    
+    Using `glm_rotate_at()` with center of bounding shape ( AABB, Sphere ... )
+    would be an easy option to rotate around object if object is not at origin.
+    
+    One another option to rotate around itself at any point is `glm_spin()`
+    which is perfect if only rotating around model position is desired e.g. not
+    specific point on model for instance center of geometry or center of mass,
+    again if geometry is not perfectly centered at origin at identity transform,
+    rotation may not be around geometry.
 
     Parameters:
       | *[in, out]* **m**     affine transform

docs/source/api_struct.rst

@@ -91,6 +91,9 @@ To configure the Struct API namespace, you can define the following macros befor
 - **CGLM_STRUCT_API_NS**:           define name space for struct api, DEFAULT is **glms**
 - **CGLM_STRUCT_API_NAME_SUFFIX**:  define name suffix, DEFAULT is **empty** e.g defining it as #define CGLM_STRUCT_API_NAME_SUFFIX  s will add s suffix to mat4_mul -> mat4s_mul
 
+❗️ IMPORTANT ❗️
+
+It's a good idea to set up your config macros in build settings like CMake, Xcode, or Visual Studio. This is especially important if you're using features like Modules in Xcode, where adding macros directly before the **cglm** headers might not work.
 
 Detailed documentation for Struct API:
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

docs/source/conf.py

@@ -62,9 +62,9 @@ author = u'Recep Aslantas'
 # built documents.
 #
 # The short X.Y version.
-version = u'0.9.2'
+version = u'0.9.4'
 # The full version, including alpha/beta/rc tags.
-release = u'0.9.2'
+release = u'0.9.4'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/source/mat2.rst

@@ -180,7 +180,7 @@ Functions documentation
     Returns:
         scalar value e.g. Matrix1x1
 
-.. c:function:: void glm_mat2_make(float * __restrict src, mat2 dest)
+.. c:function:: void glm_mat2_make(const float * __restrict src, mat2 dest)
 
     Create mat2 matrix from pointer
 

docs/source/mat2x3.rst

@@ -23,6 +23,16 @@ Functions:
 #. :c:func:`glm_mat2x3_transpose`
 #. :c:func:`glm_mat2x3_scale`
 
+Represented
+~~~~~~~~~~~
+
+.. csv-table:: mat2x3
+   :header: "", "column 1", "column 2"
+
+   "row 1", "m00", "m10"
+   "row 2", "m01", "m11"
+   "row 3", "m02", "m12"
+
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -41,7 +51,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat2x3_make(float * __restrict src, mat2x3 dest)
+.. c:function:: void glm_mat2x3_make(const float * __restrict src, mat2x3 dest)
 
     Create mat2x3 matrix from pointer
 
@@ -51,28 +61,68 @@ Functions documentation
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination matrix2x3
 
-.. c:function:: void glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest)
+.. c:function:: void glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat3 dest)
 
     multiply m1 and m2 to dest
 
     .. code-block:: c
 
-       glm_mat2x3_mul(mat2x3, mat3x2, mat2);
+       glm_mat2x3_mul(mat2x3, mat3x2, mat3);
 
     Parameters:
       | *[in]*  **m1**    left matrix (mat2x3)
       | *[in]*  **m2**    right matrix (mat3x2)
-      | *[out]* **dest**  destination matrix (mat2)
+      | *[out]* **dest**  destination matrix (mat3)
 
-.. c:function:: void glm_mat2x3_mulv(mat2x3 m, vec3 v, vec2 dest)
+    .. csv-table:: mat2x3
+        :header: "", "column 1", "column 2"
 
-    multiply mat2x3 with vec3 (column vector) and store in dest vector
+        "row 1", "a00", "a10"
+        "row 2", "a01", "a11"
+        "row 3", "a02", "a12"
+
+    .. csv-table:: mat3x2
+        :header: "", "column 1", "column 2", "column 3"
+
+        "row 1", "b00", "b10", "b20"
+        "row 2", "b01", "b11", "b21"
+
+    .. csv-table:: mat3x3
+        :header: "", "column 1", "column 2", "column 3"
+
+        "row 1", "a00 * b00 + a10 * b01", "a00 * b10 + a10 * b11", "a00 * b20 + a10 * b21"
+        "row 2", "a01 * b00 + a11 * b01", "a01 * b10 + a11 * b11", "a01 * b20 + a11 * b21"
+        "row 3", "a02 * b00 + a12 * b01", "a02 * b10 + a12 * b11", "a02 * b20 + a12 * b21"
+
+.. c:function:: void glm_mat2x3_mulv(mat2x3 m, vec2 v, vec3 dest)
+
+    multiply mat2x3 with vec2 (column vector) and store in dest column vector
 
     Parameters:
       | *[in]*  **m**     mat2x3 (left)
       | *[in]*  **v**     vec3 (right, column vector)
       | *[out]* **dest**  destination (result, column vector)
 
+    .. csv-table:: mat2x3
+        :header: "", "column 1", "column 2"
+
+        "row 1", "m00", "m10"
+        "row 2", "m01", "m11"
+        "row 3", "m02", "m12"
+
+    .. csv-table:: column vec2 (1x2)
+        :header: "", "column 1"
+
+        "row 1", "v0"
+        "row 2", "v1"
+
+    .. csv-table:: column vec3 (1x3)
+        :header: "", "column 1"
+
+        "row 1", "m00 * v0 + m10 * v1"
+        "row 2", "m01 * v0 + m11 * v1"
+        "row 3", "m02 * v0 + m12 * v1"
+
 .. c:function:: void glm_mat2x3_transpose(mat2x3 m, mat3x2 dest)
 
     transpose matrix and store in dest

docs/source/mat2x4.rst

@@ -23,6 +23,17 @@ Functions:
 #. :c:func:`glm_mat2x4_transpose`
 #. :c:func:`glm_mat2x4_scale`
 
+Represented
+~~~~~~~~~~~
+
+.. csv-table:: mat2x4
+   :header: "", "column 1", "column 2"
+
+   "row 1", "m00", "m10"
+   "row 2", "m01", "m11"
+   "row 3", "m02", "m12"
+   "row 4", "m03", "m13"
+
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -41,7 +52,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat2x4_make(float * __restrict src, mat2x4 dest)
+.. c:function:: void glm_mat2x4_make(const float * __restrict src, mat2x4 dest)
 
     Create mat2x4 matrix from pointer
 
@@ -51,28 +62,72 @@ Functions documentation
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination matrix2x4
 
-.. c:function:: void glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest)
+.. c:function:: void glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat4 dest)
 
     multiply m1 and m2 to dest
 
     .. code-block:: c
 
-       glm_mat2x4_mul(mat2x4, mat4x2, mat2);
+       glm_mat2x4_mul(mat2x4, mat4x2, mat4);
 
     Parameters:
       | *[in]*  **m1**    left matrix (mat2x4)
       | *[in]*  **m2**    right matrix (mat4x2)
-      | *[out]* **dest**  destination matrix (mat2)
+      | *[out]* **dest**  destination matrix (mat4)
 
-.. c:function:: void glm_mat2x4_mulv(mat2x4 m, vec4 v, vec2 dest)
+    .. csv-table:: mat2x4
+        :header: "", "column 1", "column 2"
 
-    multiply mat2x4 with vec4 (column vector) and store in dest vector
+        "row 1", "a00", "a10"
+        "row 2", "a01", "a11"
+        "row 3", "a02", "a12"
+        "row 4", "a03", "a13"
+
+    .. csv-table:: mat4x2
+        :header: "", "column 1", "column 2", "column 3", "column 4"
+
+        "row 1", "b00", "b10", "b20", "b30"
+        "row 2", "b01", "b11", "b21", "b31"
+
+    .. csv-table:: mat4x4
+        :header: "", "column 1", "column 2", "column 3", "column 4"
+
+        "row 1", "a00 * b00 + a10 * b01", "a00 * b10 + a10 * b11", "a00 * b20 + a10 * b21", "a00 * b30 + a10 * b31"
+        "row 2", "a01 * b00 + a11 * b01", "a01 * b10 + a11 * b11", "a01 * b20 + a11 * b21", "a01 * b30 + a11 * b31"
+        "row 3", "a02 * b00 + a12 * b01", "a02 * b10 + a12 * b11", "a02 * b20 + a12 * b21", "a02 * b30 + a12 * b31"
+        "row 4", "a03 * b00 + a13 * b01", "a03 * b10 + a13 * b11", "a03 * b20 + a13 * b21", "a03 * b30 + a13 * b31"
+
+.. c:function:: void glm_mat2x4_mulv(mat2x4 m, vec2 v, vec4 dest)
+
+    multiply mat2x4 with vec2 (column vector) and store in dest column vector
 
     Parameters:
       | *[in]*  **m**     mat2x4 (left)
-      | *[in]*  **v**     vec4 (right, column vector)
+      | *[in]*  **v**     vec2 (right, column vector)
       | *[out]* **dest**  destination (result, column vector)
 
+    .. csv-table:: mat2x4
+        :header: "", "column 1", "column 2"
+
+        "row 1", "m00", "m10"
+        "row 2", "m01", "m11"
+        "row 3", "m02", "m12"
+        "row 4", "m03", "m13"
+
+    .. csv-table:: column vec2 (1x2)
+        :header: "", "column 1"
+
+        "row 1", "v0"
+        "row 2", "v1"
+
+    .. csv-table:: column vec4 (1x4)
+        :header: "", "column 1"
+
+        "row 1", "m00 * v0 + m10 * v1"
+        "row 2", "m01 * v0 + m11 * v1"
+        "row 3", "m02 * v0 + m12 * v1"
+        "row 4", "m03 * v0 + m13 * v1"
+
 .. c:function:: void glm_mat2x4_transpose(mat2x4 m, mat4x2 dest)
 
     transpose matrix and store in dest

docs/source/mat3.rst

@@ -190,7 +190,7 @@ Functions documentation
     Returns:
         scalar value e.g. Matrix1x1
 
-.. c:function:: void glm_mat3_make(float * __restrict src, mat3 dest)
+.. c:function:: void glm_mat3_make(const float * __restrict src, mat3 dest)
 
     Create mat3 matrix from pointer
 

docs/source/mat3x2.rst

@@ -23,6 +23,15 @@ Functions:
 #. :c:func:`glm_mat3x2_transpose`
 #. :c:func:`glm_mat3x2_scale`
 
+Represented
+~~~~~~~~~~~
+
+.. csv-table:: mat3x2
+   :header: "", "column 1", "column 2", "column 3"
+
+   "row 1", "m00", "m10", "m20"
+   "row 2", "m01", "m11", "m21"
+
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -41,7 +50,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat3x2_make(float * __restrict src, mat3x2 dest)
+.. c:function:: void glm_mat3x2_make(const float * __restrict src, mat3x2 dest)
 
     Create mat3x2 matrix from pointer
 
@@ -51,28 +60,66 @@ Functions documentation
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination matrix3x2
 
-.. c:function:: void glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest)
+.. c:function:: void glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat2 dest)
 
     multiply m1 and m2 to dest
 
     .. code-block:: c
 
-       glm_mat3x2_mul(mat3x2, mat2x3, mat3);
+       glm_mat3x2_mul(mat3x2, mat2x3, mat2);
 
     Parameters:
       | *[in]*  **m1**    left matrix (mat3x2)
       | *[in]*  **m2**    right matrix (mat2x3)
-      | *[out]* **dest**  destination matrix (mat3)
+      | *[out]* **dest**  destination matrix (mat2)
 
-.. c:function:: void glm_mat3x2_mulv(mat3x2 m, vec2 v, vec3 dest)
+    .. csv-table:: mat3x2
+        :header: "", "column 1", "column 2", "column 3"
 
-    multiply mat3x2 with vec2 (column vector) and store in dest vector
+        "row 1", "a00", "a10", "a20"
+        "row 2", "a01", "a11", "a21"
+
+    .. csv-table:: mat2x3
+        :header: "", "column 1", "column 2"
+
+        "row 1", "b00", "b10"
+        "row 2", "b01", "b11"
+        "row 3", "b02", "b12"
+
+    .. csv-table:: mat2x2
+        :header: "", "column 1", "column 2"
+
+        "row 1", "a00 * b00 + a10 * b01 + a20 * b02", "a00 * b10 + a10 * b11 + a20 * b12"
+        "row 2", "a01 * b00 + a11 * b01 + a21 * b02", "a01 * b10 + a11 * b11 + a21 * b12"
+
+.. c:function:: void glm_mat3x2_mulv(mat3x2 m, vec3 v, vec2 dest)
+
+    multiply mat3x2 with vec3 (column vector) and store in dest vector
 
     Parameters:
       | *[in]*  **m**     mat3x2 (left)
       | *[in]*  **v**     vec3 (right, column vector)
       | *[out]* **dest**  destination (result, column vector)
 
+    .. csv-table:: mat3x2
+        :header: "", "column 1", "column 2", "column 3"
+
+        "row 1", "m00", "m10", "m20"
+        "row 2", "m01", "m11", "m21"
+
+    .. csv-table:: column vec3 (1x3)
+        :header: "", "column 1"
+
+        "row 1", "v0"
+        "row 2", "v1"
+        "row 3", "v2"
+
+    .. csv-table:: column vec2 (1x2)
+        :header: "", "column 1"
+
+        "row 1", "m00 * v0 + m10 * v1 + m20 * v2"
+        "row 2", "m01 * v0 + m11 * v1 + m21 * v2"
+
 .. c:function:: void glm_mat3x2_transpose(mat3x2 m, mat2x3 dest)
 
     transpose matrix and store in dest

docs/source/mat3x4.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat3x4_make(float * __restrict src, mat3x4 dest)
+.. c:function:: void glm_mat3x4_make(const float * __restrict src, mat3x4 dest)
 
     Create mat3x4 matrix from pointer
 

docs/source/mat4.rst

@@ -304,7 +304,7 @@ Functions documentation
     Returns:
         scalar value e.g. Matrix1x1
 
-.. c:function:: void glm_mat4_make(float * __restrict src, mat4 dest)
+.. c:function:: void glm_mat4_make(const float * __restrict src, mat4 dest)
 
     Create mat4 matrix from pointer
 

docs/source/mat4x2.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat4x2_make(float * __restrict src, mat4x2 dest)
+.. c:function:: void glm_mat4x2_make(const float * __restrict src, mat4x2 dest)
 
     Create mat4x2 matrix from pointer
 

docs/source/mat4x3.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat4x3_make(float * __restrict src, mat4x3 dest)
+.. c:function:: void glm_mat4x3_make(const float * __restrict src, mat4x3 dest)
 
     Create mat4x3 matrix from pointer
 

docs/source/opt.rst

@@ -5,35 +5,33 @@
 
 A few options are provided via macros.
 
+❗️ IMPORTANT ❗️
+
+It's a good idea to set up your config macros in build settings like CMake, Xcode, or Visual Studio. This is especially important if you're using features like Modules in Xcode, where adding macros directly before the **cglm** headers might not work.
+
 Alignment Option
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-As default, cglm requires types to be aligned. Alignment requirements:
+By default, **cglm** requires types to be aligned with specific byte requirements:
 
-vec3:   8 byte
+- vec3: 8 bytes
-vec4:   16 byte
+- vec4: 16 bytes
-mat4:   16 byte
+- mat4: 16 bytes (32 on AVX)
-versor: 16 byte
+- versor: 16 bytes
 
-By starting **v0.4.5** cglm provides an option to disable alignment requirement.
+Starting with **v0.4.5**, **cglm** offers an option to relax these alignment requirements. To use this option, define the **CGLM_ALL_UNALIGNED** macro before including any headers. This definition can be made within Xcode, Visual Studio, other IDEs, or directly in your build system. If using pre-compiled versions of **cglm**, you'll need to compile them with the **CGLM_ALL_UNALIGNED** macro.
-To enable this option define **CGLM_ALL_UNALIGNED** macro before all headers.
-You can define it in Xcode, Visual Studio (or other IDEs) or you can also prefer
-to define it in build system. If you use pre-compiled versions then you
-have to compile cglm with **CGLM_ALL_UNALIGNED** macro.
 
-**VERY VERY IMPORTANT:** If you use cglm in multiple projects and
+**❗️NOTE:❗️** If you're using **cglm** across multiple interdependent projects:
- those projects are depends on each other, then
 
- | *ALWAYS* or *NEVER USE* **CGLM_ALL_UNALIGNED** macro in linked projects
+- Always or never use the **CGLM_ALL_UNALIGNED** macro in all linked projects to avoid configuration conflicts. A **cglm** header from one project could require alignment, while a header from another might not, leading to **cglm** functions accessing invalid memory locations.
 
- if you do not know what you are doing. Because a cglm header included
+- **Key Point:** Maintain the same **cglm** configuration across all your projects. For example, if you activate **CGLM_ALL_UNALIGNED** in one project, ensure it's set in the others too.
- via 'project A' may force types to be aligned and another cglm header
- included via 'project B' may not require alignment. In this case
- cglm functions will read from and write to **INVALID MEMORY LOCATIONSNs**.
 
- ALWAYS USE SAME CONFIGURATION / OPTION for **cglm** if you have multiple projects.
+**❗️NOTE:❗️**
 
- For instance if you set CGLM_ALL_UNALIGNED in a project then set it in other projects too
+While **CGLM_ALL_UNALIGNED** allows for flexibility in alignment, it doesn't override C's fundamental alignment rules. For example, an array like *vec4* decays to a pointer (float*) in functions, which must adhere to the alignment requirements of a float pointer (4 bytes). This adherence is crucial because **cglm** directly dereferences these pointers instead of copying data, and failing to meet alignment requirements can lead to unpredictable errors, such as crashes.
+
+You can use `CGLM_ALIGN` and `CGLM_ALIGN_MAT` macros for aligning local variables or struct members. However, when dealing with dynamic memory allocation or custom memory locations, you'll need to ensure alignment requirements are met appropriately for those cases
 
 Clipspace Option[s]
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

docs/source/quat.rst

@@ -422,7 +422,7 @@ Functions documentation
       | *[in]*       **q**      quaternion
       | *[in]*       **pivot**  pivot
 
-.. c:function:: void glm_quat_make(float * __restrict src, versor dest)
+.. c:function:: void glm_quat_make(const float * __restrict src, versor dest)
 
     Create quaternion from pointer
 

docs/source/ray.rst

@@ -13,6 +13,8 @@ Table of contents (click to go):
 Functions:
 
 1. :c:func:`glm_ray_triangle`
+#. :c:func:`glm_ray_sphere`
+#. :c:func:`glm_ray_at`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -29,3 +31,39 @@ Functions documentation
       | *[in]*       **v2**            third vertex of triangle
       | *[in, out]*  **d**             float pointer to save distance to intersection
       | *[out]*      **intersection**  whether there is intersection
+
+.. c:function:: bool glm_ray_sphere(vec3 origin, vec3 dir, vec4 s, float * __restrict t1, float * __restrict t2)
+
+    ray sphere intersection
+
+    returns false if there is no intersection if true:
+
+    - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+    - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+    - t1 < 0, t2 < 0: no intersection ahead of the ray ( returns false )
+    - the caller can check if the intersection points (t1 and t2) fall within a
+      specific range (for example, tmin < t1, t2 < tmax) to determine if the
+      intersections are within a desired segment of the ray
+
+    Parameters:
+      | *[in]*   **origin**  ray origin
+      | *[in]*   **dir**     normalized ray direction
+      | *[in]*   **s**       sphere  [center.x, center.y, center.z, radii]
+      | *[out]*  **t1**      near point1 (closer to origin)
+      | *[out]*  **t2**      far point2 (farther from origin)
+
+    Return:
+      | whether there is intersection
+
+.. c:function:: bool  glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point)
+
+    point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+
+    Parameters:
+      | *[in]*   **origin**  ray origin
+      | *[in]*   **dir**     ray direction
+      | *[out]*  **t**       parameter
+      | *[out]*  **point**   point at t
+
+    Return:
+      | point at t

docs/source/vec2.rst

@@ -53,6 +53,8 @@ Functions:
 #. :c:func:`glm_vec2_clamp`
 #. :c:func:`glm_vec2_lerp`
 #. :c:func:`glm_vec2_make`
+#. :c:func:`glm_vec2_reflect`
+#. :c:func:`glm_vec2_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -385,7 +387,7 @@ Functions documentation
       | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
       | *[out]* **dest**   destination
 
-.. c:function:: void glm_vec2_make(float * __restrict src, vec2 dest)
+.. c:function:: void glm_vec2_make(const float * __restrict src, vec2 dest)
 
     Create two dimensional vector from pointer
 
@@ -394,3 +396,29 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec2_reflect(vec2 v, vec2 n, vec2 dest)
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **v**     incident vector
+      | *[in]*  **n**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: bool glm_vec2_refract(vec2 v, vec2 n, float eta, vec2 dest)
+
+    Computes refraction vector for an incident vector and a surface normal.
+   
+    Calculates the refraction vector based on Snell's law. If total internal reflection
+    occurs (angle too great given eta), dest is set to zero and returns false.
+    Otherwise, computes refraction vector, stores it in dest, and returns true.
+
+    Parameters:
+      | *[in]*  **v**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **n**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction (incident/transmitted)
+      | *[out]* **dest**  refraction vector if refraction occurs; zero vector otherwise
+
+    Returns:
+      returns true if refraction occurs; false if total internal reflection occurs.

docs/source/vec3.rst

@@ -80,6 +80,9 @@ Functions:
 #. :c:func:`glm_vec3_clamp`
 #. :c:func:`glm_vec3_lerp`
 #. :c:func:`glm_vec3_make`
+#. :c:func:`glm_vec3_faceforward`
+#. :c:func:`glm_vec3_reflect`
+#. :c:func:`glm_vec3_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -503,7 +506,7 @@ Functions documentation
       | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
       | *[out]* **dest**   destination
 
-.. c:function:: void glm_vec3_make(float * __restrict src, vec3 dest)
+.. c:function:: void glm_vec3_make(const float * __restrict src, vec3 dest)
 
     Create three dimensional vector from pointer
 
@@ -512,3 +515,40 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec3_faceforward(vec3 n, vec3 v, vec3 nref, vec3 dest)
+
+    A vector pointing in the same direction as another
+
+    Parameters:
+      | *[in]*  **n**     vector to orient
+      | *[in]*  **v**     incident vector
+      | *[in]*  **nref**  reference vector
+      | *[out]* **dest**  destination: oriented vector, pointing away from the surface.
+
+.. c:function:: void glm_vec3_reflect(vec3 v, vec3 n, vec3 dest)
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **v**     incident vector
+      | *[in]*  **n**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: bool glm_vec3_refract(vec3 v, vec3 n, float eta, vec3 dest)
+
+    
+    Computes refraction vector for an incident vector and a surface normal.
+   
+    Calculates the refraction vector based on Snell's law. If total internal reflection
+    occurs (angle too great given eta), dest is set to zero and returns false.
+    Otherwise, computes refraction vector, stores it in dest, and returns true.
+
+    Parameters:
+      | *[in]*  **v**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **n**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction (incident/transmitted)
+      | *[out]* **dest**  refraction vector if refraction occurs; zero vector otherwise
+
+    Returns:
+      returns true if refraction occurs; false if total internal reflection occurs.

docs/source/vec4.rst

@@ -60,6 +60,8 @@ Functions:
 #. :c:func:`glm_vec4_lerp`
 #. :c:func:`glm_vec4_cubic`
 #. :c:func:`glm_vec4_make`
+#. :c:func:`glm_vec4_reflect`
+#. :c:func:`glm_vec4_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -415,7 +417,7 @@ Functions documentation
       | *[in]*  **s**      parameter
       | *[out]* **dest**   destination
 
-.. c:function:: void glm_vec4_make(float * __restrict src, vec4 dest)
+.. c:function:: void glm_vec4_make(const float * __restrict src, vec4 dest)
 
     Create four dimensional vector from pointer
 
@@ -424,3 +426,33 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: bool glm_vec4_reflect(vec4 v, vec4 n, vec4 dest)
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **v**     incident vector
+      | *[in]*  **n**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: bool glm_vec4_refract(vec4 v, vec4 n, float eta, vec4 dest)
+
+    computes refraction vector for an incident vector and a surface normal.
+   
+    Calculates the refraction vector based on Snell's law. If total internal reflection
+    occurs (angle too great given eta), dest is set to zero and returns false.
+    Otherwise, computes refraction vector, stores it in dest, and returns true.
+   
+    This implementation does not explicitly preserve the 'w' component of the
+    incident vector 'I' in the output 'dest', users requiring the preservation of
+    the 'w' component should manually adjust 'dest' after calling this function.
+
+    Parameters:
+      | *[in]*  **v**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **n**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction (incident/transmitted)
+      | *[out]* **dest**  refraction vector if refraction occurs; zero vector otherwise
+
+    Returns:
+      returns true if refraction occurs; false if total internal reflection occurs.

include/cglm/affine-pre.h

@@ -207,7 +207,25 @@ glm_rotate_z(mat4 m, float angle, mat4 dest) {
 }
 
 /*!
- * @brief rotate existing transform matrix around given axis by angle
+ * @brief rotate existing transform matrix 
+ *        around given axis by angle at ORIGIN (0,0,0)
+ *
+ *   **❗️IMPORTANT ❗️**
+ *
+ *   If you need to rotate object around itself e.g. center of object or at
+ *   some point [of object] then `glm_rotate_at()` would be better choice to do so.
+ *
+ *   Even if object's model transform is identiy, rotation may not be around
+ *   center of object if object does not lay out at ORIGIN perfectly.
+ *
+ *   Using `glm_rotate_at()` with center of bounding shape ( AABB, Sphere ... )
+ *   would be an easy option to rotate around object if object is not at origin.
+ *
+ *   One another option to rotate around itself at any point is `glm_spin()`
+ *   which is perfect if only rotating around model position is desired e.g. not
+ *   specific point on model for instance center of geometry or center of mass,
+ *   again if geometry is not perfectly centered at origin at identity transform,
+ *   rotation may not be around geometry.
  *
  * @param[in, out]  m      affine transform
  * @param[in]       angle  angle (radians)
@@ -268,7 +286,8 @@ glm_rotate_atm(mat4 m, vec3 pivot, float angle, vec3 axis) {
 }
 
 /*!
- * @brief rotate existing transform matrix around given axis by angle around self (doesn't affected by position)
+ * @brief rotate existing transform matrix 
+ *        around given axis by angle around self (doesn't affected by position)
  *
  * @param[in, out]  m      affine transform
  * @param[in]       angle  angle (radians)

include/cglm/call/mat2.h

@@ -75,7 +75,7 @@ glmc_mat2_rmc(vec2 r, mat2 m, vec2 c);
 
 CGLM_EXPORT
 void
-glmc_mat2_make(float * __restrict src, mat2 dest);
+glmc_mat2_make(const float * __restrict src, mat2 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/mat2x3.h

@@ -23,15 +23,15 @@ glmc_mat2x3_zero(mat2x3 mat);
 
 CGLM_EXPORT
 void
-glmc_mat2x3_make(float * __restrict src, mat2x3 dest);
+glmc_mat2x3_make(const float * __restrict src, mat2x3 dest);
 
 CGLM_EXPORT
 void
-glmc_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest);
+glmc_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat3 dest);
 
 CGLM_EXPORT
 void
-glmc_mat2x3_mulv(mat2x3 m, vec3 v, vec2 dest);
+glmc_mat2x3_mulv(mat2x3 m, vec2 v, vec3 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat2x4.h

@@ -23,15 +23,15 @@ glmc_mat2x4_zero(mat2x4 mat);
 
 CGLM_EXPORT
 void
-glmc_mat2x4_make(float * __restrict src, mat2x4 dest);
+glmc_mat2x4_make(const float * __restrict src, mat2x4 dest);
 
 CGLM_EXPORT
 void
-glmc_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest);
+glmc_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat4 dest);
 
 CGLM_EXPORT
 void
-glmc_mat2x4_mulv(mat2x4 m, vec4 v, vec2 dest);
+glmc_mat2x4_mulv(mat2x4 m, vec2 v, vec4 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat3.h

@@ -82,7 +82,7 @@ glmc_mat3_rmc(vec3 r, mat3 m, vec3 c);
 
 CGLM_EXPORT
 void
-glmc_mat3_make(float * __restrict src, mat3 dest);
+glmc_mat3_make(const float * __restrict src, mat3 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/mat3x2.h

@@ -23,15 +23,15 @@ glmc_mat3x2_zero(mat3x2 mat);
 
 CGLM_EXPORT
 void
-glmc_mat3x2_make(float * __restrict src, mat3x2 dest);
+glmc_mat3x2_make(const float * __restrict src, mat3x2 dest);
 
 CGLM_EXPORT
 void
-glmc_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest);
+glmc_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat2 dest);
 
 CGLM_EXPORT
 void
-glmc_mat3x2_mulv(mat3x2 m, vec2 v, vec3 dest);
+glmc_mat3x2_mulv(mat3x2 m, vec3 v, vec2 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat3x4.h

@@ -23,7 +23,7 @@ glmc_mat3x4_zero(mat3x4 mat);
 
 CGLM_EXPORT
 void
-glmc_mat3x4_make(float * __restrict src, mat3x4 dest);
+glmc_mat3x4_make(const float * __restrict src, mat3x4 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat4.h

@@ -123,7 +123,7 @@ glmc_mat4_rmc(vec4 r, mat4 m, vec4 c);
 
 CGLM_EXPORT
 void
-glmc_mat4_make(float * __restrict src, mat4 dest);
+glmc_mat4_make(const float * __restrict src, mat4 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/mat4x2.h

@@ -23,7 +23,7 @@ glmc_mat4x2_zero(mat4x2 mat);
 
 CGLM_EXPORT
 void
-glmc_mat4x2_make(float * __restrict src, mat4x2 dest);
+glmc_mat4x2_make(const float * __restrict src, mat4x2 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat4x3.h

@@ -23,7 +23,7 @@ glmc_mat4x3_zero(mat4x3 mat);
 
 CGLM_EXPORT
 void
-glmc_mat4x3_make(float * __restrict src, mat4x3 dest);
+glmc_mat4x3_make(const float * __restrict src, mat4x3 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/quat.h

@@ -163,7 +163,7 @@ glmc_quat_rotate_atm(mat4 m, versor q, vec3 pivot);
 
 CGLM_EXPORT
 void
-glmc_quat_make(float * __restrict src, versor dest);
+glmc_quat_make(const float * __restrict src, versor dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/ray.h

@@ -21,6 +21,18 @@ glmc_ray_triangle(vec3   origin,
                   vec3   v2,
                   float *d);
 
+CGLM_EXPORT
+bool
+glmc_ray_sphere(vec3 origin,
+                vec3 dir,
+                vec4 s,
+                float * __restrict t1,
+                float * __restrict t2);
+
+CGLM_EXPORT
+void
+glmc_ray_at(vec3 orig, vec3 dir, float t, vec3 point);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/call/vec2.h

@@ -195,7 +195,15 @@ glmc_vec2_complex_conjugate(vec2 a, vec2 dest);
 
 CGLM_EXPORT
 void
-glmc_vec2_make(float * __restrict src, vec2 dest);
+glmc_vec2_make(const float * __restrict src, vec2 dest);
+
+CGLM_EXPORT
+void
+glmc_vec2_reflect(vec2 v, vec2 n, vec2 dest);
+
+CGLM_EXPORT
+bool
+glmc_vec2_refract(vec2 v, vec2 n, float eta, vec2 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/vec3.h

@@ -332,7 +332,19 @@ glmc_vec3_sqrt(vec3 v, vec3 dest);
 
 CGLM_EXPORT
 void
-glmc_vec3_make(float * __restrict src, vec3 dest);
+glmc_vec3_make(const float * __restrict src, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_faceforward(vec3 n, vec3 v, vec3 nref, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_reflect(vec3 v, vec3 n, vec3 dest);
+
+CGLM_EXPORT
+bool
+glmc_vec3_refract(vec3 v, vec3 n, float eta, vec3 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/vec4.h

@@ -309,7 +309,15 @@ glmc_vec4_sqrt(vec4 v, vec4 dest);
 
 CGLM_EXPORT
 void
-glmc_vec4_make(float * __restrict src, vec4 dest);
+glmc_vec4_make(const float * __restrict src, vec4 dest);
+
+CGLM_EXPORT
+void
+glmc_vec4_reflect(vec4 v, vec4 n, vec4 dest);
+
+CGLM_EXPORT
+bool
+glmc_vec4_refract(vec4 v, vec4 n, float eta, vec4 dest);
 
 #ifdef __cplusplus
 }

include/cglm/common.h

@@ -37,6 +37,18 @@
 #  define CGLM_INLINE static inline __attribute((always_inline))
 #endif
 
+#if defined(__GNUC__) || defined(__clang__)
+#  define CGLM_UNLIKELY(expr) __builtin_expect(!!(expr), 0)
+#  define CGLM_LIKELY(expr)   __builtin_expect(!!(expr), 1)
+#else
+#  define CGLM_UNLIKELY(expr) (expr)
+#  define CGLM_LIKELY(expr)   (expr)
+#endif
+
+#if defined(_M_FP_FAST) || defined(__FAST_MATH__)
+#  define CGLM_FAST_MATH
+#endif
+
 #define GLM_SHUFFLE4(z, y, x, w) (((z) << 6) | ((y) << 4) | ((x) << 2) | (w))
 #define GLM_SHUFFLE3(z, y, x)    (((z) << 4) | ((y) << 2) | (x))
 

include/cglm/mat2.h

@@ -354,7 +354,7 @@ glm_mat2_rmc(vec2 r, mat2 m, vec2 c) {
  */
 CGLM_INLINE
 void
-glm_mat2_make(float * __restrict src, mat2 dest) {
+glm_mat2_make(const float * __restrict src, mat2 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[1][0] = src[2];

include/cglm/mat2x3.h

@@ -13,9 +13,9 @@
  Functions:
    CGLM_INLINE void glm_mat2x3_copy(mat2x3 mat, mat2x3 dest);
    CGLM_INLINE void glm_mat2x3_zero(mat2x3 mat);
-   CGLM_INLINE void glm_mat2x3_make(float * __restrict src, mat2x3 dest);
+   CGLM_INLINE void glm_mat2x3_make(const float * __restrict src, mat2x3 dest);
-   CGLM_INLINE void glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest);
+   CGLM_INLINE void glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat3 dest);
-   CGLM_INLINE void glm_mat2x3_mulv(mat2x3 m, vec3 v, vec2 dest);
+   CGLM_INLINE void glm_mat2x3_mulv(mat2x3 m, vec2 v, vec3 dest);
    CGLM_INLINE void glm_mat2x3_transpose(mat2x3 m, mat3x2 dest);
    CGLM_INLINE void glm_mat2x3_scale(mat2x3 m, float s);
  */
@@ -68,7 +68,7 @@ glm_mat2x3_zero(mat2x3 mat) {
  */
 CGLM_INLINE
 void
-glm_mat2x3_make(float * __restrict src, mat2x3 dest) {
+glm_mat2x3_make(const float * __restrict src, mat2x3 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];
@@ -82,7 +82,7 @@ glm_mat2x3_make(float * __restrict src, mat2x3 dest) {
  * @brief multiply m1 and m2 to dest
  *
  * @code
- * glm_mat2x3_mul(mat2x3, mat3x2, mat2);
+ * glm_mat2x3_mul(mat2x3, mat3x2, mat3);
  * @endcode
  *
  * @param[in]  m1   left matrix (mat2x3)
@@ -91,7 +91,7 @@ glm_mat2x3_make(float * __restrict src, mat2x3 dest) {
  */
 CGLM_INLINE
 void
-glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest) {
+glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat3 dest) {
   float a00 = m1[0][0], a01 = m1[0][1], a02 = m1[0][2], 
         a10 = m1[1][0], a11 = m1[1][1], a12 = m1[1][2],
 
@@ -99,11 +99,17 @@ glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest) {
         b10 = m2[1][0], b11 = m2[1][1],
         b20 = m2[2][0], b21 = m2[2][1];
 
-  dest[0][0] = a00 * b00 + a01 * b10 + a02 * b20;
+  dest[0][0] = a00 * b00 + a10 * b01;
-  dest[0][1] = a00 * b01 + a01 * b11 + a02 * b21;
+  dest[0][1] = a01 * b00 + a11 * b01;
+  dest[0][2] = a02 * b00 + a12 * b01;
 
-  dest[1][0] = a10 * b00 + a11 * b10 + a12 * b20;
+  dest[1][0] = a00 * b10 + a10 * b11;
-  dest[1][1] = a10 * b01 + a11 * b11 + a12 * b21;
+  dest[1][1] = a01 * b10 + a11 * b11;
+  dest[1][2] = a02 * b10 + a12 * b11;
+
+  dest[2][0] = a00 * b20 + a10 * b21;
+  dest[2][1] = a01 * b20 + a11 * b21;
+  dest[2][2] = a02 * b20 + a12 * b21;
 }
 
 /*!
@@ -115,11 +121,12 @@ glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest) {
  */
 CGLM_INLINE
 void
-glm_mat2x3_mulv(mat2x3 m, vec3 v, vec2 dest) {
+glm_mat2x3_mulv(mat2x3 m, vec2 v, vec3 dest) {
-  float v0 = v[0], v1 = v[1], v2 = v[2];
+  float v0 = v[0], v1 = v[1];
 
-  dest[0] = m[0][0] * v0 + m[0][1] * v1 + m[0][2] * v2;
+  dest[0] = m[0][0] * v0 + m[1][0] * v1;
-  dest[1] = m[1][0] * v0 + m[1][1] * v1 + m[1][2] * v2;
+  dest[1] = m[0][1] * v0 + m[1][1] * v1;
+  dest[2] = m[0][2] * v0 + m[1][2] * v1;
 }
 
 /*!

include/cglm/mat2x4.h

@@ -13,9 +13,9 @@
  Functions:
    CGLM_INLINE void glm_mat2x4_copy(mat2x4 mat, mat2x4 dest);
    CGLM_INLINE void glm_mat2x4_zero(mat2x4 mat);
-   CGLM_INLINE void glm_mat2x4_make(float * __restrict src, mat2x4 dest);
+   CGLM_INLINE void glm_mat2x4_make(const float * __restrict src, mat2x4 dest);
-   CGLM_INLINE void glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest);
+   CGLM_INLINE void glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat4 dest);
-   CGLM_INLINE void glm_mat2x4_mulv(mat2x4 m, vec4 v, vec2 dest);
+   CGLM_INLINE void glm_mat2x4_mulv(mat2x4 m, vec2 v, vec4 dest);
    CGLM_INLINE void glm_mat2x4_transpose(mat2x4 m, mat4x2 dest);
    CGLM_INLINE void glm_mat2x4_scale(mat2x4 m, float s);
  */
@@ -64,7 +64,7 @@ glm_mat2x4_zero(mat2x4 mat) {
  */
 CGLM_INLINE
 void
-glm_mat2x4_make(float * __restrict src, mat2x4 dest) {
+glm_mat2x4_make(const float * __restrict src, mat2x4 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];
@@ -80,16 +80,16 @@ glm_mat2x4_make(float * __restrict src, mat2x4 dest) {
  * @brief multiply m1 and m2 to dest
  *
  * @code
- * glm_mat2x4_mul(mat2x4, mat4x2, mat2);
+ * glm_mat2x4_mul(mat2x4, mat4x2, mat4);
  * @endcode
  *
  * @param[in]  m1   left matrix (mat2x4)
  * @param[in]  m2   right matrix (mat4x2)
- * @param[out] dest destination matrix (mat2)
+ * @param[out] dest destination matrix (mat4)
  */
 CGLM_INLINE
 void
-glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest) {
+glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat4 dest) {
   float a00 = m1[0][0], a01 = m1[0][1], a02 = m1[0][2], a03 = m1[0][3],
         a10 = m1[1][0], a11 = m1[1][1], a12 = m1[1][2], a13 = m1[1][3],
 
@@ -98,15 +98,29 @@ glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest) {
         b20 = m2[2][0], b21 = m2[2][1],
         b30 = m2[3][0], b31 = m2[3][1];
 
-  dest[0][0] = a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30;
+  dest[0][0] = a00 * b00 + a10 * b01;
-  dest[0][1] = a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31;
+  dest[0][1] = a01 * b00 + a11 * b01;
+  dest[0][2] = a02 * b00 + a12 * b01;
+  dest[0][3] = a03 * b00 + a13 * b01;
 
-  dest[1][0] = a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30;
+  dest[1][0] = a00 * b10 + a10 * b11;
-  dest[1][1] = a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31;
+  dest[1][1] = a01 * b10 + a11 * b11;
+  dest[1][2] = a02 * b10 + a12 * b11;
+  dest[1][3] = a03 * b10 + a13 * b11;
+
+  dest[2][0] = a00 * b20 + a10 * b21;
+  dest[2][1] = a01 * b20 + a11 * b21;
+  dest[2][2] = a02 * b20 + a12 * b21;
+  dest[2][3] = a03 * b20 + a13 * b21;
+
+  dest[3][0] = a00 * b30 + a10 * b31;
+  dest[3][1] = a01 * b30 + a11 * b31;
+  dest[3][2] = a02 * b30 + a12 * b31;
+  dest[3][3] = a03 * b30 + a13 * b31;
 }
 
 /*!
- * @brief multiply matrix with column vector and store in dest vector
+ * @brief multiply matrix with column vector and store in dest column vector
  *
  * @param[in]  m    matrix (left)
  * @param[in]  v    vector (right, column vector)
@@ -114,11 +128,13 @@ glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest) {
  */
 CGLM_INLINE
 void
-glm_mat2x4_mulv(mat2x4 m, vec4 v, vec2 dest) {
+glm_mat2x4_mulv(mat2x4 m, vec2 v, vec4 dest) {
-  float v0 = v[0], v1 = v[1], v2 = v[2], v3 = v[3];
+  float v0 = v[0], v1 = v[1];
 
-  dest[0] = m[0][0] * v0 + m[0][1] * v1 + m[0][2] * v2 + m[0][3] * v3;
+  dest[0] = m[0][0] * v0 + m[1][0] * v1;
-  dest[1] = m[1][0] * v0 + m[1][1] * v1 + m[1][2] * v2 + m[1][3] * v3;
+  dest[1] = m[0][1] * v0 + m[1][1] * v1;
+  dest[2] = m[0][2] * v0 + m[1][2] * v1;
+  dest[3] = m[0][3] * v0 + m[1][3] * v1;
 }
 
 /*!

include/cglm/mat3.h

@@ -436,7 +436,7 @@ glm_mat3_rmc(vec3 r, mat3 m, vec3 c) {
  */
 CGLM_INLINE
 void
-glm_mat3_make(float * __restrict src, mat3 dest) {
+glm_mat3_make(const float * __restrict src, mat3 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/mat3x2.h

@@ -13,9 +13,9 @@
  Functions:
    CGLM_INLINE void glm_mat3x2_copy(mat3x2 mat, mat3x2 dest);
    CGLM_INLINE void glm_mat3x2_zero(mat3x2 mat);
-   CGLM_INLINE void glm_mat3x2_make(float * __restrict src, mat3x2 dest);
+   CGLM_INLINE void glm_mat3x2_make(const float * __restrict src, mat3x2 dest);
-   CGLM_INLINE void glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest);
+   CGLM_INLINE void glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat2 dest);
-   CGLM_INLINE void glm_mat3x2_mulv(mat3x2 m, vec2 v, vec3 dest);
+   CGLM_INLINE void glm_mat3x2_mulv(mat3x2 m, vec3 v, vec2 dest);
    CGLM_INLINE void glm_mat3x2_transpose(mat3x2 m, mat2x3 dest);
    CGLM_INLINE void glm_mat3x2_scale(mat3x2 m, float s);
  */
@@ -69,7 +69,7 @@ glm_mat3x2_zero(mat3x2 mat) {
  */
 CGLM_INLINE
 void
-glm_mat3x2_make(float * __restrict src, mat3x2 dest) {
+glm_mat3x2_make(const float * __restrict src, mat3x2 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
 
@@ -84,16 +84,16 @@ glm_mat3x2_make(float * __restrict src, mat3x2 dest) {
  * @brief multiply m1 and m2 to dest
  *
  * @code
- * glm_mat3x2_mul(mat3x2, mat2x3, mat3);
+ * glm_mat3x2_mul(mat3x2, mat2x3, mat2);
  * @endcode
  *
  * @param[in]  m1   left matrix (mat3x2)
  * @param[in]  m2   right matrix (mat2x3)
- * @param[out] dest destination matrix (mat3)
+ * @param[out] dest destination matrix (mat2)
  */
 CGLM_INLINE
 void
-glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest) {
+glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat2 dest) {
   float a00 = m1[0][0], a01 = m1[0][1],
         a10 = m1[1][0], a11 = m1[1][1],
         a20 = m1[2][0], a21 = m1[2][1],
@@ -101,21 +101,15 @@ glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest) {
         b00 = m2[0][0], b01 = m2[0][1], b02 = m2[0][2],
         b10 = m2[1][0], b11 = m2[1][1], b12 = m2[1][2];
 
-  dest[0][0] = a00 * b00 + a01 * b10;
+  dest[0][0] = a00 * b00 + a10 * b01 + a20 * b02;
-  dest[0][1] = a00 * b01 + a01 * b11;
+  dest[0][1] = a01 * b00 + a11 * b01 + a21 * b02;
-  dest[0][2] = a00 * b02 + a01 * b12;
 
-  dest[1][0] = a10 * b00 + a11 * b10;
+  dest[1][0] = a00 * b10 + a10 * b11 + a20 * b12;
-  dest[1][1] = a10 * b01 + a11 * b11;
+  dest[1][1] = a01 * b10 + a11 * b11 + a21 * b12;
-  dest[1][2] = a10 * b02 + a11 * b12;
-
-  dest[2][0] = a20 * b00 + a21 * b10;
-  dest[2][1] = a20 * b01 + a21 * b11;
-  dest[2][2] = a20 * b02 + a21 * b12;
 }
 
 /*!
- * @brief multiply matrix with column vector and store in dest vector
+ * @brief multiply matrix with column vector and store in dest column vector
  *
  * @param[in]  m    matrix (left)
  * @param[in]  v    vector (right, column vector)
@@ -123,12 +117,11 @@ glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest) {
  */
 CGLM_INLINE
 void
-glm_mat3x2_mulv(mat3x2 m, vec2 v, vec3 dest) {
+glm_mat3x2_mulv(mat3x2 m, vec3 v, vec2 dest) {
-  float v0 = v[0], v1 = v[1];
+  float v0 = v[0], v1 = v[1], v2 = v[2];
 
-  dest[0] = m[0][0] * v0 + m[0][1] * v1;
+  dest[0] = m[0][0] * v0 + m[1][0] * v1 + m[2][0] * v2;
-  dest[1] = m[1][0] * v0 + m[1][1] * v1;
+  dest[1] = m[0][1] * v0 + m[1][1] * v1 + m[2][1] * v2;
-  dest[2] = m[2][0] * v0 + m[2][1] * v1;
 }
 
 /*!

include/cglm/mat3x4.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat3x4_copy(mat3x4 mat, mat3x4 dest);
    CGLM_INLINE void glm_mat3x4_zero(mat3x4 mat);
-   CGLM_INLINE void glm_mat3x4_make(float * __restrict src, mat3x4 dest);
+   CGLM_INLINE void glm_mat3x4_make(const float * __restrict src, mat3x4 dest);
    CGLM_INLINE void glm_mat3x4_mul(mat3x4 m1, mat4x3 m2, mat3 dest);
    CGLM_INLINE void glm_mat3x4_mulv(mat3x4 m, vec4 v, vec3 dest);
    CGLM_INLINE void glm_mat3x4_transpose(mat3x4 m, mat4x3 dest);
@@ -66,7 +66,7 @@ glm_mat3x4_zero(mat3x4 mat) {
  */
 CGLM_INLINE
 void
-glm_mat3x4_make(float * __restrict src, mat3x4 dest) {
+glm_mat3x4_make(const float * __restrict src, mat3x4 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/mat4.h

@@ -790,7 +790,7 @@ glm_mat4_rmc(vec4 r, mat4 m, vec4 c) {
  */
 CGLM_INLINE
 void
-glm_mat4_make(float * __restrict src, mat4 dest) {
+glm_mat4_make(const float * __restrict src, mat4 dest) {
   dest[0][0] = src[0];   dest[1][0] = src[4];
   dest[0][1] = src[1];   dest[1][1] = src[5];
   dest[0][2] = src[2];   dest[1][2] = src[6];

include/cglm/mat4x2.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat4x2_copy(mat4x2 mat, mat4x2 dest);
    CGLM_INLINE void glm_mat4x2_zero(mat4x2 mat);
-   CGLM_INLINE void glm_mat4x2_make(float * __restrict src, mat4x2 dest);
+   CGLM_INLINE void glm_mat4x2_make(const float * __restrict src, mat4x2 dest);
    CGLM_INLINE void glm_mat4x2_mul(mat4x2 m1, mat2x4 m2, mat4 dest);
    CGLM_INLINE void glm_mat4x2_mulv(mat4x2 m, vec2 v, vec4 dest);
    CGLM_INLINE void glm_mat4x2_transpose(mat4x2 m, mat2x4 dest);
@@ -72,7 +72,7 @@ glm_mat4x2_zero(mat4x2 mat) {
  */
 CGLM_INLINE
 void
-glm_mat4x2_make(float * __restrict src, mat4x2 dest) {
+glm_mat4x2_make(const float * __restrict src, mat4x2 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
 

include/cglm/mat4x3.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat4x3_copy(mat4x3 mat, mat4x3 dest);
    CGLM_INLINE void glm_mat4x3_zero(mat4x3 mat);
-   CGLM_INLINE void glm_mat4x3_make(float * __restrict src, mat4x3 dest);
+   CGLM_INLINE void glm_mat4x3_make(const float * __restrict src, mat4x3 dest);
    CGLM_INLINE void glm_mat4x3_mul(mat4x3 m1, mat3x4 m2, mat4 dest);
    CGLM_INLINE void glm_mat4x3_mulv(mat4x3 m, vec3 v, vec4 dest);
    CGLM_INLINE void glm_mat4x3_transpose(mat4x3 m, mat3x4 dest);
@@ -77,7 +77,7 @@ glm_mat4x3_zero(mat4x3 mat) {
  */
 CGLM_INLINE
 void
-glm_mat4x3_make(float * __restrict src, mat4x3 dest) {
+glm_mat4x3_make(const float * __restrict src, mat4x3 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/plane.h

@@ -33,7 +33,7 @@ void
 glm_plane_normalize(vec4 plane) {
   float norm;
   
-  if ((norm = glm_vec3_norm(plane)) == 0.0f) {
+  if (CGLM_UNLIKELY((norm = glm_vec3_norm(plane)) < FLT_EPSILON)) {
     glm_vec4_zero(plane);
     return;
   }

include/cglm/quat.h

@@ -894,7 +894,7 @@ glm_quat_rotate_atm(mat4 m, versor q, vec3 pivot) {
  */
 CGLM_INLINE
 void
-glm_quat_make(float * __restrict src, versor dest) {
+glm_quat_make(const float * __restrict src, versor dest) {
   dest[0] = src[0]; dest[1] = src[1];
   dest[2] = src[2]; dest[3] = src[3];
 }

include/cglm/ray.h

@@ -7,12 +7,18 @@
 
 /*
  Functions:
-   CGLM_INLINE bool glm_line_triangle_intersect(vec3   origin,
+   CGLM_INLINE bool glm_ray_triangle(vec3   origin,
-                                                vec3   direction,
+                                     vec3   direction,
-                                                vec3   v0,
+                                     vec3   v0,
-                                                vec3   v1,
+                                     vec3   v1,
-                                                vec3   v2,
+                                     vec3   v2,
-                                                float *d);
+                                     float *d);
+ CGLM_INLINE bool glm_ray_sphere(vec3 origin,
+                                 vec3 dir,
+                                 vec4 s,
+                                 float * __restrict t1,
+                                 float * __restrict t2)
+ CGLM_INLINE void glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point);
 */
 
 #ifndef cglm_ray_h
@@ -31,7 +37,6 @@
  * @param[in, out] d         distance to intersection
  * @return whether there is intersection
  */
-
 CGLM_INLINE
 bool
 glm_ray_triangle(vec3   origin,
@@ -74,4 +79,96 @@ glm_ray_triangle(vec3   origin,
   return dist > epsilon;
 }
 
+/*!
+ * @brief ray sphere intersection
+ *
+ * returns false if there is no intersection if true:
+ *
+ * - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+ * - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+ * - t1 < 0, t2 < 0: no intersection ahead of the ray ( returns false )
+ * - the caller can check if the intersection points (t1 and t2) fall within a
+ *   specific range (for example, tmin < t1, t2 < tmax) to determine if the
+ *   intersections are within a desired segment of the ray
+ *
+ * @param[in]  origin ray origin
+ * @param[out] dir    normalized ray direction
+ * @param[in]  s      sphere  [center.x, center.y, center.z, radii]
+ * @param[in]  t1     near point1 (closer to origin)
+ * @param[in]  t2     far point2 (farther from origin)
+ *
+ * @returns whether there is intersection
+ */
+CGLM_INLINE
+bool 
+glm_ray_sphere(vec3 origin,
+               vec3 dir,
+               vec4 s,
+               float * __restrict t1,
+               float * __restrict t2) {
+  vec3  dp;
+  float r2, ddp, dpp, dscr, q, tmp, _t1, _t2;
+
+  glm_vec3_sub(s, origin, dp);
+
+  ddp = glm_vec3_dot(dir, dp);
+  dpp = glm_vec3_norm2(dp);
+
+  /* compute the remedy term for numerical stability */
+  glm_vec3_mulsubs(dir, ddp, dp); /* dp: remedy term */
+
+  r2   = s[3] * s[3];
+  dscr = r2 - glm_vec3_norm2(dp);
+
+  if (dscr < 0.0f) {
+    /* no intersection */
+    return false;
+  }
+
+  dscr = sqrtf(dscr);
+  q    = (ddp >= 0.0f) ? (ddp + dscr) : (ddp - dscr);
+
+  /*
+     include Press, William H., Saul A. Teukolsky,
+     William T. Vetterling, and Brian P. Flannery,
+     "Numerical Recipes in C," Cambridge University Press, 1992.
+   */
+  _t1 = q;
+  _t2 = (dpp - r2) / q;
+
+  /* adjust t1 and t2 to ensure t1 is the closer intersection */
+  if (_t1 > _t2) {
+    tmp = _t1;
+    _t1 = _t2;
+    _t2 = tmp;
+  }
+
+  *t1 = _t1;
+  *t2 = _t2;
+
+  /* check if the closest intersection (t1) is behind the ray's origin */
+  if (_t1 < 0.0f && _t2 < 0.0f) {
+    /* both intersections are behind the ray, no visible intersection */
+    return false;
+  }
+
+  return true;
+}
+
+/*!
+ * @brief point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+ *
+ * @param[in]  orig  origin of ray
+ * @param[in]  dir   direction of ray
+ * @param[in]  t     parameter
+ * @param[out] point point at t
+ */
+CGLM_INLINE
+void
+glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point) {
+  vec3 dst;
+  glm_vec3_scale(dir, t, dst);
+  glm_vec3_add(orig, dst, point);
+}
+
 #endif

include/cglm/simd/intrin.h

@@ -10,6 +10,9 @@
 
 #if defined( _MSC_VER )
 #  if (defined(_M_AMD64) || defined(_M_X64)) || _M_IX86_FP == 2
+#    ifndef __SSE__
+#      define __SSE__
+#    endif
 #    ifndef __SSE2__
 #      define __SSE2__
 #    endif
@@ -24,15 +27,22 @@
 #  endif
 #endif
 
-#if defined( __SSE__ ) || defined( __SSE2__ )
+#if defined(__SSE__)
 #  include <xmmintrin.h>
-#  include <emmintrin.h>
 #  define CGLM_SSE_FP 1
 #  ifndef CGLM_SIMD_x86
 #    define CGLM_SIMD_x86
 #  endif
 #endif
 
+#if defined(__SSE2__)
+#  include <emmintrin.h>
+#  define CGLM_SSE2_FP 1
+#  ifndef CGLM_SIMD_x86
+#    define CGLM_SIMD_x86
+#  endif
+#endif
+
 #if defined(__SSE3__)
 #  include <pmmintrin.h>
 #  ifndef CGLM_SIMD_x86

include/cglm/simd/x86.h

@@ -21,7 +21,7 @@
 #define glmm_set1(x) _mm_set1_ps(x)
 #define glmm_128     __m128
 
-#ifdef CGLM_USE_INT_DOMAIN
+#if defined(CGLM_USE_INT_DOMAIN) && defined(__SSE2__)
 #  define glmm_shuff1(xmm, z, y, x, w)                                        \
      _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(xmm),                \
                                         _MM_SHUFFLE(z, y, x, w)))
@@ -55,17 +55,40 @@
 #endif
 
 /* Note that `0x80000000` corresponds to `INT_MIN` for a 32-bit int. */
-#define GLMM_NEGZEROf ((int)0x80000000) /*  0x80000000 ---> -0.0f  */
 
-#define GLMM__SIGNMASKf(X, Y, Z, W)                                           \
+#if defined(__SSE2__)
+#  define GLMM_NEGZEROf ((int)0x80000000) /*  0x80000000 ---> -0.0f  */
+#  define GLMM_POSZEROf ((int)0x00000000) /*  0x00000000 ---> +0.0f  */
+#else
+#  ifdef CGLM_FAST_MATH
+     union { int i; float f; } static GLMM_NEGZEROf_TU = { .i = (int)0x80000000 };
+#    define GLMM_NEGZEROf GLMM_NEGZEROf_TU.f
+#    define GLMM_POSZEROf 0.0f
+#  else
+#    define GLMM_NEGZEROf -0.0f
+#    define GLMM_POSZEROf  0.0f
+#  endif
+#endif
+
+#if defined(__SSE2__)
+#  define GLMM__SIGNMASKf(X, Y, Z, W)                                         \
    _mm_castsi128_ps(_mm_set_epi32(X, Y, Z, W))
   /* _mm_set_ps(X, Y, Z, W); */
+#else
+#  define GLMM__SIGNMASKf(X, Y, Z, W)  _mm_set_ps(X, Y, Z, W)
+#endif
 
-#define glmm_float32x4_SIGNMASK_PNPN GLMM__SIGNMASKf(0, GLMM_NEGZEROf, 0, GLMM_NEGZEROf)
+#define glmm_float32x4_SIGNMASK_PNPN GLMM__SIGNMASKf(GLMM_POSZEROf, GLMM_NEGZEROf, GLMM_POSZEROf, GLMM_NEGZEROf)
-#define glmm_float32x4_SIGNMASK_NPNP GLMM__SIGNMASKf(GLMM_NEGZEROf, 0, GLMM_NEGZEROf, 0)
+#define glmm_float32x4_SIGNMASK_NPNP GLMM__SIGNMASKf(GLMM_NEGZEROf, GLMM_POSZEROf, GLMM_NEGZEROf, GLMM_POSZEROf)
-#define glmm_float32x4_SIGNMASK_NPPN GLMM__SIGNMASKf(GLMM_NEGZEROf, 0, 0, GLMM_NEGZEROf)
+#define glmm_float32x4_SIGNMASK_NPPN GLMM__SIGNMASKf(GLMM_NEGZEROf, GLMM_POSZEROf, GLMM_POSZEROf, GLMM_NEGZEROf)
+
+/* fasth math prevents -0.0f to work */
+#if defined(__SSE2__)
+#  define glmm_float32x4_SIGNMASK_NEG _mm_castsi128_ps(_mm_set1_epi32(GLMM_NEGZEROf)) /* _mm_set1_ps(-0.0f) */
+#else
+#  define glmm_float32x4_SIGNMASK_NEG _mm_set1_ps(GLMM_NEGZEROf)
+#endif
 
-#define glmm_float32x4_SIGNMASK_NEG _mm_castsi128_ps(_mm_set1_epi32(GLMM_NEGZEROf)) /* _mm_set1_ps(-0.0f) */
 #define glmm_float32x8_SIGNMASK_NEG _mm256_castsi256_ps(_mm256_set1_epi32(GLMM_NEGZEROf))
 
 static inline
@@ -207,6 +230,7 @@ glmm_norm_inf(__m128 a) {
   return _mm_cvtss_f32(glmm_vhmax(glmm_abs(a)));
 }
 
+#if defined(__SSE2__)
 static inline
 __m128
 glmm_load3(float v[3]) {
@@ -225,6 +249,7 @@ glmm_store3(float v[3], __m128 vx) {
   _mm_storel_pi(CGLM_CASTPTR_ASSUME_ALIGNED(v, __m64), vx);
   _mm_store_ss(&v[2], glmm_shuff1(vx, 2, 2, 2, 2));
 }
+#endif
 
 static inline
 __m128

include/cglm/struct.h

@@ -41,6 +41,7 @@ extern "C" {
 #include "struct/sphere.h"
 #include "struct/curve.h"
 #include "struct/affine2d.h"
+#include "struct/ray.h"
 
 #ifdef __cplusplus
 }

include/cglm/struct/mat2.h

@@ -27,7 +27,7 @@
    CGLM_INLINE void  glms_mat2_swap_col(mat2 mat, int col1, int col2)
    CGLM_INLINE void  glms_mat2_swap_row(mat2 mat, int row1, int row2)
    CGLM_INLINE float glms_mat2_rmc(vec2 r, mat2 m, vec2 c)
-   CGLM_INLINE mat2s glms_mat2_make(float * __restrict src);
+   CGLM_INLINE mat2s glms_mat2_make(const float * __restrict src);
  */
 
 #ifndef cglms_mat2_h
@@ -267,7 +267,7 @@ glms_mat2_(rmc)(vec2s r, mat2s m, vec2s c) {
  */
 CGLM_INLINE
 mat2s
-glms_mat2_(make)(float * __restrict src) {
+glms_mat2_(make)(const float * __restrict src) {
   mat2s r;
   glm_mat2_make(src, r.raw);
   return r;

include/cglm/struct/mat2x3.h

@@ -12,9 +12,9 @@
 
  Functions:
    CGLM_INLINE mat2x3s glms_mat2x3_zero(void);
-   CGLM_INLINE mat2x3s glms_mat2x3_make(float * __restrict src);
+   CGLM_INLINE mat2x3s glms_mat2x3_make(const float * __restrict src);
    CGLM_INLINE mat2s   glms_mat2x3_mul(mat2x3s m1, mat3x2s m2);
-   CGLM_INLINE vec2s   glms_mat2x3_mulv(mat2x3s m, vec3s v);
+   CGLM_INLINE vec3s   glms_mat2x3_mulv(mat2x3s m, vec2s v);
    CGLM_INLINE mat3x2s glms_mat2x3_transpose(mat2x3s m);
    CGLM_INLINE mat2x3s glms_mat2x3_scale(mat2x3s m, float s);
  */
@@ -55,7 +55,7 @@ glms_mat2x3_(zero)(void) {
  */
 CGLM_INLINE
 mat2x3s
-glms_mat2x3_(make)(float * __restrict src) {
+glms_mat2x3_(make)(const float * __restrict src) {
   mat2x3s r;
   glm_mat2x3_make(src, r.raw);
   return r;
@@ -73,9 +73,9 @@ glms_mat2x3_(make)(float * __restrict src) {
  * @returns destination matrix (mat2s)
  */
 CGLM_INLINE
-mat2s
+mat3s
 glms_mat2x3_(mul)(mat2x3s m1, mat3x2s m2) {
-  mat2s r;
+  mat3s r;
   glm_mat2x3_mul(m1.raw, m2.raw, r.raw);
   return r;
 }
@@ -85,12 +85,12 @@ glms_mat2x3_(mul)(mat2x3s m1, mat3x2s m2) {
  *
  * @param[in]  m    matrix (left)
  * @param[in]  v    vector (right, column vector)
- * @param[out] dest result vector
+ * @returns destination vector (vec3s)
  */
 CGLM_INLINE
-vec2s
+vec3s
-glms_mat2x3_(mulv)(mat2x3s m, vec3s v) {
+glms_mat2x3_(mulv)(mat2x3s m, vec2s v) {
-  vec2s r;
+  vec3s r;
   glm_mat2x3_mulv(m.raw, v.raw, r.raw);
   return r;
 }

include/cglm/struct/mat2x4.h

@@ -12,9 +12,9 @@
 
  Functions:
    CGLM_INLINE mat2x4s glms_mat2x4_zero(void);
-   CGLM_INLINE mat2x4s glms_mat2x4_make(float * __restrict src);
+   CGLM_INLINE mat2x4s glms_mat2x4_make(const float * __restrict src);
    CGLM_INLINE mat2s   glms_mat2x4_mul(mat2x4s m1, mat4x2s m2);
-   CGLM_INLINE vec2s   glms_mat2x4_mulv(mat2x4s m, vec4s v);
+   CGLM_INLINE vec4s   glms_mat2x4_mulv(mat2x4s m, vec2s v);
    CGLM_INLINE mat4x2s glms_mat2x4_transpose(mat2x4s m);
    CGLM_INLINE mat2x4s glms_mat2x4_scale(mat2x4s m, float s);
  */
@@ -55,7 +55,7 @@ glms_mat2x4_(zero)(void) {
  */
 CGLM_INLINE
 mat2x4s
-glms_mat2x4_(make)(float * __restrict src) {
+glms_mat2x4_(make)(const float * __restrict src) {
   mat2x4s r;
   glm_mat2x4_make(src, r.raw);
   return r;
@@ -73,24 +73,24 @@ glms_mat2x4_(make)(float * __restrict src) {
  * @returns destination matrix (mat2s)
  */
 CGLM_INLINE
-mat2s
+mat4s
 glms_mat2x4_(mul)(mat2x4s m1, mat4x2s m2) {
-  mat2s r;
+  mat4s r;
   glm_mat2x4_mul(m1.raw, m2.raw, r.raw);
   return r;
 }
 
 /*!
- * @brief multiply matrix with column vector and store in dest vector
+ * @brief multiply matrix with column vector and store in dest column vector
  *
  * @param[in]  m    matrix (left)
  * @param[in]  v    vector (right, column vector)
  * @param[out] dest result vector
  */
 CGLM_INLINE
-vec2s
+vec4s
-glms_mat2x4_(mulv)(mat2x4s m, vec4s v) {
+glms_mat2x4_(mulv)(mat2x4s m, vec2s v) {
-  vec2s r;
+  vec4s r;
   glm_mat2x4_mulv(m.raw, v.raw, r.raw);
   return r;
 }

include/cglm/struct/mat3.h

@@ -28,7 +28,7 @@
    CGLM_INLINE mat3s  glms_mat3_swap_col(mat3s mat, int col1, int col2);
    CGLM_INLINE mat3s  glms_mat3_swap_row(mat3s mat, int row1, int row2);
    CGLM_INLINE float  glms_mat3_rmc(vec3s r, mat3s m, vec3s c);
-   CGLM_INLINE mat3s  glms_mat3_make(float * __restrict src);
+   CGLM_INLINE mat3s  glms_mat3_make(const float * __restrict src);
  */
 
 #ifndef cglms_mat3s_h
@@ -294,7 +294,7 @@ glms_mat3_(rmc)(vec3s r, mat3s m, vec3s c) {
  */
 CGLM_INLINE
 mat3s
-glms_mat3_(make)(float * __restrict src) {
+glms_mat3_(make)(const float * __restrict src) {
   mat3s r;
   glm_mat3_make(src, r.raw);
   return r;

include/cglm/struct/mat3x2.h

@@ -12,9 +12,9 @@
 
  Functions:
    CGLM_INLINE mat3x2s glms_mat3x2_zero(void);
-   CGLM_INLINE mat3x2s glms_mat3x2_make(float * __restrict src);
+   CGLM_INLINE mat3x2s glms_mat3x2_make(const float * __restrict src);
-   CGLM_INLINE mat3s   glms_mat3x2_mul(mat3x2s m1, mat2x3s m2);
+   CGLM_INLINE mat2s   glms_mat3x2_mul(mat3x2s m1, mat2x3s m2);
-   CGLM_INLINE vec3s   glms_mat3x2_mulv(mat3x2s m, vec2s v);
+   CGLM_INLINE vec2s   glms_mat3x2_mulv(mat3x2s m, vec3s v);
    CGLM_INLINE mat2x3s glms_mat3x2_transpose(mat3x2s m);
    CGLM_INLINE mat3x2s glms_mat3x2_scale(mat3x2s m, float s);
  */
@@ -55,7 +55,7 @@ glms_mat3x2_(zero)(void) {
  */
 CGLM_INLINE
 mat3x2s
-glms_mat3x2_(make)(float * __restrict src) {
+glms_mat3x2_(make)(const float * __restrict src) {
   mat3x2s r;
   glm_mat3x2_make(src, r.raw);
   return r;
@@ -73,9 +73,9 @@ glms_mat3x2_(make)(float * __restrict src) {
  * @returns destination matrix (mat3s)
  */
 CGLM_INLINE
-mat3s
+mat2s
 glms_mat3x2_(mul)(mat3x2s m1, mat2x3s m2) {
-  mat3s r;
+  mat2s r;
   glm_mat3x2_mul(m1.raw, m2.raw, r.raw);
   return r;
 }
@@ -88,9 +88,9 @@ glms_mat3x2_(mul)(mat3x2s m1, mat2x3s m2) {
  * @param[out] dest result vector
  */
 CGLM_INLINE
-vec3s
+vec2s
-glms_mat3x2_(mulv)(mat3x2s m, vec2s v) {
+glms_mat3x2_(mulv)(mat3x2s m, vec3s v) {
-  vec3s r;
+  vec2s r;
   glm_mat3x2_mulv(m.raw, v.raw, r.raw);
   return r;
 }

include/cglm/struct/mat3x4.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat3x4s glms_mat3x4_zero(void);
-   CGLM_INLINE mat3x4s glms_mat3x4_make(float * __restrict src);
+   CGLM_INLINE mat3x4s glms_mat3x4_make(const float * __restrict src);
    CGLM_INLINE mat3s   glms_mat3x4_mul(mat3x4s m1, mat4x3s m2);
    CGLM_INLINE vec3s   glms_mat3x4_mulv(mat3x4s m, vec4s v);
    CGLM_INLINE mat4x3s glms_mat3x4_transpose(mat3x4s m);
@@ -55,7 +55,7 @@ glms_mat3x4_(zero)(void) {
  */
 CGLM_INLINE
 mat3x4s
-glms_mat3x4_(make)(float * __restrict src) {
+glms_mat3x4_(make)(const float * __restrict src) {
   mat3x4s r;
   glm_mat3x4_make(src, r.raw);
   return r;

include/cglm/struct/mat4.h

@@ -42,7 +42,7 @@
    CGLM_INLINE mat4s   glms_mat4_swap_col(mat4s mat, int col1, int col2);
    CGLM_INLINE mat4s   glms_mat4_swap_row(mat4s mat, int row1, int row2);
    CGLM_INLINE float   glms_mat4_rmc(vec4s r, mat4s m, vec4s c);
-   CGLM_INLINE mat4s   glms_mat4_make(float * __restrict src);
+   CGLM_INLINE mat4s   glms_mat4_make(const float * __restrict src);
  */
 
 #ifndef cglms_mat4s_h
@@ -468,7 +468,7 @@ glms_mat4_(rmc)(vec4s r, mat4s m, vec4s c) {
  */
 CGLM_INLINE
 mat4s
-glms_mat4_(make)(float * __restrict src) {
+glms_mat4_(make)(const float * __restrict src) {
   mat4s r;
   glm_mat4_make(src, r.raw);
   return r;

include/cglm/struct/mat4x2.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat4x2s glms_mat4x2_zero(void);
-   CGLM_INLINE mat4x2s glms_mat4x2_make(float * __restrict src);
+   CGLM_INLINE mat4x2s glms_mat4x2_make(const float * __restrict src);
    CGLM_INLINE mat4s   glms_mat4x2_mul(mat4x2s m1, mat2x4s m2);
    CGLM_INLINE vec4s   glms_mat4x2_mulv(mat4x2s m, vec2s v);
    CGLM_INLINE mat2x4s glms_mat4x2_transpose(mat4x2s m);
@@ -56,7 +56,7 @@ glms_mat4x2_(zero)(void) {
  */
 CGLM_INLINE
 mat4x2s
-glms_mat4x2_(make)(float * __restrict src) {
+glms_mat4x2_(make)(const float * __restrict src) {
   mat4x2s r;
   glm_mat4x2_make(src, r.raw);
   return r;

include/cglm/struct/mat4x3.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat4x3s glms_mat4x3_zero(void);
-   CGLM_INLINE mat4x3s glms_mat4x3_make(float * __restrict src);
+   CGLM_INLINE mat4x3s glms_mat4x3_make(const float * __restrict src);
    CGLM_INLINE mat4s   glms_mat4x3_mul(mat4x3s m1, mat3x4s m2);
    CGLM_INLINE vec4s   glms_mat4x3_mulv(mat4x3s m, vec3s v);
    CGLM_INLINE mat3x4s glms_mat4x3_transpose(mat4x3s m);
@@ -55,7 +55,7 @@ glms_mat4x3_(zero)(void) {
  */
 CGLM_INLINE
 mat4x3s
-glms_mat4x3_(make)(float * __restrict src) {
+glms_mat4x3_(make)(const float * __restrict src) {
   mat4x3s r;
   glm_mat4x3_make(src, r.raw);
   return r;

include/cglm/struct/quat.h

@@ -574,7 +574,7 @@ glms_quat_(rotate_atm)(versors q, vec3s pivot) {
  */
 CGLM_INLINE
 versors
-glms_quat_(make)(float * __restrict src) {
+glms_quat_(make)(const float * __restrict src) {
   versors dest;
   glm_quat_make(src, dest.raw);
   return dest;

include/cglm/struct/ray.h

@@ -0,0 +1,86 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglms_ray_h
+#define cglms_ray_h
+
+#include "../common.h"
+#include "../types-struct.h"
+#include "../ray.h"
+
+/* api definition */
+#define glms_ray_(NAME) CGLM_STRUCTAPI(ray, NAME)
+
+/*!
+ * @brief Möller–Trumbore ray-triangle intersection algorithm
+ *
+ * @param[in] origin         origin of ray
+ * @param[in] direction      direction of ray
+ * @param[in] v0             first vertex of triangle
+ * @param[in] v1             second vertex of triangle
+ * @param[in] v2             third vertex of triangle
+ * @param[in, out] d         distance to intersection
+ * @return whether there is intersection
+ */
+CGLM_INLINE
+bool
+glms_ray_(triangle)(vec3s  origin,
+                    vec3s  direction,
+                    vec3s  v0,
+                    vec3s  v1,
+                    vec3s  v2,
+                    float *d) {
+  return glm_ray_triangle(origin.raw, direction.raw, v0.raw, v1.raw, v2.raw, d);
+}
+
+/*!
+ * @brief ray sphere intersection
+ *
+ * returns false if there is no intersection if true:
+ *
+ * - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+ * - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+ * - t1 < 0, t2 < 0: no intersection ahead of the ray ( returns false )
+ * - the caller can check if the intersection points (t1 and t2) fall within a
+ *   specific range (for example, tmin < t1, t2 < tmax) to determine if the
+ *   intersections are within a desired segment of the ray
+ *
+ * @param[in]  origin ray origin
+ * @param[out] dir    normalized ray direction
+ * @param[in]  s      sphere  [center.x, center.y, center.z, radii]
+ * @param[in]  t1     near point1 (closer to origin)
+ * @param[in]  t2     far point2 (farther from origin)
+ *
+ * @returns whether there is intersection
+ */
+CGLM_INLINE
+bool
+glms_ray_(sphere)(vec3s origin,
+                  vec3s dir,
+                  vec4s s,
+                  float * __restrict t1,
+                  float * __restrict t2) {
+  return glm_ray_sphere(origin.raw, dir.raw, s.raw, t1, t2);
+}
+
+/*!
+ * @brief point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+ *
+ * @param[in]  orig  origin of ray
+ * @param[in]  dir   direction of ray
+ * @param[in]  t     parameter
+ * @returns point point at t
+ */
+CGLM_INLINE
+vec3s
+glms_ray_(at)(vec3s orig, vec3s dir, float t) {
+  vec3s r;
+  glm_ray_at(orig.raw, dir.raw, t, r.raw);
+  return r;
+}
+
+#endif /* cglms_ray_h */

include/cglm/struct/vec2.h

@@ -54,6 +54,8 @@
    CGLM_INLINE vec2s glms_vec2_clamp(vec2s v, float minVal, float maxVal)
    CGLM_INLINE vec2s glms_vec2_lerp(vec2s from, vec2s to, float t)
    CGLM_INLINE vec2s glms_vec2_make(float * restrict src)
+   CGLM_INLINE vec2s glms_vec2_reflect(vec2s v, vec2s n)
+   CGLM_INLINE bool  glms_vec2_refract(vec2s v, vec2s n, float eta, vec2s *dest)
  */
 
 #ifndef cglms_vec2s_h
@@ -685,10 +687,45 @@ glms_vec2_(lerp)(vec2s from, vec2s to, float t) {
  */
 CGLM_INLINE
 vec2s
-glms_vec2_(make)(float * __restrict src) {
+glms_vec2_(make)(const float * __restrict src) {
   vec2s dest;
   glm_vec2_make(src, dest.raw);
   return dest;
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @returns reflection result
+ */
+CGLM_INLINE
+vec2s
+glms_vec2_(reflect)(vec2s v, vec2s n) {
+  vec2s dest;
+  glm_vec2_reflect(v.raw, n.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief computes refraction vector for an incident vector and a surface normal.
+ *
+ * calculates the refraction vector based on Snell's law. If total internal reflection
+ * occurs (angle too great given eta), dest is set to zero and returns false.
+ * Otherwise, computes refraction vector, stores it in dest, and returns true.
+ *
+ * @param[in]  v    normalized incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction (incident/transmitted)
+ * @param[out] dest refraction vector if refraction occurs; zero vector otherwise
+ *
+ * @returns true if refraction occurs; false if total internal reflection occurs.
+ */
+CGLM_INLINE
+bool
+glms_vec2_(refract)(vec2s v, vec2s n, float eta, vec2s * __restrict dest) {
+  return glm_vec2_refract(v.raw, n.raw, eta, dest->raw);
+}
+
 #endif /* cglms_vec2s_h */

include/cglm/struct/vec3.h

@@ -76,6 +76,9 @@
    CGLM_INLINE vec3s glms_vec3_smoothinterpc(vec3s from, vec3s to, float t);
    CGLM_INLINE vec3s glms_vec3_swizzle(vec3s v, int mask);
    CGLM_INLINE vec3s glms_vec3_make(float * restrict src);
+   CGLM_INLINE vec3s glms_vec3_faceforward(vec3s n, vec3s v, vec3s nref);
+   CGLM_INLINE vec3s glms_vec3_reflect(vec3s v, vec3s n);
+   CGLM_INLINE bool  glms_vec3_refract(vec3s v, vec3s n, float eta, vec3s *dest)
 
  Convenient:
    CGLM_INLINE vec3s glms_cross(vec3s a, vec3s b);
@@ -1077,10 +1080,63 @@ glms_vec3_(swizzle)(vec3s v, int mask) {
  */
 CGLM_INLINE
 vec3s
-glms_vec3_(make)(float * __restrict src) {
+glms_vec3_(make)(const float * __restrict src) {
   vec3s dest;
   glm_vec3_make(src, dest.raw);
   return dest;
 }
 
+/*!
+ * @brief a vector pointing in the same direction as another
+ *
+ * orients a vector to point away from a surface as defined by its normal
+ *
+ * @param[in] n      vector to orient.
+ * @param[in] v      incident vector
+ * @param[in] nref   reference vector
+ * @returns oriented vector, pointing away from the surface.
+ */
+CGLM_INLINE
+vec3s
+glms_vec3_(faceforward)(vec3s n, vec3s v, vec3s nref) {
+  vec3s dest;
+  glm_vec3_faceforward(n.raw, v.raw, nref.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @returns reflection result
+ */
+CGLM_INLINE
+vec3s
+glms_vec3_(reflect)(vec3s v, vec3s n) {
+  vec3s dest;
+  glm_vec3_reflect(v.raw, n.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief computes refraction vector for an incident vector and a surface normal.
+ *
+ * calculates the refraction vector based on Snell's law. If total internal reflection
+ * occurs (angle too great given eta), dest is set to zero and returns false.
+ * Otherwise, computes refraction vector, stores it in dest, and returns true.
+ *
+ * @param[in]  v    normalized incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction (incident/transmitted)
+ * @param[out] dest refraction vector if refraction occurs; zero vector otherwise
+ *
+ * @returns true if refraction occurs; false if total internal reflection occurs.
+ */
+CGLM_INLINE
+bool
+glms_vec3_(refract)(vec3s v, vec3s n, float eta, vec3s * __restrict dest) {
+  return glm_vec3_refract(v.raw, n.raw, eta, dest->raw);
+}
+
 #endif /* cglms_vec3s_h */

include/cglm/struct/vec4.h

@@ -67,6 +67,8 @@
    CGLM_INLINE vec4s glms_vec4_cubic(float s);
    CGLM_INLINE vec4s glms_vec4_swizzle(vec4s v, int mask);
    CGLM_INLINE vec4s glms_vec4_make(float * restrict src);
+   CGLM_INLINE vec4s glms_vec4_reflect(vec4s v, vec4s n);
+   CGLM_INLINE bool  glms_vec4_refract(vec4s v, vec4s n, float eta, vec4s *dest)
  */
 
 #ifndef cglms_vec4s_h
@@ -921,10 +923,49 @@ glms_vec4_(swizzle)(vec4s v, int mask) {
  */
 CGLM_INLINE
 vec4s
-glms_vec4_(make)(float * __restrict src) {
+glms_vec4_(make)(const float * __restrict src) {
   vec4s dest;
   glm_vec4_make(src, dest.raw);
   return dest;
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  v    incident vector
+ * @param[in]  n    normalized normal vector
+ * @returns reflection result
+ */
+CGLM_INLINE
+vec4s
+glms_vec4_(reflect)(vec4s v, vec4s n) {
+  vec4s dest;
+  glm_vec4_reflect(v.raw, n.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief computes refraction vector for an incident vector and a surface normal.
+ *
+ * calculates the refraction vector based on Snell's law. If total internal reflection
+ * occurs (angle too great given eta), dest is set to zero and returns false.
+ * Otherwise, computes refraction vector, stores it in dest, and returns true.
+ *
+ * this implementation does not explicitly preserve the 'w' component of the
+ * incident vector 'I' in the output 'dest', users requiring the preservation of
+ * the 'w' component should manually adjust 'dest' after calling this function.
+ *
+ * @param[in]  v    normalized incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction (incident/transmitted)
+ * @param[out] dest refraction vector if refraction occurs; zero vector otherwise
+ *
+ * @returns true if refraction occurs; false if total internal reflection occurs.
+ */
+CGLM_INLINE
+bool
+glms_vec4_(refract)(vec4s v, vec4s n, float eta, vec4s * __restrict dest) {
+  return glm_vec4_refract(v.raw, n.raw, eta, dest->raw);
+}
+
 #endif /* cglms_vec4s_h */

include/cglm/vec2.h

@@ -55,7 +55,8 @@
    CGLM_INLINE void  glm_vec2_clamp(vec2 v, float minVal, float maxVal)
    CGLM_INLINE void  glm_vec2_lerp(vec2 from, vec2 to, float t, vec2 dest)
    CGLM_INLINE void  glm_vec2_make(float * restrict src, vec2 dest)
-
+   CGLM_INLINE void  glm_vec2_reflect(vec2 v, vec2 n, vec2 dest)
+   CGLM_INLINE void  glm_vec2_refract(vec2 v, vec2 n, float eta, vec2 dest)
  */
 
 #ifndef cglm_vec2_h
@@ -277,7 +278,7 @@ glm_vec2_scale_as(vec2 v, float s, vec2 dest) {
   float norm;
   norm = glm_vec2_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     glm_vec2_zero(dest);
     return;
   }
@@ -541,7 +542,7 @@ glm_vec2_normalize(vec2 v) {
 
   norm = glm_vec2_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     v[0] = v[1] = 0.0f;
     return;
   }
@@ -562,7 +563,7 @@ glm_vec2_normalize_to(vec2 v, vec2 dest) {
 
   norm = glm_vec2_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     glm_vec2_zero(dest);
     return;
   }
@@ -708,8 +709,56 @@ glm_vec2_lerp(vec2 from, vec2 to, float t, vec2 dest) {
  */
 CGLM_INLINE
 void
-glm_vec2_make(float * __restrict src, vec2 dest) {
+glm_vec2_make(const float * __restrict src, vec2 dest) {
   dest[0] = src[0]; dest[1] = src[1];
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  v    incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[out] dest destination vector for the reflection result
+ */
+CGLM_INLINE
+void
+glm_vec2_reflect(vec2 v, vec2 n, vec2 dest) {
+  vec2 temp;
+  glm_vec2_scale(n, 2.0f * glm_vec2_dot(v, n), temp);
+  glm_vec2_sub(v, temp, dest);
+}
+
+/*!
+ * @brief computes refraction vector for an incident vector and a surface normal.
+ *
+ * calculates the refraction vector based on Snell's law. If total internal reflection
+ * occurs (angle too great given eta), dest is set to zero and returns false.
+ * Otherwise, computes refraction vector, stores it in dest, and returns true.
+ *
+ * @param[in]  v    normalized incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction (incident/transmitted)
+ * @param[out] dest refraction vector if refraction occurs; zero vector otherwise
+ *
+ * @returns true if refraction occurs; false if total internal reflection occurs.
+ */
+CGLM_INLINE
+bool
+glm_vec2_refract(vec2 v, vec2 n, float eta, vec2 dest) {
+  float ndi, eni, k;
+
+  ndi = glm_vec2_dot(n, v);
+  eni = eta * ndi;
+  k   = 1.0f + eta * eta - eni * eni;
+
+  if (k < 0.0f) {
+    glm_vec2_zero(dest);
+    return false;
+  }
+
+  glm_vec2_scale(v, eta, dest);
+  glm_vec2_mulsubs(n, eni + sqrtf(k), dest);
+  return true;
+}
+
 #endif /* cglm_vec2_h */

include/cglm/vec3.h

@@ -80,6 +80,9 @@
    CGLM_INLINE void  glm_vec3_smoothinterpc(vec3 from, vec3 to, float t, vec3 dest);
    CGLM_INLINE void  glm_vec3_swizzle(vec3 v, int mask, vec3 dest);
    CGLM_INLINE void  glm_vec3_make(float * restrict src, vec3 dest);
+   CGLM_INLINE void  glm_vec3_faceforward(vec3 n, vec3 v, vec3 nref, vec3 dest);
+   CGLM_INLINE void  glm_vec3_reflect(vec3 v, vec3 n, vec3 dest);
+   CGLM_INLINE void  glm_vec3_refract(vec3 v, vec3 n, float eta, vec3 dest);
 
  Convenient:
    CGLM_INLINE void  glm_cross(vec3 a, vec3 b, vec3 d);
@@ -369,7 +372,7 @@ glm_vec3_scale_as(vec3 v, float s, vec3 dest) {
   float norm;
   norm = glm_vec3_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     glm_vec3_zero(dest);
     return;
   }
@@ -648,7 +651,7 @@ glm_vec3_normalize(vec3 v) {
 
   norm = glm_vec3_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     v[0] = v[1] = v[2] = 0.0f;
     return;
   }
@@ -669,7 +672,7 @@ glm_vec3_normalize_to(vec3 v, vec3 dest) {
 
   norm = glm_vec3_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     glm_vec3_zero(dest);
     return;
   }
@@ -1196,10 +1199,80 @@ glm_normalize_to(vec3 v, vec3 dest) {
  */
 CGLM_INLINE
 void
-glm_vec3_make(float * __restrict src, vec3 dest) {
+glm_vec3_make(const float * __restrict src, vec3 dest) {
   dest[0] = src[0];
   dest[1] = src[1];
   dest[2] = src[2];
 }
 
+/*!
+ * @brief a vector pointing in the same direction as another
+ *
+ * orients a vector to point away from a surface as defined by its normal
+ *
+ * @param[in] n      vector to orient
+ * @param[in] v      incident vector
+ * @param[in] nref   reference vector
+ * @param[out] dest  oriented vector, pointing away from the surface
+ */
+CGLM_INLINE
+void
+glm_vec3_faceforward(vec3 n, vec3 v, vec3 nref, vec3 dest) {
+  if (glm_vec3_dot(v, nref) < 0.0f) {
+    /* N is facing away from I */
+    glm_vec3_copy(n, dest);
+  } else {
+    /* N is facing towards I, negate it */
+    glm_vec3_negate_to(n, dest);
+  }
+}
+
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  v    incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[out] dest reflection result
+ */
+CGLM_INLINE
+void
+glm_vec3_reflect(vec3 v, vec3 n, vec3 dest) {
+  vec3 temp;
+  glm_vec3_scale(n, 2.0f * glm_vec3_dot(v, n), temp);
+  glm_vec3_sub(v, temp, dest);
+}
+
+/*!
+ * @brief computes refraction vector for an incident vector and a surface normal.
+ *
+ * calculates the refraction vector based on Snell's law. If total internal reflection
+ * occurs (angle too great given eta), dest is set to zero and returns false.
+ * Otherwise, computes refraction vector, stores it in dest, and returns true.
+ *
+ * @param[in]  v    normalized incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction (incident/transmitted)
+ * @param[out] dest refraction vector if refraction occurs; zero vector otherwise
+ *
+ * @returns true if refraction occurs; false if total internal reflection occurs.
+ */
+CGLM_INLINE
+bool
+glm_vec3_refract(vec3 v, vec3 n, float eta, vec3 dest) {
+  float ndi, eni, k;
+
+  ndi = glm_vec3_dot(n, v);
+  eni = eta * ndi;
+  k   = 1.0f + eta * eta - eni * eni;
+
+  if (k < 0.0f) {
+    glm_vec3_zero(dest);
+    return false;
+  }
+
+  glm_vec3_scale(v, eta, dest);
+  glm_vec3_mulsubs(n, eni + sqrtf(k), dest);
+  return true;
+}
+
 #endif /* cglm_vec3_h */

include/cglm/vec4.h

@@ -65,6 +65,8 @@
    CGLM_INLINE void  glm_vec4_smoothinterpc(vec4 from, vec4 to, float t, vec4 dest);
    CGLM_INLINE void  glm_vec4_swizzle(vec4 v, int mask, vec4 dest);
    CGLM_INLINE void  glm_vec4_make(float * restrict src, vec4 dest);
+   CGLM_INLINE void  glm_vec4_reflect(vec4 v, vec4 n, vec4 dest);
+   CGLM_INLINE void  glm_vec4_refract(vec4 v, vec4 n, float eta, vec4 dest);
 
  DEPRECATED:
    glm_vec4_dup
@@ -485,7 +487,7 @@ glm_vec4_scale_as(vec4 v, float s, vec4 dest) {
   float norm;
   norm = glm_vec4_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     glm_vec4_zero(dest);
     return;
   }
@@ -916,7 +918,7 @@ glm_vec4_normalize_to(vec4 v, vec4 dest) {
   /* dot  = _mm_cvtss_f32(xdot); */
   dot  = wasm_f32x4_extract_lane(xdot, 0);
 
-  if (dot == 0.0f) {
+  if (CGLM_UNLIKELY(dot < FLT_EPSILON)) {
     glmm_store(dest, wasm_f32x4_const_splat(0.f));
     return;
   }
@@ -930,7 +932,7 @@ glm_vec4_normalize_to(vec4 v, vec4 dest) {
   xdot = glmm_vdot(x0, x0);
   dot  = _mm_cvtss_f32(xdot);
 
-  if (dot == 0.0f) {
+  if (CGLM_UNLIKELY(dot < FLT_EPSILON)) {
     glmm_store(dest, _mm_setzero_ps());
     return;
   }
@@ -941,7 +943,7 @@ glm_vec4_normalize_to(vec4 v, vec4 dest) {
 
   norm = glm_vec4_norm(v);
 
-  if (norm == 0.0f) {
+  if (CGLM_UNLIKELY(norm < FLT_EPSILON)) {
     glm_vec4_zero(dest);
     return;
   }
@@ -1299,9 +1301,65 @@ glm_vec4_swizzle(vec4 v, int mask, vec4 dest) {
  */
 CGLM_INLINE
 void
-glm_vec4_make(float * __restrict src, vec4 dest) {
+glm_vec4_make(const float * __restrict src, vec4 dest) {
   dest[0] = src[0]; dest[1] = src[1];
   dest[2] = src[2]; dest[3] = src[3];
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  v    incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[out] dest destination vector for the reflection result
+ */
+CGLM_INLINE
+void
+glm_vec4_reflect(vec4 v, vec4 n, vec4 dest) {
+  vec4 temp;
+
+  /* TODO: direct simd touch */
+  glm_vec4_scale(n, 2.0f * glm_vec4_dot(v, n), temp);
+  glm_vec4_sub(v, temp, dest);
+
+  dest[3] = v[3];
+}
+
+/*!
+ * @brief computes refraction vector for an incident vector and a surface normal.
+ *
+ * calculates the refraction vector based on Snell's law. If total internal reflection
+ * occurs (angle too great given eta), dest is set to zero and returns false.
+ * Otherwise, computes refraction vector, stores it in dest, and returns true.
+ *
+ * this implementation does not explicitly preserve the 'w' component of the
+ * incident vector 'I' in the output 'dest', users requiring the preservation of
+ * the 'w' component should manually adjust 'dest' after calling this function.
+ *
+ * @param[in]  v    normalized incident vector
+ * @param[in]  n    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction (incident/transmitted)
+ * @param[out] dest refraction vector if refraction occurs; zero vector otherwise
+ *
+ * @returns true if refraction occurs; false if total internal reflection occurs.
+ */
+CGLM_INLINE
+bool
+glm_vec4_refract(vec4 v, vec4 n, float eta, vec4 dest) {
+  float ndi, eni, k;
+
+  ndi = glm_vec4_dot(n, v);
+  eni = eta * ndi;
+  k   = 1.0f + eta * eta - eni * eni;
+
+  if (k < 0.0f) {
+    glm_vec4_zero(dest);
+    return false;
+  }
+
+  glm_vec4_scale(v, eta, dest);
+  glm_vec4_mulsubs(n, eni + sqrtf(k), dest);
+  return true;
+}
+
 #endif /* cglm_vec4_h */

include/cglm/version.h

@@ -10,6 +10,6 @@
 
 #define CGLM_VERSION_MAJOR 0
 #define CGLM_VERSION_MINOR 9
-#define CGLM_VERSION_PATCH 3
+#define CGLM_VERSION_PATCH 4
 
 #endif /* cglm_version_h */

meson.build

@@ -1,5 +1,5 @@
 project('cglm', 'c',
-    version : '0.9.3',
+    version : '0.9.4',
     license : 'mit',
     default_options : [
         'c_std=c11',

src/mat2.c

@@ -100,6 +100,6 @@ glmc_mat2_rmc(vec2 r, mat2 m, vec2 c) {
 
 CGLM_EXPORT
 void
-glmc_mat2_make(float * __restrict src, mat2 dest) {
+glmc_mat2_make(const float * __restrict src, mat2 dest) {
   glm_mat2_make(src, dest);
 }

src/mat2x3.c

@@ -22,19 +22,19 @@ glmc_mat2x3_zero(mat2x3 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat2x3_make(float * __restrict src, mat2x3 dest) {
+glmc_mat2x3_make(const float * __restrict src, mat2x3 dest) {
   glm_mat2x3_make(src, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest) {
+glmc_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat3 dest) {
   glm_mat2x3_mul(m1, m2, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_mat2x3_mulv(mat2x3 m, vec3 v, vec2 dest) {
+glmc_mat2x3_mulv(mat2x3 m, vec2 v, vec3 dest) {
   glm_mat2x3_mulv(m, v, dest);
 }
 

src/mat2x4.c

@@ -22,19 +22,19 @@ glmc_mat2x4_zero(mat2x4 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat2x4_make(float * __restrict src, mat2x4 dest) {
+glmc_mat2x4_make(const float * __restrict src, mat2x4 dest) {
   glm_mat2x4_make(src, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest) {
+glmc_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat4 dest) {
   glm_mat2x4_mul(m1, m2, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_mat2x4_mulv(mat2x4 m, vec4 v, vec2 dest) {
+glmc_mat2x4_mulv(mat2x4 m, vec2 v, vec4 dest) {
   glm_mat2x4_mulv(m, v, dest);
 }
 

src/mat3.c

@@ -106,6 +106,6 @@ glmc_mat3_rmc(vec3 r, mat3 m, vec3 c) {
 
 CGLM_EXPORT
 void
-glmc_mat3_make(float * __restrict src, mat3 dest) {
+glmc_mat3_make(const float * __restrict src, mat3 dest) {
   glm_mat3_make(src, dest);
 }

src/mat3x2.c

@@ -22,19 +22,19 @@ glmc_mat3x2_zero(mat3x2 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat3x2_make(float * __restrict src, mat3x2 dest) {
+glmc_mat3x2_make(const float * __restrict src, mat3x2 dest) {
   glm_mat3x2_make(src, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest) {
+glmc_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat2 dest) {
   glm_mat3x2_mul(m1, m2, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_mat3x2_mulv(mat3x2 m, vec2 v, vec3 dest) {
+glmc_mat3x2_mulv(mat3x2 m, vec3 v, vec2 dest) {
   glm_mat3x2_mulv(m, v, dest);
 }
 

src/mat3x4.c

@@ -22,7 +22,7 @@ glmc_mat3x4_zero(mat3x4 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat3x4_make(float * __restrict src, mat3x4 dest) {
+glmc_mat3x4_make(const float * __restrict src, mat3x4 dest) {
   glm_mat3x4_make(src, dest);
 }
 

src/mat4.c

@@ -166,6 +166,6 @@ glmc_mat4_rmc(vec4 r, mat4 m, vec4 c) {
 
 CGLM_EXPORT
 void
-glmc_mat4_make(float * __restrict src, mat4 dest) {
+glmc_mat4_make(const float * __restrict src, mat4 dest) {
   glm_mat4_make(src, dest);
 }

src/mat4x2.c

@@ -22,7 +22,7 @@ glmc_mat4x2_zero(mat4x2 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat4x2_make(float * __restrict src, mat4x2 dest) {
+glmc_mat4x2_make(const float * __restrict src, mat4x2 dest) {
   glm_mat4x2_make(src, dest);
 }
 

src/mat4x3.c

@@ -22,7 +22,7 @@ glmc_mat4x3_zero(mat4x3 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat4x3_make(float * __restrict src, mat4x3 dest) {
+glmc_mat4x3_make(const float * __restrict src, mat4x3 dest) {
   glm_mat4x3_make(src, dest);
 }
 

src/quat.c

@@ -232,6 +232,6 @@ glmc_quat_rotate_atm(mat4 m, versor q, vec3 pivot) {
 
 CGLM_EXPORT
 void
-glmc_quat_make(float * __restrict src, versor dest) {
+glmc_quat_make(const float * __restrict src, versor dest) {
   glm_quat_make(src, dest);
 }

src/ray.c

@@ -9,5 +9,21 @@ glmc_ray_triangle(vec3   origin,
                   vec3   v1,
                   vec3   v2,
                   float *d) {
-    return glm_ray_triangle(origin, direction, v0, v1, v2, d);
+  return glm_ray_triangle(origin, direction, v0, v1, v2, d);
+}
+
+CGLM_EXPORT
+bool
+glmc_ray_sphere(vec3 origin,
+                vec3 dir,
+                vec4 s,
+                float * __restrict t1,
+                float * __restrict t2) {
+  return glm_ray_sphere(origin, dir, s, t1, t2);
+}
+
+CGLM_EXPORT
+void
+glmc_ray_at(vec3 orig, vec3 dir, float t, vec3 point) {
+  glm_ray_at(orig, dir, t, point);
 }

src/vec2.c

@@ -299,6 +299,18 @@ glmc_vec2_complex_conjugate(vec2 a, vec2 dest) {
 
 CGLM_EXPORT
 void
-glmc_vec2_make(float * __restrict src, vec2 dest) {
+glmc_vec2_make(const float * __restrict src, vec2 dest) {
   glm_vec2_make(src, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec2_reflect(vec2 v, vec2 n, vec2 dest) {
+  glm_vec2_reflect(v, n, dest);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec2_refract(vec2 v, vec2 n, float eta, vec2 dest) {
+  return glm_vec2_refract(v, n, eta, dest);
+}

src/vec3.c

@@ -456,6 +456,24 @@ glmc_vec3_sqrt(vec3 v, vec3 dest) {
 
 CGLM_EXPORT
 void
-glmc_vec3_make(float * __restrict src, vec3 dest) {
+glmc_vec3_make(const float * __restrict src, vec3 dest) {
   glm_vec3_make(src, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec3_faceforward(vec3 n, vec3 v, vec3 nref, vec3 dest) {
+  glm_vec3_faceforward(n, v, nref, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec3_reflect(vec3 v, vec3 n, vec3 dest) {
+  glm_vec3_reflect(v, n, dest);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec3_refract(vec3 v, vec3 n, float eta, vec3 dest) {
+  return glm_vec3_refract(v, n, eta, dest);
+}

src/vec4.c

@@ -420,6 +420,18 @@ glmc_vec4_sqrt(vec4 v, vec4 dest) {
 
 CGLM_EXPORT
 void
-glmc_vec4_make(float * __restrict src, vec4 dest) {
+glmc_vec4_make(const float * __restrict src, vec4 dest) {
   glm_vec4_make(src, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec4_reflect(vec4 v, vec4 n, vec4 dest) {
+  glm_vec4_reflect(v, n, dest);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_refract(vec4 v, vec4 n, float eta, vec4 dest) {
+  return glm_vec4_refract(v, n, eta, dest);
+}

test/src/test_mat2x3.h

@@ -83,42 +83,42 @@ TEST_IMPL(GLM_PREFIX, mat2x3_mul) {
   mat2x3 m1 = GLM_MAT2X3_ZERO_INIT;
   mat3x2 m2 = GLM_MAT3X2_ZERO_INIT;
 
-  mat2 m3 = GLM_MAT2_ZERO_INIT;
+  mat3 m3 = GLM_MAT3_ZERO_INIT;
-  mat2 m4 = GLM_MAT2_ZERO_INIT;
+  mat3 m4 = GLM_MAT3_ZERO_INIT;
 
-  int i, j, k;
+  int c, r, k;
 
   /* test random matrices */
   /* random matrices */
   test_rand_mat2x3(m1);
   test_rand_mat3x2(m2);
 
-  for (i = 0; i < 2; i++) {
+  for (r = 0; r < 3; r++) {
-    for (j = 0; j < 2; j++) {
+    for (c = 0; c < 3; c++) {
-      for (k = 0; k < 3; k++) {
+      for (k = 0; k < 2; k++) {
-        m4[i][j] += m1[i][k] * m2[k][j];
+        m4[c][r] += m1[k][r] * m2[c][k];
       }
     }
   }
 
   GLM(mat2x3_mul)(m1, m2, m3);
-  ASSERTIFY(test_assert_mat2_eq(m3, m4))
+  ASSERTIFY(test_assert_mat3_eq(m3, m4))
 
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, mat2x3_mulv) {
   mat2x3 mat = A_MATRIX2X3;
-  vec3 v = {11.0f, 21.0f, 31.0f};
+  vec2 v = {11.0f, 21.0f};
 
   int  i;
-  vec2 dest;
+  vec3 dest;
   float res = 0.0;
 
   GLM(mat2x3_mulv)(mat, v, dest);
 
-  for (i = 0; i < 2; i++) {
+  for (i = 0; i < 3; i++) {
-    res = mat[i][0] * v[0] + mat[i][1] * v[1] + mat[i][2] * v[2];
+    res = mat[0][i] * v[0] + mat[1][i] * v[1];
     ASSERT(test_eq(dest[i], res))
   }
 

test/src/test_mat2x4.h

@@ -86,42 +86,42 @@ TEST_IMPL(GLM_PREFIX, mat2x4_mul) {
   mat2x4 m1 = GLM_MAT2X4_ZERO_INIT;
   mat4x2 m2 = GLM_MAT4X2_ZERO_INIT;
 
-  mat2 m3 = GLM_MAT2_ZERO_INIT;
+  mat4 m3 = GLM_MAT4_ZERO_INIT;
-  mat2 m4 = GLM_MAT2_ZERO_INIT;
+  mat4 m4 = GLM_MAT4_ZERO_INIT;
 
-  int i, j, k;
+  int c, r, k;
 
   /* test random matrices */
   /* random matrices */
   test_rand_mat2x4(m1);
   test_rand_mat4x2(m2);
 
-  for (i = 0; i < 2; i++) {
+  for (r = 0; r < 4; r++) {
-    for (j = 0; j < 2; j++) {
+    for (c = 0; c < 4; c++) {
-      for (k = 0; k < 4; k++) {
+      for (k = 0; k < 2; k++) {
-        m4[i][j] += m1[i][k] * m2[k][j];
+        m4[c][r] += m1[k][r] * m2[c][k];
       }
     }
   }
 
   GLM(mat2x4_mul)(m1, m2, m3);
-  ASSERTIFY(test_assert_mat2_eq(m3, m4))
+  ASSERTIFY(test_assert_mat4_eq(m3, m4))
 
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, mat2x4_mulv) {
   mat2x4 mat = A_MATRIX2X4;
-  vec4 v = {11.0f, 21.0f, 31.0f, 41.0f};
+  vec2 v = {11.0f, 21.0f};
 
   int  i;
-  vec2 dest;
+  vec4 dest;
   float res = 0.0;
 
   GLM(mat2x4_mulv)(mat, v, dest);
 
-  for (i = 0; i < 2; i++) {
+  for (i = 0; i < 4; i++) {
-    res = mat[i][0] * v[0] + mat[i][1] * v[1] + mat[i][2] * v[2] + mat[i][3] * v[3];
+    res = mat[0][i] * v[0] + mat[1][i] * v[1];
     ASSERT(test_eq(dest[i], res))
   }
 

test/src/test_mat3x2.h

@@ -84,40 +84,40 @@ TEST_IMPL(GLM_PREFIX, mat3x2_mul) {
   mat3x2 m1 = GLM_MAT3X2_ZERO_INIT;
   mat2x3 m2 = GLM_MAT2X3_ZERO_INIT;
 
-  mat3 m3 = GLM_MAT3_ZERO_INIT;
+  mat2 m3 = GLM_MAT2_ZERO_INIT;
-  mat3 m4 = GLM_MAT3_ZERO_INIT;
+  mat2 m4 = GLM_MAT2_ZERO_INIT;
 
-  int i, j, k;
+  int c, r, k;
 
   test_rand_mat3x2(m1);
   test_rand_mat2x3(m2);
 
-  for (i = 0; i < 3; i++) {
+  for (r = 0; r < 2; r++) {
-    for (j = 0; j < 3; j++) {
+    for (c = 0; c < 2; c++) {
-      for (k = 0; k < 2; k++) {
+      for (k = 0; k < 3; k++) {
-        m4[i][j] += m1[i][k] * m2[k][j];
+        m4[c][r] += m1[k][r] * m2[c][k];
       }
     }
   }
 
   GLM(mat3x2_mul)(m1, m2, m3);
-  ASSERTIFY(test_assert_mat3_eq(m3, m4))
+  ASSERTIFY(test_assert_mat2_eq(m3, m4))
 
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, mat3x2_mulv) {
   mat3x2 mat = A_MATRIX3X2;
-  vec2 v = {11.0f, 21.0f};
+  vec3 v = {11.0f, 21.0f, 31.0f};
 
   int  i;
-  vec3 dest;
+  vec2 dest;
   float res = 0.0;
 
   GLM(mat3x2_mulv)(mat, v, dest);
 
-  for (i = 0; i < 3; i++) {
+  for (i = 0; i < 2; i++) {
-    res = mat[i][0] * v[0] + mat[i][1] * v[1];
+    res = mat[0][i] * v[0] + mat[1][i] * v[1] + mat[2][i] * v[2];
     ASSERT(test_eq(dest[i], res))
   }
 

test/src/test_project.h

@@ -26,9 +26,15 @@ TEST_IMPL(GLM_PREFIX, unprojecti) {
 
   /* unprojected of projected vector must be same as original one */
   /* we used 0.01 because of projection floating point errors */
+#ifndef CGLM_FAST_MATH
   ASSERT(fabsf(pos[0] - unprojected[0]) < 0.01)
   ASSERT(fabsf(pos[1] - unprojected[1]) < 0.01)
   ASSERT(fabsf(pos[2] - unprojected[2]) < 0.01)
+#else
+  ASSERT(fabsf(pos[0] - unprojected[0]) < 0.1)
+  ASSERT(fabsf(pos[1] - unprojected[1]) < 0.1)
+  ASSERT(fabsf(pos[2] - unprojected[2]) < 0.1)
+#endif
 
   TEST_SUCCESS
 }
@@ -50,9 +56,16 @@ TEST_IMPL(GLM_PREFIX, unproject) {
 
   /* unprojected of projected vector must be same as original one */
   /* we used 0.01 because of projection floating point errors */
+
+#ifndef CGLM_FAST_MATH
   ASSERT(fabsf(pos[0] - unprojected[0]) < 0.01)
   ASSERT(fabsf(pos[1] - unprojected[1]) < 0.01)
   ASSERT(fabsf(pos[2] - unprojected[2]) < 0.01)
+#else
+  ASSERT(fabsf(pos[0] - unprojected[0]) < 0.1)
+  ASSERT(fabsf(pos[1] - unprojected[1]) < 0.1)
+  ASSERT(fabsf(pos[2] - unprojected[2]) < 0.1)
+#endif
 
   TEST_SUCCESS
 }
@@ -74,9 +87,16 @@ TEST_IMPL(GLM_PREFIX, project) {
 
   /* unprojected of projected vector must be same as original one */
   /* we used 0.01 because of projection floating point errors */
+
+#ifndef CGLM_FAST_MATH
   ASSERT(fabsf(pos[0] - unprojected[0]) < 0.01)
   ASSERT(fabsf(pos[1] - unprojected[1]) < 0.01)
   ASSERT(fabsf(pos[2] - unprojected[2]) < 0.01)
+#else
+  ASSERT(fabsf(pos[0] - unprojected[0]) < 0.1)
+  ASSERT(fabsf(pos[1] - unprojected[1]) < 0.1)
+  ASSERT(fabsf(pos[2] - unprojected[2]) < 0.1)
+#endif
 
   /* test with no projection */
   glm_mat4_identity(mvp);

test/src/test_ray.h

@@ -32,3 +32,43 @@ TEST_IMPL(GLM_PREFIX, ray_triangle) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, ray_sphere) {
+  vec4  sphere = {5.0f, 0.0f, 0.0f, 1.0f}; /* Sphere: center at (5, 0, 0) with radius 1 */
+  float t1, t2;
+  bool  hit;
+
+  /* Case 1: Ray misses the sphere */
+  hit = GLM(ray_sphere)((vec3){10.0f, 3.0f, 0.0f}, (vec3){1.0f, 0.0f, 0.0f}, sphere, &t1, &t2);
+  ASSERT(!hit); /* Expect no intersection */
+
+  /* Case 2: Ray starts inside the sphere */
+  hit = GLM(ray_sphere)((vec3){5.0f, 0.5f, 0.0f}, (vec3){1.0f, 0.0f, 0.0f}, sphere, &t1, &t2);
+  ASSERT(hit); /* Expect an intersection */
+  ASSERT(t1 < 0 && t2 > 0); /* Ray exits at t2 */
+
+  /* Case 3: Ray intersects the sphere from outside */
+  hit = GLM(ray_sphere)((vec3){0.0f, 0.0f, 0.0f}, (vec3){1.0f, 0.0f, 0.0f}, sphere, &t1, &t2);
+  ASSERT(hit); /* Expect an intersection */
+  ASSERT(t1 > 0 && t2 > 0); /* Intersections at t1 and t2 */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, ray_at) {
+  vec3 origin    = {0.0f, 0.0f, 0.0f};
+  vec3 direction = {1.0f, 1.0f, 1.0f};  /* Diagonal direction */
+  float distance = sqrtf(3.0f);         /* Distance along the ray; sqrt(3) for unit length due to direction normalization */
+  vec3 result;
+
+  /* Normalize the direction to ensure accurate distance measurement */
+  glm_vec3_normalize(direction);
+
+  GLM(ray_at)(origin, direction, distance, result);
+  ASSERT(fabsf(result[0] - 1.0f) <= 0.0000009);  /* Expecting to be 1 unit along the x-axis */
+  ASSERT(fabsf(result[1] - 1.0f) <= 0.0000009);  /* Expecting to be 1 unit along the y-axis */
+  ASSERT(fabsf(result[2] - 1.0f) <= 0.0000009);  /* Expecting to be 1 unit along the z-axis */
+
+  TEST_SUCCESS
+}
+

test/src/test_vec2.h

@@ -241,7 +241,7 @@ TEST_IMPL(GLM_PREFIX, vec2_scale_as) {
   GLM(vec2_scale_as)(v1, s, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
 
@@ -492,7 +492,7 @@ TEST_IMPL(GLM_PREFIX, vec2_normalize) {
   GLM(vec2_normalize)(v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
 
@@ -519,7 +519,7 @@ TEST_IMPL(GLM_PREFIX, vec2_normalize_to) {
   GLM(vec2_normalize_to)(v1, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
 
@@ -752,3 +752,77 @@ TEST_IMPL(GLM_PREFIX, vec2_make) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, vec2_reflect) {
+  vec2 dest;
+
+  /* Reflecting off a "horizontal" surface in 2D */
+  vec2 I1 = {1.0f, -1.0f}; /* Incoming vector */
+  vec2 N1 = {0.0f, 1.0f};  /* Normal vector */
+  GLM(vec2_reflect)(I1, N1, dest);
+  ASSERT(fabsf(dest[0] - 1.0f) < 0.00001f &&
+         fabsf(dest[1] - 1.0f) < 0.00001f); /* Expect reflection upwards */
+
+  /* Reflecting at an angle in 2D */
+  vec2 I2 = {sqrtf(2)/2, -sqrtf(2)/2}; /* Incoming vector at 45 degrees */
+  vec2 N2 = {0.0f, 1.0f}; /* Upwards normal vector */
+  GLM(vec2_reflect)(I2, N2, dest);
+  ASSERT(fabsf(dest[0] - sqrtf(2)/2) < 0.00001f &&
+         fabsf(dest[1] - sqrtf(2)/2) < 0.00001f); /* Expect reflection upwards */
+
+  /* Reflecting off a line in 2D representing a "vertical" surface analogy */
+  vec2 I3 = {1.0f, 0.0f};  /* Incoming vector */
+  vec2 N3 = {-1.0f, 0.0f}; /* Normal vector representing a "vertical" line */
+  GLM(vec2_reflect)(I3, N3, dest);
+  ASSERT(fabsf(dest[0] + 1.0f) < 0.00001f &&
+         fabsf(dest[1]) < 0.00001f); /* Expect reflection to the left */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec2_refract) {
+  vec2  v = {sqrtf(0.5f), -sqrtf(0.5f)}; /* Incoming vector at 45 degrees to normal */
+  vec2  N = {0.0f, 1.0f};                /* Surface normal */
+  vec2  dest;
+  float eta;
+  float r;
+
+  /* Water to Air (eta = 1.33/1.0) */
+  eta = 1.33f / 1.0f;
+  r = GLM(vec2_refract)(v, N, eta, dest);
+  // In 2D, we expect a similar bending behavior as in 3D, so we check dest[1]
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f)) {
+    ASSERT(dest[1] < -sqrtf(0.5f)); // Refracted ray bends away from the normal
+    ASSERT(r == true);
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f); // Total internal reflection
+    ASSERT(r == false);
+  }
+
+  /* Air to Glass (eta = 1.0 / 1.5) */
+  eta = 1.0f / 1.5f;
+  r = GLM(vec2_refract)(v, N, eta, dest);
+  ASSERT(r == true);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal
+
+  /* Glass to Water (eta = 1.5 / 1.33) */
+  eta = 1.5f / 1.33f;
+  r = GLM(vec2_refract)(v, N, eta, dest);
+  ASSERT(r == true);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal, less bending than air to glass
+
+  /* Diamond to Air (eta = 2.42 / 1.0) */
+  eta = 2.42f / 1.0f;
+  r = GLM(vec2_refract)(v, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f)) {
+    /* High potential for total internal reflection, but if it occurs, expect significant bending */
+    ASSERT(dest[1] < -sqrtf(0.5f));
+    ASSERT(r == true);
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f); // Total internal reflection
+    ASSERT(r == false);
+  }
+
+  TEST_SUCCESS
+}
+

test/src/test_vec3.h

@@ -433,7 +433,7 @@ TEST_IMPL(GLM_PREFIX, vec3_scale_as) {
   GLM(vec3_scale_as)(v1, s, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -704,7 +704,7 @@ TEST_IMPL(GLM_PREFIX, vec3_normalize) {
   GLM(vec3_normalize)(v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -733,7 +733,7 @@ TEST_IMPL(GLM_PREFIX, vec3_normalize_to) {
   GLM(vec3_normalize_to)(v1, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -764,7 +764,7 @@ TEST_IMPL(GLM_PREFIX, normalize) {
   GLM(vec3_normalize)(v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -795,7 +795,7 @@ TEST_IMPL(GLM_PREFIX, normalize_to) {
   GLM(vec3_normalize_to)(v1, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -1673,14 +1673,16 @@ TEST_IMPL(GLM_PREFIX, vec3_eqv_eps) {
 
 TEST_IMPL(GLM_PREFIX, vec3_max) {
   vec3 v1 = {2.104f, -3.012f, -4.10f}, v2 = {-12.35f, -31.140f, -43.502f};
-  vec3 v3 = {INFINITY, 0.0f, 0.0f}, v4 = {NAN, INFINITY, 2.0f};
+  vec3 v3 = {INFINITY, 0.0f, 0.0f}/*, v4 = {NAN, INFINITY, 2.0f}*/;
-  vec3 v5 = {NAN, -1.0f, -1.0f}, v6 = {-1.0f, -11.0f, 11.0f};
+  vec3 /*v5 = {NAN, -1.0f, -1.0f}, */v6 = {-1.0f, -11.0f, 11.0f};
 
   ASSERT(test_eq(GLM(vec3_max)(v1),  2.104f))
   ASSERT(test_eq(GLM(vec3_max)(v2), -12.35f))
+#ifndef CGLM_FAST_MATH
   ASSERT(isinf(GLM(vec3_max)(v3)))
-  ASSERT(isnan(GLM(vec3_max)(v4)))
+#endif
-  ASSERT(isnan(GLM(vec3_max)(v5)))
+//  ASSERT(isnan(GLM(vec3_max)(v4)))
+//  ASSERT(isnan(GLM(vec3_max)(v5)))
   ASSERT(test_eq(GLM(vec3_max)(v6),  11.0f))
 
   TEST_SUCCESS
@@ -1688,20 +1690,21 @@ TEST_IMPL(GLM_PREFIX, vec3_max) {
 
 TEST_IMPL(GLM_PREFIX, vec3_min) {
   vec3  v1 = {2.104f, -3.012f, -4.10f}, v2 = {-12.35f, -31.140f, -43.502f};
-  vec3  v3 = {INFINITY, 0.0f, 0.0f}, v4 = {NAN, INFINITY, 2.0f};
+  vec3  v3 = {INFINITY, 0.0f, 0.0f}/*, v4 = {NAN, INFINITY, 2.0f}*/;
-  vec3  v5 = {NAN, -1.0f, -1.0f}, v6 = {-1.0f, -11.0f, 11.0f};
+  vec3  /*v5 = {NAN, -1.0f, -1.0f},*/ v6 = {-1.0f, -11.0f, 11.0f};
 
   ASSERT(test_eq(GLM(vec3_min)(v1), -4.10f))
   ASSERT(test_eq(GLM(vec3_min)(v2), -43.502f))
   ASSERT(test_eq(GLM(vec3_min)(v3),  0.0f))
-  ASSERT(isnan(GLM(vec3_min)(v4)))
+//  ASSERT(isnan(GLM(vec3_min)(v4)))
-  ASSERT(isnan(GLM(vec3_min)(v5)))
+//  ASSERT(isnan(GLM(vec3_min)(v5)))
   ASSERT(test_eq(GLM(vec3_min)(v6), -11.0f))
 
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, vec3_isnan) {
+#ifndef CGLM_FAST_MATH
   vec3  v1 = {2.104f, -3.012f, -4.10f}, v2 = {-12.35f, -31.140f, -43.502f};
   vec3  v3 = {INFINITY, 0.0f, 0.0f}, v4 = {NAN, INFINITY, 2.0f};
   vec3  v5 = {NAN, -1.0f, -1.0f}, v6 = {-1.0f, -1.0f, 11.0f};
@@ -1712,11 +1715,12 @@ TEST_IMPL(GLM_PREFIX, vec3_isnan) {
   ASSERT(GLM(vec3_isnan)(v4))
   ASSERT(GLM(vec3_isnan)(v5))
   ASSERT(!GLM(vec3_isnan)(v6))
-
+#endif
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, vec3_isinf) {
+#ifndef CGLM_FAST_MATH
   vec3  v1 = {2.104f, -3.012f, -4.10f}, v2 = {-12.35f, -31.140f, -43.502f};
   vec3  v3 = {INFINITY, 0.0f, 0.0f}, v4 = {NAN, INFINITY, 2.0f};
   vec3  v5 = {NAN, -1.0f, -1.0f}, v6 = {-1.0f, -1.0f, 11.0f};
@@ -1727,11 +1731,12 @@ TEST_IMPL(GLM_PREFIX, vec3_isinf) {
   ASSERT(GLM(vec3_isinf)(v4))
   ASSERT(!GLM(vec3_isinf)(v5))
   ASSERT(!GLM(vec3_isinf)(v6))
-
+#endif
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, vec3_isvalid) {
+#ifndef CGLM_FAST_MATH
   vec3  v1 = {2.104f, -3.012f, -4.10f}, v2 = {-12.35f, -31.140f, -43.502f};
   vec3  v3 = {INFINITY, 0.0f, 0.0f}, v4 = {NAN, INFINITY, 2.0f};
   vec3  v5 = {NAN, -1.0f, -1.0f}, v6 = {-1.0f, -1.0f, 11.0f};
@@ -1742,7 +1747,7 @@ TEST_IMPL(GLM_PREFIX, vec3_isvalid) {
   ASSERT(!GLM(vec3_isvalid)(v4))
   ASSERT(!GLM(vec3_isvalid)(v5))
   ASSERT(GLM(vec3_isvalid)(v6))
-
+#endif
   TEST_SUCCESS
 }
 
@@ -1840,3 +1845,96 @@ TEST_IMPL(GLM_PREFIX, vec3_make) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, vec3_faceforward) {
+  vec3 N = {0.0f, 1.0f, 0.0f};
+  vec3 v = {1.0f, -1.0f, 0.0f};
+  vec3 Nref = {0.0f, -1.0f, 0.0f};
+  vec3 dest;
+
+  GLM(vec3_faceforward)(N, v, Nref, dest);
+  ASSERT(dest[0] == 0.0f
+         && dest[1] == -1.0f
+         && dest[2] == 0.0f); /* Expect N flipped */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec3_reflect) {
+  vec3 dest;
+
+  /* Original test: Reflecting off a horizontal surface */
+  vec3 I1 = {1.0f, -1.0f, 0.0f}; /* Incoming vector */
+  vec3 N1 = {0.0f, 1.0f, 0.0f};  /* Normal vector   */
+  GLM(vec3_reflect)(I1, N1, dest);
+  ASSERT(fabsf(dest[0] - 1.0f) < 0.00001f
+         && fabsf(dest[1] - 1.0f) < 0.00001f
+         && fabsf(dest[2] - 0.0f) < 0.00001f); /* Expect reflection */
+
+  /* Scenario 2: Reflecting off a vertical surface  */
+  vec3 I2 = {1.0f, 0.0f, 0.0f};  /* Incoming vector */
+  vec3 N2 = {-1.0f, 0.0f, 0.0f}; /* Normal vector   */
+  GLM(vec3_reflect)(I2, N2, dest);
+  ASSERT(fabsf(dest[0] + 1.0f) < 0.00001f
+         && fabsf(dest[1]) < 0.00001f
+         && fabsf(dest[2]) < 0.00001f); /* Expect reflection to the left */
+
+  /* Scenario 3: Reflecting at an angle */
+  vec3 I3 = {sqrtf(2)/2, -sqrtf(2)/2, 0.0f}; /* Incoming vector at 45 degrees */
+  vec3 N3 = {0.0f, 1.0f, 0.0f}; /* Upwards normal vector */
+  GLM(vec3_reflect)(I3, N3, dest);
+  ASSERT(fabsf(dest[0] - sqrtf(2)/2) < 0.00001f
+         && fabsf(dest[1] - sqrtf(2)/2) < 0.00001f
+         && fabsf(dest[2]) < 0.00001f); /* Expect reflection upwards */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec3_refract) {
+  vec3  v = {sqrtf(0.5f), -sqrtf(0.5f), 0.0f}; /* Incoming vector at 45 degrees to normal */
+  vec3  N = {0.0f, 1.0f, 0.0f};                /* Surface normal */
+  vec3  dest;
+  float eta;
+  bool  r;
+
+  /* Water to Air (eta = 1.33/1.0) */
+  eta = 1.33f / 1.0f;
+  r = GLM(vec3_refract)(v, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f)) {
+    ASSERT(dest[1] < -sqrtf(0.5f));
+    ASSERT(r == true);
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f);
+    ASSERT(r == false);
+  }
+
+  /* Air to Glass (eta = 1.0 / 1.5) */
+  eta = 1.0f / 1.5f;
+  r = GLM(vec3_refract)(v, N, eta, dest);
+
+  /* Expect bending towards the normal */
+  ASSERT(r == true);
+  ASSERT(dest[1] < -sqrtf(0.5f));
+
+  /* Glass to Water (eta = 1.5 / 1.33) */
+  eta = 1.5f / 1.33f;
+  r = GLM(vec3_refract)(v, N, eta, dest);
+
+  /* Expect bending towards the normal, less bending than air to glass */
+  ASSERT(r == true);
+  ASSERT(dest[1] < -sqrtf(0.5f));
+
+  /* Diamond to Air (eta = 2.42 / 1.0) */
+  eta = 2.42f / 1.0f;
+  r = GLM(vec3_refract)(v, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f)) {
+    /* High potential for total internal reflection, but if it occurs, expect significant bending */
+    ASSERT(dest[1] < -sqrtf(0.5f));
+    ASSERT(r == true);
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f);
+    ASSERT(r == false);
+  }
+
+  TEST_SUCCESS
+}

test/src/test_vec4.h

@@ -410,7 +410,7 @@ TEST_IMPL(GLM_PREFIX, vec4_scale_as) {
   GLM(vec4_scale_as)(v1, s, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2] + v1[3] * v1[3]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -701,7 +701,7 @@ TEST_IMPL(GLM_PREFIX, vec4_normalize) {
   GLM(vec4_normalize)(v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2] + v1[3] * v1[3]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -732,7 +732,7 @@ TEST_IMPL(GLM_PREFIX, vec4_normalize_to) {
   GLM(vec4_normalize_to)(v1, v2);
 
   norm = sqrtf(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2] + v1[3] * v1[3]);
-  if (norm == 0.0f) {
+  if (norm < FLT_EPSILON) {
     ASSERT(test_eq(v1[0], 0.0f))
     ASSERT(test_eq(v1[1], 0.0f))
     ASSERT(test_eq(v1[2], 0.0f))
@@ -1345,15 +1345,17 @@ TEST_IMPL(GLM_PREFIX, vec4_max) {
   vec4 v1 = {2.104f, -3.012f, -4.10f, -4.10f};
   vec4 v2 = {-12.35f, -31.140f, -43.502f, -43.502f};
   vec4 v3 = {INFINITY, 0.0f, 0.0f, 0.0f};
-  vec4 v4 = {NAN, INFINITY, 2.0f, 2.0f};
+//  vec4 v4 = {NAN, INFINITY, 2.0f, 2.0f};
-  vec4 v5 = {NAN, -1.0f, -1.0f, -1.0f};
+//  vec4 v5 = {NAN, -1.0f, -1.0f, -1.0f};
   vec4 v6 = {-1.0f, -11.0f, 11.0f, 11.0f};
 
   ASSERT(test_eq(GLM(vec4_max)(v1),  2.104f))
   ASSERT(test_eq(GLM(vec4_max)(v2), -12.35f))
+#ifndef CGLM_FAST_MATH
   ASSERT(isinf(GLM(vec4_max)(v3)))
-  ASSERT(isnan(GLM(vec4_max)(v4)))
+#endif
-  ASSERT(isnan(GLM(vec4_max)(v5)))
+//  ASSERT(isnan(GLM(vec4_max)(v4)))
+//  ASSERT(isnan(GLM(vec4_max)(v5)))
   ASSERT(test_eq(GLM(vec4_max)(v6),  11.0f))
 
   TEST_SUCCESS
@@ -1363,21 +1365,22 @@ TEST_IMPL(GLM_PREFIX, vec4_min) {
   vec4 v1 = {2.104f, -3.012f, -4.10f, -4.10f};
   vec4 v2 = {-12.35f, -31.140f, -43.502f, -43.502f};
   vec4 v3 = {INFINITY, 0.0f, 0.0f, 0.0f};
-  vec4 v4 = {NAN, INFINITY, 2.0f, 2.0f};
+//  vec4 v4 = {NAN, INFINITY, 2.0f, 2.0f};
-  vec4 v5 = {NAN, -1.0f, -1.0f, -1.0f};
+//  vec4 v5 = {NAN, -1.0f, -1.0f, -1.0f};
   vec4 v6 = {-1.0f, -11.0f, 11.0f, 11.0f};
 
   ASSERT(test_eq(GLM(vec4_min)(v1), -4.10f))
   ASSERT(test_eq(GLM(vec4_min)(v2), -43.502f))
   ASSERT(test_eq(GLM(vec4_min)(v3),  0.0f))
-  ASSERT(isnan(GLM(vec4_min)(v4)))
+//  ASSERT(isnan(GLM(vec4_min)(v4)))
-  ASSERT(isnan(GLM(vec4_min)(v5)))
+//  ASSERT(isnan(GLM(vec4_min)(v5)))
   ASSERT(test_eq(GLM(vec4_min)(v6), -11.0f))
 
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, vec4_isnan) {
+#ifndef CGLM_FAST_MATH
   vec4 v1 = {2.104f, -3.012f, -4.10f, -4.10f};
   vec4 v2 = {-12.35f, -31.140f, -43.502f, -43.502f};
   vec4 v3 = {INFINITY, 0.0f, 0.0f, 0.0f};
@@ -1391,11 +1394,12 @@ TEST_IMPL(GLM_PREFIX, vec4_isnan) {
   ASSERT(GLM(vec4_isnan)(v4))
   ASSERT(GLM(vec4_isnan)(v5))
   ASSERT(!GLM(vec4_isnan)(v6))
-
+#endif
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, vec4_isinf) {
+#ifndef CGLM_FAST_MATH
   vec4 v1 = {2.104f, -3.012f, -4.10f, -4.10f};
   vec4 v2 = {-12.35f, -31.140f, -43.502f, -43.502f};
   vec4 v3 = {INFINITY, 0.0f, 0.0f, 0.0f};
@@ -1409,11 +1413,12 @@ TEST_IMPL(GLM_PREFIX, vec4_isinf) {
   ASSERT(GLM(vec4_isinf)(v4))
   ASSERT(!GLM(vec4_isinf)(v5))
   ASSERT(!GLM(vec4_isinf)(v6))
-
+#endif
   TEST_SUCCESS
 }
 
 TEST_IMPL(GLM_PREFIX, vec4_isvalid) {
+#ifndef CGLM_FAST_MATH
   vec4 v1 = {2.104f, -3.012f, -4.10f, -4.10f};
   vec4 v2 = {-12.35f, -31.140f, -43.502f, -43.502f};
   vec4 v3 = {INFINITY, 0.0f, 0.0f, 0.0f};
@@ -1427,7 +1432,7 @@ TEST_IMPL(GLM_PREFIX, vec4_isvalid) {
   ASSERT(!GLM(vec4_isvalid)(v4))
   ASSERT(!GLM(vec4_isvalid)(v5))
   ASSERT(GLM(vec4_isvalid)(v6))
-
+#endif
   TEST_SUCCESS
 }
 
@@ -1536,3 +1541,82 @@ TEST_IMPL(GLM_PREFIX, vec4_make) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, vec4_reflect) {
+  vec4 dest;
+
+  /* Original test: Reflecting off a horizontal surface */
+  vec4 I1 = {1.0f, -1.0f, 0.0f, 0.0f}; /* Incoming vector */
+  vec4 N1 = {0.0f, 1.0f, 0.0f, 0.0f};  /* Normal vector   */
+  GLM(vec4_reflect)(I1, N1, dest);
+  ASSERT(fabsf(dest[0] - 1.0f) < 0.00001f &&
+         fabsf(dest[1] - 1.0f) < 0.00001f &&
+         fabsf(dest[2] - 0.0f) < 0.00001f &&
+         fabsf(dest[3] - 0.0f) < 0.00001f); /* Expect reflection */
+
+  /* Scenario 2: Reflecting off a vertical surface */
+  vec4 I2 = {1.0f, 0.0f, 0.0f, 0.0f};  /* Incoming vector */
+  vec4 N2 = {-1.0f, 0.0f, 0.0f, 0.0f}; /* Normal vector   */
+  GLM(vec4_reflect)(I2, N2, dest);
+  ASSERT(fabsf(dest[0] + 1.0f) < 0.00001f &&
+         fabsf(dest[1]) < 0.00001f &&
+         fabsf(dest[2]) < 0.00001f &&
+         fabsf(dest[3] - 0.0f) < 0.00001f); /* Expect reflection to the left */
+
+  /* Scenario 3: Reflecting at an angle */
+  vec4 I3 = {sqrtf(2)/2, -sqrtf(2)/2, 0.0f, 0.0f}; /* Incoming vector at 45 degrees */
+  vec4 N3 = {0.0f, 1.0f, 0.0f, 0.0f}; /* Upwards normal vector */
+  GLM(vec4_reflect)(I3, N3, dest);
+  ASSERT(fabsf(dest[0] - sqrtf(2)/2) < 0.00001f &&
+         fabsf(dest[1] - sqrtf(2)/2) < 0.00001f &&
+         fabsf(dest[2]) < 0.00001f &&
+         fabsf(dest[3] - 0.0f) < 0.00001f); /* Expect reflection upwards */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec4_refract) {
+  vec4  v = {sqrtf(0.5f), -sqrtf(0.5f), 0.0f, 0.0f}; /* Incoming vector */
+  vec4  N = {0.0f, 1.0f, 0.0f, 0.0f};                /* Surface normal */
+  vec4  dest;
+  float eta;
+  float r;
+
+  /* Water to Air (eta = 1.33/1.0) */
+  eta = 1.33f / 1.0f;
+  r = GLM(vec4_refract)(v, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f)) {
+    ASSERT(dest[1] < -sqrtf(0.5f));
+    ASSERT(r == true);
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f);
+    ASSERT(r == false);
+  }
+
+  /* Air to Glass (eta = 1.0 / 1.5) */
+  eta = 1.0f / 1.5f;
+  r = GLM(vec4_refract)(v, N, eta, dest);
+  ASSERT(r == true);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal
+
+  /* Glass to Water (eta = 1.5 / 1.33) */
+  eta = 1.5f / 1.33f;
+  r = GLM(vec4_refract)(v, N, eta, dest);
+  ASSERT(r == true);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal, less bending than air to glass
+
+  /* Diamond to Air (eta = 2.42 / 1.0) */
+  eta = 2.42f / 1.0f;
+  r = GLM(vec4_refract)(v, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f)) {
+    /* High potential for total internal reflection, but if it occurs, expect significant bending */
+    ASSERT(dest[1] < -sqrtf(0.5f));
+    ASSERT(r == true);
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f);
+    ASSERT(r == false);
+  }
+
+  TEST_SUCCESS
+}
+

test/tests.h

@@ -392,7 +392,12 @@ TEST_DECLARE(euler)
 
 /* ray */
 TEST_DECLARE(glm_ray_triangle)
+TEST_DECLARE(glm_ray_sphere)
+TEST_DECLARE(glm_ray_at)
+
 TEST_DECLARE(glmc_ray_triangle)
+TEST_DECLARE(glmc_ray_sphere)
+TEST_DECLARE(glmc_ray_at)
 
 /* quat */
 TEST_DECLARE(MACRO_GLM_QUAT_IDENTITY_INIT)
@@ -530,6 +535,8 @@ TEST_DECLARE(glm_vec2_lerp)
 TEST_DECLARE(glm_vec2_complex_mul)
 TEST_DECLARE(glm_vec2_complex_div)
 TEST_DECLARE(glm_vec2_make)
+TEST_DECLARE(glm_vec2_reflect)
+TEST_DECLARE(glm_vec2_refract)
 
 TEST_DECLARE(glmc_vec2)
 TEST_DECLARE(glmc_vec2_copy)
@@ -576,6 +583,8 @@ TEST_DECLARE(glmc_vec2_lerp)
 TEST_DECLARE(glmc_vec2_complex_mul)
 TEST_DECLARE(glmc_vec2_complex_div)
 TEST_DECLARE(glmc_vec2_make)
+TEST_DECLARE(glmc_vec2_reflect)
+TEST_DECLARE(glmc_vec2_refract)
 
 /* vec3 */
 TEST_DECLARE(MACRO_GLM_VEC3_ONE_INIT)
@@ -678,6 +687,9 @@ TEST_DECLARE(glm_vec3_fract)
 TEST_DECLARE(glm_vec3_hadd)
 TEST_DECLARE(glm_vec3_sqrt)
 TEST_DECLARE(glm_vec3_make)
+TEST_DECLARE(glm_vec3_faceforward)
+TEST_DECLARE(glm_vec3_reflect)
+TEST_DECLARE(glm_vec3_refract)
 
 TEST_DECLARE(glmc_vec3)
 TEST_DECLARE(glmc_vec3_copy)
@@ -754,6 +766,9 @@ TEST_DECLARE(glmc_vec3_fract)
 TEST_DECLARE(glmc_vec3_hadd)
 TEST_DECLARE(glmc_vec3_sqrt)
 TEST_DECLARE(glmc_vec3_make)
+TEST_DECLARE(glmc_vec3_faceforward)
+TEST_DECLARE(glmc_vec3_reflect)
+TEST_DECLARE(glmc_vec3_refract)
 
 /* vec4 */
 TEST_DECLARE(MACRO_GLM_VEC4_ONE_INIT)
@@ -842,6 +857,8 @@ TEST_DECLARE(glm_vec4_fract)
 TEST_DECLARE(glm_vec4_hadd)
 TEST_DECLARE(glm_vec4_sqrt)
 TEST_DECLARE(glm_vec4_make)
+TEST_DECLARE(glm_vec4_reflect)
+TEST_DECLARE(glm_vec4_refract)
 
 TEST_DECLARE(glmc_vec4)
 TEST_DECLARE(glmc_vec4_copy3)
@@ -914,6 +931,8 @@ TEST_DECLARE(glmc_vec4_fract)
 TEST_DECLARE(glmc_vec4_hadd)
 TEST_DECLARE(glmc_vec4_sqrt)
 TEST_DECLARE(glmc_vec4_make)
+TEST_DECLARE(glmc_vec4_reflect)
+TEST_DECLARE(glmc_vec4_refract)
 
 /* ivec2 */
 TEST_DECLARE(glm_ivec2)
@@ -1475,13 +1494,13 @@ TEST_LIST {
   TEST_ENTRY(glm_mat2x4_transpose)
   TEST_ENTRY(glm_mat2x4_scale)
 
-  TEST_ENTRY(glm_mat2x4_copy)
+  TEST_ENTRY(glmc_mat2x4_copy)
-  TEST_ENTRY(glm_mat2x4_zero)
+  TEST_ENTRY(glmc_mat2x4_zero)
-  TEST_ENTRY(glm_mat2x4_make)
+  TEST_ENTRY(glmc_mat2x4_make)
-  TEST_ENTRY(glm_mat2x4_mul)
+  TEST_ENTRY(glmc_mat2x4_mul)
-  TEST_ENTRY(glm_mat2x4_mulv)
+  TEST_ENTRY(glmc_mat2x4_mulv)
-  TEST_ENTRY(glm_mat2x4_transpose)
+  TEST_ENTRY(glmc_mat2x4_transpose)
-  TEST_ENTRY(glm_mat2x4_scale)
+  TEST_ENTRY(glmc_mat2x4_scale)
 
   /* camera (incl [LR]H cross [NZ]O) */
   TEST_ENTRY(glm_perspective_lh_zo)
@@ -1551,7 +1570,12 @@ TEST_LIST {
 
   /* ray */
   TEST_ENTRY(glm_ray_triangle)
+  TEST_ENTRY(glm_ray_sphere)
+  TEST_ENTRY(glm_ray_at)
+
   TEST_ENTRY(glmc_ray_triangle)
+  TEST_ENTRY(glmc_ray_sphere)
+  TEST_ENTRY(glmc_ray_at)
 
   /* quat */
   TEST_ENTRY(MACRO_GLM_QUAT_IDENTITY_INIT)
@@ -1688,6 +1712,8 @@ TEST_LIST {
   TEST_ENTRY(glm_vec2_complex_mul)
   TEST_ENTRY(glm_vec2_complex_div)
   TEST_ENTRY(glm_vec2_make)
+  TEST_ENTRY(glm_vec2_reflect)
+  TEST_ENTRY(glm_vec2_refract)
 
   TEST_ENTRY(glmc_vec2)
   TEST_ENTRY(glmc_vec2_copy)
@@ -1734,6 +1760,8 @@ TEST_LIST {
   TEST_ENTRY(glmc_vec2_complex_mul)
   TEST_ENTRY(glmc_vec2_complex_div)
   TEST_ENTRY(glmc_vec2_make)
+  TEST_ENTRY(glmc_vec2_reflect)
+  TEST_ENTRY(glmc_vec2_refract)
 
   /* vec3 */
   TEST_ENTRY(MACRO_GLM_VEC3_ONE_INIT)
@@ -1835,6 +1863,9 @@ TEST_LIST {
   TEST_ENTRY(glm_vec3_hadd)
   TEST_ENTRY(glm_vec3_sqrt)
   TEST_ENTRY(glm_vec3_make)
+  TEST_ENTRY(glm_vec3_faceforward)
+  TEST_ENTRY(glm_vec3_reflect)
+  TEST_ENTRY(glm_vec3_refract)
 
   TEST_ENTRY(glmc_vec3)
   TEST_ENTRY(glmc_vec3_copy)
@@ -1911,6 +1942,9 @@ TEST_LIST {
   TEST_ENTRY(glmc_vec3_hadd)
   TEST_ENTRY(glmc_vec3_sqrt)
   TEST_ENTRY(glmc_vec3_make)
+  TEST_ENTRY(glmc_vec3_faceforward)
+  TEST_ENTRY(glmc_vec3_reflect)
+  TEST_ENTRY(glmc_vec3_refract)
 
   /* vec4 */
   TEST_ENTRY(MACRO_GLM_VEC4_ONE_INIT)
@@ -1999,6 +2033,8 @@ TEST_LIST {
   TEST_ENTRY(glm_vec4_hadd)
   TEST_ENTRY(glm_vec4_sqrt)
   TEST_ENTRY(glm_vec4_make)
+  TEST_ENTRY(glm_vec4_reflect)
+  TEST_ENTRY(glm_vec4_refract)
 
   TEST_ENTRY(glmc_vec4)
   TEST_ENTRY(glmc_vec4_copy3)
@@ -2071,6 +2107,8 @@ TEST_LIST {
   TEST_ENTRY(glmc_vec4_hadd)
   TEST_ENTRY(glmc_vec4_sqrt)
   TEST_ENTRY(glmc_vec4_make)
+  TEST_ENTRY(glmc_vec4_reflect)
+  TEST_ENTRY(glmc_vec4_refract)
 
   /* ivec2 */
   TEST_ENTRY(glm_ivec2)

win/cglm.vcxproj

@@ -238,6 +238,9 @@
     <ClInclude Include="..\include\cglm\struct\handed\euler_to_quat_lh.h" />
     <ClInclude Include="..\include\cglm\struct\handed\euler_to_quat_rh.h" />
     <ClInclude Include="..\include\cglm\struct\io.h" />
+    <ClInclude Include="..\include\cglm\struct\ivec2.h" />
+    <ClInclude Include="..\include\cglm\struct\ivec3.h" />
+    <ClInclude Include="..\include\cglm\struct\ivec4.h" />
     <ClInclude Include="..\include\cglm\struct\mat2.h" />
     <ClInclude Include="..\include\cglm\struct\mat2x3.h" />
     <ClInclude Include="..\include\cglm\struct\mat2x4.h" />
@@ -250,6 +253,7 @@
     <ClInclude Include="..\include\cglm\struct\plane.h" />
     <ClInclude Include="..\include\cglm\struct\project.h" />
     <ClInclude Include="..\include\cglm\struct\quat.h" />
+    <ClInclude Include="..\include\cglm\struct\ray.h" />
     <ClInclude Include="..\include\cglm\struct\sphere.h" />
     <ClInclude Include="..\include\cglm\struct\vec2-ext.h" />
     <ClInclude Include="..\include\cglm\struct\vec2.h" />

win/cglm.vcxproj.filters

@@ -720,5 +720,17 @@
     <ClInclude Include="..\include\cglm\struct\handed\euler_to_quat_rh.h">
       <Filter>include\cglm\struct\handed</Filter>
     </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ivec2.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ivec3.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ivec4.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ray.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
   </ItemGroup>
 </Project>