quat: fix quaternion inverse and tests about it

* multiplication quaternion and its inverse must be identity
Merge pull request #43 from recp/quaternion
2026-02-17 03:39:05 +00:00 · 2018-04-11 16:50:37 +03:00 · 2018-04-11 12:43:48 +03:00 · 2018-04-11 12:36:39 +03:00 · 2018-04-11 12:34:20 +03:00 · 2018-04-11 12:31:29 +03:00
58 changed files with 3370 additions and 509 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -60,3 +60,4 @@ cglm_test_iosTests/*
 docs/build/*
 win/cglm_test_*
 * copy.*
+*.o
--- a/7
+++ b/7
@@ -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
--- a/README.md
+++ b/README.md
@@ -4,6 +4,8 @@
 [![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&amp;utm_medium=referral&amp;utm_content=recp/cglm&amp;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:
--- a/build-deps.sh
+++ b/build-deps.sh
@@ -16,14 +16,18 @@ cd $(dirname "$0")
 if [ "$(uname)" = "Darwin" ]; then
  libtoolBin=$(which glibtoolize)
  libtoolBinDir=$(dirname "${libtoolBin}")
+
+  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
--- a/configure.ac
+++ b/configure.ac
@@ -7,7 +7,7 @@
 #*****************************************************************************

 AC_PREREQ([2.69])
-AC_INIT([cglm], [0.3.5], [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])
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -41,6 +41,7 @@ Follow the :doc:`build` documentation for this
   vec4-ext
   color
   plane
+   project
   util
   io
   call
--- a/docs/source/box.rst
+++ b/docs/source/box.rst
@@ -24,6 +24,10 @@ Functions:
 #. :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
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -91,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
--- a/docs/source/euler.rst
+++ b/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
--- a/docs/source/io.rst
+++ b/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
--- a/docs/source/project.rst
+++ b/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
--- a/docs/source/quat.rst
+++ b/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
--- a/docs/source/util.rst
+++ b/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
--- a/docs/source/vec3-ext.rst
+++ b/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))
--- a/docs/source/vec3.rst
+++ b/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
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -158,6 +161,14 @@ Functions documentation
    Parameters:
      | *[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)

    make vector as inverse/opposite of itself
@@ -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
--- a/docs/source/vec4-ext.rst
+++ b/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))
--- a/docs/source/vec4.rst
+++ b/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
--- a/include/cglm/box.h
+++ b/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 */
--- a/include/cglm/call.h
+++ b/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
 }
--- a/include/cglm/call/euler.h
+++ b/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);
--- a/include/cglm/call/mat4.h
+++ b/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);
--- a/include/cglm/call/project.h
+++ b/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 */
+
+
--- a/include/cglm/call/quat.h
+++ b/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
 }
--- a/include/cglm/call/vec3.h
+++ b/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
--- a/include/cglm/call/vec4.h
+++ b/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
--- a/include/cglm/cglm.h
+++ b/include/cglm/cglm.h
@@ -23,5 +23,6 @@
 #include "color.h"
 #include "util.h"
 #include "io.h"
+#include "project.h"

 #endif /* cglm_h */
--- a/include/cglm/euler.h
+++ b/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,97 +70,61 @@ 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(vec3 angles, mat4 dest) {
+glm_euler_xyz(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]);

-  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;
-}
-
-/*!
- * @brief build rotation matrix from euler angles
- *
- * @param[in]  angles angles as vector [Ez, Ey, Ex]
- * @param[out] dest   rotation matrix
- */
-CGLM_INLINE
-void
-glm_euler_zyx(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]);
+  czsx = cz * sx;
+  cxcz = cx * cz;
+  sysz = sy * sz;

  dest[0][0] =  cy * cz;
-  dest[0][1] = cz * sx * sy + cx * sz;
-  dest[0][2] =-cx * cz * sy + sx * sz;
+  dest[0][1] =  czsx * sy + cx * sz;
+  dest[0][2] = -cxcz * 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[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;
@@ -167,26 +140,82 @@ glm_euler_zyx(vec3 angles,
 /*!
 * @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(vec3 angles, mat4 dest) {
+  glm_euler_xyz(angles, dest);
+}
+
+/*!
+ * @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_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]);

-  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;
+  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] = cy * cz * sx + sy * sz;
+  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 [Xangle, Yangle, Zangle]
+ * @param[out] dest   rotation matrix
+ */
+CGLM_INLINE
+void
+glm_euler_yxz(vec3 angles, mat4 dest) {
+  float cx, cy, cz,
+        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]);
+
+  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;
@@ -202,29 +231,33 @@ glm_euler_zxy(vec3 angles,
 /*!
 * @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]);

+  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] = sx * sy - cx * cy * sz;
+  dest[1][0] =  sxsy - cxcy * sz;
  dest[1][1] =  cx * cz;
-  dest[1][2] = cy * sx + cx * sy * sz;
-  dest[2][0] = cx * sy + cy * sx * sz;
+  dest[1][2] =  cysx + cxsy * sz;
+  dest[2][0] =  cxsy + cysx * sz;
  dest[2][1] = -cz * sx;
-  dest[2][2] = cx * cy - sx * sy * sz;
+  dest[2][2] =  cxcy - sxsy * sz;
  dest[0][3] =  0.0f;
  dest[1][3] =  0.0f;
  dest[2][3] =  0.0f;
@@ -237,63 +270,31 @@ glm_euler_xzy(vec3 angles,
 /*!
 * @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]);

-  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;

-/*!
- * @brief build rotation matrix from euler angles
- *
- * @param[in]  angles angles as vector [Ey, Ex, Ez]
- * @param[out] dest   rotation matrix
- */
-CGLM_INLINE
-void
-glm_euler_yxz(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 - sx * sy * sz;
-  dest[0][1] = cz * sx * sy + 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 + cy * sx * sz;
-  dest[2][1] =-cy * cz * sx + sy * sz;
+  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;
@@ -307,7 +308,45 @@ glm_euler_yxz(vec3 angles,
 /*!
 * @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[out] dest   rotation matrix
+ */
+CGLM_INLINE
+void
+glm_euler_zyx(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] =  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 [Xangle, Yangle, Zangle]
 * @param[in]  ord    euler order
 * @param[out] dest   rotation matrix
 */
@@ -332,40 +371,18 @@ 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[0][1] =  sx * sy + cx * cysz;
+      dest[0][2] = -cx * sy + cysx * sz;
+      dest[1][0] = -sz;
      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[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:
+    case GLM_EULER_XYZ:
      dest[0][0] =  cycz;
      dest[0][1] =  czsx * sy + cxsz;
      dest[0][2] = -cx * czsy + sx * sz;
@@ -377,6 +394,28 @@ glm_euler_by_order(vec3 angles, glm_euler_sq ord, mat4 dest) {
      dest[2][2] =  cxcy;
      break;
    case GLM_EULER_YXZ:
+      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] =  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;
@@ -387,16 +426,16 @@ glm_euler_by_order(vec3 angles, glm_euler_sq ord, mat4 dest) {
      dest[2][1] = -cycz * sx + sysz;
      dest[2][2] =  cxcy;
      break;
-    case GLM_EULER_YZX:
+    case GLM_EULER_ZYX:
      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][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;
  }

--- a/include/cglm/io.h
+++ b/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 */
--- a/include/cglm/mat3.h
+++ b/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
 *
--- a/include/cglm/mat4.h
+++ b/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

+#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 */
--- a/include/cglm/project.h
+++ b/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 */
--- a/include/cglm/quat.h
+++ b/include/cglm/quat.h
@@ -12,14 +12,40 @@

 Functions:
   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 v);
+   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_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 multiplies two quaternion and stores result in dest
+ * @brief conjugate of quaternion
 *
- * @param[in]   q1    quaternion 1
- * @param[in]   q2    quaternion 2
+ * @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_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_add(versor p, versor q, versor dest) {
+  glm_vec4_add(p, q, dest);
+}

-  glm_quat_normalize(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
+ *
+ * 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_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 */
--- a/include/cglm/simd/intrin.h
+++ b/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 */
--- a/include/cglm/simd/sse2/quat.h
+++ b/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)));
+
+  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)));
+
+  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)));
+
+  _mm_store_ps(dest, r);
 }

-  if (cosTheta >= 1.0f) {
-    _mm_store_ps(dest, xmm_q);
-    return;
-  }
-
-  sinTheta = sqrtf(1.0f - cosTheta * cosTheta);
-
-  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)));
-}

 #endif
 #endif /* cglm_quat_simd_h */
--- a/include/cglm/types.h
+++ b/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];
--- a/include/cglm/util.h
+++ b/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 */
--- a/include/cglm/vec3-ext.h
+++ b/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 */
--- a/include/cglm/vec3.h
+++ b/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
 *
--- a/include/cglm/vec4-ext.h
+++ b/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 */
--- a/include/cglm/vec4.h
+++ b/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 */
--- a/include/cglm/version.h
+++ b/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 5
+#define CGLM_VERSION_MINOR 4
+#define CGLM_VERSION_PATCH 0

 #endif /* cglm_version_h */
--- a/makefile.am
+++ b/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
--- a/src/euler.c
+++ b/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) {
--- a/src/mat4.c
+++ b/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) {
--- a/src/project.c
+++ b/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);
+}
--- a/src/quat.c
+++ b/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);
 }
--- a/src/vec3.c
+++ b/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);
+}
--- a/src/vec4.c
+++ b/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);
+}
--- a/test/src/test_clamp.c
+++ b/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);
+}
--- a/test/src/test_common.c
+++ b/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);
+}
+
--- a/test/src/test_common.h
+++ b/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 */
--- a/test/src/test_euler.c
+++ b/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);
+}
--- a/test/src/test_main.c
+++ b/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);
--- a/test/src/test_project.c
+++ b/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);
+}
--- a/test/src/test_quat.c
+++ b/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 */
+}
--- a/test/src/test_tests.h
+++ b/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 */
--- a/test/src/test_vec4.c
+++ b/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);
+  }
+}
--- a/win/cglm.vcxproj
+++ b/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" />
--- a/win/cglm.vcxproj.filters
+++ b/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>
Author	SHA1	Message	Date
Recep Aslantas	2d77123999	quat: fix quaternion inverse and tests about it * multiplication quaternion and its inverse must be identity	2018-04-11 16:50:37 +03:00
Recep Aslantas	462067cfdc	Merge pull request #43 from recp/quaternion quaternion improvements and new features	2018-04-11 12:43:48 +03:00
Recep Aslantas	9ae8da3e0a	update version to v0.4.0	2018-04-11 12:36:39 +03:00
Recep Aslantas	0e63c245d4	update docs	2018-04-11 12:34:20 +03:00
Recep Aslantas	de55850136	add call version of vector extensions	2018-04-11 12:31:29 +03:00
Recep Aslantas	51278b26b4	quat: update call versions of quaternion	2018-04-11 11:19:13 +03:00
Recep Aslantas	fdea13507b	replace mat4_mulq with glm_quat_rotate * glm_quat_rotate is better name to rotate transform matrix using quaternion. * we may use mat4_mulq in the future for another purpose e.g. left multiplication quat with matrix	2018-04-11 10:49:53 +03:00
Recep Aslantas	80d255e6d9	rotate vector using quaternion	2018-04-11 00:47:11 +03:00
Recep Aslantas	d447876c70	improve glm_vec_rotate	2018-04-11 00:46:23 +03:00
Recep Aslantas	b1fa7ff597	normalize axis quaternion axis-angle constructor	2018-04-11 00:36:39 +03:00
Recep Aslantas	010dcc9837	optimize normalize quaternion with SIMD * provide _to version for storing into another quat	2018-04-11 00:17:41 +03:00
Recep Aslantas	5dec68823c	add additional tests and comments to quat tests	2018-04-10 17:41:25 +03:00
Recep Aslantas	4c79fee5d3	quat: additional tests for angle, axis, mul (hamilton product)	2018-04-10 17:16:31 +03:00
Recep Aslantas	18ef0d7af1	quat: quaternion for look rotation ( from source point to dest point )	2018-04-10 16:52:52 +03:00
Recep Aslantas	9466182c10	quat: create view wmatrix with quaternion helper	2018-04-10 16:01:23 +03:00
Recep Aslantas	f0a51b35ad	quat: transposed/inverted version of quat2mat	2018-04-10 15:41:09 +03:00
Recep Aslantas	290bcf134c	quat: add lerp and improve slerp	2018-04-10 12:38:54 +03:00
Recep Aslantas	416e2f4452	vec: lerp for vec3 and vec4	2018-04-10 11:44:16 +03:00
Recep Aslantas	1fb82a1922	quat: use vector functions for available operations * provide quat_copy function	2018-04-10 10:47:55 +03:00
Recep Aslantas	591c881376	vec: extend flip sign to store result in another vector	2018-04-10 10:46:45 +03:00
Recep Aslantas	6f69da361b	quaternion multiplication * convert quaternion multiplication to xyzw * previous implementation may be wrong, wikipedia version implemented * implement SSE version	2018-04-09 23:56:09 +03:00
Recep Aslantas	93a08fce17	quat: axis angle of quaternion	2018-04-09 23:12:44 +03:00
Recep Aslantas	cc1d3b53ea	quat: implement add, sub, real and imag helpers	2018-04-09 22:32:55 +03:00
Recep Aslantas	b21df8fc37	inverse of quaternion	2018-04-09 22:26:23 +03:00
Recep Aslantas	76e9f74020	conjugate of quaternion	2018-04-09 21:54:53 +03:00
Recep Aslantas	d79e58486d	update credits file	2018-04-09 21:54:35 +03:00
Recep Aslantas	3dc93c56e8	convert quaterinon to xyzw order (part 1)	2018-04-09 18:49:12 +03:00
Recep Aslantas	7615f785ac	improve quaternion to matrix	2018-04-09 00:53:14 +03:00
Recep Aslantas	f0daaca58b	improve matrix to quaternion	2018-04-09 00:46:00 +03:00
Recep Aslantas	381b2fdcc0	fix vec4_norm2, use dot for vec3_norm2	2018-04-09 00:01:56 +03:00
Recep Aslantas	e4e0fa623c	sse2 version of vec4 dot product * use this for normalizing vector	2018-04-08 18:27:54 +03:00
Recep Aslantas	932f638d5a	optimize mat4 to quaternion * add SSE2 version and optimize scalar version	2018-04-08 12:31:32 +03:00
Recep Aslantas	81bda7439d	vector square root	2018-04-08 12:30:15 +03:00
Recep Aslantas	b27603c268	normalize quaternion before converting to matrix * because it must be unit quaternion and didn't specified this in docs. * we must provide alternative func for unit quat	2018-04-08 00:09:40 +03:00
Recep Aslantas	12c5307447	vec3 and vec4 sign helper	2018-04-07 21:53:22 +03:00
Recep Aslantas	257c57d41f	mat4 to quaternion	2018-04-07 19:46:46 +03:00
Recep Aslantas	f5140ea005	quat: mat4_mul_quat helper * the quaternion is used as right matrix	2018-04-07 13:47:20 +03:00
Recep Aslantas	619ecdc5a4	quat: improve normalize	2018-04-07 13:46:46 +03:00
Recep Aslantas	9b8748acc4	quat: quaternion to mat3	2018-04-07 13:27:40 +03:00
Recep Aslantas	ae06c51746	improve glm_mat4_mulN for non-DEBUG environment	2018-04-07 13:22:44 +03:00
Recep Aslantas	11430559b4	fix isnan and isinf	2018-04-07 08:28:37 +03:00
Recep Aslantas	58f0043417	vector utils: isnan and isinf * a vector which has least one NaN or INF member, is assumed not valid vector.	2018-04-06 22:57:24 +03:00
Recep Aslantas	967fb1afad	Update README.md	2018-04-03 17:32:10 +03:00
Recep Aslantas	7411ac36c1	update docs for euler angles	2018-04-03 17:05:45 +03:00
Recep Aslantas	238609f2c0	Merge pull request #31 from recp/proj add project / unproject functions	2018-04-03 16:51:55 +03:00
Recep Aslantas	ea0a10ade9	suppress warnings	2018-04-03 16:47:59 +03:00
Recep Aslantas	429fdfd5c5	update build scripts	2018-04-03 16:47:51 +03:00
Recep Aslantas	024412f00e	add docs for project/unproject	2018-04-03 16:41:13 +03:00
Recep Aslantas	e8615ea14c	fix tests list	2018-04-03 12:35:30 +03:00
Recep Aslantas	be81d73895	Update test_main.c	2018-04-03 12:32:21 +03:00
Recep Aslantas	b16f0ded85	Merge branch 'master' into proj	2018-04-03 12:30:03 +03:00
Recep Aslantas	63acfd681e	fix unproject, add tests to project/unproject	2018-04-03 12:27:20 +03:00
Recep Aslantas	eb527e39b4	optimize project	2018-04-03 11:25:33 +03:00
Recep Aslantas	9f389ab8ec	project function	2018-04-03 11:09:13 +03:00
Recep Aslantas	3399595dc2	add vec2 type	2018-04-03 10:46:46 +03:00
Recep Aslantas	2513d46102	Merge pull request #41 from winduptoy/patch-1 Fix small typo.	2018-04-02 20:26:52 +03:00
Matt Reyer	c298f4a4d7	Fix small typo.	2018-04-02 11:33:58 -04:00
Recep Aslantas	84cdbd5072	Merge pull request #40 from recp/aabb-ext Axis-Aligned Bounding Box (AABB) Extensions	2018-04-02 16:40:23 +03:00
Recep Aslantas	74f9865884	add docs for new aabb functions	2018-04-02 16:36:55 +03:00
Recep Aslantas	dbd1e334ea	aabb box size and radius	2018-04-02 16:26:14 +03:00
Recep Aslantas	acda316c12	get sign of float helper as -1, +1 and 0 * add clarification for zero input	2018-04-02 16:18:50 +03:00
Recep Aslantas	86efe64b8e	helper for check aabb is valid or not	2018-04-02 12:35:22 +03:00
Recep Aslantas	b0991342a6	aabb printer function	2018-04-02 12:08:08 +03:00
Recep Aslantas	984916d520	invalidate axis-aligned boundng box util	2018-04-02 11:50:53 +03:00
Recep Aslantas	54c44ff224	Merge pull request #39 from opencollective/opencollective Activating Open Collective	2018-04-01 22:26:33 +03:00
Jess	db4761b437	Added backers and sponsors on the README	2018-04-01 18:12:46 +09:00
Recep Aslantas	ca504f7058	now working on v0.3.6	2018-03-29 00:12:16 +03:00
Recep Aslantas	5a7b9caf16	Update README.md	2018-03-29 00:02:37 +03:00
Recep Aslantas	43b3df992d	Merge pull request #37 from recp/euler fix euler angles (extrinsic -> intrinsic)	2018-03-28 23:58:15 +03:00
Recep Aslantas	26110f83d1	euler: fix thetaY in extracting angles	2018-03-27 12:35:19 +03:00
Recep Aslantas	d1f3feeb6e	test: add tests for euler XYZ	2018-03-27 12:14:46 +03:00
Recep Aslantas	4298211795	euler: fix extracting XYZ angles	2018-03-27 12:14:12 +03:00
Recep Aslantas	c244b68e73	build: improve build-deps	2018-03-27 11:22:05 +03:00
Recep Aslantas	205d13aa93	fix euler angles (extrinsic -> intrinsic) because cglm uses intrinsics for these rotations	2018-03-27 11:14:26 +03:00
Recep Aslantas	45f13217c3	Merge pull request #35 from recp/clamp clamp functions	2018-03-22 21:28:19 +03:00
Recep Aslantas	21ec45b2af	add tests for clamp	2018-03-22 21:24:41 +03:00
Recep Aslantas	71b48b530e	add documentation for clamp	2018-03-22 21:24:26 +03:00
Recep Aslantas	48b7b30e42	add call version for clamp	2018-03-22 21:18:08 +03:00
Recep Aslantas	86055097e1	clamp functions	2018-03-22 18:10:10 +03:00
Recep Aslantas	08be94a89b	Merge pull request #34 from NoxNode/typofix typo fixes	2018-03-20 10:41:11 +03:00
mcsquizzy123	91b2a989e2	typo fixes - heaer and haeder	2018-03-19 18:37:49 -07:00
Recep Aslantas	780179ff0d	fix unproject	2018-03-08 22:29:10 +03:00
Recep Aslantas	c148eacdc2	fix unproject's parameters	2018-03-08 13:12:08 +03:00
Recep Aslantas	29996d0bdd	add unproject function	2018-03-08 13:02:33 +03:00