diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,14 @@
 # h-raylib changelog
 
+## Version 4.6.0.1
+_2 April, 2023_
+
+- Created the `Raylib.Util.Math` and `Raylib.Util.Camera` modules. They are Haskell implementations of `raymath` and `rcamera`.
+
+\[[#15](https://github.com/Anut-py/h-raylib/pull/15)\]
+
+- Fixed a memory issue with `getFontDefault`
+
 ## Version 4.5.3.4
 _19 March, 2023_
 
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -12,7 +12,7 @@
 
 _This section only contains h-raylib specific information. For information about raylib in general, view the [raylib wiki](https://github.com/raysan5/raylib/wiki)._
 
-This project is split into 10 public modules. `Raylib.Types` contains all of raylib's types and low-level code to convert them to and from raw bytes. `Raylib.Util` contains miscellaneous utility functions. `Raylib.Util.Colors` contains some colors defined by raylib. The other 8 public modules, `Raylib.Core`, `Raylib.Core.Camera`, `Raylib.Core.Shapes`, `Raylib.Core.Textures`, `Raylib.Core.Text`, `Raylib.Core.Models`, `Raylib.Core.Audio`, and `Raylib.Util.RLGL`, correspond to their respective raylib modules.
+This project is split into 12 public modules. `Raylib.Types` contains all of raylib's types and low-level code to convert them to and from raw bytes. `Raylib.Util` contains miscellaneous utility functions. `Raylib.Util.Colors` contains some colors defined by raylib. The other 10 public modules, `Raylib.Core`, `Raylib.Core.Camera`, `Raylib.Core.Shapes`, `Raylib.Core.Textures`, `Raylib.Core.Text`, `Raylib.Core.Models`, `Raylib.Core.Audio`, `Raylib.Util.Camera`, `Raylib.Util.Math`, and `Raylib.Util.RLGL`, correspond to their respective raylib modules.
 
 The functions in h-raylib are an almost one-to-one mapping to their corresponding raylib functions. The types are, in some cases, slightly modified if it is possible to utilize Haskell features.
 
@@ -20,12 +20,10 @@
 
 ### Raylib.Types
 
-`Raylib.Types` has 4 sections: one for enumerations, one for typeclasses, one for structures, and one for callbacks.
+`Raylib.Types` has 3 sections: one for enumerations, one for structures, and one for callbacks.
 
 The enumerations section contains Haskell sum types that are instances of `Enum`. Each of these types corresponds to a raylib `enum` or set of `define`s. The `fromEnum` and `toEnum` functions for these types use the numbers associated with these values in the C `enum`s. Most of these types are instances of `Storable` so they can be converted to raw bytes and passed to a C function. _NOTE: Some of these Haskell types correspond to C `enum`s that are in C source files, rather than header files._
 
-The typeclasses section contains Haskell typeclasses that are derived by some of the types in the structures section.
-
 The structures section contains Haskell types that correspond to each of raylib's `structs`. Each field in these types is named `typeName'fieldName` (e.g. the C struct `Vector2`'s `x` field is called `vector2'x` in Haskell). These structs also all derive the typeclass `Freeable` (declared in the internal `Raylib.ForeignUtil` module). This typeclass allows types to describe how to properly free all the data associated with a pointer to that type. For example, `Image`'s implementation of `Freeable` also frees the pointer stored in the `Image.data` field in C. Finally, all of these types derive `Storable`, obviously, to convert them to and from pointers.
 
 The callbacks section contains `FunPtr` types that are passed to some functions. _NOTE: These callbacks are very unlikely to be used, so they may be removed in the future._
@@ -38,7 +36,7 @@
 
 `Raylib.Util.Colors` is very simple: it declares 26 colors defined in `raylib.h`, namely `lightGray`, `gray`, `darkGray`, `yellow`, `gold`, `orange`, `pink`, `red`, `maroon`, `green`, `lime`, `darkGreen`, `skyBlue`, `blue`, `darkBlue`, `purple`, `violet`, `darkPurple`, `beige`, `brown`, `darkBrown`, `white`, `black`, `blank`, `magenta`, and `rayWhite`.
 
-### The other 8 modules
+### The other 10 modules
 
 These modules contain only functions. Each of these functions corresponds to a C function.
 
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -78,6 +78,11 @@
 
 To contribute or learn more about the project, see [CONTRIBUTING.md](https://github.com/Anut-py/h-raylib/blob/master/CONTRIBUTING.md).
 
+## FAQ
+
+- When I try to run an h-raylib program I get an error saying "*The code execution cannot proceed because libwinpthread-1.dll was not found. Reinstalling the program may fix this problem.*"
+  - See [#14](https://github.com/Anut-py/h-raylib/issues/14)
+
 ## License
 
 This project is licensed under the Apache License 2.0. See more in `LICENSE`.
diff --git a/ROADMAP.md b/ROADMAP.md
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -3,12 +3,12 @@
 Items higher on the list have higher priority.
 
 ## Uncompleted
-- Bind `rcamera`
-- Bind `raymath`
 - Bind `raygui`
 - Add web build support
 
 ## Completed
 - Make it easier to pass shader parameters (Added in `4.5.3.0`)
 - Allow manual unloading of assets for larger projects (Added in `4.5.3.1`)
-- Bind `rlgl` (Not published to Hackage)
+- Bind `rlgl` (Added in `4.5.3.2`)
+- Bind `rcamera` (Not published to hackage)
+- Bind `raymath` (Not published to hackage)
diff --git a/h-raylib.cabal b/h-raylib.cabal
--- a/h-raylib.cabal
+++ b/h-raylib.cabal
@@ -1,6 +1,6 @@
 cabal-version:      2.4
 name:               h-raylib
-version:            4.5.3.4
+version:            4.6.0.1
 synopsis:           Raylib bindings for Haskell
 category:           graphics
 description:
@@ -153,7 +153,9 @@
     Raylib.Core.Textures
     Raylib.Types
     Raylib.Util
+    Raylib.Util.Camera
     Raylib.Util.Colors
+    Raylib.Util.Math
     Raylib.Util.RLGL
 
   other-modules:
diff --git a/lib/rl_bindings.c b/lib/rl_bindings.c
--- a/lib/rl_bindings.c
+++ b/lib/rl_bindings.c
@@ -976,9 +976,24 @@
 
 Font *GetFontDefault_()
 {
-    Font *ptr = (Font *)malloc(sizeof(Font));
-    *ptr = GetFontDefault();
-    return ptr;
+    Font defaultFont      = GetFontDefault();
+    Font *defaultFontCopy = (Font *)malloc(sizeof(Font));
+
+    defaultFontCopy->baseSize     = defaultFont.baseSize;
+    defaultFontCopy->glyphCount   = defaultFont.glyphCount;
+    defaultFontCopy->glyphPadding = defaultFont.glyphPadding; 
+    defaultFontCopy->texture      = defaultFont.texture;
+    
+    defaultFontCopy->glyphs = malloc(sizeof(GlyphInfo) * defaultFont.glyphCount);
+    defaultFontCopy->recs = malloc(sizeof(Rectangle) * defaultFont.glyphCount);
+
+    for (int i = 0; i < defaultFont.glyphCount; i++) {
+        defaultFontCopy->glyphs[i] = defaultFont.glyphs[i];
+        defaultFontCopy->glyphs[i].image = ImageCopy(defaultFont.glyphs[i].image);
+        defaultFontCopy->recs[i]   = defaultFont.recs[i];
+    }
+    
+    return defaultFontCopy;
 }
 
 Font *LoadFont_(char *a)
diff --git a/raylib/examples/core/core_2d_camera_platformer.c b/raylib/examples/core/core_2d_camera_platformer.c
--- a/raylib/examples/core/core_2d_camera_platformer.c
+++ b/raylib/examples/core/core_2d_camera_platformer.c
@@ -90,7 +90,7 @@
         "Follow player center",
         "Follow player center, but clamp to map edges",
         "Follow player center; smoothed",
-        "Follow player center horizontally; updateplayer center vertically after landing",
+        "Follow player center horizontally; update player center vertically after landing",
         "Player push camera on getting too close to screen edge"
     };
 
diff --git a/raylib/src/external/qoaplay.c b/raylib/src/external/qoaplay.c
--- a/raylib/src/external/qoaplay.c
+++ b/raylib/src/external/qoaplay.c
@@ -62,7 +62,7 @@
 extern "C" {            // Prevents name mangling of functions
 #endif
 
-qoaplay_desc *qoaplay_open(char *path);
+qoaplay_desc *qoaplay_open(const char *path);
 qoaplay_desc *qoaplay_open_memory(const unsigned char *data, int data_size);
 void qoaplay_close(qoaplay_desc *qoa_ctx);
 
@@ -83,7 +83,7 @@
 //----------------------------------------------------------------------------------
 
 // Open QOA file, keep FILE pointer to keep reading from file
-qoaplay_desc *qoaplay_open(char *path)
+qoaplay_desc *qoaplay_open(const char *path)
 {
     FILE *file = fopen(path, "rb");
     if (!file) return NULL;
diff --git a/raylib/src/external/rl_gputex.h b/raylib/src/external/rl_gputex.h
--- a/raylib/src/external/rl_gputex.h
+++ b/raylib/src/external/rl_gputex.h
@@ -814,5 +814,4 @@
 
     return data_size;
 }
-
 #endif // RL_GPUTEX_IMPLEMENTATION
diff --git a/raylib/src/raudio.c b/raylib/src/raudio.c
--- a/raylib/src/raudio.c
+++ b/raylib/src/raudio.c
@@ -1287,7 +1287,7 @@
 #if defined(SUPPORT_FILEFORMAT_QOA)
     else if (IsFileExtension(fileName, ".qoa"))
     {
-        qoaplay_desc *ctxQoa = qoaplay_open((char *)fileName);
+        qoaplay_desc *ctxQoa = qoaplay_open(fileName);
         music.ctxType = MUSIC_AUDIO_QOA;
         music.ctxData = ctxQoa;
 
diff --git a/raylib/src/raylib.h b/raylib/src/raylib.h
--- a/raylib/src/raylib.h
+++ b/raylib/src/raylib.h
@@ -1,6 +1,6 @@
 /**********************************************************************************************
 *
-*   raylib v4.5 - A simple and easy-to-use library to enjoy videogames programming (www.raylib.com)
+*   raylib v4.6-dev - A simple and easy-to-use library to enjoy videogames programming (www.raylib.com)
 *
 *   FEATURES:
 *       - NO external dependencies, all required libraries included with raylib
@@ -82,9 +82,9 @@
 #include <stdarg.h>     // Required for: va_list - Only used by TraceLogCallback
 
 #define RAYLIB_VERSION_MAJOR 4
-#define RAYLIB_VERSION_MINOR 5
+#define RAYLIB_VERSION_MINOR 6
 #define RAYLIB_VERSION_PATCH 0
-#define RAYLIB_VERSION  "4.5"
+#define RAYLIB_VERSION  "4.6-dev"
 
 // Function specifiers in case library is build/used as a shared library (Windows)
 // NOTE: Microsoft specifiers to tell compiler that symbols are imported/exported from a .dll
diff --git a/raylib/src/raymath.h b/raylib/src/raymath.h
--- a/raylib/src/raymath.h
+++ b/raylib/src/raymath.h
@@ -703,12 +703,14 @@
     Vector3 result = v;
 
     float length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
-    if (length == 0.0f) length = 1.0f;
-    float ilength = 1.0f/length;
+    if (length != 0.0f) 
+    {
+        float ilength = 1.0f/length;
 
-    result.x *= ilength;
-    result.y *= ilength;
-    result.z *= ilength;
+        result.x *= ilength;
+        result.y *= ilength;
+        result.z *= ilength;
+    }
 
     return result;
 }
diff --git a/raylib/src/rcore.c b/raylib/src/rcore.c
--- a/raylib/src/rcore.c
+++ b/raylib/src/rcore.c
@@ -595,14 +595,14 @@
     "ACTION_SETTARGETFPS"
 };
 
-// Automation Event (20 bytes)
+// Automation Event (24 bytes)
 typedef struct AutomationEvent {
     unsigned int frame;                 // Event frame
-    unsigned int type;                  // Event type (AutoEventType)
-    int params[3];                      // Event parameters (if required)
+    unsigned int type;                  // Event type (AutomationEventType)
+    int params[4];                      // Event parameters (if required)
 } AutomationEvent;
 
-static AutomationEvent *events = NULL;        // Events array
+static AutomationEvent *events = NULL;  // Events array
 static unsigned int eventCount = 0;     // Events count
 static bool eventsPlaying = false;      // Play events
 static bool eventsRecording = false;    // Record events
@@ -925,7 +925,7 @@
 #endif
 
 #if defined(SUPPORT_EVENTS_AUTOMATION)
-    events = (AutomationEvent *)malloc(MAX_CODE_AUTOMATION_EVENTS*sizeof(AutomationEvent));
+    events = (AutomationEvent *)RL_CALLOC(MAX_CODE_AUTOMATION_EVENTS, sizeof(AutomationEvent));
     CORE.Time.frameCounter = 0;
 #endif
 
@@ -1073,7 +1073,7 @@
 #endif
 
 #if defined(SUPPORT_EVENTS_AUTOMATION)
-    free(events);
+    RL_FREE(events);
 #endif
 
     CORE.Window.ready = false;
@@ -6945,11 +6945,11 @@
 // TODO: This system should probably be redesigned
 static void LoadAutomationEvents(const char *fileName)
 {
-    //unsigned char fileId[4] = { 0 };
-
-    // Load binary
+    // Load events file (binary)
     /*
     FILE *repFile = fopen(fileName, "rb");
+    unsigned char fileId[4] = { 0 };
+    
     fread(fileId, 1, 4, repFile);
 
     if ((fileId[0] == 'r') && (fileId[1] == 'E') && (fileId[2] == 'P') && (fileId[1] == ' '))
@@ -6962,7 +6962,7 @@
     fclose(repFile);
     */
 
-    // Load events (text file)
+    // Load events file (text)
     FILE *repFile = fopen(fileName, "rt");
 
     if (repFile != NULL)
diff --git a/raylib/src/rmodels.c b/raylib/src/rmodels.c
--- a/raylib/src/rmodels.c
+++ b/raylib/src/rmodels.c
@@ -3400,7 +3400,7 @@
     {
         if (mesh->vboId[SHADER_LOC_VERTEX_TANGENT] != 0)
         {
-            // Upate existing vertex buffer
+            // Update existing vertex buffer
             rlUpdateVertexBuffer(mesh->vboId[SHADER_LOC_VERTEX_TANGENT], mesh->tangents, mesh->vertexCount*4*sizeof(float), 0);
         }
         else
diff --git a/raylib/src/rtext.c b/raylib/src/rtext.c
--- a/raylib/src/rtext.c
+++ b/raylib/src/rtext.c
@@ -694,14 +694,25 @@
     // NOTE: Rectangles memory is loaded here!
     Rectangle *recs = (Rectangle *)RL_MALLOC(glyphCount*sizeof(Rectangle));
 
-    // Calculate image size based on required pixel area
-    // NOTE 1: Image is forced to be squared and POT... very conservative!
-    // NOTE 2: SDF font characters already contain an internal padding,
-    // so image size would result bigger than default font type
-    float requiredArea = 0;
-    for (int i = 0; i < glyphCount; i++) requiredArea += ((chars[i].image.width + 2*padding)*(fontSize + 2*padding));
-    float guessSize = sqrtf(requiredArea)*1.4f;
-    int imageSize = (int)powf(2, ceilf(logf((float)guessSize)/logf(2)));    // Calculate next POT
+    // Calculate image size based on total glyph width and glyph row count
+    int totalWidth = 0;
+    int maxGlyphWidth = 0;
+    
+    for (int i = 0; i < glyphCount; i++)
+    {
+        if (chars[i].image.width > maxGlyphWidth) maxGlyphWidth = chars[i].image.width;
+        totalWidth += chars[i].image.width + 2*padding;
+    }
+    
+    int rowCount = 0;
+    int imageSize = 64;  // Define minimum starting value to avoid unnecessary calculation steps for very small images
+    
+    // NOTE: maxGlyphWidth is maximum possible space left at the end of row
+    while (totalWidth > (imageSize - maxGlyphWidth)*rowCount) 
+    {
+        imageSize *= 2;                                 // Double the size of image (to keep POT)
+        rowCount = imageSize/(fontSize + 2*padding);    // Calculate new row count for the new image size
+    }
 
     atlas.width = imageSize;   // Atlas bitmap width
     atlas.height = imageSize;  // Atlas bitmap height
@@ -720,24 +731,7 @@
         // NOTE: Using simple packaging, one char after another
         for (int i = 0; i < glyphCount; i++)
         {
-            // Copy pixel data from fc.data to atlas
-            for (int y = 0; y < chars[i].image.height; y++)
-            {
-                for (int x = 0; x < chars[i].image.width; x++)
-                {
-                    ((unsigned char *)atlas.data)[(offsetY + y)*atlas.width + (offsetX + x)] = ((unsigned char *)chars[i].image.data)[y*chars[i].image.width + x];
-                }
-            }
-
-            // Fill chars rectangles in atlas info
-            recs[i].x = (float)offsetX;
-            recs[i].y = (float)offsetY;
-            recs[i].width = (float)chars[i].image.width;
-            recs[i].height = (float)chars[i].image.height;
-
-            // Move atlas position X for next character drawing
-            offsetX += (chars[i].image.width + 2*padding);
-
+            // Check remaining space for glyph
             if (offsetX >= (atlas.width - chars[i].image.width - 2*padding))
             {
                 offsetX = padding;
@@ -761,6 +755,24 @@
                     break;
                 }
             }
+
+            // Copy pixel data from fc.data to atlas
+            for (int y = 0; y < chars[i].image.height; y++)
+            {
+                for (int x = 0; x < chars[i].image.width; x++)
+                {
+                    ((unsigned char *)atlas.data)[(offsetY + y)*atlas.width + (offsetX + x)] = ((unsigned char *)chars[i].image.data)[y*chars[i].image.width + x];
+                }
+            }
+
+            // Fill chars rectangles in atlas info
+            recs[i].x = (float)offsetX;
+            recs[i].y = (float)offsetY;
+            recs[i].width = (float)chars[i].image.width;
+            recs[i].height = (float)chars[i].image.height;
+
+            // Move atlas position X for next character drawing
+            offsetX += (chars[i].image.width + 2*padding);
         }
     }
     else if (packMethod == 1)  // Use Skyline rect packing algorithm (stb_pack_rect)
diff --git a/raylib/src/rtextures.c b/raylib/src/rtextures.c
--- a/raylib/src/rtextures.c
+++ b/raylib/src/rtextures.c
@@ -219,7 +219,7 @@
 //----------------------------------------------------------------------------------
 // Other Modules Functions Declaration (required by text)
 //----------------------------------------------------------------------------------
-// ...
+extern void LoadFontDefault(void);          // [Module: text] Loads default font, required by ImageDrawText()
 
 //----------------------------------------------------------------------------------
 // Module specific Functions Declaration
@@ -3152,6 +3152,9 @@
 void ImageDrawText(Image *dst, const char *text, int posX, int posY, int fontSize, Color color)
 {
 #if defined(SUPPORT_MODULE_RTEXT)
+    // Make sure default font is loaded to be used on image text drawing
+    if (GetFontDefault().texture.id == 0) LoadFontDefault();
+
     Vector2 position = { (float)posX, (float)posY };
     // NOTE: For default font, spacing is set to desired font size / default font size (10)
     ImageDrawTextEx(dst, GetFontDefault(), text, position, (float)fontSize, (float)fontSize/10, color);   // WARNING: Module required: rtext
diff --git a/src/Raylib/Types.hs b/src/Raylib/Types.hs
--- a/src/Raylib/Types.hs
+++ b/src/Raylib/Types.hs
@@ -577,24 +577,24 @@
 
 -- I don't know if there's a cleaner way to do this
 unpackShaderUniformData :: ShaderUniformData -> IO (ShaderUniformDataType, Ptr ())
-unpackShaderUniformData x = do
-  case x of
+unpackShaderUniformData u = do
+  case u of
     (ShaderUniformFloat f) ->
       do
         ptr <- malloc
         poke ptr (realToFrac f :: CFloat)
         return (ShaderUniformFloatType, castPtr ptr)
-    (ShaderUniformVec2 v) ->
+    (ShaderUniformVec2 (Vector2 x y)) ->
       do
-        ptr <- newArray (map realToFrac (asList v) :: [CFloat])
+        ptr <- newArray (map realToFrac [x, y] :: [CFloat])
         return (ShaderUniformVec2Type, castPtr ptr)
-    (ShaderUniformVec3 v) ->
+    (ShaderUniformVec3 (Vector3 x y z)) ->
       do
-        ptr <- newArray (map realToFrac (asList v) :: [CFloat])
+        ptr <- newArray (map realToFrac [x, y, z] :: [CFloat])
         return (ShaderUniformVec3Type, castPtr ptr)
-    (ShaderUniformVec4 v) ->
+    (ShaderUniformVec4 (Vector4 x y z w)) ->
       do
-        ptr <- newArray (map realToFrac (asList v) :: [CFloat])
+        ptr <- newArray (map realToFrac [x, y, z, w] :: [CFloat])
         return (ShaderUniformVec4Type, castPtr ptr)
     (ShaderUniformInt i) ->
       do
@@ -628,15 +628,15 @@
         return (ShaderUniformFloatType, castPtr ptr, length fs)
     (ShaderUniformVec2V vs) ->
       do
-        ptr <- newArray (map realToFrac $ concatMap asList vs :: [CFloat])
+        ptr <- newArray (map realToFrac $ concatMap (\(Vector2 x y) -> [x, y]) vs :: [CFloat])
         return (ShaderUniformVec2Type, castPtr ptr, length vs)
     (ShaderUniformVec3V vs) ->
       do
-        ptr <- newArray (map realToFrac $ concatMap asList vs :: [CFloat])
+        ptr <- newArray (map realToFrac $ concatMap (\(Vector3 x y z) -> [x, y, z]) vs :: [CFloat])
         return (ShaderUniformVec3Type, castPtr ptr, length vs)
     (ShaderUniformVec4V vs) ->
       do
-        ptr <- newArray (map realToFrac $ concatMap asList vs :: [CFloat])
+        ptr <- newArray (map realToFrac $ concatMap (\(Vector4 x y z w) -> [x, y, z, w]) vs :: [CFloat])
         return (ShaderUniformVec4Type, castPtr ptr, length vs)
     (ShaderUniformIntV is) ->
       do
@@ -1551,45 +1551,6 @@
   poke ptr v = poke (castPtr ptr) (fromIntegral (fromEnum v) :: CInt)
 
 ------------------------------------------------
--- Raylib typeclasses --------------------------
-------------------------------------------------
-
-class Vector a where
-  -- List representation of the vector
-  asList :: a -> [Float]
-
-  -- Vector-vector addition
-  (|+|) :: a -> a -> a
-
-  -- Vector-vector subtraction
-  (|-|) :: a -> a -> a
-  v1 |-| v2 = v1 |+| inverse v2
-
-  -- Vector-scalar multiplication
-  (|*|) :: a -> Float -> a
-
-  -- Vector-scalar division
-  (|/|) :: a -> Float -> a
-  v |/| num = v |*| (1 / num)
-
-  -- Vector-vector dot product
-  (|.|) :: a -> a -> Float
-
-  -- Zero vector
-  zero :: a
-
-  -- Vector additive inverse
-  inverse :: a -> a
-
-  -- Normalize vector (same direction, magnitude 1)
-  normalize :: a -> a
-  normalize v = v |/| magnitude v
-
-  -- Vector magnitude
-  magnitude :: a -> Float
-  magnitude x = sqrt $ x |.| x
-
-------------------------------------------------
 -- Raylib structures ---------------------------
 ------------------------------------------------
 
@@ -1613,17 +1574,6 @@
     pokeByteOff _p 4 (realToFrac y :: CFloat)
     return ()
 
-instance Vector Vector2 where
-  asList (Vector2 x y) = [x, y]
-
-  (Vector2 x1 y1) |+| (Vector2 x2 y2) = Vector2 (x1 + x2) (y1 + y2)
-  (Vector2 x y) |*| num = Vector2 (x * num) (y * num)
-
-  (Vector2 x1 y1) |.| (Vector2 x2 y2) = (x1 * x2) + (y1 * y2)
-
-  zero = Vector2 0 0
-  inverse (Vector2 x y) = Vector2 (- x) (- y)
-
 data Vector3 = Vector3
   { vector3'x :: Float,
     vector3'y :: Float,
@@ -1645,21 +1595,6 @@
     pokeByteOff _p 8 (realToFrac z :: CFloat)
     return ()
 
--- Vector cross-product
-cross :: Vector3 -> Vector3 -> Vector3
-(Vector3 x1 y1 z1) `cross` (Vector3 x2 y2 z2) = Vector3 (y1 * z2 - z1 * y2) (z1 * x2 - x1 * z2) (x1 * y2 - y1 * x2)
-
-instance Vector Vector3 where
-  asList (Vector3 x y z) = [x, y, z]
-
-  (Vector3 x1 y1 z1) |+| (Vector3 x2 y2 z2) = Vector3 (x1 + x2) (y1 + y2) (z1 + z2)
-  (Vector3 x y z) |*| num = Vector3 (x * num) (y * num) (z * num)
-
-  (Vector3 x1 y1 z1) |.| (Vector3 x2 y2 z2) = (x1 * x2) + (y1 * y2) + (z1 * z2)
-
-  zero = Vector3 0 0 0
-  inverse (Vector3 x y z) = Vector3 (- x) (- y) (- z)
-
 data Vector4 = Vector4
   { vector4'x :: Float,
     vector4'y :: Float,
@@ -1683,17 +1618,6 @@
     pokeByteOff _p 8 (realToFrac z :: CFloat)
     pokeByteOff _p 12 (realToFrac w :: CFloat)
     return ()
-
-instance Vector Vector4 where
-  asList (Vector4 x y z w) = [x, y, z, w]
-
-  (Vector4 x1 y1 z1 w1) |+| (Vector4 x2 y2 z2 w2) = Vector4 (x1 + x2) (y1 + y2) (z1 + z2) (w1 + w2)
-  (Vector4 x y z w) |*| num = Vector4 (x * num) (y * num) (z * num) (w * num)
-
-  (Vector4 x1 y1 z1 w1) |.| (Vector4 x2 y2 z2 w2) = (x1 * x2) + (y1 * y2) + (z1 * z2) + (w1 * w2)
-
-  zero = Vector4 0 0 0 0
-  inverse (Vector4 x y z w) = Vector4 (- x) (- y) (- z) (- w)
 
 type Quaternion = Vector4
 
diff --git a/src/Raylib/Util.hs b/src/Raylib/Util.hs
--- a/src/Raylib/Util.hs
+++ b/src/Raylib/Util.hs
@@ -1,10 +1,11 @@
 {-# OPTIONS -Wall #-}
 {-# LANGUAGE CPP #-}
 
-module Raylib.Util (WindowResources, cameraDirectionRay, whileWindowOpen, whileWindowOpen_, whileWindowOpen0, setMaterialShader, inGHCi) where
+module Raylib.Util (WindowResources, cameraDirectionRay, whileWindowOpen, whileWindowOpen_, whileWindowOpen0, setMaterialShader, inGHCi, Freeable (..)) where
 
 import Control.Monad (void)
 import Raylib.Core (windowShouldClose)
+import Raylib.ForeignUtil (Freeable (..))
 import Raylib.Internal (WindowResources)
 import Raylib.Types
   ( Camera3D (camera3D'position, camera3D'target),
@@ -12,12 +13,12 @@
     Model (model'materials),
     Ray (Ray),
     Shader,
-    Vector (normalize, (|-|)),
   )
+import Raylib.Util.Math (Vector (vectorNormalize, (|-|)))
 
 -- | Gets the direction of a camera as a ray.
 cameraDirectionRay :: Camera3D -> Ray
-cameraDirectionRay camera = Ray (camera3D'position camera) (normalize $ camera3D'target camera |-| camera3D'position camera)
+cameraDirectionRay camera = Ray (camera3D'position camera) (vectorNormalize $ camera3D'target camera |-| camera3D'position camera)
 
 -- | Calls the game loop every frame as long as the window is open.
 --  For larger projects, instead of using this function, consider
diff --git a/src/Raylib/Util/Camera.hs b/src/Raylib/Util/Camera.hs
new file mode 100644
--- /dev/null
+++ b/src/Raylib/Util/Camera.hs
@@ -0,0 +1,187 @@
+{-# OPTIONS -Wall #-}
+
+module Raylib.Util.Camera
+  ( getCameraForward,
+    getCameraUp,
+    getCameraRight,
+    cameraMove,
+    cameraMoveForward,
+    cameraMoveUp,
+    cameraMoveRight,
+    cameraRotate,
+    cameraYaw,
+    cameraPitch,
+    cameraRoll,
+    getCameraViewMatrix,
+    getCameraProjectionMatrix,
+  )
+where
+
+import Raylib.Types (Camera, Camera3D (..), CameraProjection (CameraOrthographic, CameraPerspective), Matrix, Vector3 (..))
+import Raylib.Util.Math (Vector (..), clamp, deg2Rad, matrixLookAt, matrixOrtho, matrixPerspective, vector3Angle, vector3CrossProduct, vector3RotateByAxisAngle)
+
+-- | The camera's forward vector (normalized)
+getCameraForward :: Camera -> Vector3
+getCameraForward cam = vectorNormalize $ camera3D'target cam |-| camera3D'position cam
+
+-- | The camera's up vector (normalized)
+getCameraUp :: Camera -> Vector3
+getCameraUp cam = vectorNormalize $ camera3D'up cam
+
+-- | The camera's right vector (normalized)
+getCameraRight :: Camera -> Vector3
+getCameraRight cam = vector3CrossProduct (getCameraForward cam) (getCameraUp cam)
+
+-- | Move the camera by a specific vector
+cameraMove :: Camera -> Vector3 -> Camera
+cameraMove cam dir =
+  cam {camera3D'position = camera3D'position cam |+| dir, camera3D'target = camera3D'target cam |+| dir}
+
+-- | Move the camera in its forward direction
+cameraMoveForward ::
+  Camera ->
+  -- | Distance to move
+  Float ->
+  -- | Move in world plane (i.e. no vertical movement if enabled)
+  Bool ->
+  Camera
+cameraMoveForward cam distance moveInWorldPlane =
+  cameraMove cam (forward |* distance)
+  where
+    forward = if moveInWorldPlane then camForward {vector3'y = 0} else camForward
+    camForward = getCameraForward cam
+
+-- | Move the camera in its up direction
+cameraMoveUp ::
+  Camera ->
+  -- | Distance to move
+  Float ->
+  Camera
+cameraMoveUp cam distance =
+  cameraMove cam (up |* distance)
+  where
+    up = getCameraUp cam
+
+-- | Move the camera in its right direction
+cameraMoveRight ::
+  Camera ->
+  -- | Distance to move
+  Float ->
+  -- | Move in world plane (i.e. no vertical movement if enabled)
+  Bool ->
+  Camera
+cameraMoveRight cam distance moveInWorldPlane =
+  cameraMove cam (right |* distance)
+  where
+    right = if moveInWorldPlane then camRight {vector3'y = 0} else camRight
+    camRight = getCameraRight cam
+
+-- | Rotate the camera using an axis and angle
+cameraRotate ::
+  Camera ->
+  -- | Axis of rotation
+  Vector3 ->
+  -- | Angle to rotate by
+  Float ->
+  -- | Rotate around target (if false, the camera rotates around its position)
+  Bool ->
+  Camera
+cameraRotate cam axis angle rotateAroundTarget =
+  cam
+    { camera3D'position = if rotateAroundTarget then target |-| viewRot else pos,
+      camera3D'target = if rotateAroundTarget then target else pos |+| viewRot
+    }
+  where
+    viewVec = target |-| pos
+    viewRot = vector3RotateByAxisAngle viewVec axis angle
+    pos = camera3D'position cam
+    target = camera3D'target cam
+
+-- | Rotate the camera around its up vector.
+--   Yaw is "looking left and right".
+cameraYaw ::
+  Camera ->
+  -- | Angle in radians
+  Float ->
+  -- | Rotate around target (if false, the camera rotates around its position)
+  Bool ->
+  Camera
+cameraYaw cam angle rotateAroundTarget =
+  cam
+    { camera3D'position = if rotateAroundTarget then target |-| viewRot else pos,
+      camera3D'target = if rotateAroundTarget then target else pos |+| viewRot
+    }
+  where
+    viewVec = target |-| pos
+    viewRot = vector3RotateByAxisAngle viewVec (getCameraUp cam) angle
+    pos = camera3D'position cam
+    target = camera3D'target cam
+
+-- | Rotate the camera around its right vector.
+--   Pitch is "looking up and down".
+cameraPitch ::
+  Camera ->
+  -- | Angle in radians
+  Float ->
+  -- | Lock view (prevents camera overrotation, aka "somersaults")
+  Bool ->
+  -- | Rotate around target (if false, the camera rotates around its position)
+  Bool ->
+  -- | Rotate the camera's up vector to match the new pitch
+  Bool ->
+  Camera
+cameraPitch cam angle lockView rotateAroundTarget rotateUp =
+  cam
+    { camera3D'position = if rotateAroundTarget then target |-| viewRot else pos,
+      camera3D'target = if rotateAroundTarget then target else pos |+| viewRot,
+      camera3D'up = if not rotateUp then up else vector3RotateByAxisAngle up right angle'
+    }
+  where
+    angle' = if not lockView then angle else clamp angle maxAngleDown maxAngleUp
+    maxAngleUp = vector3Angle up viewVec - 0.001
+    maxAngleDown = (- vector3Angle (additiveInverse up) viewVec) - 0.001
+
+    viewVec = target |-| pos
+    viewRot = vector3RotateByAxisAngle viewVec right angle'
+
+    pos = camera3D'position cam
+    target = camera3D'target cam
+    up = getCameraUp cam
+    right = getCameraRight cam
+
+-- | Rotates the camera around its forward vector.
+--   Roll is "turning your head sideways to the left or right".
+cameraRoll ::
+  Camera ->
+  -- | Angle in radians
+  Float ->
+  Camera
+cameraRoll cam angle =
+  cam
+    { camera3D'up = vector3RotateByAxisAngle up forward angle
+    }
+  where
+    forward = getCameraForward cam
+    up = getCameraUp cam
+
+-- | View matrix from camera
+getCameraViewMatrix :: Camera -> Matrix
+getCameraViewMatrix cam = matrixLookAt (camera3D'position cam) (camera3D'target cam) (camera3D'up cam)
+
+-- | Projection matrix from camera
+getCameraProjectionMatrix ::
+  Camera ->
+  -- | Aspect ratio
+  Float ->
+  -- | Near clipping plane distance (recommended: 0.01)
+  Float ->
+  -- | Far clipping plane distance (recommended: 1000)
+  Float ->
+  Matrix
+getCameraProjectionMatrix cam aspect near far =
+  case camera3D'projection cam of
+    CameraPerspective -> matrixPerspective (camera3D'fovy cam * deg2Rad) aspect near far
+    CameraOrthographic -> matrixOrtho (- right) right (- top) top near far
+      where
+        top = camera3D'fovy cam / 2
+        right = top * aspect
diff --git a/src/Raylib/Util/Math.hs b/src/Raylib/Util/Math.hs
new file mode 100644
--- /dev/null
+++ b/src/Raylib/Util/Math.hs
@@ -0,0 +1,1005 @@
+{-# OPTIONS -Wall #-}
+
+module Raylib.Util.Math
+  ( -- * Utility constants
+    epsilon,
+    deg2Rad,
+    rad2Deg,
+
+    -- * Float math
+    clamp,
+    lerp,
+    normalize,
+    remap,
+    wrap,
+    floatEquals,
+
+    -- * General vector math
+    Vector (..),
+
+    -- * Vector2 math
+    vector2Angle,
+    vector2LineAngle,
+    vector2Transform,
+    vector2Reflect,
+    vector2Rotate,
+
+    -- * Vector3 math
+    vector3CrossProduct,
+    vector3Perpendicular,
+    vector3Angle,
+    vector3OrthoNormalize,
+    vector3Transform,
+    vector3RotateByQuaternion,
+    vector3RotateByAxisAngle,
+    vector3Reflect,
+    vector3Barycenter,
+    vector3Unproject,
+    vector3Refract,
+
+    -- * Matrix math
+    matrixToList,
+    matrixFromList,
+    matrixConstant,
+    matrixDeterminant,
+    matrixTrace,
+    matrixTranspose,
+    matrixInvert,
+    matrixIdentity,
+    matrixAdd,
+    (/+/),
+    matrixSubtract,
+    (/-/),
+    matrixMultiply,
+    (/*/),
+    matrixTranslate,
+    matrixRotate,
+    matrixRotateX,
+    matrixRotateY,
+    matrixRotateZ,
+    matrixRotateXYZ,
+    matrixRotateZYX,
+    matrixScale,
+    matrixFrustum,
+    matrixPerspective,
+    matrixOrtho,
+    matrixLookAt,
+
+    -- * Quaternion math
+    quaternionIdentity,
+    quaternionInvert,
+    quaternionMultiply,
+    quaternionNormalize,
+    quaternionLerp,
+    quaternionNLerp,
+    quaternionSLerp,
+    quaternionFromVector3ToVector3,
+    quaternionFromMatrix,
+    quaternionToMatrix,
+    quaternionFromAxisAngle,
+    quaternionToAxisAngle,
+    quaternionFromEuler,
+    quaternionToEuler,
+    quaternionTransform,
+  )
+where
+
+import Raylib.Types (Matrix (..), Quaternion, Vector2 (Vector2), Vector3 (Vector3), Vector4 (Vector4))
+
+epsilon :: Float
+epsilon = 0.000001
+
+deg2Rad :: Float
+deg2Rad = pi / 180
+
+rad2Deg :: Float
+rad2Deg = 180 / pi
+
+------------------------------------------------
+-- Float math ----------------------------------
+------------------------------------------------
+
+-- | Clamp float to range
+clamp ::
+  -- | Value to clamp
+  Float ->
+  -- | Lower bound
+  Float ->
+  -- | Upper bound
+  Float ->
+  Float
+clamp value low high
+  | value < low = low
+  | value > high = high
+  | otherwise = 0
+
+-- | Calculate linear interpolation between two floats
+lerp ::
+  -- | Starting value
+  Float ->
+  -- | Ending value
+  Float ->
+  -- | Lerp amount
+  Float ->
+  Float
+lerp start end amount = start + amount * (end - start)
+
+-- | Normalize input value within input range
+normalize ::
+  -- | Value to normalize
+  Float ->
+  -- | Starting value of range
+  Float ->
+  -- | Ending value of range
+  Float ->
+  Float
+normalize value start end = (value - start) / (end - start)
+
+-- | Remap input value within input range to output range
+remap ::
+  -- | Input value
+  Float ->
+  -- | Input range start
+  Float ->
+  -- | Input range end
+  Float ->
+  -- | Output range start
+  Float ->
+  -- | Output range end
+  Float ->
+  Float
+remap value inputStart inputEnd outputStart outputEnd = (value - inputStart) / (inputEnd - inputStart) * (outputEnd - outputStart) + outputStart
+
+-- | Wrap input value from min to max
+wrap ::
+  -- | Input value
+  Float ->
+  -- | Min value
+  Float ->
+  -- | Max value
+  Float ->
+  Float
+wrap value low high = value - (high - low) * fromIntegral (floor ((value - low) / (high - low)) :: Integer)
+
+-- | Check if two floats are close to equal
+floatEquals :: Float -> Float -> Bool
+floatEquals x y = abs (x - y) <= (epsilon * max 1 (max (abs x) (abs y)))
+
+------------------------------------------------
+-- Vector math ---------------------------------
+------------------------------------------------
+
+class Vector a where
+  -- | List representation of a vector
+  asList :: a -> [Float]
+
+  -- | Vector representation of a list
+  fromList :: [Float] -> a
+
+  -- | Vector-vector addition
+  (|+|) :: a -> a -> a
+  a |+| b = fromList $ zipWith (+) (asList a) (asList b)
+
+  -- | Vector-vector subtraction
+  (|-|) :: a -> a -> a
+  a |-| b = fromList $ zipWith (-) (asList a) (asList b)
+
+  -- | Vector-scalar addition
+  (|+) :: a -> Float -> a
+  a |+ b = a |+| constant b
+
+  -- | Vector-scalar subtraction
+  (|-) :: a -> Float -> a
+  a |- b = a |-| constant b
+
+  -- | Vector-vector multiplication
+  (|*|) :: a -> a -> a
+  a |*| b = fromList $ zipWith (*) (asList a) (asList b)
+
+  -- | Vector-vector division
+  (|/|) :: a -> a -> a
+  a |/| b = fromList $ zipWith (/) (asList a) (asList b)
+
+  -- | Vector-scalar multiplication
+  (|*) :: a -> Float -> a
+  a |* b = a |*| constant b
+
+  -- | Vector-scalar division
+  (|/) :: a -> Float -> a
+  a |/ b = a |/| constant b
+
+  -- | Vector-vector dot product
+  (|.|) :: a -> a -> Float
+  a |.| b = sum . asList $ a |*| b
+
+  -- | Zero vector
+  zero :: a
+  zero = constant 0
+
+  -- | One vector
+  one :: a
+  one = constant 1
+
+  -- | Scalar to vector (all elements are set to the scalar)
+  constant :: Float -> a
+  constant val = fromList $ repeat val
+
+  -- | Vector additive inverse
+  additiveInverse :: a -> a
+  additiveInverse v = fromList $ map negate (asList v)
+
+  -- | Vector multiplicative inverse
+  multiplicativeInverse :: a -> a
+  multiplicativeInverse v = fromList $ map (1 /) (asList v)
+
+  -- | Squared magnitude of a vector
+  magnitudeSqr :: a -> Float
+  magnitudeSqr v = v |.| v
+
+  -- | Vector magnitude
+  magnitude :: a -> Float
+  magnitude v = if m == 1 then m else sqrt m
+    where
+      m = magnitudeSqr v
+
+  -- | Squared distance between two vectors
+  vectorDistanceSqr :: a -> a -> Float
+  vectorDistanceSqr a b = magnitudeSqr $ a |-| b
+
+  -- | Distance between two vectors
+  vectorDistance :: a -> a -> Float
+  vectorDistance a b = sqrt $ vectorDistanceSqr a b
+
+  -- | Normalize vector (same direction, magnitude 1)
+  vectorNormalize :: a -> a
+  vectorNormalize v = v |/ magnitude v
+
+  -- | Lerp between two vectors
+  vectorLerp :: a -> a -> Float -> a
+  vectorLerp a b amount = fromList $ zipWith (\v1 v2 -> lerp v1 v2 amount) (asList a) (asList b)
+
+  -- | Move vector towards target
+  vectorMoveTowards ::
+    -- | Vector to move
+    a ->
+    -- | Target vector
+    a ->
+    -- | Max distance to move by
+    Float ->
+    a
+  vectorMoveTowards v target maxDistance =
+    if distSquared <= maxDistance * maxDistance
+      then target
+      else v |+| fromList (map (* (maxDistance / dist)) (asList diff))
+    where
+      diff = target |-| v
+      distSquared = magnitudeSqr diff
+      dist = sqrt distSquared
+
+  -- | Clamp vector to range
+  vectorClamp ::
+    -- | Vector to clamp
+    a ->
+    -- | Lower bound
+    a ->
+    -- | Upper bound
+    a ->
+    a
+  vectorClamp v low high = fromList $ zipWith3 clamp (asList v) (asList low) (asList high)
+
+  -- | Clamp the magnitude of a vector to a range
+  vectorClampValue ::
+    -- | Vector to clamp
+    a ->
+    -- | Lower bound
+    Float ->
+    -- | Upper bound
+    Float ->
+    a
+  vectorClampValue v low high = v |* (clamp size low high / size)
+    where
+      size = magnitude v
+
+  -- | Min value for each pair of components
+  vectorMin :: a -> a -> a
+  vectorMin v1 v2 = fromList $ zipWith min (asList v1) (asList v2)
+
+  -- | Max value for each pair of components
+  vectorMax :: a -> a -> a
+  vectorMax v1 v2 = fromList $ zipWith max (asList v1) (asList v2)
+
+instance Vector Vector2 where
+  asList (Vector2 x y) = [x, y]
+  fromList (x : y : _) = Vector2 x y
+  fromList _ = error "(Vector2.fromList) Input list must have at least two elements!"
+
+instance Vector Vector3 where
+  asList (Vector3 x y z) = [x, y, z]
+  fromList (x : y : z : _) = Vector3 x y z
+  fromList _ = error "(Vector3.fromList) Input list must have at least three elements!"
+
+instance Vector Vector4 where
+  asList (Vector4 x y z w) = [x, y, z, w]
+  fromList (x : y : z : w : _) = Vector4 x y z w
+  fromList _ = error "(Vector4.fromList) Input list must have at least four elements!"
+
+------------------------------------------------
+-- Vector2 math --------------------------------
+------------------------------------------------
+
+-- | Angle between two 2D vectors
+vector2Angle :: Vector2 -> Vector2 -> Float
+vector2Angle (Vector2 x1 y1) (Vector2 x2 y2) = atan2 (y2 - y1) (x2 - x1)
+
+-- | Angle created by the line between two 2D vectors (parameters must be normalized)
+vector2LineAngle :: Vector2 -> Vector2 -> Float
+vector2LineAngle a b = acos $ clamp (a |.| b) (-1) 1
+
+-- | Transform a 2D vector by the given matrix
+vector2Transform :: Vector2 -> Matrix -> Vector2
+vector2Transform (Vector2 x y) mat = Vector2 ((x * matrix'm0 mat) + (y * matrix'm4 mat) + matrix'm12 mat) ((x * matrix'm1 mat) + (y * matrix'm5 mat) + matrix'm13 mat)
+
+-- | Reflect 2D vector to normal
+vector2Reflect ::
+  -- | Input vector
+  Vector2 ->
+  -- | Normal vector
+  Vector2 ->
+  Vector2
+vector2Reflect v normal = v |-| (normal |* (2 * (v |.| normal)))
+
+-- | Rotate 2D vector by angle
+vector2Rotate :: Vector2 -> Float -> Vector2
+vector2Rotate (Vector2 x y) angle = Vector2 (x * c - y * s) (x * s + y * c)
+  where
+    c = cos angle
+    s = sin angle
+
+------------------------------------------------
+-- Vector3 math --------------------------------
+------------------------------------------------
+
+-- | 3D vector cross-product
+vector3CrossProduct :: Vector3 -> Vector3 -> Vector3
+vector3CrossProduct (Vector3 x1 y1 z1) (Vector3 x2 y2 z2) = Vector3 (y1 * z2 - z1 * y2) (z1 * x2 - x1 * z2) (x1 * y2 - y1 * x2)
+
+-- | Perpendicular to given 3D vector
+vector3Perpendicular :: Vector3 -> Vector3
+vector3Perpendicular v@(Vector3 x y z) = vector3CrossProduct v cardinalAxis
+  where
+    cardinalAxis
+      | abs y < abs x = Vector3 0 1 0
+      | abs z < abs x = Vector3 0 0 1
+      | otherwise = Vector3 1 0 0
+
+-- | Angle between two 3D vectors
+vector3Angle :: Vector3 -> Vector3 -> Float
+vector3Angle v1 v2 = atan2 len dot
+  where
+    cross = vector3CrossProduct v1 v2
+    len = magnitude cross
+    dot = v1 |.| v2
+
+-- | Orthonormalize provided vectors.
+--   Makes vectors normalized and orthogonal to each other.
+--   Gram-Schmidt function implementation.
+vector3OrthoNormalize :: Vector3 -> Vector3 -> (Vector3, Vector3)
+vector3OrthoNormalize v1 v2 = (v1', v2')
+  where
+    v1' = vectorNormalize v1
+    v2' = vector3CrossProduct vn1 v1'
+    vn1 = vectorNormalize $ vector3CrossProduct v1' v2
+
+-- | Transform a 3D vector by a matrix
+vector3Transform :: Vector3 -> Matrix -> Vector3
+vector3Transform
+  (Vector3 x y z)
+  (Matrix m0 m4 m8 m12 m1 m5 m9 m13 m2 m6 m10 m14 _ _ _ _) =
+    Vector3
+      (x * m0 + y * m4 + z * m8 + m12)
+      (x * m1 + y * m5 + z * m9 + m13)
+      (x * m2 + y * m6 + z * m10 + m14)
+
+-- | Transform a 3D vector by quaternion rotation
+vector3RotateByQuaternion :: Vector3 -> Quaternion -> Vector3
+vector3RotateByQuaternion (Vector3 a b c) (Vector4 x y z w) =
+  Vector3
+    (a * (x * x + w * w - y * y - z * z) + b * (2 * x * y - 2 * w * z) + c * (2 * x * z + 2 * w * y))
+    (a * (2 * w * z + 2 * x * y) + b * (w * w - x * x + y * y - z * z) + c * ((-2) * w * x + 2 * y * z))
+    (a * ((-2) * w * y + 2 * x * z) + b * (2 * w * x + 2 * y * z) + c * (w * w - x * x - y * y + z * z))
+
+-- | Rotate a 3D vector around an axis
+vector3RotateByAxisAngle ::
+  -- | Vector to rotate
+  Vector3 ->
+  -- | Axis to rotate around
+  Vector3 ->
+  -- | Angle to rotate by
+  Float ->
+  Vector3
+vector3RotateByAxisAngle v axis angle = v |+| (wv |* (2 * a)) |+| (wwv |* 2)
+  where
+    (Vector3 ax ay az) = vectorNormalize axis
+    s = sin (angle / 2)
+    a = cos (angle / 2)
+    b = ax * s
+    c = ay * s
+    d = az * s
+    w = Vector3 b c d
+    wv = vector3CrossProduct w v
+    wwv = vector3CrossProduct w wv
+
+-- | Reflect 3D vector to normal
+vector3Reflect ::
+  -- | Input vector
+  Vector3 ->
+  -- | Normal vector
+  Vector3 ->
+  Vector3
+vector3Reflect v normal = v |-| (normal |* (2 * v |.| normal))
+
+-- | Compute barycenter coordinates (u, v, w) for a point with respect to a triangle.
+--   NOTE: Assumes the point is on the plane of the triangle.
+vector3Barycenter ::
+  -- | Input point
+  Vector3 ->
+  -- | Triangle vertices
+  (Vector3, Vector3, Vector3) ->
+  Vector3
+vector3Barycenter p (a, b, c) = Vector3 (1 - y - z) y z
+  where
+    v0 = b |-| a
+    v1 = c |-| a
+    v2 = p |-| a
+    d00 = v0 |.| v0
+    d01 = v0 |.| v1
+    d11 = v1 |.| v1
+    d20 = v2 |.| v0
+    d21 = v2 |.| v1
+    denom = d00 * d11 - d01 * d01
+    y = (d11 * d20 - d01 * d21) / denom
+    z = (d00 * d21 - d01 * d20) / denom
+
+-- | Project a Vector3 from screen space into object space
+vector3Unproject ::
+  -- | Vector to unproject
+  Vector3 ->
+  -- | Projection matrix
+  Matrix ->
+  -- | View matrix
+  Matrix ->
+  Vector3
+vector3Unproject (Vector3 x y z) projection view = Vector3 (rx / rw) (ry / rw) (rz / rw)
+  where
+    matViewProj = view /*/ projection
+    a00 = matrix'm0 matViewProj
+    a01 = matrix'm1 matViewProj
+    a02 = matrix'm2 matViewProj
+    a03 = matrix'm3 matViewProj
+    a10 = matrix'm4 matViewProj
+    a11 = matrix'm5 matViewProj
+    a12 = matrix'm6 matViewProj
+    a13 = matrix'm7 matViewProj
+    a20 = matrix'm8 matViewProj
+    a21 = matrix'm9 matViewProj
+    a22 = matrix'm10 matViewProj
+    a23 = matrix'm11 matViewProj
+    a30 = matrix'm12 matViewProj
+    a31 = matrix'm13 matViewProj
+    a32 = matrix'm14 matViewProj
+    a33 = matrix'm15 matViewProj
+    b00 = a00 * a11 - a01 * a10
+    b01 = a00 * a12 - a02 * a10
+    b02 = a00 * a13 - a03 * a10
+    b03 = a01 * a12 - a02 * a11
+    b04 = a01 * a13 - a03 * a11
+    b05 = a02 * a13 - a03 * a12
+    b06 = a20 * a31 - a21 * a30
+    b07 = a20 * a32 - a22 * a30
+    b08 = a20 * a33 - a23 * a30
+    b09 = a21 * a32 - a22 * a31
+    b10 = a21 * a33 - a23 * a31
+    b11 = a22 * a33 - a23 * a32
+    invDet = 1 / (b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06)
+    matViewProjInv =
+      Matrix
+        ((a11 * b11 - a12 * b10 + a13 * b09) * invDet)
+        ((- a01 * b11 + a02 * b10 - a03 * b09) * invDet)
+        ((a31 * b05 - a32 * b04 + a33 * b03) * invDet)
+        ((- a21 * b05 + a22 * b04 - a23 * b03) * invDet)
+        ((- a10 * b11 + a12 * b08 - a13 * b07) * invDet)
+        ((a00 * b11 - a02 * b08 + a03 * b07) * invDet)
+        ((- a30 * b05 + a32 * b02 - a33 * b01) * invDet)
+        ((a20 * b05 - a22 * b02 + a23 * b01) * invDet)
+        ((a10 * b10 - a11 * b08 + a13 * b06) * invDet)
+        ((- a00 * b10 + a01 * b08 - a03 * b06) * invDet)
+        ((a30 * b04 - a31 * b02 + a33 * b00) * invDet)
+        ((- a20 * b04 + a21 * b02 - a23 * b00) * invDet)
+        ((- a10 * b09 + a11 * b07 - a12 * b06) * invDet)
+        ((a00 * b09 - a01 * b07 + a02 * b06) * invDet)
+        ((- a30 * b03 + a31 * b01 - a32 * b00) * invDet)
+        ((a20 * b03 - a21 * b01 + a22 * b00) * invDet)
+    (Vector4 rx ry rz rw) = quaternionTransform (Vector4 x y z 1) matViewProjInv
+
+-- | Compute the direction of a refracted ray
+vector3Refract ::
+  -- | Normalized direction of the incoming ray
+  Vector3 ->
+  -- | Normalized normal vector of the interface of two optical media
+  Vector3 ->
+  -- | Ratio of the refractive index of the medium from where the ray
+  --   comes to the refractive index of the medium on the other side of
+  --   the surface
+  Float ->
+  Vector3
+vector3Refract (Vector3 x y z) (Vector3 nx ny nz) r = Vector3 (r * x - (r * dot + d) * nx) (r * y - (r * dot + d) * ny) (r * z - (r * dot + d) * nz)
+  where
+    dot = x * nx + y * ny + z * nz
+    d = sqrt $ 1 - r * r * (1 - dot * dot)
+
+------------------------------------------------
+-- Matrix math ---------------------------------
+------------------------------------------------
+
+-- | Utility function
+matrixToList :: Matrix -> [Float]
+matrixToList
+  (Matrix a00 a10 a20 a30 a01 a11 a21 a31 a02 a12 a22 a32 a03 a13 a23 a33) =
+    [a00, a10, a20, a30, a01, a11, a21, a31, a02, a12, a22, a32, a03, a13, a23, a33]
+
+-- | Utility function
+matrixFromList :: [Float] -> Matrix
+matrixFromList
+  (a00 : a10 : a20 : a30 : a01 : a11 : a21 : a31 : a02 : a12 : a22 : a32 : a03 : a13 : a23 : a33 : _) =
+    Matrix a00 a10 a20 a30 a01 a11 a21 a31 a02 a12 a22 a32 a03 a13 a23 a33
+matrixFromList _ = error "matrixFromList expects a list of at least 16 elements!"
+
+-- | Scalar to matrix (all elements are set to the scalar)
+matrixConstant :: Float -> Matrix
+matrixConstant n = matrixFromList $ repeat n
+
+-- | Compute matrix determinant
+matrixDeterminant :: Matrix -> Float
+matrixDeterminant
+  (Matrix a00 a10 a20 a30 a01 a11 a21 a31 a02 a12 a22 a32 a03 a13 a23 a33) =
+    a30 * a21 * a12 * a03 - a20 * a31 * a12 * a03 - a30 * a11 * a22 * a03 + a10 * a31 * a22 * a03
+      + a20 * a11 * a32 * a03 - a10 * a21 * a32 * a03 - a30 * a21 * a02 * a13
+      + a20 * a31 * a02 * a13
+      + a30 * a01 * a22 * a13 - a00 * a31 * a22 * a13 - a20 * a01 * a32 * a13
+      + a00 * a21 * a32 * a13
+      + a30 * a11 * a02 * a23 - a10 * a31 * a02 * a23 - a30 * a01 * a12 * a23
+      + a00 * a31 * a12 * a23
+      + a10 * a01 * a32 * a23 - a00 * a11 * a32 * a23 - a20 * a11 * a02 * a33
+      + a10 * a21 * a02 * a33
+      + a20 * a01 * a12 * a33 - a00 * a21 * a12 * a33 - a10 * a01 * a22 * a33
+      + a00 * a11 * a22 * a33
+
+-- | Trace of a matrix (sum of the values along the diagonal)
+matrixTrace :: Matrix -> Float
+matrixTrace mat = matrix'm0 mat + matrix'm5 mat + matrix'm10 mat + matrix'm15 mat
+
+-- | Transpose a matrix
+matrixTranspose :: Matrix -> Matrix
+matrixTranspose
+  (Matrix m0 m4 m8 m12 m1 m5 m9 m13 m2 m6 m10 m14 m3 m7 m11 m15) =
+    Matrix m0 m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15
+
+-- | Invert a matrix
+matrixInvert :: Matrix -> Matrix
+matrixInvert
+  (Matrix a00 a10 a20 a30 a01 a11 a21 a31 a02 a12 a22 a32 a03 a13 a23 a33) =
+    Matrix
+      ((a11 * b11 - a12 * b10 + a13 * b09) * invDet)
+      ((- a10 * b11 + a12 * b08 - a13 * b07) * invDet)
+      ((a10 * b10 - a11 * b08 + a13 * b06) * invDet)
+      ((- a10 * b09 + a11 * b07 - a12 * b06) * invDet)
+      ((- a01 * b11 + a02 * b10 - a03 * b09) * invDet)
+      ((a00 * b11 - a02 * b08 + a03 * b07) * invDet)
+      ((- a00 * b10 + a01 * b08 - a03 * b06) * invDet)
+      ((a00 * b09 - a01 * b07 + a02 * b06) * invDet)
+      ((a31 * b05 - a32 * b04 + a33 * b03) * invDet)
+      ((- a30 * b05 + a32 * b02 - a33 * b01) * invDet)
+      ((a30 * b04 - a31 * b02 + a33 * b00) * invDet)
+      ((- a30 * b03 + a31 * b01 - a32 * b00) * invDet)
+      ((- a21 * b05 + a22 * b04 - a23 * b03) * invDet)
+      ((a20 * b05 - a22 * b02 + a23 * b01) * invDet)
+      ((- a20 * b04 + a21 * b02 - a23 * b00) * invDet)
+      ((a20 * b03 - a21 * b01 + a22 * b00) * invDet)
+    where
+      b00 = a00 * a11 - a01 * a10
+      b01 = a00 * a12 - a02 * a10
+      b02 = a00 * a13 - a03 * a10
+      b03 = a01 * a12 - a02 * a11
+      b04 = a01 * a13 - a03 * a11
+      b05 = a02 * a13 - a03 * a12
+      b06 = a20 * a31 - a21 * a30
+      b07 = a20 * a32 - a22 * a30
+      b08 = a20 * a33 - a23 * a30
+      b09 = a21 * a32 - a22 * a31
+      b10 = a21 * a33 - a23 * a31
+      b11 = a22 * a33 - a23 * a32
+      invDet = 1 / (b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06)
+
+-- | Identity matrix
+matrixIdentity :: Matrix
+matrixIdentity = Matrix 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
+
+-- | Add two matrices
+matrixAdd :: Matrix -> Matrix -> Matrix
+matrixAdd left right = matrixFromList $ zipWith (+) (matrixToList left) (matrixToList right)
+
+-- | Alias for 'matrixAdd'
+(/+/) :: Matrix -> Matrix -> Matrix
+(/+/) = matrixAdd
+
+-- | Subtract two matrices
+matrixSubtract :: Matrix -> Matrix -> Matrix
+matrixSubtract left right = matrixFromList $ zipWith (-) (matrixToList left) (matrixToList right)
+
+-- | Alias for 'matrixSubtract'
+(/-/) :: Matrix -> Matrix -> Matrix
+(/-/) = matrixSubtract
+
+-- | Multiply two matrices (order matters!)
+matrixMultiply :: Matrix -> Matrix -> Matrix
+matrixMultiply
+  (Matrix l0 l4 l8 l12 l1 l5 l9 l13 l2 l6 l10 l14 l3 l7 l11 l15)
+  (Matrix r0 r4 r8 r12 r1 r5 r9 r13 r2 r6 r10 r14 r3 r7 r11 r15) =
+    Matrix
+      (l0 * r0 + l1 * r4 + l2 * r8 + l3 * r12)
+      (l4 * r0 + l5 * r4 + l6 * r8 + l7 * r12)
+      (l8 * r0 + l9 * r4 + l10 * r8 + l11 * r12)
+      (l12 * r0 + l13 * r4 + l14 * r8 + l15 * r12)
+      (l0 * r1 + l1 * r5 + l2 * r9 + l3 * r13)
+      (l4 * r1 + l5 * r5 + l6 * r9 + l7 * r13)
+      (l8 * r1 + l9 * r5 + l10 * r9 + l11 * r13)
+      (l12 * r1 + l13 * r5 + l14 * r9 + l15 * r13)
+      (l0 * r2 + l1 * r6 + l2 * r10 + l3 * r14)
+      (l4 * r2 + l5 * r6 + l6 * r10 + l7 * r14)
+      (l8 * r2 + l9 * r6 + l10 * r10 + l11 * r14)
+      (l12 * r2 + l13 * r6 + l14 * r10 + l15 * r14)
+      (l0 * r3 + l1 * r7 + l2 * r11 + l3 * r15)
+      (l4 * r3 + l5 * r7 + l6 * r11 + l7 * r15)
+      (l8 * r3 + l9 * r7 + l10 * r11 + l11 * r15)
+      (l12 * r3 + l13 * r7 + l14 * r11 + l15 * r15)
+
+-- | Alias for 'matrixMultiply'
+(/*/) :: Matrix -> Matrix -> Matrix
+(/*/) = matrixMultiply
+
+-- | Translation matrix
+matrixTranslate ::
+  -- | x translation
+  Float ->
+  -- | y translation
+  Float ->
+  -- | z translation
+  Float ->
+  Matrix
+matrixTranslate x y z = Matrix 1 0 0 x 0 1 0 y 0 0 1 z 0 0 0 1
+
+-- | Axis-angle rotation matrix (angle should be in radians)
+matrixRotate ::
+  -- | Axis to rotate around
+  Vector3 ->
+  -- | Angle to rotate by
+  Float ->
+  Matrix
+matrixRotate axis angle = Matrix (x * x * t + c) (x * y * t - z * s) (x * z * t + y * s) 0 (y * x * t + z * s) (y * y * t + c) (y * z * t - x * s) 0 (z * x * t - y * s) (z * y * t + x * s) (z * z * t + c) 0 0 0 0 1
+  where
+    (Vector3 x y z) = vectorNormalize axis
+    s = sin angle
+    c = cos angle
+    t = 1 - c
+
+-- | x-rotation matrix (angle should be in radians)
+matrixRotateX :: Float -> Matrix
+matrixRotateX angle = Matrix 1 0 0 0 0 c (- s) 0 0 s c 0 0 0 0 1
+  where
+    s = sin angle
+    c = cos angle
+
+-- | y-rotation matrix (angle should be in radians)
+matrixRotateY :: Float -> Matrix
+matrixRotateY angle = Matrix c 0 s 0 0 1 0 0 (- s) 0 c 0 0 0 0 1
+  where
+    s = sin angle
+    c = cos angle
+
+-- | z-rotation matrix (angle should be in radians)
+matrixRotateZ :: Float -> Matrix
+matrixRotateZ angle = Matrix c (- s) 0 0 s c 0 0 0 0 1 0 0 0 0 1
+  where
+    s = sin angle
+    c = cos angle
+
+-- | Euler angle xyz rotation matrix (angles should be in radians)
+matrixRotateXYZ :: Vector3 -> Matrix
+matrixRotateXYZ (Vector3 x y z) = Matrix (cz - cy) (sz * cy) (- sy) 0 (cz * sy * sx - sz * cx) (sz * sy * sx + cz * cx) (cy * sx) 0 (cz * sy * cx + sz * sx) (sz * sy * cx - cz * sx) (cy * cx) 0 0 0 0 1
+  where
+    cx = cos (- x)
+    sx = sin (- x)
+    cy = cos (- y)
+    sy = sin (- y)
+    cz = cos (- z)
+    sz = sin (- z)
+
+-- | Euler angle zyx rotation matrix (angles should be in radians)
+matrixRotateZYX :: Vector3 -> Matrix
+matrixRotateZYX (Vector3 x y z) = Matrix (cz * cy) (cz * sy * sx - cx * sz) (sz * sx + cz * cx * sy) 0 (cy * sz) (cz * cx + sz * sy * sx) (cx * sz * sy - cz * sx) 0 (- sy) (cy * sx) (cy * cx) 0 0 0 0 1
+  where
+    cz = cos z
+    sz = sin z
+    cy = cos y
+    sy = sin y
+    cx = cos x
+    sx = sin x
+
+-- | Scaling matrix
+matrixScale :: Vector3 -> Matrix
+matrixScale (Vector3 x y z) = Matrix x 0 0 0 0 y 0 0 0 0 z 0 0 0 0 1
+
+-- | Frustum projection matrix
+matrixFrustum ::
+  -- | Left edge distance
+  Float ->
+  -- | Right edge distance
+  Float ->
+  -- | Bottom edge distance
+  Float ->
+  -- | Top edge distance
+  Float ->
+  -- | Near clipping plane distance
+  Float ->
+  -- | Far clipping plane distance
+  Float ->
+  Matrix
+matrixFrustum left right bottom top near far =
+  Matrix
+    (near * 2 / x)
+    0
+    ((right + left) / x)
+    0
+    0
+    (near * 2 / y)
+    ((top + bottom) / y)
+    0
+    0
+    0
+    (- (far + near) / z)
+    (- far * near * 2 / z)
+    0
+    0
+    (-1)
+    0
+  where
+    x = right - left
+    y = top - bottom
+    z = far - near
+
+-- | Perspective projection matrix
+matrixPerspective ::
+  -- | y-fov angle (should be in radians)
+  Float ->
+  -- | Aspect ratio
+  Float ->
+  -- | Near clipping plane distance
+  Float ->
+  -- | Far clipping plane distance
+  Float ->
+  Matrix
+matrixPerspective fovy aspect near far = matrixFrustum left right bottom top near far
+  where
+    top = near * tan (fovy / 2)
+    bottom = - top
+    right = top * aspect
+    left = - right
+
+-- | Orthographic projection matrix
+matrixOrtho ::
+  -- | Left edge distance
+  Float ->
+  -- | Right edge distance
+  Float ->
+  -- | Bottom edge distance
+  Float ->
+  -- | Top edge distance
+  Float ->
+  -- | Near clipping plane distance
+  Float ->
+  -- | Far clipping plane distance
+  Float ->
+  Matrix
+matrixOrtho left right bottom top near far =
+  Matrix (2 / x) 0 0 (- (left + right) / x) 0 (2 / y) 0 (- (top + bottom) / y) 0 0 (-2 / z) (- (far + near) / z) 0 0 0 1
+  where
+    x = right - left
+    y = top - bottom
+    z = far - near
+
+-- | Camera look-at matrix (view matrix)
+matrixLookAt ::
+  -- | Camera position
+  Vector3 ->
+  -- | Camera target
+  Vector3 ->
+  -- | World up vector
+  Vector3 ->
+  Matrix
+matrixLookAt eye target up = Matrix xx xy xz (- vx |.| eye) yx yy yz (- vy |.| eye) zx zy zz (- vz |.| eye) 0 0 0 1
+  where
+    vz@(Vector3 zx zy zz) = vectorNormalize $ eye |-| target
+    vx@(Vector3 xx xy xz) = vectorNormalize $ vector3CrossProduct up vz
+    vy@(Vector3 yx yy yz) = vector3CrossProduct vz vx
+
+------------------------------------------------
+-- Quaternion math -----------------------------
+------------------------------------------------
+
+-- | Identity quaternion
+quaternionIdentity :: Quaternion
+quaternionIdentity = Vector4 0 0 0 1
+
+-- | Invert a quaternion
+quaternionInvert :: Quaternion -> Quaternion
+quaternionInvert q@(Vector4 x y z w) = Vector4 (- x * invLength) (- y * invLength) (- z * invLength) (w * invLength)
+  where
+    invLength = 1 / magnitudeSqr q
+
+-- | Multiply two quaternions
+quaternionMultiply :: Quaternion -> Quaternion -> Quaternion
+quaternionMultiply (Vector4 qax qay qaz qaw) (Vector4 qbx qby qbz qbw) = Vector4 (qax * qbw + qaw * qbx + qay * qbz - qaz * qby) (qay * qbw + qaw * qby + qaz * qbx - qax * qbz) (qaz * qbw + qaw * qbz + qax * qby - qay * qbx) (qaw * qbw - qax * qbx - qay * qby - qaz * qbz)
+
+-- | Normalize a quaternion (alias for 'vectorNormalize')
+quaternionNormalize :: Quaternion -> Quaternion
+quaternionNormalize = vectorNormalize
+
+-- | Lerp between two quaternions (alias for 'quaternionLerp')
+quaternionLerp ::
+  -- | Lerp start value
+  Quaternion ->
+  -- | Lerp end value
+  Quaternion ->
+  -- | Lerp amount
+  Float ->
+  Quaternion
+quaternionLerp = vectorLerp
+
+-- | Slerp-optimized interpolation between two quaternions
+quaternionNLerp ::
+  -- | Lerp start value
+  Quaternion ->
+  -- | Lerp end value
+  Quaternion ->
+  -- | Lerp amount
+  Float ->
+  Quaternion
+quaternionNLerp q1 q2 amount = vectorNormalize $ quaternionLerp q1 q2 amount
+
+-- | Spherical linear interpolation between two quaternions
+quaternionSLerp ::
+  -- | Lerp start value
+  Quaternion ->
+  -- | Lerp end value
+  Quaternion ->
+  -- | Lerp amount
+  Float ->
+  Quaternion
+quaternionSLerp q1 q2 amount
+  | cosHalfTheta >= 1 = q1
+  | cosHalfTheta > 0.95 = quaternionNLerp q1 q2' amount
+  | abs sinHalfTheta < 0.001 = (q1 |+| q2') |/ 2
+  | otherwise = (q1 |* ratioA) |+| (q2 |* ratioB)
+  where
+    cosHalfTheta = if dot < 0 then - dot else dot
+    sinHalfTheta = sqrt (1 - cosHalfTheta * cosHalfTheta)
+    halfTheta = acos cosHalfTheta
+
+    ratioA = sin ((1 - amount) * halfTheta) / sinHalfTheta
+    ratioB = sin (amount * halfTheta) / sinHalfTheta
+
+    q2' = if dot < 0 then additiveInverse q2 else q2
+    dot = q1 |.| q2
+
+-- | Quaternion based on the rotation between two vectors
+quaternionFromVector3ToVector3 :: Vector3 -> Vector3 -> Quaternion
+quaternionFromVector3ToVector3 from to = quaternionNormalize (Vector4 x y z (1 + cos2Theta))
+  where
+    cos2Theta = from |.| to
+    (Vector3 x y z) = vector3CrossProduct from to
+
+-- | Create a quaternion from a rotation matrix
+quaternionFromMatrix :: Matrix -> Quaternion
+quaternionFromMatrix mat
+  | fourBiggestSquaredMinus1 == fourWSquaredMinus1 = Vector4 ((matrix'm6 mat - matrix'm9 mat) * mult) ((matrix'm8 mat - matrix'm2 mat) * mult) ((matrix'm1 mat - matrix'm4 mat) * mult) biggestVal
+  | fourBiggestSquaredMinus1 == fourXSquaredMinus1 = Vector4 biggestVal ((matrix'm1 mat + matrix'm4 mat) * mult) ((matrix'm8 mat + matrix'm2 mat) * mult) ((matrix'm6 mat - matrix'm9 mat) * mult)
+  | fourBiggestSquaredMinus1 == fourYSquaredMinus1 = Vector4 ((matrix'm1 mat + matrix'm4 mat) * mult) biggestVal ((matrix'm6 mat + matrix'm9 mat) * mult) ((matrix'm8 mat - matrix'm2 mat) * mult)
+  | fourBiggestSquaredMinus1 == fourZSquaredMinus1 = Vector4 ((matrix'm8 mat + matrix'm2 mat) * mult) ((matrix'm6 mat + matrix'm9 mat) * mult) biggestVal ((matrix'm1 mat - matrix'm4 mat) * mult)
+  | otherwise = error "(quaternionFromMatrix) This error should never happen"
+  where
+    fourWSquaredMinus1 = matrix'm0 mat + matrix'm5 mat + matrix'm10 mat
+    fourXSquaredMinus1 = matrix'm0 mat - matrix'm5 mat - matrix'm10 mat
+    fourYSquaredMinus1 = matrix'm5 mat - matrix'm0 mat - matrix'm10 mat
+    fourZSquaredMinus1 = matrix'm10 mat - matrix'm0 mat - matrix'm5 mat
+    fourBiggestSquaredMinus1 =
+      maximum
+        [ fourWSquaredMinus1,
+          fourXSquaredMinus1,
+          fourYSquaredMinus1,
+          fourZSquaredMinus1
+        ]
+    biggestVal = sqrt (fourBiggestSquaredMinus1 + 1) / 2
+    mult = 0.5 / biggestVal
+
+-- | Create a rotation matrix from a quaternion
+quaternionToMatrix :: Quaternion -> Matrix
+quaternionToMatrix (Vector4 x y z w) = Matrix (1 - 2 * (b2 + c2)) (2 * (ab - cd)) (2 * (ac + bd)) 0 (2 * (ab + cd)) (1 - 2 * (a2 + c2)) (2 * (bc - ad)) 0 (2 * (ac - bd)) (2 * (bc + ad)) (1 - 2 * (a2 + b2)) 0 0 0 0 1
+  where
+    a2 = x * x
+    b2 = y * y
+    c2 = z * z
+    ac = x * z
+    ab = x * y
+    bc = y * z
+    ad = w * x
+    bd = w * y
+    cd = w * z
+
+-- | Create a quaternion for an angle and axis
+quaternionFromAxisAngle ::
+  -- | Rotation axis
+  Vector3 ->
+  -- | Angle in radians
+  Float ->
+  Quaternion
+quaternionFromAxisAngle axis angle = quaternionNormalize $ Vector4 (x * s) (y * s) (z * s) c -- I'm pretty sure normalizing the quaternion here is redundant
+  where
+    (Vector3 x y z) = vectorNormalize axis
+    s = sin (angle / 2)
+    c = cos (angle / 2)
+
+-- | Convert a quaternion to axis-angle representation
+quaternionToAxisAngle :: Quaternion -> (Vector3, Float)
+quaternionToAxisAngle q = (axis, angle)
+  where
+    (Vector4 x y z w) = quaternionNormalize q
+    s = sqrt (1 - w * w)
+    axis = if w > 0.9999 then Vector3 1 0 0 else Vector3 (x / s) (y / s) (z / s)
+    angle = acos w * 2
+
+-- | Create a quaternion from Euler angles (ZYX rotation order, angles should be in radians)
+quaternionFromEuler ::
+  -- | Pitch
+  Float ->
+  -- | Yaw
+  Float ->
+  -- | Roll
+  Float ->
+  Quaternion
+quaternionFromEuler pitch yaw roll = Vector4 (x1 * y0 * z0 - x0 * y1 * z1) (x0 * y1 * z0 + x1 * y0 * z1) (x0 * y0 * z1 - x1 * y1 * z0) (x0 * y0 * z0 + x1 * y1 * z1)
+  where
+    x0 = cos (pitch * 0.5)
+    x1 = sin (pitch * 0.5)
+    y0 = cos (yaw * 0.5)
+    y1 = sin (yaw * 0.5)
+    z0 = cos (roll * 0.5)
+    z1 = sin (roll * 0.5)
+
+-- | Convert a quaternion to Euler angle representation (Vector3 roll pitch yaw, all angles in radians)
+quaternionToEuler :: Quaternion -> Vector3
+quaternionToEuler (Vector4 x y z w) = Vector3 (atan2 x0 x1) (asin y0) (atan2 z0 z1)
+  where
+    x0 = 2 * (w * x + y * z)
+    x1 = 1 - 2 * (x * x + y * y)
+    y0 = clamp (2 * (w * y - z * x)) (-1) 1
+    z0 = 2 * (w * z + x * y)
+    z1 = 1 - 2 * (y * y + z * z)
+
+-- | Transform a quaternion given a transformation matrix
+quaternionTransform :: Quaternion -> Matrix -> Quaternion
+quaternionTransform (Vector4 x y z w) mat =
+  Vector4
+    (matrix'm0 mat * x + matrix'm4 mat * y + matrix'm8 mat * z + matrix'm12 mat * w)
+    (matrix'm1 mat * x + matrix'm5 mat * y + matrix'm9 mat * z + matrix'm13 mat * w)
+    (matrix'm2 mat * x + matrix'm6 mat * y + matrix'm10 mat * z + matrix'm14 mat * w)
+    (matrix'm3 mat * x + matrix'm7 mat * y + matrix'm11 mat * z + matrix'm15 mat * w)