quat: fix quaternion inverse and tests about it

* multiplication quaternion and its inverse must be identity

.gitignore

@@ -59,4 +59,5 @@ cglm_test_ios/*
 cglm_test_iosTests/*
 docs/build/*
 win/cglm_test_*
-* copy.*
+* copy.*
+*.o

CREDITS

@@ -43,3 +43,10 @@ https://github.com/erich666/GraphicsGems/blob/master/gems/TransBox.c
 6. Cull frustum
 http://www.txutxi.com/?p=584
 http://old.cescg.org/CESCG-2002/DSykoraJJelinek/
+
+7. Quaternions
+Initial mat4_quat is borrowed from Apple's simd library
+
+
+8. Vector Rotation using Quaternion
+https://gamedev.stackexchange.com/questions/28395/rotating-vector3-by-a-quaternion

README.md

@@ -1,9 +1,11 @@
 # 🎥 OpenGL Mathematics (glm) for `C`
 [![Build Status](https://travis-ci.org/recp/cglm.svg?branch=master)](https://travis-ci.org/recp/cglm)
-[![Build status](https://ci.appveyor.com/api/projects/status/av7l3gc0yhfex8y4/branch/master?svg=true)](https://ci.appveyor.com/project/recp/cglm/branch/master)
+ [![Build status](https://ci.appveyor.com/api/projects/status/av7l3gc0yhfex8y4/branch/master?svg=true)](https://ci.appveyor.com/project/recp/cglm/branch/master)
 [![Documentation Status](https://readthedocs.org/projects/cglm/badge/?version=latest)](http://cglm.readthedocs.io/en/latest/?badge=latest)
 [![Coverage Status](https://coveralls.io/repos/github/recp/cglm/badge.svg?branch=master)](https://coveralls.io/github/recp/cglm?branch=master)
 [![Codacy Badge](https://api.codacy.com/project/badge/Grade/6a62b37d5f214f178ebef269dc4a6bf1)](https://www.codacy.com/app/recp/cglm?utm_source=github.com&utm_medium=referral&utm_content=recp/cglm&utm_campaign=Badge_Grade)
+[![Backers on Open Collective](https://opencollective.com/cglm/backers/badge.svg)](#backers) 
+[![Sponsors on Open Collective](https://opencollective.com/cglm/sponsors/badge.svg)](#sponsors)
 
 The original glm library is for C++ only (templates, namespaces, classes...), this library targeted to C99 but currently you can use it for C89 safely by language extensions e.g `__restrict`
 
@@ -17,6 +19,7 @@ Complete documentation: http://cglm.readthedocs.io
 - _dup (duplicate) is changed to _copy. For instance `glm_vec_dup -> glm_vec_copy`
 - OpenGL related functions are dropped to make this lib platform/third-party independent
 - make sure you have latest version and feel free to report bugs, troubles
+- **[bugfix]** euler angles was implemented in reverse order (extrinsic) it was fixed, now they are intrinsic. Make sure that you have the latest version
 
 #### Note for C++ developers:
 If you don't aware about original GLM library yet, you may also want to look at:
@@ -73,6 +76,7 @@ Currently *cglm* uses default clip space configuration (-1, 1) for camera functi
 - inline or pre-compiled function call
 - frustum (extract view frustum planes, corners...)
 - bounding box  (AABB in Frustum (culling), crop, merge...)
+- project, unproject
 
 <hr />
 
@@ -114,6 +118,36 @@ glm_mul(T, R, modelMat);
 glm_inv_tr(modelMat);
 ```
 
+## Contributors
+
+This project exists thanks to all the people who contribute. [[Contribute](CONTRIBUTING.md)].
+<a href="graphs/contributors"><img src="https://opencollective.com/cglm/contributors.svg?width=890&button=false" /></a>
+
+
+## Backers
+
+Thank you to all our backers! 🙏 [[Become a backer](https://opencollective.com/cglm#backer)]
+
+<a href="https://opencollective.com/cglm#backers" target="_blank"><img src="https://opencollective.com/cglm/backers.svg?width=890"></a>
+
+
+## Sponsors
+
+Support this project by becoming a sponsor. Your logo will show up here with a link to your website. [[Become a sponsor](https://opencollective.com/cglm#sponsor)]
+
+<a href="https://opencollective.com/cglm/sponsor/0/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/0/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/1/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/1/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/2/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/2/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/3/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/3/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/4/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/4/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/5/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/5/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/6/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/6/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/7/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/7/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/8/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/8/avatar.svg"></a>
+<a href="https://opencollective.com/cglm/sponsor/9/website" target="_blank"><img src="https://opencollective.com/cglm/sponsor/9/avatar.svg"></a>
+
+
+
 ## License
 MIT. check the LICENSE file
 
@@ -161,7 +195,7 @@ If you want to use inline versions of funcstions then; include main header
 ```C
 #include <cglm/cglm.h>
 ```
-the haeder will include all headers. Then call func you want e.g. rotate vector by axis:
+the header will include all headers. Then call func you want e.g. rotate vector by axis:
 ```C
 glm_vec_rotate(v1, glm_rad(45), (vec3){1.0f, 0.0f, 0.0f});
 ```
@@ -180,7 +214,7 @@ to call pre-compiled versions include header with `c` postfix, c means call. Pre
 ```C
 #include <cglm/call.h>
 ```
-this header will include all heaers with c postfix. You need to call functions with c posfix:
+this header will include all headers with c postfix. You need to call functions with c posfix:
 ```C
 glmc_vec_normalize(vec);
 ```
@@ -194,6 +228,27 @@ glm_mat4_mul(m1, m1, m1);
 ```
 the first two parameter are **[in]** and the last one is **[out]** parameter. After multiplied *m1* and *m2* the result is stored in *m1*. This is why we send *m1* twice. You may store result in different matrix, this just an example.
 
+### Example: Computing MVP matrix
+
+#### Option 1
+```C
+mat4 proj, view, model, mvp;
+
+/* init proj, view and model ... */
+
+glm_mat4_mul(proj, view, viewProj);
+glm_mat4_mul(viewProj, model, mvp);
+```
+
+#### Option 2
+```C
+mat4 proj, view, model, mvp;
+
+/* init proj, view and model ... */
+
+glm_mat4_mulN((mat4 *[]){&proj, &view, &model}, 3, mvp);
+```
+
 ## How to send matrix to OpenGL
 
 mat4 is array of vec4 and vec4 is array of floats. `glUniformMatrix4fv` functions accecpts `float*` as `value` (last param), so you can cast mat4 to float* or you can pass first column of matrix as beginning of memory of matrix:

build-deps.sh

@@ -16,14 +16,18 @@ cd $(dirname "$0")
 if [ "$(uname)" = "Darwin" ]; then
   libtoolBin=$(which glibtoolize)
   libtoolBinDir=$(dirname "${libtoolBin}")
-  ln -s $libtoolBin "${libtoolBinDir}/libtoolize"
+
+  if [ ! -f "${libtoolBinDir}/libtoolize" ]; then
+    ln -s $libtoolBin "${libtoolBinDir}/libtoolize"
+  fi
 fi
 
 # general deps: gcc make autoconf automake libtool cmake
 
 # test - cmocka
 cd ./test/lib/cmocka
-mkdir build
+rm -rf build
+mkdir -p build
 cd build
 cmake -DCMAKE_INSTALL_PREFIX=/usr -DCMAKE_BUILD_TYPE=Debug ..
 make -j8

configure.ac

@@ -7,7 +7,7 @@
 #*****************************************************************************
 
 AC_PREREQ([2.69])
-AC_INIT([cglm], [0.3.4], [info@recp.me])
+AC_INIT([cglm], [0.4.0], [info@recp.me])
 AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects])
 
 AC_CONFIG_MACRO_DIR([m4])

docs/source/api.rst

@@ -41,6 +41,7 @@ Follow the :doc:`build` documentation for this
    vec4-ext
    color
    plane
+   project
    util
    io
    call

docs/source/box.rst

@@ -23,6 +23,11 @@ Functions:
 #. :c:func:`glm_aabb_merge`
 #. :c:func:`glm_aabb_crop`
 #. :c:func:`glm_aabb_crop_until`
+#. :c:func:`glm_aabb_frustum`
+#. :c:func:`glm_aabb_invalidate`
+#. :c:func:`glm_aabb_isvalid`
+#. :c:func:`glm_aabb_size`
+#. :c:func:`glm_aabb_radius`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -61,7 +66,7 @@ Functions documentation
       | *[in]*  **cropBox**  crop box
       | *[out]* **dest**     cropped bounding box
 
-.. c:function:: void  glm_frustum_box(vec4 corners[8], mat4 m, vec3 box[2])
+.. c:function:: void  glm_aabb_crop_until(vec3 box[2], vec3 cropBox[2], vec3 clampBox[2], vec3 dest[2])
 
     | crops a bounding box with another one.
 
@@ -90,3 +95,39 @@ Functions documentation
     Parameters:
       | *[in]*   **box**     bounding box
       | *[out]*  **planes**  frustum planes
+
+.. c:function:: void  glm_aabb_invalidate(vec3 box[2])
+
+    | invalidate AABB min and max values
+
+    | It fills *max* values with -FLT_MAX and *min* values with +FLT_MAX
+
+    Parameters:
+      | *[in, out]*   **box**     bounding box
+
+.. c:function:: bool  glm_aabb_isvalid(vec3 box[2])
+
+    | check if AABB is valid or not
+
+    Parameters:
+      | *[in]*   **box**     bounding box
+
+    Returns:
+      returns true if aabb is valid otherwise false
+
+.. c:function:: float  glm_aabb_size(vec3 box[2])
+
+    | distance between of min and max
+
+    Parameters:
+      | *[in]*   **box**     bounding box
+
+    Returns:
+      distance between min - max
+
+.. c:function:: float  glm_aabb_radius(vec3 box[2])
+
+    | radius of sphere which surrounds AABB
+
+    Parameters:
+      | *[in]*   **box**     bounding box

docs/source/euler.rst

@@ -70,6 +70,7 @@ Functions:
 1. :c:func:`glm_euler_order`
 #. :c:func:`glm_euler_angles`
 #. :c:func:`glm_euler`
+#. :c:func:`glm_euler_xyz`
 #. :c:func:`glm_euler_zyx`
 #. :c:func:`glm_euler_zxy`
 #. :c:func:`glm_euler_xzy`
@@ -115,8 +116,18 @@ Functions documentation
 
     | build rotation matrix from euler angles
 
+    this is alias of glm_euler_xyz function
+
     Parameters:
-      | *[in]*  **angles**  angles as vector [Ex, Ey, Ez]
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
+      | *[in]*  **dest**    rotation matrix
+
+.. c:function:: void  glm_euler_xyz(vec3 angles, mat4 dest)
+
+    | build rotation matrix from euler angles
+
+    Parameters:
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **dest**    rotation matrix
 
 .. c:function:: void  glm_euler_zyx(vec3 angles, mat4 dest)
@@ -124,7 +135,7 @@ Functions documentation
     | build rotation matrix from euler angles
 
     Parameters:
-      | *[in]*  **angles**  angles as vector [Ez, Ey, Ex]
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **dest**    rotation matrix
 
 .. c:function:: void  glm_euler_zxy(vec3 angles, mat4 dest)
@@ -132,7 +143,7 @@ Functions documentation
     | build rotation matrix from euler angles
 
     Parameters:
-      | *[in]*  **angles**  angles as vector [Ez, Ex, Ey]
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **dest**    rotation matrix
 
 .. c:function:: void  glm_euler_xzy(vec3 angles, mat4 dest)
@@ -140,7 +151,7 @@ Functions documentation
     | build rotation matrix from euler angles
 
     Parameters:
-      | *[in]*  **angles**  angles as vector [Ex, Ez, Ey]
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **dest**    rotation matrix
 
 .. c:function:: void  glm_euler_yzx(vec3 angles, mat4 dest)
@@ -148,7 +159,7 @@ Functions documentation
     build rotation matrix from euler angles
 
     Parameters:
-      | *[in]*  **angles**  angles as vector [Ey, Ez, Ex]
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **dest**    rotation matrix
 
 .. c:function:: void  glm_euler_yxz(vec3 angles, mat4 dest)
@@ -156,7 +167,7 @@ Functions documentation
     | build rotation matrix from euler angles
 
     Parameters:
-      | *[in]*  **angles**  angles as vector [Ey, Ex, Ez]
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **dest**    rotation matrix
 
 .. c:function:: void glm_euler_by_order(vec3 angles, glm_euler_sq ord, mat4 dest)
@@ -166,6 +177,6 @@ Functions documentation
     Use :c:func:`glm_euler_order` function to build *ord* parameter
 
     Parameters:
-      | *[in]*  **angles**  angles as vector (ord parameter spceifies angles order)
+      | *[in]*  **angles**  angles as vector [Xangle, Yangle, Zangle]
       | *[in]*  **ord**     euler order
       | *[in]*  **dest**    rotation matrix

docs/source/io.rst

@@ -39,6 +39,7 @@ Functions:
 #. :c:func:`glm_vec3_print`
 #. :c:func:`glm_ivec3_print`
 #. :c:func:`glm_versor_print`
+#. :c:func:`glm_aabb_print`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -90,3 +91,12 @@ Functions documentation
     Parameters:
       | *[in]*  **vec**      quaternion
       | *[in]*  **ostream**  FILE to write
+
+.. c:function:: void  glm_aabb_print(versor vec, const char * __restrict tag, FILE * __restrict ostream)
+
+    | print aabb to given stream
+
+    Parameters:
+      | *[in]*  **vec**      aabb (axis-aligned bounding box)
+      | *[in]*  **tag**      tag to find it more easly in logs
+      | *[in]*  **ostream**  FILE to write

docs/source/project.rst

@@ -0,0 +1,102 @@
+.. default-domain:: C
+
+Project / UnProject
+================================================================================
+
+Header: cglm/project.h
+
+Viewport is required as *vec4* **[X, Y, Width, Height]** but this doesn't mean
+that you should store it as **vec4**. You can convert your data representation
+to vec4 before passing it to related functions.
+
+Table of contents (click to go):
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Functions:
+
+1. :c:func:`glm_unprojecti`
+#. :c:func:`glm_unproject`
+#. :c:func:`glm_project`
+
+Functions documentation
+~~~~~~~~~~~~~~~~~~~~~~~
+
+.. c:function:: void  glm_unprojecti(vec3 pos, mat4 invMat, vec4 vp, vec3 dest)
+
+    | maps the specified viewport coordinates into specified space [1]
+      the matrix should contain projection matrix.
+
+    if you don't have ( and don't want to have ) an inverse matrix then use
+    glm_unproject version. You may use existing inverse of matrix in somewhere
+    else, this is why glm_unprojecti exists to save save inversion cost
+
+    [1] space:
+      - if m = invProj:     View Space
+      - if m = invViewProj: World Space
+      - if m = invMVP:      Object Space
+
+    You probably want to map the coordinates into object space
+    so use invMVP as m
+
+    Computing viewProj:
+
+    .. code-block:: c
+
+       glm_mat4_mul(proj, view, viewProj);
+       glm_mat4_mul(viewProj, model, MVP);
+       glm_mat4_inv(viewProj, invMVP);
+
+    Parameters:
+      | *[in]*  **pos**     point/position in viewport coordinates
+      | *[in]*  **invMat**  matrix (see brief)
+      | *[in]*  **vp**      viewport as [x, y, width, height]
+      | *[out]* **dest**    unprojected coordinates
+
+.. c:function:: void  glm_unproject(vec3 pos, mat4 m, vec4 vp, vec3 dest)
+
+    | maps the specified viewport coordinates into specified space [1]
+      the matrix should contain projection matrix.
+
+    this is same as glm_unprojecti except this function get inverse matrix for
+    you.
+
+    [1] space:
+      - if m = proj:     View Space
+      - if m = viewProj: World Space
+      - if m = MVP:      Object Space
+
+    You probably want to map the coordinates into object space so use MVP as m
+
+    Computing viewProj and MVP:
+
+    .. code-block:: c
+
+       glm_mat4_mul(proj, view, viewProj);
+       glm_mat4_mul(viewProj, model, MVP);
+
+    Parameters:
+      | *[in]*  **pos**  point/position in viewport coordinates
+      | *[in]*  **m**    matrix (see brief)
+      | *[in]*  **vp**   viewport as [x, y, width, height]
+      | *[out]* **dest** unprojected coordinates
+
+.. c:function:: void  glm_project(vec3 pos, mat4 m, vec4 vp, vec3 dest)
+
+    | map object coordinates to window coordinates
+
+    Computing MVP:
+
+    .. code-block:: c
+
+       glm_mat4_mul(proj, view, viewProj);
+       glm_mat4_mul(viewProj, model, MVP);
+
+    this could be useful for gettng a bbox which fits with view frustum and
+    object bounding boxes. In this case you crop view frustum box with objects
+    box
+
+    Parameters:
+      | *[in]*  **pos**      object coordinates
+      | *[in]*  **m**        MVP matrix
+      | *[in]*  **vp**       viewport as [x, y, width, height]
+      | *[out]* **dest**     projected coordinates

docs/source/quat.rst

@@ -5,17 +5,16 @@ quaternions
 
 Header: cglm/quat.h
 
- **Important:** *cglm* stores quaternion as [w, x, y, z] in memory, don't
- forget that when changing quaternion items manually. For instance *quat[3]*
- is *quat.z* and *quat[0*] is *quat.w*. This may change in the future if *cglm*
- will got enough request to do that. Probably it will not be changed in near
- future
+ **Important:** *cglm* stores quaternion as **[x, y, z, w]** in memory
+ since **v0.4.0** it was **[w, x, y, z]**
+ before v0.4.0 ( **v0.3.5 and earlier** ). w is real part.
 
-There are some TODOs for quaternions check TODO list to see them.
+What you can do with quaternions with existing functions is (Some of them):
 
-Also **versor** is identity quaternion so the type may change to **vec4** or
-something else. This will not affect existing functions for your engine because
-*versor* is alias of *vec4*
+- You can rotate transform matrix using quaterion
+- You can rotate vector using quaterion
+- You can create view matrix using quaterion
+- You can create a lookrotation (from source point to dest)
 
 Table of contents (click to go):
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -28,14 +27,35 @@ Macros:
 Functions:
 
 1. :c:func:`glm_quat_identity`
+#. :c:func:`glm_quat_init`
 #. :c:func:`glm_quat`
 #. :c:func:`glm_quatv`
+#. :c:func:`glm_quat_copy`
 #. :c:func:`glm_quat_norm`
 #. :c:func:`glm_quat_normalize`
+#. :c:func:`glm_quat_normalize_to`
 #. :c:func:`glm_quat_dot`
-#. :c:func:`glm_quat_mulv`
+#. :c:func:`glm_quat_conjugate`
+#. :c:func:`glm_quat_inv`
+#. :c:func:`glm_quat_add`
+#. :c:func:`glm_quat_sub`
+#. :c:func:`glm_quat_real`
+#. :c:func:`glm_quat_imag`
+#. :c:func:`glm_quat_imagn`
+#. :c:func:`glm_quat_imaglen`
+#. :c:func:`glm_quat_angle`
+#. :c:func:`glm_quat_axis`
+#. :c:func:`glm_quat_mul`
 #. :c:func:`glm_quat_mat4`
+#. :c:func:`glm_quat_mat4t`
+#. :c:func:`glm_quat_mat3`
+#. :c:func:`glm_quat_mat3t`
+#. :c:func:`glm_quat_lerp`
 #. :c:func:`glm_quat_slerp`
+#. :c:func:`glm_quat_look`
+#. :c:func:`glm_quat_for`
+#. :c:func:`glm_quat_forp`
+#. :c:func:`glm_quat_rotatev`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -47,10 +67,23 @@ Functions documentation
     Parameters:
       | *[in, out]*  **q**    quaternion
 
+.. c:function:: void  glm_quat_init(versor q, float x, float y, float z, float w)
+
+    | inits quaternion with given values
+
+    Parameters:
+      | *[out]* **q**      quaternion
+      | *[in]*  **x**      imag.x
+      | *[in]*  **y**      imag.y
+      | *[in]*  **z**      imag.z
+      | *[in]*  **w**      w (real part)
+
 .. c:function:: void  glm_quat(versor q, float  angle, float  x, float  y, float  z)
 
     | creates NEW quaternion with individual axis components
 
+    | given axis will be normalized
+
     Parameters:
       | *[out]* **q**      quaternion
       | *[in]*  **angle**  angle (radians)
@@ -58,14 +91,24 @@ Functions documentation
       | *[in]*  **y**      axis.y
       | *[in]*  **z**      axis.z
 
-.. c:function:: void  glm_quatv(versor q, float  angle, vec3   v)
+.. c:function:: void  glm_quatv(versor q, float  angle, vec3   axis)
 
     | creates NEW quaternion with axis vector
 
+    | given axis will be normalized
+
     Parameters:
       | *[out]* **q**      quaternion
       | *[in]*  **angle**  angle (radians)
-      | *[in]*  **v**      axis
+      | *[in]*  **axis**   axis (will be normalized)
+
+.. c:function:: void  glm_quat_copy(versor q, versor dest)
+
+    | copy quaternion to another one
+
+    Parameters:
+      | *[in]*  **q**     source quaternion
+      | *[out]* **dest**  destination quaternion
 
 .. c:function:: float  glm_quat_norm(versor q)
 
@@ -77,6 +120,14 @@ Functions documentation
     Returns:
       norm (magnitude)
 
+.. c:function:: void  glm_quat_normalize_to(versor q, versor dest)
+
+    | normalize quaternion and store result in dest, original one will not be normalized
+
+    Parameters:
+      | *[in]*  **q**    quaternion to normalize into
+      | *[out]* **dest** destination quaternion
+
 .. c:function:: void  glm_quat_normalize(versor q)
 
     | normalize quaternion
@@ -84,24 +135,118 @@ Functions documentation
     Parameters:
       | *[in, out]*  **q** quaternion
 
-.. c:function:: float  glm_quat_dot(versor q, versor r)
+.. c:function:: float  glm_quat_dot(versor p, versor q)
 
     dot product of two quaternion
 
     Parameters:
-      | *[in]*  **q1**   quaternion 1
-      | *[in]*  **q2**   quaternion 2
+      | *[in]*  **p**   quaternion 1
+      | *[in]*  **q**   quaternion 2
 
     Returns:
       dot product
 
-.. c:function:: void  glm_quat_mulv(versor q1, versor q2, versor dest)
+.. c:function:: void  glm_quat_conjugate(versor q, versor dest)
+
+    conjugate of quaternion
+
+    Parameters:
+      | *[in]*  **q**      quaternion
+      | *[in]*  **dest**   conjugate
+
+.. c:function:: void  glm_quat_inv(versor q, versor dest)
+
+    inverse of non-zero quaternion
+
+    Parameters:
+      | *[in]*  **q**      quaternion
+      | *[in]*  **dest**   inverse quaternion
+
+.. c:function:: void  glm_quat_add(versor p, versor q, versor dest)
+
+    add (componentwise) two quaternions and store result in dest
+
+    Parameters:
+      | *[in]*  **p**      quaternion 1
+      | *[in]*  **q**      quaternion 2
+      | *[in]*  **dest**   result quaternion
+
+.. c:function:: void  glm_quat_sub(versor p, versor q, versor dest)
+
+    subtract (componentwise) two quaternions and store result in dest
+
+    Parameters:
+      | *[in]*  **p**      quaternion 1
+      | *[in]*  **q**      quaternion 2
+      | *[in]*  **dest**   result quaternion
+
+.. c:function:: float  glm_quat_real(versor q)
+
+    returns real part of quaternion
+
+    Parameters:
+      | *[in]*  **q**   quaternion
+
+    Returns:
+      real part (quat.w)
+
+.. c:function:: void  glm_quat_imag(versor q, vec3 dest)
+
+    returns imaginary part of quaternion
+
+    Parameters:
+      | *[in]*   **q**      quaternion
+      | *[out]*  **dest**   imag
+
+.. c:function:: void  glm_quat_imagn(versor q, vec3 dest)
+
+    returns normalized imaginary part of quaternion
+
+    Parameters:
+      | *[in]*   **q**      quaternion
+      | *[out]*  **dest**   imag
+
+.. c:function:: float  glm_quat_imaglen(versor q)
+
+    returns length of imaginary part of quaternion
+
+    Parameters:
+      | *[in]*   **q**      quaternion
+
+    Returns:
+      norm of imaginary part
+
+.. c:function:: float  glm_quat_angle(versor q)
+
+    returns angle of quaternion
+
+    Parameters:
+      | *[in]*  **q**   quaternion
+
+    Returns:
+      angles of quat (radians)
+
+.. c:function:: void  glm_quat_axis(versor q, versor dest)
+
+    axis of quaternion
+
+    Parameters:
+      | *[in]*   **p**      quaternion
+      | *[out]*  **dest**   axis of quaternion
+
+.. c:function:: void  glm_quat_mul(versor p, versor q, versor dest)
 
     | multiplies two quaternion and stores result in dest
 
+    | this is also called Hamilton Product
+
+    | According to WikiPedia:
+    | The product of two rotation quaternions [clarification needed] will be
+      equivalent to the rotation q followed by the rotation p
+
     Parameters:
-      | *[in]*  **q1**    quaternion 1
-      | *[in]*  **q2**    quaternion 2
+      | *[in]*  **p**     quaternion 1 (first rotation)
+      | *[in]*  **q**     quaternion 2 (second rotation)
       | *[out]* **dest**  result quaternion
 
 .. c:function:: void  glm_quat_mat4(versor q, mat4 dest)
@@ -112,13 +257,100 @@ Functions documentation
       | *[in]*  **q**     quaternion
       | *[out]* **dest**  result matrix
 
+.. c:function:: void  glm_quat_mat4t(versor q, mat4 dest)
+
+    | convert quaternion to mat4 (transposed). This is transposed version of glm_quat_mat4
+
+    Parameters:
+      | *[in]*  **q**     quaternion
+      | *[out]* **dest**  result matrix
+
+.. c:function:: void  glm_quat_mat3(versor q, mat3 dest)
+
+    | convert quaternion to mat3
+
+    Parameters:
+      | *[in]*  **q**     quaternion
+      | *[out]* **dest**  result matrix
+
+.. c:function:: void  glm_quat_mat3t(versor q, mat3 dest)
+
+    | convert quaternion to mat3 (transposed). This is transposed version of glm_quat_mat3
+
+    Parameters:
+      | *[in]*  **q**     quaternion
+      | *[out]* **dest**  result matrix
+
+.. c:function:: void  glm_quat_lerp(versor from, versor to, float t, versor dest)
+
+    | interpolates between two quaternions
+    | using spherical linear interpolation (LERP)
+
+    Parameters:
+      | *[in]*  **from**  from
+      | *[in]*  **to**    to
+      | *[in]*  **t**     interpolant (amount) clamped between 0 and 1
+      | *[out]* **dest**  result quaternion
+
 .. c:function:: void glm_quat_slerp(versor q, versor r, float  t, versor dest)
 
     | interpolates between two quaternions
     | using spherical linear interpolation (SLERP)
 
     Parameters:
-      | *[in]*  **q**     from
-      | *[in]*  **r**     to
-      | *[in]*  **t**     amout
+      | *[in]*  **from**  from
+      | *[in]*  **to**    to
+      | *[in]*  **t**     interpolant (amount) clamped between 0 and 1
       | *[out]* **dest**  result quaternion
+
+.. c:function:: void  glm_quat_look(vec3 eye, versor ori, mat4 dest)
+
+    | creates view matrix using quaternion as camera orientation
+
+    Parameters:
+      | *[in]*  **eye**   eye
+      | *[in]*  **ori**   orientation in world space as quaternion
+      | *[out]* **dest**  result matrix
+
+.. c:function:: void  glm_quat_for(vec3 dir, vec3 fwd, vec3 up, versor dest)
+
+    | creates look rotation quaternion
+
+    Parameters:
+      | *[in]*  **dir**   direction to look
+      | *[in]*  **fwd**   forward vector
+      | *[in]*  **up**    up vector
+      | *[out]* **dest**  result matrix
+
+.. c:function:: void  glm_quat_forp(vec3 from, vec3 to, vec3 fwd, vec3 up, versor dest)
+
+    | creates look rotation quaternion using source and  destination positions p suffix stands for position
+
+    | this is similar to glm_quat_for except this computes direction for glm_quat_for for you.
+
+    Parameters:
+      | *[in]*  **from**  source point
+      | *[in]*  **to**    destination point
+      | *[in]*  **fwd**   forward vector
+      | *[in]*  **up**    up vector
+      | *[out]* **dest**  result matrix
+
+.. c:function:: void  glm_quat_rotatev(versor q, vec3 v, vec3 dest)
+
+    | crotate vector using using quaternion
+
+    Parameters:
+      | *[in]*  **q**     quaternion
+      | *[in]*  **v**     vector to rotate
+      | *[out]* **dest**  rotated vector
+
+.. c:function:: void glm_quat_rotate(mat4 m, versor q, mat4 dest)
+
+    | rotate existing transform matrix using quaternion
+
+    instead of passing identity matrix, consider to use quat_mat4 functions
+
+    Parameters:
+      | *[in]*  **m**     existing transform matrix to rotate
+      | *[in]*  **q**     quaternion
+      | *[out]* **dest**  rotated matrix/transform

docs/source/util.rst

@@ -13,18 +13,25 @@ Table of contents (click to go):
 Functions:
 
 1. :c:func:`glm_sign`
+#. :c:func:`glm_signf`
 #. :c:func:`glm_rad`
 #. :c:func:`glm_deg`
 #. :c:func:`glm_make_rad`
 #. :c:func:`glm_make_deg`
 #. :c:func:`glm_pow2`
+#. :c:func:`glm_min`
+#. :c:func:`glm_max`
+#. :c:func:`glm_clamp`
+#. :c:func:`glm_lerp`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
 
 .. c:function:: int  glm_sign(int val)
 
-    | returns sign of 32 bit integer as +1 or -1
+    | returns sign of 32 bit integer as +1, -1, 0
+
+    | **Important**: It returns 0 for zero input
 
     Parameters:
       | *[in]*  **val**   an integer
@@ -32,6 +39,18 @@ Functions documentation
     Returns:
       sign of given number
 
+.. c:function:: float  glm_signf(float val)
+
+    | returns sign of 32 bit integer as +1.0, -1.0, 0.0
+
+    | **Important**: It returns 0.0f for zero input
+
+    Parameters:
+      | *[in]*  **val**   a float
+
+    Returns:
+      sign of given number
+
 .. c:function:: float  glm_rad(float deg)
 
     | convert degree to radians
@@ -91,3 +110,29 @@ Functions documentation
 
     Returns:
       maximum value
+
+.. c:function:: void  glm_clamp(float val, float minVal, float maxVal)
+
+    constrain a value to lie between two further values
+
+    Parameters:
+      | *[in]*  **val**     input value
+      | *[in]*  **minVal**  minimum value
+      | *[in]*  **maxVal**  maximum value
+
+    Returns:
+      clamped value
+
+.. c:function:: float  glm_lerp(float from, float to, float t)
+
+    linear interpolation between two number
+
+    | formula:  from + s * (to - from)
+
+    Parameters:
+      | *[in]*  **from**   from value
+      | *[in]*  **to**     to value
+      | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
+
+    Returns:
+       interpolated value

docs/source/vec3-ext.rst

@@ -23,6 +23,11 @@ Functions:
 #. :c:func:`glm_vec_eqv_eps`
 #. :c:func:`glm_vec_max`
 #. :c:func:`glm_vec_min`
+#. :c:func:`glm_vec_isnan`
+#. :c:func:`glm_vec_isinf`
+#. :c:func:`glm_vec_isvalid`
+#. :c:func:`glm_vec_sign`
+#. :c:func:`glm_vec_sqrt`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -96,3 +101,43 @@ Functions documentation
 
     Parameters:
       | *[in]*  **v**  vector
+
+.. c:function:: bool glm_vec_isnan(vec3 v)
+
+    | check if one of items is NaN (not a number)
+    | you should only use this in DEBUG mode or very critical asserts
+
+    Parameters:
+      | *[in]*  **v**  vector
+
+.. c:function:: bool glm_vec_isinf(vec3 v)
+
+    | check if one of items is INFINITY
+    | you should only use this in DEBUG mode or very critical asserts
+
+    Parameters:
+      | *[in]*  **v**  vector
+
+.. c:function:: bool glm_vec_isvalid(vec3 v)
+
+    | check if all items are valid number
+    | you should only use this in DEBUG mode or very critical asserts
+
+    Parameters:
+      | *[in]*  **v**  vector
+
+.. c:function:: void glm_vec_sign(vec3 v, vec3 dest)
+
+    get sign of 32 bit float as +1, -1, 0
+
+    Parameters:
+      | *[in]*   **v**     vector
+      | *[out]*  **dest**  sign vector (only keeps signs as -1, 0, -1)
+
+.. c:function:: void glm_vec_sqrt(vec3 v, vec3 dest)
+
+    square root of each vector item
+
+    Parameters:
+      | *[in]*   **v**     vector
+      | *[out]*  **dest**  destination vector (sqrt(v))

docs/source/vec3.rst

@@ -40,6 +40,7 @@ Functions:
 #. :c:func:`glm_vec_scale`
 #. :c:func:`glm_vec_scale_as`
 #. :c:func:`glm_vec_flipsign`
+#. :c:func:`glm_vec_flipsign_to`
 #. :c:func:`glm_vec_inv`
 #. :c:func:`glm_vec_inv_to`
 #. :c:func:`glm_vec_normalize`
@@ -53,6 +54,8 @@ Functions:
 #. :c:func:`glm_vec_maxv`
 #. :c:func:`glm_vec_minv`
 #. :c:func:`glm_vec_ortho`
+#. :c:func:`glm_vec_clamp`
+#. :c:func:`glm_vec_lerp`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -156,7 +159,15 @@ Functions documentation
     flip sign of all vec3 members
 
     Parameters:
-    | *[in, out]*  **v**    vector
+      | *[in, out]*  **v**    vector
+
+.. c:function:: void  glm_vec_flipsign_to(vec3 v, vec3 dest)
+
+    flip sign of all vec3 members and store result in dest
+
+    Parameters:
+      | *[in]*  **v**       vector
+      | *[out]* **dest**    negated vector
 
 .. c:function:: void  glm_vec_inv(vec3 v)
 
@@ -205,7 +216,7 @@ Functions documentation
 
     Parameters:
       | *[in, out]*  **v**      vector
-      | *[in]*       **axis**   axis vector (must be unit vector)
+      | *[in]*       **axis**   axis vector (will be normalized)
       | *[out]*      **angle**  angle (radians)
 
 .. c:function:: void  glm_vec_rotate_m4(mat4 m, vec3 v, vec3 dest)
@@ -271,3 +282,24 @@ Functions documentation
     Parameters:
       | *[in]*  **mat**   vector
       | *[out]* **dest**  orthogonal/perpendicular vector
+
+.. c:function:: void  glm_vec_clamp(vec3 v, float minVal, float maxVal)
+
+    constrain a value to lie between two further values
+
+    Parameters:
+      | *[in, out]*  **v**       vector
+      | *[in]*       **minVal**  minimum value
+      | *[in]*       **maxVal**  maximum value
+
+.. c:function:: void  glm_vec_lerp(vec3 from, vec3 to, float t, vec3 dest)
+
+    linear interpolation between two vector
+
+    | formula:  from + s * (to - from)
+
+    Parameters:
+      | *[in]*  **from**   from value
+      | *[in]*  **to**     to value
+      | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
+      | *[out]* **dest**   destination

docs/source/vec4-ext.rst

@@ -96,3 +96,43 @@ Functions documentation
 
     Parameters:
       | *[in]*  **v**  vector
+
+.. c:function:: bool glm_vec4_isnan(vec4 v)
+
+    | check if one of items is NaN (not a number)
+    | you should only use this in DEBUG mode or very critical asserts
+
+    Parameters:
+      | *[in]*  **v**  vector
+
+.. c:function:: bool glm_vec4_isinf(vec4 v)
+
+    | check if one of items is INFINITY
+    | you should only use this in DEBUG mode or very critical asserts
+
+    Parameters:
+      | *[in]*  **v**  vector
+
+.. c:function:: bool glm_vec4_isvalid(vec4 v)
+
+    | check if all items are valid number
+    | you should only use this in DEBUG mode or very critical asserts
+
+    Parameters:
+      | *[in]*  **v**  vector
+
+.. c:function:: void glm_vec4_sign(vec4 v, vec4 dest)
+
+    get sign of 32 bit float as +1, -1, 0
+
+    Parameters:
+      | *[in]*   **v**     vector
+      | *[out]*  **dest**  sign vector (only keeps signs as -1, 0, -1)
+
+.. c:function:: void glm_vec4_sqrt(vec4 v, vec4 dest)
+
+    square root of each vector item
+
+    Parameters:
+      | *[in]*   **v**     vector
+      | *[out]*  **dest**  destination vector (sqrt(v))

docs/source/vec4.rst

@@ -32,6 +32,7 @@ Functions:
 #. :c:func:`glm_vec4_scale`
 #. :c:func:`glm_vec4_scale_as`
 #. :c:func:`glm_vec4_flipsign`
+#. :c:func:`glm_vec_flipsign_to`
 #. :c:func:`glm_vec4_inv`
 #. :c:func:`glm_vec4_inv_to`
 #. :c:func:`glm_vec4_normalize`
@@ -39,6 +40,13 @@ Functions:
 #. :c:func:`glm_vec4_distance`
 #. :c:func:`glm_vec4_maxv`
 #. :c:func:`glm_vec4_minv`
+#. :c:func:`glm_vec4_clamp`
+#. :c:func:`glm_vec4_lerp`
+#. :c:func:`glm_vec4_isnan`
+#. :c:func:`glm_vec4_isinf`
+#. :c:func:`glm_vec4_isvalid`
+#. :c:func:`glm_vec4_sign`
+#. :c:func:`glm_vec4_sqrt`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -145,6 +153,14 @@ Functions documentation
     Parameters:
     | *[in, out]*  **v**    vector
 
+.. c:function:: void  glm_vec4_flipsign_to(vec4 v, vec4 dest)
+
+    flip sign of all vec4 members and store result in dest
+
+    Parameters:
+      | *[in]*  **v**       vector
+      | *[out]* **dest**    negated vector
+
 .. c:function:: void  glm_vec4_inv(vec4 v)
 
     make vector as inverse/opposite of itself
@@ -203,3 +219,24 @@ Functions documentation
       | *[in]*  **v1**    vector1
       | *[in]*  **v2**    vector2
       | *[out]* **dest**  destination
+
+.. c:function:: void  glm_vec4_clamp(vec4 v, float minVal, float maxVal)
+
+    constrain a value to lie between two further values
+
+    Parameters:
+      | *[in, out]*  **v**       vector
+      | *[in]*       **minVal**  minimum value
+      | *[in]*       **maxVal**  maximum value
+
+.. c:function:: void  glm_vec4_lerp(vec4 from, vec4 to, float t, vec4 dest)
+
+    linear interpolation between two vector
+
+    | formula:  from + s * (to - from)
+
+    Parameters:
+      | *[in]*  **from**   from value
+      | *[in]*  **to**     to value
+      | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
+      | *[out]* **dest**   destination

include/cglm/box.h

@@ -153,4 +153,50 @@ glm_aabb_frustum(vec3 box[2], vec4 planes[6]) {
   return true;
 }
 
+/*!
+ * @brief invalidate AABB min and max values
+ *
+ * @param[in, out]  box bounding box
+ */
+CGLM_INLINE
+void
+glm_aabb_invalidate(vec3 box[2]) {
+  glm_vec_broadcast(FLT_MAX,  box[0]);
+  glm_vec_broadcast(-FLT_MAX, box[1]);
+}
+
+/*!
+ * @brief check if AABB is valid or not
+ *
+ * @param[in]  box bounding box
+ */
+CGLM_INLINE
+bool
+glm_aabb_isvalid(vec3 box[2]) {
+  return glm_vec_max(box[0]) != FLT_MAX
+         && glm_vec_min(box[1]) != -FLT_MAX;
+}
+
+/*!
+ * @brief distance between of min and max
+ *
+ * @param[in]  box bounding box
+ */
+CGLM_INLINE
+float
+glm_aabb_size(vec3 box[2]) {
+  return glm_vec_distance(box[0], box[1]);
+}
+
+/*!
+ * @brief radius of sphere which surrounds AABB
+ *
+ * @param[in]  box bounding box
+ */
+CGLM_INLINE
+float
+glm_aabb_radius(vec3 box[2]) {
+  return glm_aabb_size(box) * 0.5f;
+}
+
 #endif /* cglm_box_h */

include/cglm/call.h

@@ -24,6 +24,7 @@ extern "C" {
 #include "call/frustum.h"
 #include "call/box.h"
 #include "call/io.h"
+#include "call/project.h"
 
 #ifdef __cplusplus
 }

include/cglm/call/euler.h

@@ -21,6 +21,10 @@ CGLM_EXPORT
 void
 glmc_euler(vec3 angles, mat4 dest);
 
+CGLM_EXPORT
+void
+glmc_euler_xyz(vec3 angles,  mat4 dest);
+
 CGLM_EXPORT
 void
 glmc_euler_zyx(vec3 angles,  mat4 dest);

include/cglm/call/mat4.h

@@ -47,12 +47,16 @@ glmc_mat4_mul(mat4 m1, mat4 m2, mat4 dest);
 
 CGLM_EXPORT
 void
-glmc_mat4_mulN(mat4 * __restrict matrices[], int len, mat4 dest);
+glmc_mat4_mulN(mat4 * __restrict matrices[], uint32_t len, mat4 dest);
 
 CGLM_EXPORT
 void
 glmc_mat4_mulv(mat4 m, vec4 v, vec4 dest);
 
+CGLM_EXPORT
+void
+glmc_mat4_quat(mat4 m, versor dest);
+
 CGLM_EXPORT
 void
 glmc_mat4_transpose_to(mat4 m, mat4 dest);

include/cglm/call/project.h

@@ -0,0 +1,33 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglmc_project_h
+#define cglmc_project_h
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "../cglm.h"
+
+CGLM_EXPORT
+void
+glmc_unprojecti(vec3 pos, mat4 invMat, vec4 vp, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_unproject(vec3 pos, mat4 m, vec4 vp, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_project(vec3 pos, mat4 m, vec4 vp, vec3 dest);
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* cglmc_project_h */
+
+

include/cglm/call/quat.h

@@ -19,33 +19,79 @@ glmc_quat_identity(versor q);
 
 CGLM_EXPORT
 void
-glmc_quat(versor q,
-          float angle,
-          float x,
-          float y,
-          float z);
+glmc_quat_init(versor q, float x, float y, float z, float w);
 
 CGLM_EXPORT
 void
-glmc_quatv(versor q,
-          float  angle,
-          vec3   v);
+glmc_quat(versor q, float angle, float x, float y, float z);
+
+CGLM_EXPORT
+void
+glmc_quatv(versor q, float angle, vec3 axis);
+
+CGLM_EXPORT
+void
+glmc_quat_copy(versor q, versor dest);
 
 CGLM_EXPORT
 float
 glmc_quat_norm(versor q);
 
+CGLM_EXPORT
+void
+glmc_quat_normalize_to(versor q, versor dest);
+
 CGLM_EXPORT
 void
 glmc_quat_normalize(versor q);
 
 CGLM_EXPORT
 float
-glmc_quat_dot(versor q, versor r);
+glmc_quat_dot(versor p, versor q);
 
 CGLM_EXPORT
 void
-glmc_quat_mulv(versor q1, versor q2, versor dest);
+glmc_quat_conjugate(versor q, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_inv(versor q, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_add(versor p, versor q, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_sub(versor p, versor q, versor dest);
+
+CGLM_EXPORT
+float
+glmc_quat_real(versor q);
+
+CGLM_EXPORT
+void
+glmc_quat_imag(versor q, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_quat_imagn(versor q, vec3 dest);
+
+CGLM_EXPORT
+float
+glmc_quat_imaglen(versor q);
+
+CGLM_EXPORT
+float
+glmc_quat_angle(versor q);
+
+CGLM_EXPORT
+void
+glmc_quat_axis(versor q, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_mul(versor p, versor q, versor dest);
 
 CGLM_EXPORT
 void
@@ -53,10 +99,43 @@ glmc_quat_mat4(versor q, mat4 dest);
 
 CGLM_EXPORT
 void
-glmc_quat_slerp(versor q,
-                versor r,
-                float  t,
-                versor dest);
+glmc_quat_mat4t(versor q, mat4 dest);
+
+CGLM_EXPORT
+void
+glmc_quat_mat3(versor q, mat3 dest);
+
+CGLM_EXPORT
+void
+glmc_quat_mat3t(versor q, mat3 dest);
+
+CGLM_EXPORT
+void
+glmc_quat_lerp(versor from, versor to, float t, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_slerp(versor q, versor r, float t, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_look(vec3 eye, versor ori, mat4 dest);
+
+CGLM_EXPORT
+void
+glmc_quat_for(vec3 dir, vec3 fwd, vec3 up, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_forp(vec3 from, vec3 to, vec3 fwd, vec3 up, versor dest);
+
+CGLM_EXPORT
+void
+glmc_quat_rotatev(versor from, vec3 to, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_quat_rotate(mat4 m, versor q, mat4 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/vec3.h

@@ -16,6 +16,10 @@ extern "C" {
 /* DEPRECATED! use _copy, _ucopy versions */
 #define glmc_vec_dup(v, dest) glmc_vec_copy(v, dest)
 
+CGLM_EXPORT
+void
+glmc_vec3(vec4 v4, vec3 dest);
+
 CGLM_EXPORT
 void
 glmc_vec_copy(vec3 a, vec3 dest);
@@ -64,6 +68,10 @@ CGLM_EXPORT
 void
 glmc_vec_flipsign(vec3 v);
 
+CGLM_EXPORT
+void
+glmc_vec_flipsign_to(vec3 v, vec3 dest);
+
 CGLM_EXPORT
 void
 glmc_vec_inv(vec3 v);
@@ -104,6 +112,76 @@ CGLM_EXPORT
 void
 glmc_vec_minv(vec3 v1, vec3 v2, vec3 dest);
 
+CGLM_EXPORT
+void
+glmc_vec_clamp(vec3 v, float minVal, float maxVal);
+
+CGLM_EXPORT
+void
+glmc_vec_ortho(vec3 v, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec_lerp(vec3 from, vec3 to, float t, vec3 dest);
+
+/* ext */
+
+CGLM_EXPORT
+void
+glmc_vec_mulv(vec3 a, vec3 b, vec3 d);
+
+CGLM_EXPORT
+void
+glmc_vec_broadcast(float val, vec3 d);
+
+CGLM_EXPORT
+bool
+glmc_vec_eq(vec3 v, float val);
+
+CGLM_EXPORT
+bool
+glmc_vec_eq_eps(vec3 v, float val);
+
+CGLM_EXPORT
+bool
+glmc_vec_eq_all(vec3 v);
+
+CGLM_EXPORT
+bool
+glmc_vec_eqv(vec3 v1, vec3 v2);
+
+CGLM_EXPORT
+bool
+glmc_vec_eqv_eps(vec3 v1, vec3 v2);
+
+CGLM_EXPORT
+float
+glmc_vec_max(vec3 v);
+
+CGLM_EXPORT
+float
+glmc_vec_min(vec3 v);
+
+CGLM_EXPORT
+bool
+glmc_vec_isnan(vec3 v);
+
+CGLM_EXPORT
+bool
+glmc_vec_isinf(vec3 v);
+
+CGLM_EXPORT
+bool
+glmc_vec_isvalid(vec3 v);
+
+CGLM_EXPORT
+void
+glmc_vec_sign(vec3 v, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec_sqrt(vec3 v, vec3 dest);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/call/vec4.h

@@ -17,6 +17,10 @@ extern "C" {
 #define glmc_vec4_dup3(v, dest) glmc_vec4_copy3(v, dest)
 #define glmc_vec4_dup(v, dest)  glmc_vec4_copy(v, dest)
 
+CGLM_EXPORT
+void
+glmc_vec4(vec3 v3, float last, vec4 dest);
+
 CGLM_EXPORT
 void
 glmc_vec4_copy3(vec4 a, vec3 dest);
@@ -65,6 +69,10 @@ CGLM_EXPORT
 void
 glmc_vec4_flipsign(vec4 v);
 
+CGLM_EXPORT
+void
+glmc_vec4_flipsign_to(vec4 v, vec4 dest);
+
 CGLM_EXPORT
 void
 glmc_vec4_inv(vec4 v);
@@ -85,6 +93,72 @@ CGLM_EXPORT
 void
 glmc_vec4_minv(vec4 v1, vec4 v2, vec4 dest);
 
+CGLM_EXPORT
+void
+glmc_vec4_clamp(vec4 v, float minVal, float maxVal);
+
+CGLM_EXPORT
+void
+glmc_vec4_lerp(vec4 from, vec4 to, float t, vec4 dest);
+
+/* ext */
+
+CGLM_EXPORT
+void
+glmc_vec4_mulv(vec4 a, vec4 b, vec4 d);
+
+CGLM_EXPORT
+void
+glmc_vec4_broadcast(float val, vec4 d);
+
+CGLM_EXPORT
+bool
+glmc_vec4_eq(vec4 v, float val);
+
+CGLM_EXPORT
+bool
+glmc_vec4_eq_eps(vec4 v, float val);
+
+CGLM_EXPORT
+bool
+glmc_vec4_eq_all(vec4 v);
+
+CGLM_EXPORT
+bool
+glmc_vec4_eqv(vec4 v1, vec4 v2);
+
+CGLM_EXPORT
+bool
+glmc_vec4_eqv_eps(vec4 v1, vec4 v2);
+
+CGLM_EXPORT
+float
+glmc_vec4_max(vec4 v);
+
+CGLM_EXPORT
+float
+glmc_vec4_min(vec4 v);
+
+CGLM_EXPORT
+bool
+glmc_vec4_isnan(vec4 v);
+
+CGLM_EXPORT
+bool
+glmc_vec4_isinf(vec4 v);
+
+CGLM_EXPORT
+bool
+glmc_vec4_isvalid(vec4 v);
+
+CGLM_EXPORT
+void
+glmc_vec4_sign(vec4 v, vec4 dest);
+
+CGLM_EXPORT
+void
+glmc_vec4_sqrt(vec4 v, vec4 dest);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/cglm.h

@@ -23,5 +23,6 @@
 #include "color.h"
 #include "util.h"
 #include "io.h"
+#include "project.h"
 
 #endif /* cglm_h */

include/cglm/euler.h

@@ -5,6 +5,14 @@
  * Full license can be found in the LICENSE file
  */
 
+/*
+ NOTE:
+  angles must be passed as [X-Angle, Y-Angle, Z-angle] order
+  For instance you don't pass angles as [Z-Angle, X-Angle, Y-angle] to
+  glm_euler_zxy funciton, All RELATED functions accept angles same order
+  which is [X, Y, Z].
+ */
+
 /*
  Types:
    enum glm_euler_sq
@@ -13,6 +21,7 @@
    CGLM_INLINE glm_euler_sq glm_euler_order(int newOrder[3]);
    CGLM_INLINE void glm_euler_angles(mat4 m, vec3 dest);
    CGLM_INLINE void glm_euler(vec3 angles, mat4 dest);
+   CGLM_INLINE void glm_euler_xyz(vec3 angles, mat4 dest);
    CGLM_INLINE void glm_euler_zyx(vec3 angles, mat4 dest);
    CGLM_INLINE void glm_euler_zxy(vec3 angles, mat4 dest);
    CGLM_INLINE void glm_euler_xzy(vec3 angles, mat4 dest);
@@ -61,253 +70,283 @@ glm_euler_order(int ord[3]) {
 CGLM_INLINE
 void
 glm_euler_angles(mat4 m, vec3 dest) {
-  if (m[0][2] < 1.0f) {
-    if (m[0][2] > -1.0f) {
-      vec3  a[2];
-      float cy1, cy2;
-      int   path;
+  float m00, m01, m10, m11, m20, m21, m22;
+  float thetaX, thetaY, thetaZ;
 
-      a[0][1] = asinf(-m[0][2]);
-      a[1][1] = CGLM_PI - a[0][1];
+  m00 = m[0][0];  m10 = m[1][0];  m20 = m[2][0];
+  m01 = m[0][1];  m11 = m[1][1];  m21 = m[2][1];
+                                  m22 = m[2][2];
 
-      cy1 = cosf(a[0][1]);
-      cy2 = cosf(a[1][1]);
-
-      a[0][0] = atan2f(m[1][2] / cy1, m[2][2] / cy1);
-      a[1][0] = atan2f(m[1][2] / cy2, m[2][2] / cy2);
-
-      a[0][2] = atan2f(m[0][1] / cy1, m[0][0] / cy1);
-      a[1][2] = atan2f(m[0][1] / cy2, m[0][0] / cy2);
-
-      path = (fabsf(a[0][0]) + fabsf(a[0][1]) + fabsf(a[0][2])) >=
-               (fabsf(a[1][0]) + fabsf(a[1][1]) + fabsf(a[1][2]));
-
-      glm_vec_copy(a[path], dest);
-    } else {
-      dest[0] = atan2f(m[1][0], m[2][0]);
-      dest[1] = CGLM_PI_2;
-      dest[2] = 0.0f;
+  if (m20 < 1.0f) {
+    if (m20 > -1.0f) {
+      thetaY = asinf(m20);
+      thetaX = atan2f(-m21, m22);
+      thetaZ = atan2f(-m10, m00);
+    } else { /* m20 == -1 */
+      /* Not a unique solution */
+      thetaY = -CGLM_PI_2;
+      thetaX = -atan2f(m01, m11);
+      thetaZ =  0.0f;
     }
-  } else {
-    dest[0] = atan2f(-m[1][0], -m[2][0]);
-    dest[1] =-CGLM_PI_2;
-    dest[2] = 0.0f;
+  } else { /* m20 == +1 */
+    thetaY = CGLM_PI_2;
+    thetaX = atan2f(m01, m11);
+    thetaZ = 0.0f;
   }
+
+  dest[0] = thetaX;
+  dest[1] = thetaY;
+  dest[2] = thetaZ;
 }
 
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector [Ex, Ey, Ez]
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
+ * @param[out] dest   rotation matrix
+ */
+CGLM_INLINE
+void
+glm_euler_xyz(vec3 angles, mat4 dest) {
+  float cx, cy, cz,
+        sx, sy, sz, czsx, cxcz, sysz;
+
+  sx   = sinf(angles[0]); cx = cosf(angles[0]);
+  sy   = sinf(angles[1]); cy = cosf(angles[1]);
+  sz   = sinf(angles[2]); cz = cosf(angles[2]);
+
+  czsx = cz * sx;
+  cxcz = cx * cz;
+  sysz = sy * sz;
+
+  dest[0][0] =  cy * cz;
+  dest[0][1] =  czsx * sy + cx * sz;
+  dest[0][2] = -cxcz * sy + sx * sz;
+  dest[1][0] = -cy * sz;
+  dest[1][1] =  cxcz - sx * sysz;
+  dest[1][2] =  czsx + cx * sysz;
+  dest[2][0] =  sy;
+  dest[2][1] = -cy * sx;
+  dest[2][2] =  cx * cy;
+  dest[0][3] =  0.0f;
+  dest[1][3] =  0.0f;
+  dest[2][3] =  0.0f;
+  dest[3][0] =  0.0f;
+  dest[3][1] =  0.0f;
+  dest[3][2] =  0.0f;
+  dest[3][3] =  1.0f;
+}
+
+/*!
+ * @brief build rotation matrix from euler angles
+ *
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[out] dest   rotation matrix
  */
 CGLM_INLINE
 void
 glm_euler(vec3 angles, mat4 dest) {
-  float cx, cy, cz,
-        sx, sy, sz;
-
-  sx = sinf(angles[0]); cx = cosf(angles[0]);
-  sy = sinf(angles[1]); cy = cosf(angles[1]);
-  sz = sinf(angles[2]); cz = cosf(angles[2]);
-
-  dest[0][0] = cy * cz;
-  dest[0][1] = cy * sz;
-  dest[0][2] =-sy;
-  dest[1][0] = cz * sx * sy - cx * sz;
-  dest[1][1] = cx * cz + sx * sy * sz;
-  dest[1][2] = cy * sx;
-  dest[2][0] = cx * cz * sy + sx * sz;
-  dest[2][1] =-cz * sx + cx * sy * sz;
-  dest[2][2] = cx * cy;
-  dest[0][3] = 0.0f;
-  dest[1][3] = 0.0f;
-  dest[2][3] = 0.0f;
-  dest[3][0] = 0.0f;
-  dest[3][1] = 0.0f;
-  dest[3][2] = 0.0f;
-  dest[3][3] = 1.0f;
+  glm_euler_xyz(angles, dest);
 }
 
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector [Ez, Ey, Ex]
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[out] dest   rotation matrix
  */
 CGLM_INLINE
 void
-glm_euler_zyx(vec3 angles,
-              mat4 dest) {
+glm_euler_xzy(vec3 angles, mat4 dest) {
   float cx, cy, cz,
-        sx, sy, sz;
+  sx, sy, sz, sxsy, cysx, cxsy, cxcy;
 
-  sx = sinf(angles[0]); cx = cosf(angles[0]);
-  sy = sinf(angles[1]); cy = cosf(angles[1]);
-  sz = sinf(angles[2]); cz = cosf(angles[2]);
+  sx   = sinf(angles[0]); cx = cosf(angles[0]);
+  sy   = sinf(angles[1]); cy = cosf(angles[1]);
+  sz   = sinf(angles[2]); cz = cosf(angles[2]);
 
-  dest[0][0] = cy * cz;
-  dest[0][1] = cz * sx * sy + cx * sz;
-  dest[0][2] =-cx * cz * sy + sx * sz;
-  dest[1][0] =-cy * sz;
-  dest[1][1] = cx * cz - sx * sy * sz;
-  dest[1][2] = cz * sx + cx * sy * sz;
-  dest[2][0] = sy;
-  dest[2][1] =-cy * sx;
-  dest[2][2] = cx * cy;
-  dest[0][3] = 0.0f;
-  dest[1][3] = 0.0f;
-  dest[2][3] = 0.0f;
-  dest[3][0] = 0.0f;
-  dest[3][1] = 0.0f;
-  dest[3][2] = 0.0f;
-  dest[3][3] = 1.0f;
+  sxsy = sx * sy;
+  cysx = cy * sx;
+  cxsy = cx * sy;
+  cxcy = cx * cy;
+
+  dest[0][0] =  cy * cz;
+  dest[0][1] =  sxsy + cxcy * sz;
+  dest[0][2] = -cxsy + cysx * sz;
+  dest[1][0] = -sz;
+  dest[1][1] =  cx * cz;
+  dest[1][2] =  cz * sx;
+  dest[2][0] =  cz * sy;
+  dest[2][1] = -cysx + cxsy * sz;
+  dest[2][2] =  cxcy + sxsy * sz;
+  dest[0][3] =  0.0f;
+  dest[1][3] =  0.0f;
+  dest[2][3] =  0.0f;
+  dest[3][0] =  0.0f;
+  dest[3][1] =  0.0f;
+  dest[3][2] =  0.0f;
+  dest[3][3] =  1.0f;
 }
 
+
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector [Ez, Ex, Ey]
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[out] dest   rotation matrix
  */
 CGLM_INLINE
 void
-glm_euler_zxy(vec3 angles,
-              mat4 dest) {
+glm_euler_yxz(vec3 angles, mat4 dest) {
   float cx, cy, cz,
-        sx, sy, sz;
+        sx, sy, sz, cycz, sysz, czsy, cysz;
 
-  sx = sinf(angles[0]); cx = cosf(angles[0]);
-  sy = sinf(angles[1]); cy = cosf(angles[1]);
-  sz = sinf(angles[2]); cz = cosf(angles[2]);
+  sx   = sinf(angles[0]); cx = cosf(angles[0]);
+  sy   = sinf(angles[1]); cy = cosf(angles[1]);
+  sz   = sinf(angles[2]); cz = cosf(angles[2]);
 
-  dest[0][0] = cy * cz + sx * sy * sz;
-  dest[0][1] = cx * sz;
-  dest[0][2] =-cz * sy + cy * sx * sz;
-  dest[1][0] = cz * sx * sy - cy * sz;
-  dest[1][1] = cx * cz;
-  dest[1][2] = cy * cz * sx + sy * sz;
-  dest[2][0] = cx * sy;
-  dest[2][1] =-sx;
-  dest[2][2] = cx * cy;
-  dest[0][3] = 0.0f;
-  dest[1][3] = 0.0f;
-  dest[2][3] = 0.0f;
-  dest[3][0] = 0.0f;
-  dest[3][1] = 0.0f;
-  dest[3][2] = 0.0f;
-  dest[3][3] = 1.0f;
+  cycz = cy * cz;
+  sysz = sy * sz;
+  czsy = cz * sy;
+  cysz = cy * sz;
+
+  dest[0][0] =  cycz + sx * sysz;
+  dest[0][1] =  cx * sz;
+  dest[0][2] = -czsy + cysz * sx;
+  dest[1][0] = -cysz + czsy * sx;
+  dest[1][1] =  cx * cz;
+  dest[1][2] =  cycz * sx + sysz;
+  dest[2][0] =  cx * sy;
+  dest[2][1] = -sx;
+  dest[2][2] =  cx * cy;
+  dest[0][3] =  0.0f;
+  dest[1][3] =  0.0f;
+  dest[2][3] =  0.0f;
+  dest[3][0] =  0.0f;
+  dest[3][1] =  0.0f;
+  dest[3][2] =  0.0f;
+  dest[3][3] =  1.0f;
 }
 
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector [Ex, Ez, Ey]
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[out] dest   rotation matrix
  */
 CGLM_INLINE
 void
-glm_euler_xzy(vec3 angles,
-              mat4 dest) {
+glm_euler_yzx(vec3 angles, mat4 dest) {
   float cx, cy, cz,
-        sx, sy, sz;
+        sx, sy, sz, sxsy, cxcy, cysx, cxsy;
 
-  sx = sinf(angles[0]); cx = cosf(angles[0]);
-  sy = sinf(angles[1]); cy = cosf(angles[1]);
-  sz = sinf(angles[2]); cz = cosf(angles[2]);
+  sx   = sinf(angles[0]); cx = cosf(angles[0]);
+  sy   = sinf(angles[1]); cy = cosf(angles[1]);
+  sz   = sinf(angles[2]); cz = cosf(angles[2]);
 
-  dest[0][0] = cy * cz;
-  dest[0][1] = sz;
-  dest[0][2] =-cz * sy;
-  dest[1][0] = sx * sy - cx * cy * sz;
-  dest[1][1] = cx * cz;
-  dest[1][2] = cy * sx + cx * sy * sz;
-  dest[2][0] = cx * sy + cy * sx * sz;
-  dest[2][1] =-cz * sx;
-  dest[2][2] = cx * cy - sx * sy * sz;
-  dest[0][3] = 0.0f;
-  dest[1][3] = 0.0f;
-  dest[2][3] = 0.0f;
-  dest[3][0] = 0.0f;
-  dest[3][1] = 0.0f;
-  dest[3][2] = 0.0f;
-  dest[3][3] = 1.0f;
+  sxsy = sx * sy;
+  cxcy = cx * cy;
+  cysx = cy * sx;
+  cxsy = cx * sy;
+
+  dest[0][0] =  cy * cz;
+  dest[0][1] =  sz;
+  dest[0][2] = -cz * sy;
+  dest[1][0] =  sxsy - cxcy * sz;
+  dest[1][1] =  cx * cz;
+  dest[1][2] =  cysx + cxsy * sz;
+  dest[2][0] =  cxsy + cysx * sz;
+  dest[2][1] = -cz * sx;
+  dest[2][2] =  cxcy - sxsy * sz;
+  dest[0][3] =  0.0f;
+  dest[1][3] =  0.0f;
+  dest[2][3] =  0.0f;
+  dest[3][0] =  0.0f;
+  dest[3][1] =  0.0f;
+  dest[3][2] =  0.0f;
+  dest[3][3] =  1.0f;
 }
 
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector [Ey, Ez, Ex]
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[out] dest   rotation matrix
  */
 CGLM_INLINE
 void
-glm_euler_yzx(vec3 angles,
-              mat4 dest) {
+glm_euler_zxy(vec3 angles, mat4 dest) {
   float cx, cy, cz,
-        sx, sy, sz;
+        sx, sy, sz, cycz, sxsy, cysz;
 
-  sx = sinf(angles[0]); cx = cosf(angles[0]);
-  sy = sinf(angles[1]); cy = cosf(angles[1]);
-  sz = sinf(angles[2]); cz = cosf(angles[2]);
+  sx   = sinf(angles[0]); cx = cosf(angles[0]);
+  sy   = sinf(angles[1]); cy = cosf(angles[1]);
+  sz   = sinf(angles[2]); cz = cosf(angles[2]);
 
-  dest[0][0] = cy * cz;
-  dest[0][1] = sx * sy + cx * cy * sz;
-  dest[0][2] =-cx * sy + cy * sx * sz;
-  dest[1][0] =-sz;
-  dest[1][1] = cx * cz;
-  dest[1][2] = cz * sx;
-  dest[2][0] = cz * sy;
-  dest[2][1] =-cy * sx + cx * sy * sz;
-  dest[2][2] = cx * cy + sx * sy * sz;
-  dest[0][3] = 0.0f;
-  dest[1][3] = 0.0f;
-  dest[2][3] = 0.0f;
-  dest[3][0] = 0.0f;
-  dest[3][1] = 0.0f;
-  dest[3][2] = 0.0f;
-  dest[3][3] = 1.0f;
+  cycz = cy * cz;
+  sxsy = sx * sy;
+  cysz = cy * sz;
+
+  dest[0][0] =  cycz - sxsy * sz;
+  dest[0][1] =  cz * sxsy + cysz;
+  dest[0][2] = -cx * sy;
+  dest[1][0] = -cx * sz;
+  dest[1][1] =  cx * cz;
+  dest[1][2] =  sx;
+  dest[2][0] =  cz * sy + cysz * sx;
+  dest[2][1] = -cycz * sx + sy * sz;
+  dest[2][2] =  cx * cy;
+  dest[0][3] =  0.0f;
+  dest[1][3] =  0.0f;
+  dest[2][3] =  0.0f;
+  dest[3][0] =  0.0f;
+  dest[3][1] =  0.0f;
+  dest[3][2] =  0.0f;
+  dest[3][3] =  1.0f;
 }
 
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector [Ey, Ex, Ez]
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[out] dest   rotation matrix
  */
 CGLM_INLINE
 void
-glm_euler_yxz(vec3 angles,
-              mat4 dest) {
+glm_euler_zyx(vec3 angles, mat4 dest) {
   float cx, cy, cz,
-        sx, sy, sz;
+        sx, sy, sz, czsx, cxcz, sysz;
 
-  sx = sinf(angles[0]); cx = cosf(angles[0]);
-  sy = sinf(angles[1]); cy = cosf(angles[1]);
-  sz = sinf(angles[2]); cz = cosf(angles[2]);
+  sx   = sinf(angles[0]); cx = cosf(angles[0]);
+  sy   = sinf(angles[1]); cy = cosf(angles[1]);
+  sz   = sinf(angles[2]); cz = cosf(angles[2]);
 
-  dest[0][0] = cy * cz - sx * sy * sz;
-  dest[0][1] = cz * sx * sy + cy * sz;
-  dest[0][2] =-cx * sy;
-  dest[1][0] =-cx * sz;
-  dest[1][1] = cx * cz;
-  dest[1][2] = sx;
-  dest[2][0] = cz * sy + cy * sx * sz;
-  dest[2][1] =-cy * cz * sx + sy * sz;
-  dest[2][2] = cx * cy;
-  dest[0][3] = 0.0f;
-  dest[1][3] = 0.0f;
-  dest[2][3] = 0.0f;
-  dest[3][0] = 0.0f;
-  dest[3][1] = 0.0f;
-  dest[3][2] = 0.0f;
-  dest[3][3] = 1.0f;
+  czsx = cz * sx;
+  cxcz = cx * cz;
+  sysz = sy * sz;
+
+  dest[0][0] =  cy * cz;
+  dest[0][1] =  cy * sz;
+  dest[0][2] = -sy;
+  dest[1][0] =  czsx * sy - cx * sz;
+  dest[1][1] =  cxcz + sx * sysz;
+  dest[1][2] =  cy * sx;
+  dest[2][0] =  cxcz * sy + sx * sz;
+  dest[2][1] = -czsx + cx * sysz;
+  dest[2][2] =  cx * cy;
+  dest[0][3] =  0.0f;
+  dest[1][3] =  0.0f;
+  dest[2][3] =  0.0f;
+  dest[3][0] =  0.0f;
+  dest[3][1] =  0.0f;
+  dest[3][2] =  0.0f;
+  dest[3][3] =  1.0f;
 }
 
 /*!
  * @brief build rotation matrix from euler angles
  *
- * @param[in]  angles angles as vector (ord parameter spceifies angles order)
+ * @param[in]  angles angles as vector [Xangle, Yangle, Zangle]
  * @param[in]  ord    euler order
  * @param[out] dest   rotation matrix
  */
@@ -332,71 +371,71 @@ glm_euler_by_order(vec3 angles, glm_euler_sq ord, mat4 dest) {
   sysz = sy * sz;
 
   switch (ord) {
-    case GLM_EULER_XYZ:
-      dest[0][0] = cycz;
-      dest[0][1] = cysz;
-      dest[0][2] =-sy;
-      dest[1][0] = czsx * sy - cxsz;
-      dest[1][1] = cxcz + sx * sysz;
-      dest[1][2] = cysx;
-      dest[2][0] = cx * czsy + sx * sz;
-      dest[2][1] =-czsx + cx * sysz;
-      dest[2][2] = cxcy;
-      break;
     case GLM_EULER_XZY:
-      dest[0][0] = cycz;
-      dest[0][1] = sz;
-      dest[0][2] =-czsy;
-      dest[1][0] = sx * sy - cx * cysz;
-      dest[1][1] = cxcz;
-      dest[1][2] = cysx + cx * sysz;
-      dest[2][0] = cx * sy + cysx * sz;
-      dest[2][1] =-czsx;
-      dest[2][2] = cxcy - sx * sysz;
+      dest[0][0] =  cycz;
+      dest[0][1] =  sx * sy + cx * cysz;
+      dest[0][2] = -cx * sy + cysx * sz;
+      dest[1][0] = -sz;
+      dest[1][1] =  cxcz;
+      dest[1][2] =  czsx;
+      dest[2][0] =  czsy;
+      dest[2][1] = -cysx + cx * sysz;
+      dest[2][2] =  cxcy + sx * sysz;
       break;
-    case GLM_EULER_ZXY:
-      dest[0][0] = cycz + sx * sysz;
-      dest[0][1] = cxsz;
-      dest[0][2] =-czsy + cysx * sz;
-      dest[1][0] = czsx * sy - cysz;
-      dest[1][1] = cxcz;
-      dest[1][2] = cycz * sx + sysz;
-      dest[2][0] = cx * sy;
-      dest[2][1] =-sx;
-      dest[2][2] = cxcy;
-      break;
-    case GLM_EULER_ZYX:
-      dest[0][0] = cycz;
-      dest[0][1] = czsx * sy + cxsz;
-      dest[0][2] =-cx * czsy + sx * sz;
-      dest[1][0] =-cysz;
-      dest[1][1] = cxcz - sx * sysz;
-      dest[1][2] = czsx + cx * sysz;
-      dest[2][0] = sy;
-      dest[2][1] =-cysx;
-      dest[2][2] = cxcy;
+    case GLM_EULER_XYZ:
+      dest[0][0] =  cycz;
+      dest[0][1] =  czsx * sy + cxsz;
+      dest[0][2] = -cx * czsy + sx * sz;
+      dest[1][0] = -cysz;
+      dest[1][1] =  cxcz - sx * sysz;
+      dest[1][2] =  czsx + cx * sysz;
+      dest[2][0] =  sy;
+      dest[2][1] = -cysx;
+      dest[2][2] =  cxcy;
       break;
     case GLM_EULER_YXZ:
-      dest[0][0] = cycz - sx * sysz;
-      dest[0][1] = czsx * sy + cysz;
-      dest[0][2] =-cx * sy;
-      dest[1][0] =-cxsz;
-      dest[1][1] = cxcz;
-      dest[1][2] = sx;
-      dest[2][0] = czsy + cysx * sz;
-      dest[2][1] =-cycz * sx + sysz;
-      dest[2][2] = cxcy;
+      dest[0][0] =  cycz + sx * sysz;
+      dest[0][1] =  cxsz;
+      dest[0][2] = -czsy + cysx * sz;
+      dest[1][0] =  czsx * sy - cysz;
+      dest[1][1] =  cxcz;
+      dest[1][2] =  cycz * sx + sysz;
+      dest[2][0] =  cx * sy;
+      dest[2][1] = -sx;
+      dest[2][2] =  cxcy;
       break;
     case GLM_EULER_YZX:
-      dest[0][0] = cycz;
-      dest[0][1] = sx * sy + cx * cysz;
-      dest[0][2] =-cx * sy + cysx * sz;
-      dest[1][0] =-sz;
-      dest[1][1] = cxcz;
-      dest[1][2] = czsx;
-      dest[2][0] = czsy;
-      dest[2][1] =-cysx + cx * sysz;
-      dest[2][2] = cxcy + sx * sysz;
+      dest[0][0] =  cycz;
+      dest[0][1] =  sz;
+      dest[0][2] = -czsy;
+      dest[1][0] =  sx * sy - cx * cysz;
+      dest[1][1] =  cxcz;
+      dest[1][2] =  cysx + cx * sysz;
+      dest[2][0] =  cx * sy + cysx * sz;
+      dest[2][1] = -czsx;
+      dest[2][2] =  cxcy - sx * sysz;
+      break;
+    case GLM_EULER_ZXY:
+      dest[0][0] =  cycz - sx * sysz;
+      dest[0][1] =  czsx * sy + cysz;
+      dest[0][2] = -cx * sy;
+      dest[1][0] = -cxsz;
+      dest[1][1] =  cxcz;
+      dest[1][2] =  sx;
+      dest[2][0] =  czsy + cysx * sz;
+      dest[2][1] = -cycz * sx + sysz;
+      dest[2][2] =  cxcy;
+      break;
+    case GLM_EULER_ZYX:
+      dest[0][0] =  cycz;
+      dest[0][1] =  cysz;
+      dest[0][2] = -sy;
+      dest[1][0] =  czsx * sy - cxsz;
+      dest[1][1] =  cxcz + sx * sysz;
+      dest[1][2] =  cysx;
+      dest[2][0] =  cx * czsy + sx * sz;
+      dest[2][1] = -czsx + cx * sysz;
+      dest[2][2] =  cxcy;
       break;
   }
 

include/cglm/io.h

@@ -171,4 +171,33 @@ glm_versor_print(versor vec,
 #undef m
 }
 
+CGLM_INLINE
+void
+glm_aabb_print(vec3                    bbox[2],
+               const char * __restrict tag,
+               FILE       * __restrict ostream) {
+  int i, j;
+
+#define m 3
+
+  fprintf(ostream, "AABB (%s):\n", tag ? tag: "float");
+
+  for (i = 0; i < 2; i++) {
+    fprintf(ostream, "\t|");
+
+    for (j = 0; j < m; j++) {
+      fprintf(ostream, "%0.4f", bbox[i][j]);
+
+      if (j != m - 1)
+        fprintf(ostream, "\t");
+    }
+
+    fprintf(ostream, "|\n");
+  }
+
+  fprintf(ostream, "\n");
+
+#undef m
+}
+
 #endif /* cglm_io_h */

include/cglm/mat3.h

@@ -186,6 +186,56 @@ glm_mat3_mulv(mat3 m, vec3 v, vec3 dest) {
   dest[2] = m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2];
 }
 
+
+/*!
+ * @brief convert mat4's rotation part to quaternion
+ *
+ * @param[in]  m    left matrix
+ * @param[out] dest destination quaternion
+ */
+CGLM_INLINE
+void
+glm_mat3_quat(mat3 m, versor dest) {
+  float trace, r, rinv;
+
+  /* it seems using like m12 instead of m[1][2] causes extra instructions */
+
+  trace = m[0][0] + m[1][1] + m[2][2];
+  if (trace >= 0.0f) {
+    r       = sqrtf(1.0f + trace);
+    rinv    = 0.5f / r;
+
+    dest[0] = rinv * (m[1][2] - m[2][1]);
+    dest[1] = rinv * (m[2][0] - m[0][2]);
+    dest[2] = rinv * (m[0][1] - m[1][0]);
+    dest[3] = r    * 0.5f;
+  } else if (m[0][0] >= m[1][1] && m[0][0] >= m[2][2]) {
+    r       = sqrtf(1.0f - m[1][1] - m[2][2] + m[0][0]);
+    rinv    = 0.5f / r;
+
+    dest[0] = r    * 0.5f;
+    dest[1] = rinv * (m[0][1] + m[1][0]);
+    dest[2] = rinv * (m[0][2] + m[2][0]);
+    dest[3] = rinv * (m[1][2] - m[2][1]);
+  } else if (m[1][1] >= m[2][2]) {
+    r       = sqrtf(1.0f - m[0][0] - m[2][2] + m[1][1]);
+    rinv    = 0.5f / r;
+
+    dest[0] = rinv * (m[0][1] + m[1][0]);
+    dest[1] = r    * 0.5f;
+    dest[2] = rinv * (m[1][2] + m[2][1]);
+    dest[3] = rinv * (m[2][0] - m[0][2]);
+  } else {
+    r       = sqrtf(1.0f - m[0][0] - m[1][1] + m[2][2]);
+    rinv    = 0.5f / r;
+
+    dest[0] = rinv * (m[0][2] + m[2][0]);
+    dest[1] = rinv * (m[1][2] + m[2][1]);
+    dest[2] = r    * 0.5f;
+    dest[3] = rinv * (m[0][1] - m[1][0]);
+  }
+}
+
 /*!
  * @brief scale (multiply with scalar) matrix
  *

include/cglm/mat4.h

@@ -45,6 +45,7 @@
 #define cglm_mat_h
 
 #include "common.h"
+#include "quat.h"
 
 #ifdef CGLM_SSE_FP
 #  include "simd/sse2/mat4.h"
@@ -58,7 +59,9 @@
 #  include "simd/neon/mat4.h"
 #endif
 
-#include <assert.h>
+#ifdef DEBUG
+# include <assert.h>
+#endif
 
 #define GLM_MAT4_IDENTITY_INIT  {{1.0f, 0.0f, 0.0f, 0.0f},                    \
                                  {0.0f, 1.0f, 0.0f, 0.0f},                    \
@@ -281,19 +284,17 @@ glm_mat4_mul(mat4 m1, mat4 m2, mat4 dest) {
  */
 CGLM_INLINE
 void
-glm_mat4_mulN(mat4 * __restrict matrices[], int len, mat4 dest) {
-  int i;
+glm_mat4_mulN(mat4 * __restrict matrices[], uint32_t len, mat4 dest) {
+  uint32_t i;
 
+#ifdef DEBUG
   assert(len > 1 && "there must be least 2 matrices to go!");
+#endif
 
-  glm_mat4_mul(*matrices[0],
-               *matrices[1],
-               dest);
+  glm_mat4_mul(*matrices[0], *matrices[1], dest);
 
   for (i = 2; i < len; i++)
-    glm_mat4_mul(dest,
-                 *matrices[i],
-                 dest);
+    glm_mat4_mul(dest, *matrices[i], dest);
 }
 
 /*!
@@ -318,6 +319,55 @@ glm_mat4_mulv(mat4 m, vec4 v, vec4 dest) {
 #endif
 }
 
+/*!
+ * @brief convert mat4's rotation part to quaternion
+ *
+ * @param[in]  m    left matrix
+ * @param[out] dest destination quaternion
+ */
+CGLM_INLINE
+void
+glm_mat4_quat(mat4 m, versor dest) {
+  float trace, r, rinv;
+
+  /* it seems using like m12 instead of m[1][2] causes extra instructions */
+
+  trace = m[0][0] + m[1][1] + m[2][2];
+  if (trace >= 0.0f) {
+    r       = sqrtf(1.0f + trace);
+    rinv    = 0.5f / r;
+
+    dest[0] = rinv * (m[1][2] - m[2][1]);
+    dest[1] = rinv * (m[2][0] - m[0][2]);
+    dest[2] = rinv * (m[0][1] - m[1][0]);
+    dest[3] = r    * 0.5f;
+  } else if (m[0][0] >= m[1][1] && m[0][0] >= m[2][2]) {
+    r       = sqrtf(1.0f - m[1][1] - m[2][2] + m[0][0]);
+    rinv    = 0.5f / r;
+
+    dest[0] = r    * 0.5f;
+    dest[1] = rinv * (m[0][1] + m[1][0]);
+    dest[2] = rinv * (m[0][2] + m[2][0]);
+    dest[3] = rinv * (m[1][2] - m[2][1]);
+  } else if (m[1][1] >= m[2][2]) {
+    r       = sqrtf(1.0f - m[0][0] - m[2][2] + m[1][1]);
+    rinv    = 0.5f / r;
+
+    dest[0] = rinv * (m[0][1] + m[1][0]);
+    dest[1] = r    * 0.5f;
+    dest[2] = rinv * (m[1][2] + m[2][1]);
+    dest[3] = rinv * (m[2][0] - m[0][2]);
+  } else {
+    r       = sqrtf(1.0f - m[0][0] - m[1][1] + m[2][2]);
+    rinv    = 0.5f / r;
+
+    dest[0] = rinv * (m[0][2] + m[2][0]);
+    dest[1] = rinv * (m[1][2] + m[2][1]);
+    dest[2] = r    * 0.5f;
+    dest[3] = rinv * (m[0][1] - m[1][0]);
+  }
+}
+
 /*!
  * @brief multiply vector with mat4's mat3 part(rotation)
  *
@@ -568,5 +618,4 @@ glm_mat4_swap_row(mat4 mat, int row1, int row2) {
   mat[3][row2] = tmp[3];
 }
 
-#else
 #endif /* cglm_mat_h */

include/cglm/project.h

@@ -0,0 +1,117 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglm_project_h
+#define cglm_project_h
+
+#include "mat4.h"
+#include "vec3.h"
+#include "vec4.h"
+
+/*!
+ * @brief maps the specified viewport coordinates into specified space [1]
+ *        the matrix should contain projection matrix.
+ *
+ * if you don't have ( and don't want to have ) an inverse matrix then use
+ * glm_unproject version. You may use existing inverse of matrix in somewhere
+ * else, this is why glm_unprojecti exists to save save inversion cost
+ *
+ * [1] space:
+ *  1- if m = invProj:     View Space
+ *  2- if m = invViewProj: World Space
+ *  3- if m = invMVP:      Object Space
+ *
+ * You probably want to map the coordinates into object space
+ * so use invMVP as m
+ *
+ * Computing viewProj:
+ *   glm_mat4_mul(proj, view, viewProj);
+ *   glm_mat4_mul(viewProj, model, MVP);
+ *   glm_mat4_inv(viewProj, invMVP);
+ *
+ * @param[in]  pos      point/position in viewport coordinates
+ * @param[in]  invMat   matrix (see brief)
+ * @param[in]  vp       viewport as [x, y, width, height]
+ * @param[out] dest     unprojected coordinates
+ */
+CGLM_INLINE
+void
+glm_unprojecti(vec3 pos, mat4 invMat, vec4 vp, vec3 dest) {
+  vec4 v;
+
+  v[0] = 2.0f * (pos[0] - vp[0]) / vp[2] - 1.0f;
+  v[1] = 2.0f * (pos[1] - vp[1]) / vp[3] - 1.0f;
+  v[2] = 2.0f *  pos[2]                  - 1.0f;
+  v[3] = 1.0f;
+
+  glm_mat4_mulv(invMat, v, v);
+  glm_vec4_scale(v, 1.0f / v[3], v);
+  glm_vec3(v, dest);
+}
+
+/*!
+ * @brief maps the specified viewport coordinates into specified space [1]
+ *        the matrix should contain projection matrix.
+ *
+ * this is same as glm_unprojecti except this function get inverse matrix for
+ * you.
+ *
+ * [1] space:
+ *  1- if m = proj:     View Space
+ *  2- if m = viewProj: World Space
+ *  3- if m = MVP:      Object Space
+ *
+ * You probably want to map the coordinates into object space
+ * so use MVP as m
+ *
+ * Computing viewProj and MVP:
+ *   glm_mat4_mul(proj, view, viewProj);
+ *   glm_mat4_mul(viewProj, model, MVP);
+ *
+ * @param[in]  pos      point/position in viewport coordinates
+ * @param[in]  m        matrix (see brief)
+ * @param[in]  vp       viewport as [x, y, width, height]
+ * @param[out] dest     unprojected coordinates
+ */
+CGLM_INLINE
+void
+glm_unproject(vec3 pos, mat4 m, vec4 vp, vec3 dest) {
+  mat4 inv;
+  glm_mat4_inv(m, inv);
+  glm_unprojecti(pos, inv, vp, dest);
+}
+
+/*!
+ * @brief map object coordinates to window coordinates
+ *
+ * Computing MVP:
+ *   glm_mat4_mul(proj, view, viewProj);
+ *   glm_mat4_mul(viewProj, model, MVP);
+ *
+ * @param[in]  pos      object coordinates
+ * @param[in]  m        MVP matrix
+ * @param[in]  vp       viewport as [x, y, width, height]
+ * @param[out] dest     projected coordinates
+ */
+CGLM_INLINE
+void
+glm_project(vec3 pos, mat4 m, vec4 vp, vec3 dest) {
+  vec4 pos4, vone = GLM_VEC4_ONE_INIT;
+
+  glm_vec4(pos, 1.0f, pos4);
+
+  glm_mat4_mulv(m, pos4, pos4);
+  glm_vec4_scale(pos4, 1.0f / pos4[3], pos4); /* pos = pos / pos.w */
+  glm_vec4_add(pos4, vone, pos4);
+  glm_vec4_scale(pos4, 0.5f, pos4);
+
+  dest[0] = pos4[0] * vp[2] + vp[0];
+  dest[1] = pos4[1] * vp[3] + vp[1];
+  dest[2] = pos4[2];
+}
+
+#endif /* cglm_project_h */

include/cglm/quat.h

@@ -11,15 +11,41 @@
    GLM_QUAT_IDENTITY
 
  Functions:
-   CGLM_INLINE void  glm_quat_identity(versor q);
-   CGLM_INLINE void  glm_quat(versor q, float angle, float x, float y, float z);
-   CGLM_INLINE void  glm_quatv(versor q, float angle, vec3 v);
+   CGLM_INLINE void glm_quat_identity(versor q);
+   CGLM_INLINE void glm_quat_init(versor q, float x, float y, float z, float w);
+   CGLM_INLINE void glm_quat(versor q, float angle, float x, float y, float z);
+   CGLM_INLINE void glm_quatv(versor q, float angle, vec3 axis);
+   CGLM_INLINE void glm_quat_copy(versor q, versor dest);
    CGLM_INLINE float glm_quat_norm(versor q);
-   CGLM_INLINE void  glm_quat_normalize(versor q);
-   CGLM_INLINE float glm_quat_dot(versor q, versor r);
-   CGLM_INLINE void  glm_quat_mulv(versor q1, versor q2, versor dest);
-   CGLM_INLINE void  glm_quat_mat4(versor q, mat4 dest);
-   CGLM_INLINE void  glm_quat_slerp(versor q, versor r, float t, versor dest);
+   CGLM_INLINE void glm_quat_normalize(versor q);
+   CGLM_INLINE void glm_quat_normalize_to(versor q, versor dest);
+   CGLM_INLINE float glm_quat_dot(versor q1, versor q2);
+   CGLM_INLINE void glm_quat_conjugate(versor q, versor dest);
+   CGLM_INLINE void glm_quat_inv(versor q, versor dest);
+   CGLM_INLINE void glm_quat_add(versor p, versor q, versor dest);
+   CGLM_INLINE void glm_quat_sub(versor p, versor q, versor dest);
+   CGLM_INLINE float glm_quat_real(versor q);
+   CGLM_INLINE void glm_quat_imag(versor q, vec3 dest);
+   CGLM_INLINE void glm_quat_imagn(versor q, vec3 dest);
+   CGLM_INLINE float glm_quat_imaglen(versor q);
+   CGLM_INLINE float glm_quat_angle(versor q);
+   CGLM_INLINE void glm_quat_axis(versor q, versor dest);
+   CGLM_INLINE void glm_quat_mul(versor p, versor q, versor dest);
+   CGLM_INLINE void glm_quat_mat4(versor q, mat4 dest);
+   CGLM_INLINE void glm_quat_mat4t(versor q, mat4 dest);
+   CGLM_INLINE void glm_quat_mat3(versor q, mat3 dest);
+   CGLM_INLINE void glm_quat_mat3t(versor q, mat3 dest);
+   CGLM_INLINE void glm_quat_lerp(versor from, versor to, float t, versor dest);
+   CGLM_INLINE void glm_quat_slerp(versor q, versor r, float t, versor dest);
+   CGLM_INLINE void glm_quat_look(vec3 eye, versor ori, mat4 dest);
+   CGLM_INLINE void glm_quat_for(vec3 dir, vec3 fwd, vec3 up, versor dest);
+   CGLM_INLINE void glm_quat_forp(vec3 from,
+                                  vec3 to,
+                                  vec3 fwd,
+                                  vec3 up,
+                                  versor dest);
+   CGLM_INLINE void glm_quat_rotatev(versor q, vec3 v, vec3 dest);
+   CGLM_INLINE void glm_quat_rotate(mat4 m, versor q, mat4 dest);
  */
 
 #ifndef cglm_quat_h
@@ -27,25 +53,32 @@
 
 #include "common.h"
 #include "vec4.h"
+#include "mat4.h"
+#include "mat3.h"
 
 #ifdef CGLM_SSE_FP
 #  include "simd/sse2/quat.h"
 #endif
 
+CGLM_INLINE
+void
+glm_mat4_mulv(mat4 m, vec4 v, vec4 dest);
+
+CGLM_INLINE
+void
+glm_mat4_mul(mat4 m1, mat4 m2, mat4 dest);
+
 /*
- * IMPORTANT! cglm stores quat as [w, x, y, z]
+ * IMPORTANT:
+ * ----------------------------------------------------------------------------
+ * cglm stores quat as [x, y, z, w] since v0.3.6
  *
- * Possible changes (these may be changed in the future):
- *  - versor is identity quat, we can define new type for quat.
- *    it can't be quat or quaternion becuase someone can use that name for
- *    variable name. maybe just vec4.
- *  - it stores [w, x, y, z] but it may change to [x, y, z, w] if we get enough
- *    feedback to change it.
- *  - in general we use last param as dest, but this header used first param
- *    as dest this may be changed but decided yet
+ * it was [w, x, y, z] before v0.3.6 it has been changed to [x, y, z, w]
+ * with v0.3.6 version.
+ * ----------------------------------------------------------------------------
  */
 
-#define GLM_QUAT_IDENTITY_INIT  {1.0f, 0.0f, 0.0f, 0.0f}
+#define GLM_QUAT_IDENTITY_INIT  {0.0f, 0.0f, 0.0f, 1.0f}
 #define GLM_QUAT_IDENTITY       ((versor)GLM_QUAT_IDENTITY_INIT)
 
 /*!
@@ -60,6 +93,49 @@ glm_quat_identity(versor q) {
   glm_vec4_copy(v, q);
 }
 
+/*!
+ * @brief inits quaterion with raw values
+ *
+ * @param[out]  q     quaternion
+ * @param[in]   x     x
+ * @param[in]   y     y
+ * @param[in]   z     z
+ * @param[in]   w     w (real part)
+ */
+CGLM_INLINE
+void
+glm_quat_init(versor q, float x, float y, float z, float w) {
+  q[0] = x;
+  q[1] = y;
+  q[2] = z;
+  q[3] = w;
+}
+
+/*!
+ * @brief creates NEW quaternion with axis vector
+ *
+ * @param[out]  q     quaternion
+ * @param[in]   angle angle (radians)
+ * @param[in]   axis  axis
+ */
+CGLM_INLINE
+void
+glm_quatv(versor q, float angle, vec3 axis) {
+  vec3  k;
+  float a, c, s;
+
+  a = angle * 0.5f;
+  c = cosf(a);
+  s = sinf(a);
+
+  glm_normalize_to(axis, k);
+
+  q[0] = s * k[0];
+  q[1] = s * k[1];
+  q[2] = s * k[2];
+  q[3] = c;
+}
+
 /*!
  * @brief creates NEW quaternion with individual axis components
  *
@@ -71,45 +147,21 @@ glm_quat_identity(versor q) {
  */
 CGLM_INLINE
 void
-glm_quat(versor q,
-         float  angle,
-         float  x,
-         float  y,
-         float  z) {
-  float a, c, s;
-
-  a = angle * 0.5f;
-  c = cosf(a);
-  s = sinf(a);
-
-  q[0] = c;
-  q[1] = s * x;
-  q[2] = s * y;
-  q[3] = s * z;
+glm_quat(versor q, float angle, float x, float y, float z) {
+  vec3 axis = {x, y, z};
+  glm_quatv(q, angle, axis);
 }
 
 /*!
- * @brief creates NEW quaternion with axis vector
+ * @brief copy quaternion to another one
  *
- * @param[out]  q     quaternion
- * @param[in]   angle angle (radians)
- * @param[in]   v     axis
+ * @param[in]  q     quaternion
+ * @param[out] dest  destination
  */
 CGLM_INLINE
 void
-glm_quatv(versor q,
-          float  angle,
-          vec3   v) {
-  float a, c, s;
-
-  a = angle * 0.5f;
-  c = cosf(a);
-  s = sinf(a);
-
-  q[0] = c;
-  q[1] = s * v[0];
-  q[2] = s * v[1];
-  q[3] = s * v[2];
+glm_quat_copy(versor q, versor dest) {
+  glm_vec4_copy(q, dest);
 }
 
 /*!
@@ -123,6 +175,43 @@ glm_quat_norm(versor q) {
   return glm_vec4_norm(q);
 }
 
+/*!
+ * @brief normalize quaternion and store result in dest
+ *
+ * @param[in]   q     quaternion to normalze
+ * @param[out]  dest  destination quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_normalize_to(versor q, versor dest) {
+#if defined( __SSE2__ ) || defined( __SSE2__ )
+  __m128 xdot, x0;
+  float  dot;
+
+  x0   = _mm_load_ps(q);
+  xdot = glm_simd_dot(x0, x0);
+  dot  = _mm_cvtss_f32(xdot);
+
+  if (dot <= 0.0f) {
+    glm_quat_identity(dest);
+    return;
+  }
+
+  _mm_store_ps(dest, _mm_div_ps(x0, _mm_sqrt_ps(xdot)));
+#else
+  float dot;
+
+  dot = glm_vec4_norm2(q);
+
+  if (dot <= 0.0f) {
+    glm_quat_identity(q);
+    return;
+  }
+
+  glm_vec4_scale(q, 1.0f / sqrtf(dot), dest);
+#endif
+}
+
 /*!
  * @brief normalize quaternion
  *
@@ -131,45 +220,178 @@ glm_quat_norm(versor q) {
 CGLM_INLINE
 void
 glm_quat_normalize(versor q) {
-  float sum;
-
-  sum = q[0] * q[0] + q[1] * q[1]
-          + q[2] * q[2] + q[3] * q[3];
-
-  if (fabs(1.0f - sum) < 0.0001f)
-    return;
-
-  glm_vec4_scale(q, 1.0f / sqrtf(sum), q);
+  glm_quat_normalize_to(q, q);
 }
 
 /*!
  * @brief dot product of two quaternion
  *
- * @param[in]  q  quaternion 1
- * @param[in]  r  quaternion 2
+ * @param[in]  p  quaternion 1
+ * @param[in]  q  quaternion 2
  */
 CGLM_INLINE
 float
-glm_quat_dot(versor q, versor r) {
-  return glm_vec4_dot(q, r);
+glm_quat_dot(versor p, versor q) {
+  return glm_vec4_dot(p, q);
+}
+
+/*!
+ * @brief conjugate of quaternion
+ *
+ * @param[in]   q     quaternion
+ * @param[out]  dest  conjugate
+ */
+CGLM_INLINE
+void
+glm_quat_conjugate(versor q, versor dest) {
+  glm_vec4_flipsign_to(q, dest);
+  dest[3] = -dest[3];
+}
+
+/*!
+ * @brief inverse of non-zero quaternion
+ *
+ * @param[in]   q    quaternion
+ * @param[out]  dest inverse quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_inv(versor q, versor dest) {
+  versor conj;
+  glm_quat_conjugate(q, conj);
+  glm_vec4_scale(conj, 1.0f / glm_vec4_norm2(q), dest);
+}
+
+/*!
+ * @brief add (componentwise) two quaternions and store result in dest
+ *
+ * @param[in]   p    quaternion 1
+ * @param[in]   q    quaternion 2
+ * @param[out]  dest result quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_add(versor p, versor q, versor dest) {
+  glm_vec4_add(p, q, dest);
+}
+
+/*!
+ * @brief subtract (componentwise) two quaternions and store result in dest
+ *
+ * @param[in]   p    quaternion 1
+ * @param[in]   q    quaternion 2
+ * @param[out]  dest result quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_sub(versor p, versor q, versor dest) {
+  glm_vec4_sub(p, q, dest);
+}
+
+/*!
+ * @brief returns real part of quaternion
+ *
+ * @param[in]   q    quaternion
+ */
+CGLM_INLINE
+float
+glm_quat_real(versor q) {
+  return q[3];
+}
+
+/*!
+ * @brief returns imaginary part of quaternion
+ *
+ * @param[in]   q    quaternion
+ * @param[out]  dest imag
+ */
+CGLM_INLINE
+void
+glm_quat_imag(versor q, vec3 dest) {
+  dest[0] = q[0];
+  dest[1] = q[1];
+  dest[2] = q[2];
+}
+
+/*!
+ * @brief returns normalized imaginary part of quaternion
+ *
+ * @param[in]   q    quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_imagn(versor q, vec3 dest) {
+  glm_normalize_to(q, dest);
+}
+
+/*!
+ * @brief returns length of imaginary part of quaternion
+ *
+ * @param[in]   q    quaternion
+ */
+CGLM_INLINE
+float
+glm_quat_imaglen(versor q) {
+  return glm_vec_norm(q);
+}
+
+/*!
+ * @brief returns angle of quaternion
+ *
+ * @param[in]   q    quaternion
+ */
+CGLM_INLINE
+float
+glm_quat_angle(versor q) {
+  /*
+   sin(theta / 2) = length(x*x + y*y + z*z)
+   cos(theta / 2) = w
+   theta          = 2 * atan(sin(theta / 2) / cos(theta / 2))
+   */
+  return 2.0f * atan2f(glm_quat_imaglen(q), glm_quat_real(q));
+}
+
+/*!
+ * @brief axis of quaternion
+ *
+ * @param[in]   q    quaternion
+ * @param[out]  dest axis of quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_axis(versor q, versor dest) {
+  glm_quat_imagn(q, dest);
 }
 
 /*!
  * @brief multiplies two quaternion and stores result in dest
+ *        this is also called Hamilton Product
  *
- * @param[in]   q1    quaternion 1
- * @param[in]   q2    quaternion 2
+ * According to WikiPedia:
+ * The product of two rotation quaternions [clarification needed] will be
+ * equivalent to the rotation q followed by the rotation p
+ *
+ * @param[in]   p     quaternion 1
+ * @param[in]   q     quaternion 2
  * @param[out]  dest  result quaternion
  */
 CGLM_INLINE
 void
-glm_quat_mulv(versor q1, versor q2, versor dest) {
-  dest[0] = q2[0] * q1[0] - q2[1] * q1[1] - q2[2] * q1[2] - q2[3] * q1[3];
-  dest[1] = q2[0] * q1[1] + q2[1] * q1[0] - q2[2] * q1[3] + q2[3] * q1[2];
-  dest[2] = q2[0] * q1[2] + q2[1] * q1[3] + q2[2] * q1[0] - q2[3] * q1[1];
-  dest[3] = q2[0] * q1[3] - q2[1] * q1[2] + q2[2] * q1[1] + q2[3] * q1[0];
-
-  glm_quat_normalize(dest);
+glm_quat_mul(versor p, versor q, versor dest) {
+  /*
+    + (a1 b2 + b1 a2 + c1 d2 − d1 c2)i
+    + (a1 c2 − b1 d2 + c1 a2 + d1 b2)j
+    + (a1 d2 + b1 c2 − c1 b2 + d1 a2)k
+       a1 a2 − b1 b2 − c1 c2 − d1 d2
+   */
+#if defined( __SSE__ ) || defined( __SSE2__ )
+  glm_quat_mul_sse2(p, q, dest);
+#else
+  dest[0] = p[3] * q[0] + p[0] * q[3] + p[1] * q[2] - p[2] * q[1];
+  dest[1] = p[3] * q[1] - p[0] * q[2] + p[1] * q[3] + p[2] * q[0];
+  dest[2] = p[3] * q[2] + p[0] * q[1] - p[1] * q[0] + p[2] * q[3];
+  dest[3] = p[3] * q[3] - p[0] * q[0] - p[1] * q[1] - p[2] * q[2];
+#endif
 }
 
 /*!
@@ -181,19 +403,22 @@ glm_quat_mulv(versor q1, versor q2, versor dest) {
 CGLM_INLINE
 void
 glm_quat_mat4(versor q, mat4 dest) {
-  float w, x, y, z;
-  float xx, yy, zz;
-  float xy, yz, xz;
-  float wx, wy, wz;
+  float w, x, y, z,
+        xx, yy, zz,
+        xy, yz, xz,
+        wx, wy, wz, norm, s;
 
-  w = q[0];
-  x = q[1];
-  y = q[2];
-  z = q[3];
+  norm = glm_quat_norm(q);
+  s    = norm > 0.0f ? 2.0f / norm : 0.0f;
 
-  xx = 2.0f * x * x;   xy = 2.0f * x * y;   wx = 2.0f * w * x;
-  yy = 2.0f * y * y;   yz = 2.0f * y * z;   wy = 2.0f * w * y;
-  zz = 2.0f * z * z;   xz = 2.0f * x * z;   wz = 2.0f * w * z;
+  x = q[0];
+  y = q[1];
+  z = q[2];
+  w = q[3];
+
+  xx = s * x * x;   xy = s * x * y;   wx = s * w * x;
+  yy = s * y * y;   yz = s * y * z;   wy = s * w * y;
+  zz = s * z * z;   xz = s * x * z;   wz = s * w * z;
 
   dest[0][0] = 1.0f - yy - zz;
   dest[1][1] = 1.0f - xx - zz;
@@ -207,8 +432,8 @@ glm_quat_mat4(versor q, mat4 dest) {
   dest[2][1] = yz - wx;
   dest[0][2] = xz - wy;
 
-  dest[1][3] = 0.0f;
   dest[0][3] = 0.0f;
+  dest[1][3] = 0.0f;
   dest[2][3] = 0.0f;
   dest[3][0] = 0.0f;
   dest[3][1] = 0.0f;
@@ -216,69 +441,303 @@ glm_quat_mat4(versor q, mat4 dest) {
   dest[3][3] = 1.0f;
 }
 
+/*!
+ * @brief convert quaternion to mat4 (transposed)
+ *
+ * @param[in]   q     quaternion
+ * @param[out]  dest  result matrix as transposed
+ */
+CGLM_INLINE
+void
+glm_quat_mat4t(versor q, mat4 dest) {
+  float w, x, y, z,
+        xx, yy, zz,
+        xy, yz, xz,
+        wx, wy, wz, norm, s;
+
+  norm = glm_quat_norm(q);
+  s    = norm > 0.0f ? 2.0f / norm : 0.0f;
+
+  x = q[0];
+  y = q[1];
+  z = q[2];
+  w = q[3];
+
+  xx = s * x * x;   xy = s * x * y;   wx = s * w * x;
+  yy = s * y * y;   yz = s * y * z;   wy = s * w * y;
+  zz = s * z * z;   xz = s * x * z;   wz = s * w * z;
+
+  dest[0][0] = 1.0f - yy - zz;
+  dest[1][1] = 1.0f - xx - zz;
+  dest[2][2] = 1.0f - xx - yy;
+
+  dest[1][0] = xy + wz;
+  dest[2][1] = yz + wx;
+  dest[0][2] = xz + wy;
+
+  dest[0][1] = xy - wz;
+  dest[1][2] = yz - wx;
+  dest[2][0] = xz - wy;
+
+  dest[0][3] = 0.0f;
+  dest[1][3] = 0.0f;
+  dest[2][3] = 0.0f;
+  dest[3][0] = 0.0f;
+  dest[3][1] = 0.0f;
+  dest[3][2] = 0.0f;
+  dest[3][3] = 1.0f;
+}
+
+/*!
+ * @brief convert quaternion to mat3
+ *
+ * @param[in]   q     quaternion
+ * @param[out]  dest  result matrix
+ */
+CGLM_INLINE
+void
+glm_quat_mat3(versor q, mat3 dest) {
+  float w, x, y, z,
+        xx, yy, zz,
+        xy, yz, xz,
+        wx, wy, wz, norm, s;
+
+  norm = glm_quat_norm(q);
+  s    = norm > 0.0f ? 2.0f / norm : 0.0f;
+
+  x = q[0];
+  y = q[1];
+  z = q[2];
+  w = q[3];
+
+  xx = s * x * x;   xy = s * x * y;   wx = s * w * x;
+  yy = s * y * y;   yz = s * y * z;   wy = s * w * y;
+  zz = s * z * z;   xz = s * x * z;   wz = s * w * z;
+
+  dest[0][0] = 1.0f - yy - zz;
+  dest[1][1] = 1.0f - xx - zz;
+  dest[2][2] = 1.0f - xx - yy;
+
+  dest[0][1] = xy + wz;
+  dest[1][2] = yz + wx;
+  dest[2][0] = xz + wy;
+
+  dest[1][0] = xy - wz;
+  dest[2][1] = yz - wx;
+  dest[0][2] = xz - wy;
+}
+
+/*!
+ * @brief convert quaternion to mat3 (transposed)
+ *
+ * @param[in]   q     quaternion
+ * @param[out]  dest  result matrix
+ */
+CGLM_INLINE
+void
+glm_quat_mat3t(versor q, mat3 dest) {
+  float w, x, y, z,
+        xx, yy, zz,
+        xy, yz, xz,
+        wx, wy, wz, norm, s;
+
+  norm = glm_quat_norm(q);
+  s    = norm > 0.0f ? 2.0f / norm : 0.0f;
+
+  x = q[0];
+  y = q[1];
+  z = q[2];
+  w = q[3];
+
+  xx = s * x * x;   xy = s * x * y;   wx = s * w * x;
+  yy = s * y * y;   yz = s * y * z;   wy = s * w * y;
+  zz = s * z * z;   xz = s * x * z;   wz = s * w * z;
+
+  dest[0][0] = 1.0f - yy - zz;
+  dest[1][1] = 1.0f - xx - zz;
+  dest[2][2] = 1.0f - xx - yy;
+
+  dest[1][0] = xy + wz;
+  dest[2][1] = yz + wx;
+  dest[0][2] = xz + wy;
+
+  dest[0][1] = xy - wz;
+  dest[1][2] = yz - wx;
+  dest[2][0] = xz - wy;
+}
+
+/*!
+ * @brief interpolates between two quaternions
+ *        using linear interpolation (LERP)
+ *
+ * @param[in]   from  from
+ * @param[in]   to    to
+ * @param[in]   t     interpolant (amount) clamped between 0 and 1
+ * @param[out]  dest  result quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_lerp(versor from, versor to, float t, versor dest) {
+  glm_vec4_lerp(from, to, t, dest);
+}
+
 /*!
  * @brief interpolates between two quaternions
  *        using spherical linear interpolation (SLERP)
  *
- * @param[in]   q     from
- * @param[in]   r     to
+ * @param[in]   from  from
+ * @param[in]   to    to
  * @param[in]   t     amout
  * @param[out]  dest  result quaternion
  */
 CGLM_INLINE
 void
-glm_quat_slerp(versor q,
-               versor r,
-               float  t,
-               versor dest) {
-  /* https://en.wikipedia.org/wiki/Slerp */
-#if defined( __SSE__ ) || defined( __SSE2__ )
-  glm_quat_slerp_sse2(q, r, t, dest);
-#else
-  float cosTheta, sinTheta, angle, a, b, c;
+glm_quat_slerp(versor from, versor to, float t, versor dest) {
+  vec4  q1, q2;
+  float cosTheta, sinTheta, angle;
 
-  cosTheta = glm_quat_dot(q, r);
-  if (cosTheta < 0.0f) {
-    q[0] *= -1.0f;
-    q[1] *= -1.0f;
-    q[2] *= -1.0f;
-    q[3] *= -1.0f;
+  cosTheta = glm_quat_dot(from, to);
+  glm_quat_copy(from, q1);
 
-    cosTheta = -cosTheta;
-  }
-
-  if (fabs(cosTheta) >= 1.0f) {
-    dest[0] = q[0];
-    dest[1] = q[1];
-    dest[2] = q[2];
-    dest[3] = q[3];
+  if (fabsf(cosTheta) >= 1.0f) {
+    glm_quat_copy(q1, dest);
     return;
   }
 
-  sinTheta = sqrt(1.0f - cosTheta * cosTheta);
+  if (cosTheta < 0.0f) {
+    glm_vec4_flipsign(q1);
+    cosTheta = -cosTheta;
+  }
 
-  c = 1.0f - t;
+  sinTheta = sqrtf(1.0f - cosTheta * cosTheta);
 
-  /* LERP */
-  /* TODO: FLT_EPSILON vs 0.001? */
-  if (sinTheta < 0.001f) {
-    dest[0] = c * q[0] + t * r[0];
-    dest[1] = c * q[1] + t * r[1];
-    dest[2] = c * q[2] + t * r[2];
-    dest[3] = c * q[3] + t * r[3];
+  /* LERP to avoid zero division */
+  if (fabsf(sinTheta) < 0.001f) {
+    glm_quat_lerp(from, to, t, dest);
     return;
   }
 
   /* SLERP */
   angle = acosf(cosTheta);
-  a = sinf(c * angle);
-  b = sinf(t * angle);
+  glm_vec4_scale(q1, sinf((1.0f - t) * angle), q1);
+  glm_vec4_scale(to, sinf(t * angle), q2);
 
-  dest[0] = (q[0] * a + r[0] * b) / sinTheta;
-  dest[1] = (q[1] * a + r[1] * b) / sinTheta;
-  dest[2] = (q[2] * a + r[2] * b) / sinTheta;
-  dest[3] = (q[3] * a + r[3] * b) / sinTheta;
-#endif
+  glm_vec4_add(q1, q2, q1);
+  glm_vec4_scale(q1, 1.0f / sinTheta, dest);
+}
+
+/*!
+ * @brief creates view matrix using quaternion as camera orientation
+ *
+ * @param[in]   eye   eye
+ * @param[in]   ori   orientation in world space as quaternion
+ * @param[out]  dest  view matrix
+ */
+CGLM_INLINE
+void
+glm_quat_look(vec3 eye, versor ori, mat4 dest) {
+  vec4 t;
+
+  /* orientation */
+  glm_quat_mat4t(ori, dest);
+
+  /* translate */
+  glm_vec4(eye, 1.0f, t);
+  glm_mat4_mulv(dest, t, t);
+  glm_vec_flipsign_to(t, dest[3]);
+}
+
+/*!
+ * @brief creates look rotation quaternion
+ *
+ * @param[in]   dir   direction to look
+ * @param[in]   fwd   forward vector
+ * @param[in]   up    up vector
+ * @param[out]  dest  destination quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_for(vec3 dir, vec3 fwd, vec3 up, versor dest) {
+  vec3  axis;
+  float dot, angle;
+
+  dot = glm_vec_dot(dir, fwd);
+  if (fabsf(dot + 1.0f)  < 0.000001f) {
+    glm_quat_init(dest, up[0], up[1], up[2], CGLM_PI);
+    return;
+  }
+
+  if (fabsf(dot - 1.0f) < 0.000001f) {
+    glm_quat_identity(dest);
+    return;
+  }
+
+  angle = acosf(dot);
+  glm_cross(fwd, dir, axis);
+  glm_normalize(axis);
+
+  glm_quatv(dest, angle, axis);
+}
+
+/*!
+ * @brief creates look rotation quaternion using source and
+ *        destination positions p suffix stands for position
+ *
+ * @param[in]   from  source point
+ * @param[in]   to    destination point
+ * @param[in]   fwd   forward vector
+ * @param[in]   up    up vector
+ * @param[out]  dest  destination quaternion
+ */
+CGLM_INLINE
+void
+glm_quat_forp(vec3 from, vec3 to, vec3 fwd, vec3 up, versor dest) {
+  vec3 dir;
+  glm_vec_sub(to, from, dir);
+  glm_quat_for(dir, fwd, up, dest);
+}
+
+/*!
+ * @brief rotate vector using using quaternion
+ *
+ * @param[in]   q     quaternion
+ * @param[in]   v     vector to rotate
+ * @param[out]  dest  rotated vector
+ */
+CGLM_INLINE
+void
+glm_quat_rotatev(versor q, vec3 v, vec3 dest) {
+  versor p;
+  vec3   u, v1, v2;
+  float  s;
+
+  glm_quat_normalize_to(q, p);
+  glm_quat_imag(p, u);
+  s = glm_quat_real(p);
+
+  glm_vec_scale(u, 2.0f * glm_vec_dot(u, v), v1);
+  glm_vec_scale(v, s * s - glm_vec_dot(u, u), v2);
+  glm_vec_add(v1, v2, v1);
+
+  glm_vec_cross(u, v, v2);
+  glm_vec_scale(v2, 2.0f * s, v2);
+
+  glm_vec_add(v1, v2, dest);
+}
+
+/*!
+ * @brief rotate existing transform matrix using quaternion
+ *
+ * @param[in]   m     existing transform matrix
+ * @param[in]   q     quaternion
+ * @param[out]  dest  rotated matrix/transform
+ */
+CGLM_INLINE
+void
+glm_quat_rotate(mat4 m, versor q, mat4 dest) {
+  mat4 rot;
+  glm_quat_mat4(q, rot);
+  glm_mat4_mul(m, rot, dest);
 }
 
 #endif /* cglm_quat_h */

include/cglm/simd/intrin.h

@@ -30,6 +30,16 @@
 #  define _mm_shuffle2_ps(a, b, z0, y0, x0, w0, z1, y1, x1, w1)                \
      _mm_shuffle1_ps(_mm_shuffle_ps(a, b, _MM_SHUFFLE(z0, y0, x0, w0)),        \
                                     z1, y1, x1, w1)
+
+CGLM_INLINE
+__m128
+glm_simd_dot(__m128 a, __m128 b) {
+  __m128 x0;
+  x0 = _mm_mul_ps(a, b);
+  x0 = _mm_add_ps(x0, _mm_shuffle1_ps(x0, 1, 0, 3, 2));
+  return _mm_add_ps(x0, _mm_shuffle1_ps(x0, 0, 1, 0, 1));
+}
+
 #endif
 
 /* x86, x64 */

include/cglm/simd/sse2/quat.h

@@ -14,56 +14,33 @@
 
 CGLM_INLINE
 void
-glm_quat_slerp_sse2(versor q,
-                    versor r,
-                    float  t,
-                    versor dest) {
-  /* https://en.wikipedia.org/wiki/Slerp */
-  float cosTheta, sinTheta, angle, a, b, c;
+glm_quat_mul_sse2(versor p, versor q, versor dest) {
+  /*
+   + (a1 b2 + b1 a2 + c1 d2 − d1 c2)i
+   + (a1 c2 − b1 d2 + c1 a2 + d1 b2)j
+   + (a1 d2 + b1 c2 − c1 b2 + d1 a2)k
+     a1 a2 − b1 b2 − c1 c2 − d1 d2
+   */
 
-  __m128 xmm_q;
+  __m128 xp, xq, x0, r;
 
-  xmm_q = _mm_load_ps(q);
+  xp = _mm_load_ps(p); /* 3 2 1 0 */
+  xq = _mm_load_ps(q);
 
-  cosTheta = glm_vec4_dot(q, r);
-  if (cosTheta < 0.0f) {
-    _mm_store_ps(q,
-                 _mm_xor_ps(xmm_q,
-                            _mm_set1_ps(-0.f))) ;
+  r  = _mm_mul_ps(_mm_shuffle1_ps1(xp, 3), xq);
 
-    cosTheta = -cosTheta;
-  }
+  x0 = _mm_xor_ps(_mm_shuffle1_ps1(xp, 0), _mm_set_ps(-0.f, 0.f, -0.f, 0.f));
+  r  = _mm_add_ps(r, _mm_mul_ps(x0, _mm_shuffle1_ps(xq, 0, 1, 2, 3)));
 
-  if (cosTheta >= 1.0f) {
-    _mm_store_ps(dest, xmm_q);
-    return;
-  }
+  x0 = _mm_xor_ps(_mm_shuffle1_ps1(xp, 1), _mm_set_ps(-0.f, -0.f, 0.f, 0.f));
+  r  = _mm_add_ps(r, _mm_mul_ps(x0, _mm_shuffle1_ps(xq, 1, 0, 3, 2)));
 
-  sinTheta = sqrtf(1.0f - cosTheta * cosTheta);
+  x0 = _mm_xor_ps(_mm_shuffle1_ps1(xp, 2), _mm_set_ps(-0.f, 0.f, 0.f, -0.f));
+  r  = _mm_add_ps(r, _mm_mul_ps(x0, _mm_shuffle1_ps(xq, 2, 3, 0, 1)));
 
-  c = 1.0f - t;
-
-  /* LERP */
-  if (sinTheta < 0.001f) {
-    _mm_store_ps(dest, _mm_add_ps(_mm_mul_ps(_mm_set1_ps(c),
-                                             xmm_q),
-                                  _mm_mul_ps(_mm_set1_ps(t),
-                                             _mm_load_ps(r))));
-    return;
-  }
-
-  /* SLERP */
-  angle = acosf(cosTheta);
-  a = sinf(c * angle);
-  b = sinf(t * angle);
-
-  _mm_store_ps(dest,
-               _mm_div_ps(_mm_add_ps(_mm_mul_ps(_mm_set1_ps(a),
-                                                xmm_q),
-                                     _mm_mul_ps(_mm_set1_ps(b),
-                                                _mm_load_ps(r))),
-                          _mm_set1_ps(sinTheta)));
+  _mm_store_ps(dest, r);
 }
 
+
 #endif
 #endif /* cglm_quat_simd_h */

include/cglm/types.h

@@ -14,6 +14,7 @@
 #  define CGLM_ALIGN(X) __attribute((aligned(X)))
 #endif
 
+typedef float vec2[2];
 typedef float vec3[3];
 typedef int  ivec3[3];
 typedef CGLM_ALIGN(16) float vec4[4];

include/cglm/util.h

@@ -21,7 +21,9 @@
 #include "common.h"
 
 /*!
- * @brief get sign of 32 bit integer as +1 or -1
+ * @brief get sign of 32 bit integer as +1, -1, 0
+ *
+ * Important: It returns 0 for zero input
  *
  * @param val integer value
  */
@@ -31,6 +33,19 @@ glm_sign(int val) {
   return ((val >> 31) - (-val >> 31));
 }
 
+/*!
+ * @brief get sign of 32 bit float as +1, -1, 0
+ *
+ * Important: It returns 0 for zero/NaN input
+ *
+ * @param val float value
+ */
+CGLM_INLINE
+float
+glm_signf(float val) {
+  return (float)((val > 0.0f) - (val < 0.0f));
+}
+
 /*!
  * @brief convert degree to radians
  *
@@ -115,4 +130,32 @@ glm_max(float a, float b) {
   return b;
 }
 
+/*!
+ * @brief clamp a number between min and max
+ *
+ * @param[in] val    value to clamp
+ * @param[in] minVal minimum value
+ * @param[in] maxVal maximum value
+ */
+CGLM_INLINE
+float
+glm_clamp(float val, float minVal, float maxVal) {
+  return glm_min(glm_max(val, minVal), maxVal);
+}
+
+/*!
+ * @brief linear interpolation between two number
+ *
+ * formula:  from + s * (to - from)
+ *
+ * @param[in]   from from value
+ * @param[in]   to   to value
+ * @param[in]   t    interpolant (amount) clamped between 0 and 1
+ */
+CGLM_INLINE
+float
+glm_lerp(float from, float to, float t) {
+  return from + glm_clamp(t, 0.0f, 1.0f) * (to - from);
+}
+
 #endif /* cglm_util_h */

include/cglm/vec3-ext.h

@@ -26,6 +26,7 @@
 #define cglm_vec3_ext_h
 
 #include "common.h"
+#include "util.h"
 #include <stdbool.h>
 #include <math.h>
 #include <float.h>
@@ -160,4 +161,69 @@ glm_vec_min(vec3 v) {
   return min;
 }
 
+/*!
+ * @brief check if all items are NaN (not a number)
+ *        you should only use this in DEBUG mode or very critical asserts
+ *
+ * @param[in] v vector
+ */
+CGLM_INLINE
+bool
+glm_vec_isnan(vec3 v) {
+  return isnan(v[0]) || isnan(v[1]) || isnan(v[2]);
+}
+
+/*!
+ * @brief check if all items are INFINITY
+ *        you should only use this in DEBUG mode or very critical asserts
+ *
+ * @param[in] v vector
+ */
+CGLM_INLINE
+bool
+glm_vec_isinf(vec3 v) {
+  return isinf(v[0]) || isinf(v[1]) || isinf(v[2]);
+}
+
+/*!
+ * @brief check if all items are valid number
+ *        you should only use this in DEBUG mode or very critical asserts
+ *
+ * @param[in] v vector
+ */
+CGLM_INLINE
+bool
+glm_vec_isvalid(vec3 v) {
+  return !glm_vec_isnan(v) && !glm_vec_isinf(v);
+}
+
+/*!
+ * @brief get sign of 32 bit float as +1, -1, 0
+ *
+ * Important: It returns 0 for zero/NaN input
+ *
+ * @param v vector
+ */
+CGLM_INLINE
+void
+glm_vec_sign(vec3 v, vec3 dest) {
+  dest[0] = glm_signf(v[0]);
+  dest[1] = glm_signf(v[1]);
+  dest[2] = glm_signf(v[2]);
+}
+
+/*!
+ * @brief square root of each vector item
+ *
+ * @param[in]  v    vector
+ * @param[out] dest destination vector
+ */
+CGLM_INLINE
+void
+glm_vec_sqrt(vec3 v, vec3 dest) {
+  dest[0] = sqrtf(v[0]);
+  dest[1] = sqrtf(v[1]);
+  dest[2] = sqrtf(v[2]);
+}
+
 #endif /* cglm_vec3_ext_h */

include/cglm/vec3.h

@@ -46,6 +46,7 @@
    CGLM_INLINE void  glm_vec_maxv(vec3 v1, vec3 v2, vec3 dest);
    CGLM_INLINE void  glm_vec_minv(vec3 v1, vec3 v2, vec3 dest);
    CGLM_INLINE void  glm_vec_ortho(vec3 v, vec3 dest);
+   CGLM_INLINE void  glm_vec_clamp(vec3 v, float minVal, float maxVal);
 
  Convenient:
    CGLM_INLINE void  glm_cross(vec3 a, vec3 b, vec3 d);
@@ -146,7 +147,7 @@ glm_vec_cross(vec3 a, vec3 b, vec3 d) {
 CGLM_INLINE
 float
 glm_vec_norm2(vec3 v) {
-  return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
+  return glm_vec_dot(v, v);
 }
 
 /*!
@@ -241,6 +242,20 @@ glm_vec_flipsign(vec3 v) {
   v[2] = -v[2];
 }
 
+/*!
+ * @brief flip sign of all vec3 members and store result in dest
+ *
+ * @param[in]   v     vector
+ * @param[out]  dest  result vector
+ */
+CGLM_INLINE
+void
+glm_vec_flipsign_to(vec3 v, vec3 dest) {
+  dest[0] = -v[0];
+  dest[1] = -v[1];
+  dest[2] = -v[2];
+}
+
 /*!
  * @brief make vector as inverse/opposite of itself
  *
@@ -324,12 +339,6 @@ glm_vec_angle(vec3 v1, vec3 v2) {
   return acosf(glm_vec_dot(v1, v2) * norm);
 }
 
-CGLM_INLINE
-void
-glm_quatv(versor q,
-          float  angle,
-          vec3   v);
-
 /*!
  * @brief rotate vec3 around axis by angle using Rodrigues' rotation formula
  *
@@ -340,31 +349,26 @@ glm_quatv(versor q,
 CGLM_INLINE
 void
 glm_vec_rotate(vec3 v, float angle, vec3 axis) {
-  versor q;
-  vec3   v1, v2, v3;
+  vec3   v1, v2, k;
   float  c, s;
 
   c = cosf(angle);
   s = sinf(angle);
 
+  glm_vec_normalize_to(axis, k);
+
   /* Right Hand, Rodrigues' rotation formula:
         v = v*cos(t) + (kxv)sin(t) + k*(k.v)(1 - cos(t))
    */
-
-  /* quaternion */
-  glm_quatv(q, angle, v);
-
   glm_vec_scale(v, c, v1);
 
-  glm_vec_cross(axis, v, v2);
+  glm_vec_cross(k, v, v2);
   glm_vec_scale(v2, s, v2);
 
-  glm_vec_scale(axis,
-                glm_vec_dot(axis, v) * (1.0f - c),
-                v3);
-
   glm_vec_add(v1, v2, v1);
-  glm_vec_add(v1, v3, v);
+
+  glm_vec_scale(k, glm_vec_dot(k, v) * (1.0f - c), v2);
+  glm_vec_add(v1, v2, v);
 }
 
 /*!
@@ -478,6 +482,43 @@ glm_vec_ortho(vec3 v, vec3 dest) {
   dest[2] = v[0] - v[1];
 }
 
+/*!
+ * @brief clamp vector's individual members between min and max values
+ *
+ * @param[in, out]  v      vector
+ * @param[in]       minVal minimum value
+ * @param[in]       maxVal maximum value
+ */
+CGLM_INLINE
+void
+glm_vec_clamp(vec3 v, float minVal, float maxVal) {
+  v[0] = glm_clamp(v[0], minVal, maxVal);
+  v[1] = glm_clamp(v[1], minVal, maxVal);
+  v[2] = glm_clamp(v[2], minVal, maxVal);
+}
+
+/*!
+ * @brief linear interpolation between two vector
+ *
+ * formula:  from + s * (to - from)
+ *
+ * @param[in]   from from value
+ * @param[in]   to   to value
+ * @param[in]   t    interpolant (amount) clamped between 0 and 1
+ * @param[out]  dest destination
+ */
+CGLM_INLINE
+void
+glm_vec_lerp(vec3 from, vec3 to, float t, vec3 dest) {
+  vec3 s, v;
+
+  /* from + s * (to - from) */
+  glm_vec_broadcast(glm_clamp(t, 0.0f, 1.0f), s);
+  glm_vec_sub(to, from, v);
+  glm_vec_mulv(s, v, v);
+  glm_vec_add(from, v, dest);
+}
+
 /*!
  * @brief vec3 cross product
  *

include/cglm/vec4-ext.h

@@ -174,5 +174,88 @@ glm_vec4_min(vec4 v) {
   return min;
 }
 
-#endif /* cglm_vec4_ext_h */
+/*!
+ * @brief check if one of items is NaN (not a number)
+ *        you should only use this in DEBUG mode or very critical asserts
+ *
+ * @param[in] v vector
+ */
+CGLM_INLINE
+bool
+glm_vec4_isnan(vec4 v) {
+  return isnan(v[0]) || isnan(v[1]) || isnan(v[2]) || isnan(v[3]);
+}
 
+/*!
+ * @brief check if one of items is INFINITY
+ *        you should only use this in DEBUG mode or very critical asserts
+ *
+ * @param[in] v vector
+ */
+CGLM_INLINE
+bool
+glm_vec4_isinf(vec4 v) {
+  return isinf(v[0]) || isinf(v[1]) || isinf(v[2]) || isinf(v[3]);
+}
+
+/*!
+ * @brief check if all items are valid number
+ *        you should only use this in DEBUG mode or very critical asserts
+ *
+ * @param[in] v vector
+ */
+CGLM_INLINE
+bool
+glm_vec4_isvalid(vec4 v) {
+  return !glm_vec4_isnan(v) && !glm_vec4_isinf(v);
+}
+
+/*!
+ * @brief get sign of 32 bit float as +1, -1, 0
+ *
+ * Important: It returns 0 for zero/NaN input
+ *
+ * @param v vector
+ */
+CGLM_INLINE
+void
+glm_vec4_sign(vec4 v, vec4 dest) {
+#if defined( __SSE2__ ) || defined( __SSE2__ )
+  __m128 x0, x1, x2, x3, x4;
+
+  x0 = _mm_load_ps(v);
+  x1 = _mm_set_ps(0.0f, 0.0f, 1.0f, -1.0f);
+  x2 = _mm_shuffle1_ps1(x1, 2);
+
+  x3 = _mm_and_ps(_mm_cmpgt_ps(x0, x2), _mm_shuffle1_ps1(x1, 1));
+  x4 = _mm_and_ps(_mm_cmplt_ps(x0, x2), _mm_shuffle1_ps1(x1, 0));
+
+  _mm_store_ps(dest, _mm_or_ps(x3, x4));
+#else
+  dest[0] = glm_signf(v[0]);
+  dest[1] = glm_signf(v[1]);
+  dest[2] = glm_signf(v[2]);
+  dest[3] = glm_signf(v[3]);
+#endif
+}
+
+/*!
+ * @brief square root of each vector item
+ *
+ * @param[in]  v    vector
+ * @param[out] dest destination vector
+ */
+CGLM_INLINE
+void
+glm_vec4_sqrt(vec4 v, vec4 dest) {
+#if defined( __SSE__ ) || defined( __SSE2__ )
+  _mm_store_ps(dest, _mm_sqrt_ps(_mm_load_ps(v)));
+#else
+  dest[0] = sqrtf(v[0]);
+  dest[1] = sqrtf(v[1]);
+  dest[2] = sqrtf(v[2]);
+  dest[3] = sqrtf(v[3]);
+#endif
+}
+
+#endif /* cglm_vec4_ext_h */

include/cglm/vec4.h

@@ -40,6 +40,7 @@
    CGLM_INLINE float glm_vec4_distance(vec4 v1, vec4 v2);
    CGLM_INLINE void  glm_vec4_maxv(vec4 v1, vec4 v2, vec4 dest);
    CGLM_INLINE void  glm_vec4_minv(vec4 v1, vec4 v2, vec4 dest);
+   CGLM_INLINE void  glm_vec4_clamp(vec4 v, float minVal, float maxVal);
  */
 
 #ifndef cglm_vec4_h
@@ -121,7 +122,14 @@ glm_vec4_copy(vec4 v, vec4 dest) {
 CGLM_INLINE
 float
 glm_vec4_dot(vec4 a, vec4 b) {
+#if defined( __SSE__ ) || defined( __SSE2__ )
+  __m128 x0;
+  x0 = _mm_mul_ps(_mm_load_ps(a), _mm_load_ps(b));
+  x0 = _mm_add_ps(x0, _mm_shuffle1_ps(x0, 1, 0, 3, 2));
+  return _mm_cvtss_f32(_mm_add_ss(x0, _mm_shuffle1_ps(x0, 0, 1, 0, 1)));
+#else
   return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
+#endif
 }
 
 /*!
@@ -138,7 +146,7 @@ glm_vec4_dot(vec4 a, vec4 b) {
 CGLM_INLINE
 float
 glm_vec4_norm2(vec4 v) {
-  return v[0] * v[0] + v[1] * v[1] + v[2] * v[2] + v[3] * v[3];
+  return glm_vec4_dot(v, v);
 }
 
 /*!
@@ -260,6 +268,26 @@ glm_vec4_flipsign(vec4 v) {
 #endif
 }
 
+/*!
+ * @brief flip sign of all vec4 members and store result in dest
+ *
+ * @param[in]  v     vector
+ * @param[out] dest  vector
+ */
+CGLM_INLINE
+void
+glm_vec4_flipsign_to(vec4 v, vec4 dest) {
+#if defined( __SSE__ ) || defined( __SSE2__ )
+  _mm_store_ps(dest, _mm_xor_ps(_mm_load_ps(v),
+                                _mm_set1_ps(-0.0f)));
+#else
+  dest[0] = -v[0];
+  dest[1] = -v[1];
+  dest[2] = -v[2];
+  dest[3] = -v[3];
+#endif
+}
+
 /*!
  * @brief make vector as inverse/opposite of itself
  *
@@ -373,4 +401,42 @@ glm_vec4_minv(vec4 v1, vec4 v2, vec4 dest) {
   dest[3] = glm_min(v1[3], v2[3]);
 }
 
+/*!
+ * @brief clamp vector's individual members between min and max values
+ *
+ * @param[in, out]  v      vector
+ * @param[in]       minVal minimum value
+ * @param[in]       maxVal maximum value
+ */
+CGLM_INLINE
+void
+glm_vec4_clamp(vec4 v, float minVal, float maxVal) {
+  v[0] = glm_clamp(v[0], minVal, maxVal);
+  v[1] = glm_clamp(v[1], minVal, maxVal);
+  v[2] = glm_clamp(v[2], minVal, maxVal);
+  v[3] = glm_clamp(v[3], minVal, maxVal);
+}
+
+/*!
+ * @brief linear interpolation between two vector
+ *
+ * formula:  from + s * (to - from)
+ *
+ * @param[in]   from from value
+ * @param[in]   to   to value
+ * @param[in]   t    interpolant (amount) clamped between 0 and 1
+ * @param[out]  dest destination
+ */
+CGLM_INLINE
+void
+glm_vec4_lerp(vec4 from, vec4 to, float t, vec4 dest) {
+  vec4 s, v;
+
+  /* from + s * (to - from) */
+  glm_vec4_broadcast(glm_clamp(t, 0.0f, 1.0f), s);
+  glm_vec4_sub(to, from, v);
+  glm_vec4_mulv(s, v, v);
+  glm_vec4_add(from, v, dest);
+}
+
 #endif /* cglm_vec4_h */

include/cglm/version.h

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

makefile.am

@@ -98,13 +98,19 @@ libcglm_la_SOURCES=\
     src/mat4.c \
     src/plane.c \
     src/frustum.c \
-    src/box.c
+    src/box.c \
+    src/project.c
 
 test_tests_SOURCES=\
     test/src/test_common.c \
     test/src/test_main.c \
     test/src/test_mat4.c \
-    test/src/test_cam.c
+    test/src/test_cam.c \
+    test/src/test_project.c \
+    test/src/test_clamp.c \
+    test/src/test_euler.c \
+    test/src/test_quat.c \
+    test/src/test_vec4.c
 
 all-local:
 	sh ./post-build.sh

src/euler.c

@@ -20,6 +20,12 @@ glmc_euler(vec3 angles, mat4 dest) {
   glm_euler(angles, dest);
 }
 
+CGLM_EXPORT
+void
+glmc_euler_xyz(vec3 angles,  mat4 dest) {
+  glm_euler_xyz(angles, dest);
+}
+
 CGLM_EXPORT
 void
 glmc_euler_zyx(vec3 angles,  mat4 dest) {

src/mat4.c

@@ -52,7 +52,7 @@ glmc_mat4_mul(mat4 m1, mat4 m2, mat4 dest) {
 
 CGLM_EXPORT
 void
-glmc_mat4_mulN(mat4 * __restrict matrices[], int len, mat4 dest) {
+glmc_mat4_mulN(mat4 * __restrict matrices[], uint32_t len, mat4 dest) {
   glm_mat4_mulN(matrices, len, dest);
 }
 
@@ -62,6 +62,12 @@ glmc_mat4_mulv(mat4 m, vec4 v, vec4 dest) {
   glm_mat4_mulv(m, v, dest);
 }
 
+CGLM_EXPORT
+void
+glmc_mat4_quat(mat4 m, versor dest) {
+  glm_mat4_quat(m, dest);
+}
+
 CGLM_EXPORT
 void
 glmc_mat4_transpose_to(mat4 m, mat4 dest) {

src/project.c

@@ -0,0 +1,27 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "../include/cglm/cglm.h"
+#include "../include/cglm/call.h"
+
+CGLM_EXPORT
+void
+glmc_unprojecti(vec3 pos, mat4 invMat, vec4 vp, vec3 dest) {
+  glm_unprojecti(pos, invMat, vp, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_unproject(vec3 pos, mat4 m, vec4 vp, vec3 dest) {
+  glm_unproject(pos, m, vp, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_project(vec3 pos, mat4 m, vec4 vp, vec3 dest) {
+  glm_project(pos, m, vp, dest);
+}

src/quat.c

@@ -8,6 +8,7 @@
 #include "../include/cglm/cglm.h"
 #include "../include/cglm/call.h"
 
+
 CGLM_EXPORT
 void
 glmc_quat_identity(versor q) {
@@ -16,20 +17,26 @@ glmc_quat_identity(versor q) {
 
 CGLM_EXPORT
 void
-glmc_quat(versor q,
-               float angle,
-               float x,
-               float y,
-               float z) {
+glmc_quat_init(versor q, float x, float y, float z, float w) {
+  glm_quat_init(q, x, y, z, w);
+}
+
+CGLM_EXPORT
+void
+glmc_quat(versor q, float angle, float x, float y, float z) {
   glm_quat(q, angle, x, y, z);
 }
 
 CGLM_EXPORT
 void
-glmc_quatv(versor q,
-           float  angle,
-           vec3   v) {
-  glm_quatv(q, angle, v);
+glmc_quatv(versor q, float angle, vec3 axis) {
+  glm_quatv(q, angle, axis);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_copy(versor q, versor dest) {
+  glm_quat_copy(q, dest);
 }
 
 CGLM_EXPORT
@@ -40,20 +47,86 @@ glmc_quat_norm(versor q) {
 
 CGLM_EXPORT
 void
-glmc_quat_normalize(versor q) {
-  glm_quat_normalize(q);
-}
-
-CGLM_EXPORT
-float
-glmc_quat_dot(versor q, versor r) {
-  return glm_quat_dot(q, r);
+glmc_quat_normalize_to(versor q, versor dest) {
+  glm_quat_normalize_to(q, dest);
 }
 
 CGLM_EXPORT
 void
-glmc_quat_mulv(versor q1, versor q2, versor dest) {
-  glm_quat_mulv(q1, q2, dest);
+glmc_quat_normalize(versor q) {
+  glm_quat_norm(q);
+}
+
+CGLM_EXPORT
+float
+glmc_quat_dot(versor p, versor q) {
+  return glm_quat_dot(p, q);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_conjugate(versor q, versor dest) {
+  glm_quat_conjugate(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_inv(versor q, versor dest) {
+  glm_quat_inv(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_add(versor p, versor q, versor dest) {
+  glm_quat_add(p, q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_sub(versor p, versor q, versor dest) {
+  glm_quat_sub(p, q, dest);
+}
+
+CGLM_EXPORT
+float
+glmc_quat_real(versor q) {
+  return glm_quat_real(q);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_imag(versor q, vec3 dest) {
+  glm_quat_imag(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_imagn(versor q, vec3 dest) {
+  glm_quat_imagn(q, dest);
+}
+
+CGLM_EXPORT
+float
+glmc_quat_imaglen(versor q) {
+  return glm_quat_imaglen(q);
+}
+
+CGLM_EXPORT
+float
+glmc_quat_angle(versor q) {
+  return glm_quat_angle(q);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_axis(versor q, versor dest) {
+  glm_quat_axis(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_mul(versor p, versor q, versor dest) {
+  glm_quat_mul(p, q, dest);
 }
 
 CGLM_EXPORT
@@ -64,9 +137,60 @@ glmc_quat_mat4(versor q, mat4 dest) {
 
 CGLM_EXPORT
 void
-glmc_quat_slerp(versor q,
-                versor r,
-                float  t,
-                versor dest) {
-  glm_quat_slerp(q, r, t, dest);
+glmc_quat_mat4t(versor q, mat4 dest) {
+  glm_quat_mat4t(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_mat3(versor q, mat3 dest) {
+  glm_quat_mat3(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_mat3t(versor q, mat3 dest) {
+  glm_quat_mat3t(q, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_lerp(versor from, versor to, float t, versor dest) {
+  glm_quat_lerp(from, to, t, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_slerp(versor from, versor to, float t, versor dest) {
+  glm_quat_slerp(from, to, t, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_look(vec3 eye, versor ori, mat4 dest) {
+  glm_quat_look(eye, ori, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_for(vec3 dir, vec3 fwd, vec3 up, versor dest) {
+  glm_quat_for(dir, fwd, up, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_forp(vec3 from, vec3 to, vec3 fwd, vec3 up, versor dest) {
+  glm_quat_forp(from, to, fwd, up, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_rotatev(versor q, vec3 v, vec3 dest) {
+  glm_quat_rotatev(q, v, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_quat_rotate(mat4 m, versor q, mat4 dest) {
+  glm_quat_rotate(m, q, dest);
 }

src/vec3.c

@@ -8,6 +8,12 @@
 #include "../include/cglm/cglm.h"
 #include "../include/cglm/call.h"
 
+CGLM_EXPORT
+void
+glmc_vec3(vec4 v4, vec3 dest) {
+  glm_vec3(v4, dest);
+}
+
 CGLM_EXPORT
 void
 glmc_vec_copy(vec3 a, vec3 dest) {
@@ -80,6 +86,12 @@ glmc_vec_flipsign(vec3 v) {
   glm_vec_flipsign(v);
 }
 
+CGLM_EXPORT
+void
+glmc_vec_flipsign_to(vec3 v, vec3 dest) {
+  glm_vec_flipsign_to(v, dest);
+}
+
 CGLM_EXPORT
 void
 glmc_vec_inv(vec3 v) {
@@ -139,3 +151,107 @@ void
 glmc_vec_minv(vec3 v1, vec3 v2, vec3 dest) {
   glm_vec_maxv(v1, v2, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec_clamp(vec3 v, float minVal, float maxVal) {
+  glm_vec_clamp(v, minVal, maxVal);
+}
+
+CGLM_EXPORT
+void
+glmc_vec_ortho(vec3 v, vec3 dest) {
+  glm_vec_ortho(v, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec_lerp(vec3 from, vec3 to, float t, vec3 dest) {
+  glm_vec_lerp(from, to, t, dest);
+}
+
+/* ext */
+
+CGLM_EXPORT
+void
+glmc_vec_mulv(vec3 a, vec3 b, vec3 d) {
+  glm_vec_mulv(a, b, d);
+}
+
+CGLM_EXPORT
+void
+glmc_vec_broadcast(float val, vec3 d) {
+  glm_vec_broadcast(val, d);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_eq(vec3 v, float val) {
+  return glm_vec_eq(v, val);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_eq_eps(vec3 v, float val) {
+  return glm_vec_eq_eps(v, val);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_eq_all(vec3 v) {
+  return glm_vec_eq_all(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_eqv(vec3 v1, vec3 v2) {
+  return glm_vec_eqv(v1, v2);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_eqv_eps(vec3 v1, vec3 v2) {
+  return glm_vec_eqv_eps(v1, v2);
+}
+
+CGLM_EXPORT
+float
+glmc_vec_max(vec3 v) {
+  return glm_vec_max(v);
+}
+
+CGLM_EXPORT
+float
+glmc_vec_min(vec3 v) {
+  return glm_vec_min(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_isnan(vec3 v) {
+  return glm_vec_isnan(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_isinf(vec3 v) {
+  return glm_vec_isinf(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec_isvalid(vec3 v) {
+  return glm_vec_isvalid(v);
+}
+
+CGLM_EXPORT
+void
+glmc_vec_sign(vec3 v, vec3 dest) {
+  glm_vec_sign(v, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec_sqrt(vec3 v, vec3 dest) {
+  glm_vec_sqrt(v, dest);
+}

src/vec4.c

@@ -8,6 +8,12 @@
 #include "../include/cglm/cglm.h"
 #include "../include/cglm/call.h"
 
+CGLM_EXPORT
+void
+glmc_vec4(vec3 v3, float last, vec4 dest) {
+  glm_vec4(v3, last, dest);
+}
+
 CGLM_EXPORT
 void
 glmc_vec4_copy3(vec4 a, vec3 dest) {
@@ -80,6 +86,12 @@ glmc_vec4_flipsign(vec4 v) {
   glm_vec4_flipsign(v);
 }
 
+CGLM_EXPORT
+void
+glmc_vec4_flipsign_to(vec4 v, vec4 dest) {
+  glm_vec4_flipsign_to(v, dest);
+}
+
 CGLM_EXPORT
 void
 glmc_vec4_inv(vec4 v) {
@@ -109,3 +121,101 @@ void
 glmc_vec4_minv(vec4 v1, vec4 v2, vec4 dest) {
   glm_vec4_maxv(v1, v2, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec4_clamp(vec4 v, float minVal, float maxVal) {
+  glm_vec4_clamp(v, minVal, maxVal);
+}
+
+CGLM_EXPORT
+void
+glmc_vec4_lerp(vec4 from, vec4 to, float t, vec4 dest) {
+  glm_vec4_lerp(from, to, t, dest);
+}
+
+/* ext */
+
+CGLM_EXPORT
+void
+glmc_vec4_mulv(vec4 a, vec4 b, vec4 d) {
+  glm_vec4_mulv(a, b, d);
+}
+
+CGLM_EXPORT
+void
+glmc_vec4_broadcast(float val, vec4 d) {
+  glm_vec4_broadcast(val, d);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_eq(vec4 v, float val) {
+  return glm_vec4_eq(v, val);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_eq_eps(vec4 v, float val) {
+  return glm_vec4_eq_eps(v, val);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_eq_all(vec4 v) {
+  return glm_vec4_eq_all(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_eqv(vec4 v1, vec4 v2) {
+  return glm_vec4_eqv(v1, v2);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_eqv_eps(vec4 v1, vec4 v2) {
+  return glm_vec4_eqv_eps(v1, v2);
+}
+
+CGLM_EXPORT
+float
+glmc_vec4_max(vec4 v) {
+  return glm_vec4_max(v);
+}
+
+CGLM_EXPORT
+float
+glmc_vec4_min(vec4 v) {
+  return glm_vec4_min(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_isnan(vec4 v) {
+  return glm_vec4_isnan(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_isinf(vec4 v) {
+  return glm_vec4_isinf(v);
+}
+
+CGLM_EXPORT
+bool
+glmc_vec4_isvalid(vec4 v) {
+  return glm_vec4_isvalid(v);
+}
+
+CGLM_EXPORT
+void
+glmc_vec4_sign(vec4 v, vec4 dest) {
+  glm_vec4_sign(v, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec4_sqrt(vec4 v, vec4 dest) {
+  glm_vec4_sqrt(v, dest);
+}

test/src/test_clamp.c

@@ -0,0 +1,30 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "test_common.h"
+
+void
+test_clamp(void **state) {
+  vec3 v3 = {15.07, 0.4, 17.3};
+  vec4 v4 = {5.07,  2.3, 1.3, 1.4};
+
+  assert_true(glm_clamp(1.6f, 0.0f, 1.0f)  == 1.0f);
+  assert_true(glm_clamp(-1.6f, 0.0f, 1.0f) == 0.0f);
+  assert_true(glm_clamp(0.6f, 0.0f, 1.0f)  == 0.6f);
+
+  glm_vec_clamp(v3, 0.0, 1.0);
+  glm_vec4_clamp(v4, 1.5, 3.0);
+
+  assert_true(v3[0] == 1.0f);
+  assert_true(v3[1] == 0.4f);
+  assert_true(v3[2] == 1.0f);
+
+  assert_true(v4[0] == 3.0f);
+  assert_true(v4[1] == 2.3f);
+  assert_true(v4[2] == 1.5f);
+  assert_true(v4[3] == 1.5f);
+}

test/src/test_common.c

@@ -27,6 +27,39 @@ test_rand_mat4(mat4 dest) {
   /* glm_scale(dest, (vec3){drand48(), drand48(), drand48()}); */
 }
 
+void
+test_rand_vec3(vec3 dest) {
+  srand((unsigned int)time(NULL));
+
+  dest[0] = drand48();
+  dest[1] = drand48();
+  dest[2] = drand48();
+}
+
+void
+test_rand_vec4(vec4 dest) {
+  srand((unsigned int)time(NULL));
+
+  dest[0] = drand48();
+  dest[1] = drand48();
+  dest[2] = drand48();
+  dest[3] = drand48();
+}
+
+float
+test_rand_angle(void) {
+  srand((unsigned int)time(NULL));
+  return drand48();
+}
+
+void
+test_rand_quat(versor q) {
+  srand((unsigned int)time(NULL));
+
+  glm_quat(q, drand48(), drand48(), drand48(), drand48());
+  glm_quat_normalize(q);
+}
+
 void
 test_assert_mat4_eq(mat4 m1, mat4 m2) {
   int i, j, k;
@@ -50,3 +83,27 @@ test_assert_mat4_eq2(mat4 m1, mat4 m2, float eps) {
     }
   }
 }
+
+void
+test_assert_vec3_eq(vec3 v1, vec3 v2) {
+  assert_true(fabsf(v1[0] - v2[0]) <= 0.000009); /* rounding errors */
+  assert_true(fabsf(v1[1] - v2[1]) <= 0.000009);
+  assert_true(fabsf(v1[2] - v2[2]) <= 0.000009);
+}
+
+void
+test_assert_quat_eq_abs(versor v1, versor v2) {
+  assert_true(fabsf(fabsf(v1[0]) - fabsf(v2[0])) <= 0.0009); /* rounding errors */
+  assert_true(fabsf(fabsf(v1[1]) - fabsf(v2[1])) <= 0.0009);
+  assert_true(fabsf(fabsf(v1[2]) - fabsf(v2[2])) <= 0.0009);
+  assert_true(fabsf(fabsf(v1[3]) - fabsf(v2[3])) <= 0.0009);
+}
+
+void
+test_assert_quat_eq(versor v1, versor v2) {
+  assert_true(fabsf(v1[0] - v2[0]) <= 0.000009); /* rounding errors */
+  assert_true(fabsf(v1[1] - v2[1]) <= 0.000009);
+  assert_true(fabsf(v1[2] - v2[2]) <= 0.000009);
+  assert_true(fabsf(v1[3] - v2[3]) <= 0.000009);
+}
+

test/src/test_common.h

@@ -31,4 +31,25 @@ test_assert_mat4_eq(mat4 m1, mat4 m2);
 void
 test_assert_mat4_eq2(mat4 m1, mat4 m2, float eps);
 
+void
+test_assert_vec3_eq(vec3 v1, vec3 v2);
+
+void
+test_assert_quat_eq(versor v1, versor v2);
+
+void
+test_assert_quat_eq_abs(versor v1, versor v2);
+
+void
+test_rand_vec3(vec3 dest);
+
+void
+test_rand_vec4(vec4 dest) ;
+
+float
+test_rand_angle(void);
+
+void
+test_rand_quat(versor q);
+
 #endif /* test_common_h */

test/src/test_euler.c

@@ -0,0 +1,44 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "test_common.h"
+
+void
+test_euler(void **state) {
+  mat4  rot1, rot2;
+  vec3  inAngles, outAngles;
+
+  inAngles[0] = glm_rad(-45.0f);  /* X angle */
+  inAngles[1] = glm_rad(88.0f);   /* Y angle */
+  inAngles[2] = glm_rad(18.0f);   /* Z angle */
+
+  glm_euler_xyz(inAngles, rot1);
+
+  /* extract angles */
+  glmc_euler_angles(rot1, outAngles);
+
+  /* angles must be equal in that range */
+  test_assert_vec3_eq(inAngles, outAngles);
+
+  /* matrices must be equal */
+  glmc_euler_xyz(outAngles, rot2);
+  test_assert_mat4_eq(rot1, rot2);
+
+  /* change range */
+  inAngles[0] = glm_rad(-145.0f);  /* X angle */
+  inAngles[1] = glm_rad(818.0f);   /* Y angle */
+  inAngles[2] = glm_rad(181.0f);   /* Z angle */
+
+  glm_euler_xyz(inAngles, rot1);
+  glmc_euler_angles(rot1, outAngles);
+
+  /* angles may not be equal but matrices MUST! */
+
+  /* matrices must be equal */
+  glmc_euler_xyz(outAngles, rot2);
+  test_assert_mat4_eq(rot1, rot2);
+}

test/src/test_main.c

@@ -14,7 +14,22 @@ main(int argc, const char * argv[]) {
 
     /* camera */
     cmocka_unit_test(test_camera_lookat),
-    cmocka_unit_test(test_camera_decomp)
+    cmocka_unit_test(test_camera_decomp),
+
+    /* project */
+    cmocka_unit_test(test_project),
+
+    /* vector */
+    cmocka_unit_test(test_clamp),
+
+    /* euler */
+    cmocka_unit_test(test_euler),
+
+    /* quaternion */
+    cmocka_unit_test(test_quat),
+
+    /* vec4 */
+    cmocka_unit_test(test_vec4)
   };
 
   return cmocka_run_group_tests(tests, NULL, NULL);

test/src/test_project.c

@@ -0,0 +1,31 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "test_common.h"
+
+void
+test_project(void **state) {
+  mat4 model, view, proj, mvp;
+  vec4 viewport = {0.0f, 0.0f, 800.0f, 600.0f};
+  vec3 pos      = {13.0f, 45.0f, 0.74f};
+  vec3 projected, unprojected;
+
+  glm_translate_make(model, (vec3){0.0f, 0.0f, -10.0f});
+  glm_lookat((vec3){0.0f, 0.0f, 0.0f}, pos, GLM_YUP, view);
+
+  glm_perspective_default(0.5f, proj);
+  glm_mat4_mulN((mat4 *[]){&proj, &view, &model}, 3, mvp);
+
+  glmc_project(pos, mvp, viewport, projected);
+  glmc_unproject(projected, mvp, viewport, unprojected);
+
+  /* unprojected of projected vector must be same as original one */
+  /* we used 0.01 because of projection floating point errors */
+  assert_true(fabsf(pos[0] - unprojected[0]) < 0.01);
+  assert_true(fabsf(pos[1] - unprojected[1]) < 0.01);
+  assert_true(fabsf(pos[2] - unprojected[2]) < 0.01);
+}

test/src/test_quat.c

@@ -0,0 +1,199 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "test_common.h"
+
+CGLM_INLINE
+void
+test_quat_mul_raw(versor p, versor q, versor dest) {
+  dest[0] = p[3] * q[0] + p[0] * q[3] + p[1] * q[2] - p[2] * q[1];
+  dest[1] = p[3] * q[1] - p[0] * q[2] + p[1] * q[3] + p[2] * q[0];
+  dest[2] = p[3] * q[2] + p[0] * q[1] - p[1] * q[0] + p[2] * q[3];
+  dest[3] = p[3] * q[3] - p[0] * q[0] - p[1] * q[1] - p[2] * q[2];
+}
+
+void
+test_quat(void **state) {
+  mat4   inRot, outRot, view1, view2, rot1, rot2;
+  versor inQuat, outQuat, q3, q4, q5;
+  vec3   eye, axis, imag, v1, v2;
+  int    i;
+
+  /* 0. test identiy quat */
+  glm_quat_identity(q4);
+  assert_true(glm_quat_real(q4) == cosf(glm_rad(0.0f) * 0.5f));
+  glm_quat_mat4(q4, rot1);
+  test_assert_mat4_eq2(rot1, GLM_MAT4_IDENTITY, 0.000009);
+
+  /* 1. test quat to mat and mat to quat */
+  for (i = 0; i < 1000; i++) {
+    test_rand_quat(inQuat);
+
+    glmc_quat_mat4(inQuat, inRot);
+    glmc_mat4_quat(inRot, outQuat);
+    glmc_quat_mat4(outQuat, outRot);
+
+    /* 2. test first quat and generated one equality */
+    test_assert_quat_eq_abs(inQuat, outQuat);
+
+    /* 3. test first rot and second rotation */
+    test_assert_mat4_eq2(inRot, outRot, 0.000009); /* almost equal */
+
+    /* 4. test SSE mul and raw mul */
+    test_quat_mul_raw(inQuat, outQuat, q3);
+    glm_quat_mul_sse2(inQuat, outQuat, q4);
+    test_assert_quat_eq(q3, q4);
+  }
+
+  /* 5. test lookat */
+  test_rand_vec3(eye);
+  glm_quatv(q3, glm_rad(-90.0f), GLM_YUP);
+
+  /* now X axis must be forward axis, Z must be right axis */
+  glm_look(eye, GLM_XUP, GLM_YUP, view1);
+
+  /* create view matrix with quaternion */
+  glm_quat_look(eye, q3, view2);
+
+  test_assert_mat4_eq2(view1, view2, 0.000009);
+
+  /* 6. test quaternion rotation matrix result */
+  test_rand_quat(q3);
+  glm_quat_mat4(q3, rot1);
+
+  /* 6.1 test axis and angle of quat */
+  glm_quat_axis(q3, axis);
+  glm_rotate_make(rot2, glm_quat_angle(q3), axis);
+
+  test_assert_mat4_eq2(rot1, rot2, 0.000009);
+
+  /* 7. test quaternion multiplication (hamilton product),
+        final rotation = first rotation + second = quat1 * quat2
+   */
+  test_rand_quat(q3);
+  test_rand_quat(q4);
+
+  glm_quat_mul(q3, q4, q5);
+
+  glm_quat_axis(q3, axis);
+  glm_rotate_make(rot1, glm_quat_angle(q3), axis);
+
+  glm_quat_axis(q4, axis);
+  glm_rotate(rot1, glm_quat_angle(q4), axis);
+
+  /* rot2 is combine of two rotation now test with quaternion result */
+  glm_quat_mat4(q5, rot2);
+
+  /* result must be same (almost) */
+  test_assert_mat4_eq2(rot1, rot2, 0.000009);
+
+  /* 8. test quaternion for look rotation */
+
+  /* 8.1 same direction */
+  /* look at from 0, 0, 1 to zero, direction = 0, 0, -1 */
+  glm_quat_for((vec3){0, 0, -1}, (vec3){0, 0, -1}, GLM_YUP, q3);
+
+  /* result must be identity */
+  glm_quat_identity(q4);
+  test_assert_quat_eq(q3, q4);
+
+  /* look at from 0, 0, 1 to zero, direction = 0, 0, -1 */
+  glm_quat_forp(GLM_ZUP, GLM_VEC3_ZERO, (vec3){0, 0, -1}, GLM_YUP, q3);
+
+  /* result must be identity */
+  glm_quat_identity(q4);
+  test_assert_quat_eq(q3, q4);
+
+  /* 8.2 perpendicular */
+  glm_quat_for(GLM_XUP, (vec3){0, 0, -1}, GLM_YUP, q3);
+
+  /* result must be -90 */
+  glm_quatv(q4, glm_rad(-90.0f), GLM_YUP);
+  test_assert_quat_eq(q3, q4);
+
+  /* 9. test imag, real */
+
+  /* 9.1 real */
+  assert_true(glm_quat_real(q4) == cosf(glm_rad(-90.0f) * 0.5f));
+
+  /* 9.1 imag */
+  glm_quat_imag(q4, imag);
+
+  /* axis = Y_UP * sinf(angle * 0.5), YUP = 0, 1, 0 */
+  axis[0] = 0.0f;
+  axis[1] = sinf(glm_rad(-90.0f) * 0.5f) * 1.0f;
+  axis[2] = 0.0f;
+
+  assert_true(glm_vec_eqv_eps(imag, axis));
+
+  /* 9.2 axis */
+  glm_quat_axis(q4, axis);
+  imag[0] =  0.0f;
+  imag[1] = -1.0f;
+  imag[2] =  0.0f;
+
+  test_assert_vec3_eq(imag, axis);
+
+  /* 10. test rotate vector using quat */
+  /* (0,0,-1) around (1,0,0) must give (0,1,0) */
+  v1[0] = 0.0f; v1[1] = 0.0f; v1[2] = -1.0f;
+  v2[0] = 0.0f; v2[1] = 0.0f; v2[2] = -1.0f;
+
+  glm_vec_rotate(v1, glm_rad(90.0f), (vec3){1.0f, 0.0f, 0.0f});
+  glm_quatv(q3, glm_rad(90.0f), (vec3){1.0f, 0.0f, 0.0f});
+
+  glm_vec4_scale(q3, 1.5, q3);
+  glm_quat_rotatev(q3, v2, v2);
+
+  /* result must be : (0,1,0) */
+  assert_true(fabsf(v1[0]) <= 0.00009f
+              && fabsf(v1[1] - 1.0f) <= 0.00009f
+              && fabsf(v1[2]) <= 0.00009f);
+
+  test_assert_vec3_eq(v1, v2);
+
+  /* 11. test rotate transform */
+  glm_translate_make(rot1, (vec3){-10.0, 45.0f, 8.0f});
+  glm_rotate(rot1, glm_rad(-90), GLM_ZUP);
+
+  glm_quatv(q3, glm_rad(-90.0f), GLM_ZUP);
+  glm_translate_make(rot2, (vec3){-10.0, 45.0f, 8.0f});
+  glm_quat_rotate(rot2, q3, rot2);
+
+  /* result must be same (almost) */
+  test_assert_mat4_eq2(rot1, rot2, 0.000009);
+
+  glm_rotate_make(rot1, glm_rad(-90), GLM_ZUP);
+  glm_translate(rot1, (vec3){-10.0, 45.0f, 8.0f});
+
+  glm_quatv(q3, glm_rad(-90.0f), GLM_ZUP);
+  glm_mat4_identity(rot2);
+  glm_quat_rotate(rot2, q3, rot2);
+  glm_translate(rot2, (vec3){-10.0, 45.0f, 8.0f});
+
+  /* result must be same (almost) */
+  test_assert_mat4_eq2(rot1, rot2, 0.000009);
+
+  /* reverse */
+  glm_rotate_make(rot1, glm_rad(-90), GLM_ZUP);
+  glm_quatv(q3, glm_rad(90.0f), GLM_ZUP);
+  glm_quat_rotate(rot1, q3, rot1);
+
+  /* result must be identity */
+  test_assert_mat4_eq2(rot1, GLM_MAT4_IDENTITY, 0.000009);
+
+  test_rand_quat(q3);
+
+  /* 12. inverse of quat, multiplication must be IDENTITY */
+  glm_quat_inv(q3, q4);
+  glm_quat_mul(q3, q4, q5);
+
+  glm_quat_identity(q3);
+  test_assert_quat_eq(q3, q5);
+
+  /* TODO: add tests for slerp, lerp */
+}

test/src/test_tests.h

@@ -16,4 +16,19 @@ test_camera_lookat(void **state);
 void
 test_camera_decomp(void **state);
 
+void
+test_project(void **state);
+
+void
+test_clamp(void **state);
+
+void
+test_euler(void **state);
+
+void
+test_quat(void **state);
+
+void
+test_vec4(void **state);
+
 #endif /* test_tests_h */

test/src/test_vec4.c

@@ -0,0 +1,30 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "test_common.h"
+
+CGLM_INLINE
+float
+test_vec4_dot(vec4 a, vec4 b) {
+  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
+}
+
+void
+test_vec4(void **state) {
+  vec4  v;
+  int   i;
+  float d1, d2;
+
+  /* test SSE/SIMD dot product */
+  for (i = 0; i < 100; i++) {
+    test_rand_vec4(v);
+    d1 = glm_vec4_dot(v, v);
+    d2 = test_vec4_dot(v, v);
+
+    assert_true(fabsf(d1 - d2) <= 0.000009);
+  }
+}

win/cglm.vcxproj

@@ -29,6 +29,7 @@
     <ClCompile Include="..\src\mat3.c" />
     <ClCompile Include="..\src\mat4.c" />
     <ClCompile Include="..\src\plane.c" />
+    <ClCompile Include="..\src\project.c" />
     <ClCompile Include="..\src\quat.c" />
     <ClCompile Include="..\src\vec3.c" />
     <ClCompile Include="..\src\vec4.c" />
@@ -47,6 +48,7 @@
     <ClInclude Include="..\include\cglm\call\mat3.h" />
     <ClInclude Include="..\include\cglm\call\mat4.h" />
     <ClInclude Include="..\include\cglm\call\plane.h" />
+    <ClInclude Include="..\include\cglm\call\project.h" />
     <ClInclude Include="..\include\cglm\call\quat.h" />
     <ClInclude Include="..\include\cglm\call\vec3.h" />
     <ClInclude Include="..\include\cglm\call\vec4.h" />
@@ -60,6 +62,7 @@
     <ClInclude Include="..\include\cglm\mat3.h" />
     <ClInclude Include="..\include\cglm\mat4.h" />
     <ClInclude Include="..\include\cglm\plane.h" />
+    <ClInclude Include="..\include\cglm\project.h" />
     <ClInclude Include="..\include\cglm\quat.h" />
     <ClInclude Include="..\include\cglm\simd\avx\affine.h" />
     <ClInclude Include="..\include\cglm\simd\avx\mat4.h" />

win/cglm.vcxproj.filters

@@ -75,6 +75,9 @@
     <ClCompile Include="..\src\box.c">
       <Filter>src</Filter>
     </ClCompile>
+    <ClCompile Include="..\src\project.c">
+      <Filter>src</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\src\config.h">
@@ -206,5 +209,11 @@
     <ClInclude Include="..\include\cglm\color.h">
       <Filter>include\cglm</Filter>
     </ClInclude>
+    <ClInclude Include="..\include\cglm\project.h">
+      <Filter>include\cglm</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\call\project.h">
+      <Filter>include\cglm\call</Filter>
+    </ClInclude>
   </ItemGroup>
 </Project>