diff --git a/Makefile b/Makefile
index 6ae8c96..40e962d 100644
--- a/Makefile
+++ b/Makefile
@@ -5,7 +5,7 @@ CC = gcc
 
 PROG = libsweet.so
 
-OBJS = sweet_math.o sweet_matrix.o 
+OBJS = sweet_math.o sweet_matrix.o sweet_matrix_stack.o
 
 CFLAGS = -pedantic -static -fPIC
 LD = -shared 
diff --git a/sweet.h b/sweet.h
index e56c25f..876f767 100644
--- a/sweet.h
+++ b/sweet.h
@@ -1,5 +1,5 @@
 /*
- *   Sweet is a small library for basic math and small matrix operation.
+ *   Sweet is a small library for basic math and small matrix operations.
  *   Copyright 2014 Luc Girod.
  *
  * This library is free software: you can redistribute it and/or
@@ -24,6 +24,7 @@
 #include "sweet_types.h"
 #include "sweet_math.h"
 #include "sweet_matrix.h"
+#include "sweet_matrix_stack.h"
 
 #endif
 
diff --git a/sweet_math.c b/sweet_math.c
index f0f86fd..8fe96fb 100644
--- a/sweet_math.c
+++ b/sweet_math.c
@@ -1,5 +1,5 @@
 /*
- *   Sweet is a small library for basic math and small matrix operation.
+ *   Sweet is a small library for basic math and small matrix operations.
  *   Copyright 2014 Luc Girod.
  *
  * This library is free software: you can redistribute it and/or
diff --git a/sweet_math.h b/sweet_math.h
index 97b3ff6..8e46ec5 100644
--- a/sweet_math.h
+++ b/sweet_math.h
@@ -1,5 +1,5 @@
 /*
- *   Sweet is a small library for basic math and small matrix operation.
+ *   Sweet is a small library for basic math and small matrix operations.
  *   Copyright 2014 Luc Girod.
  *
  * This library is free software: you can redistribute it and/or
diff --git a/sweet_matrix.c b/sweet_matrix.c
index e6e91d7..06d293f 100644
--- a/sweet_matrix.c
+++ b/sweet_matrix.c
@@ -1,5 +1,5 @@
 /*
- *   Sweet is a small library for basic math and small matrix operation.
+ *   Sweet is a small library for basic math and small matrix operations.
  *   Copyright 2014 Luc Girod.
  *
  * This library is free software: you can redistribute it and/or
@@ -536,44 +536,54 @@ sweet_matrix_texture_bias ()
 void
 sweet_matrix_transpose2 (mat2 * mt, mat2 * m)
 {
+    float mv1 = m->v[1];
     mt->v[0] = m->v[0];
     mt->v[1] = m->v[2];
-    mt->v[2] = m->v[1];
+    mt->v[2] = mv1;
     mt->v[3] = m->v[3];
 }
 
 void
 sweet_matrix_transpose3 (mat3 * mt, mat3 * m)
 {
+    float mv1 = m->v[1];
+    float mv2 = m->v[2];
+    float mv5 = m->v[5];
     mt->v[0] = m->v[0];
     mt->v[1] = m->v[3];
     mt->v[2] = m->v[6];
-    mt->v[3] = m->v[1];
+    mt->v[3] = mv1;
     mt->v[4] = m->v[4];
     mt->v[5] = m->v[7];
-    mt->v[6] = m->v[2];
-    mt->v[7] = m->v[5];
+    mt->v[6] = mv2;
+    mt->v[7] = mv5;
     mt->v[8] = m->v[8];
 }
 
 void
 sweet_matrix_transpose4 (mat4 * mt, mat4 * m)
 {
+    float mv2 = m->v[2];
+    float mv3 = m->v[3];
+    float mv6 = m->v[6];
+    float mv7 = m->v[7];
+    float mv11 = m->v[11];
+
     mt->v[0] = m->v[0];
     mt->v[1] = m->v[4];
     mt->v[2] = m->v[8];
     mt->v[3] = m->v[12];
-    mt->v[4] = m->v[2];
+    mt->v[4] = mv2;
     mt->v[5] = m->v[5];
     mt->v[6] = m->v[9];
     mt->v[7] = m->v[13];
-    mt->v[8] = m->v[2];
-    mt->v[9] = m->v[6];
+    mt->v[8] = mv2;
+    mt->v[9] = mv6;
     mt->v[10] = m->v[10];
     mt->v[11] = m->v[14];
-    mt->v[12] = m->v[3];
-    mt->v[13] = m->v[7];
-    mt->v[14] = m->v[11];
+    mt->v[12] = mv3;
+    mt->v[13] = mv7;
+    mt->v[14] = mv11;
     mt->v[15] = m->v[15];
 }
 
@@ -619,35 +629,55 @@ sweet_matrix_inverse2 (mat2 * inv, mat2 * m)
 {
     float invdet = 1.0 / sweet_matrix_det2 (m);
 
+    float mv0 = m->v[0];
     inv->v[0] = m->v[3] * invdet;
     inv->v[1] = m->v[1] * (-invdet);
     inv->v[2] = m->v[2] * (-invdet);
-    inv->v[3] = m->v[0] * invdet;
+    inv->v[3] = mv0 * invdet;
 }
 
 void
 sweet_matrix_inverse3 (mat3 * inv, mat3 * m)
 {
+    float mv0 = m->v[0];
+    float mv1 = m->v[1];
+    float mv2 = m->v[2];
+    float mv3 = m->v[3];
+    float mv4 = m->v[4];
+    float mv6 = m->v[6];
     float invdet;
 
     invdet = 1.0 / sweet_matrix_det3 (m);
 
     inv->v[0] = ((m->v[4] * m->v[8]) - (m->v[5] * m->v[7])) * invdet;
     inv->v[1] = (-(m->v[1] * m->v[8]) + (m->v[2] * m->v[7])) * invdet;
-    inv->v[2] = ((m->v[1] * m->v[5]) - (m->v[2] * m->v[4])) * invdet;
+    inv->v[2] = ((mv1 * m->v[5]) - (m->v[2] * m->v[4])) * invdet;
 
     inv->v[3] = (-(m->v[3] * m->v[8]) + (m->v[5] * m->v[6])) * invdet;
-    inv->v[4] = ((m->v[0] * m->v[8]) - (m->v[2] * m->v[6])) * invdet;
-    inv->v[5] = (-(m->v[0] * m->v[5]) + (m->v[2] * m->v[3])) * invdet;
+    inv->v[4] = ((mv0 * m->v[8]) - (mv2 * m->v[6])) * invdet;
+    inv->v[5] = (-(mv0 * m->v[5]) + (mv2 * mv3)) * invdet;
 
-    inv->v[6] = ((m->v[3] * m->v[7]) - (m->v[4] * m->v[6])) * invdet;
-    inv->v[7] = (-(m->v[0] * m->v[7]) + (m->v[1] * m->v[6])) * invdet;
-    inv->v[8] = ((m->v[0] * m->v[4]) - (m->v[1] * m->v[3])) * invdet;
+    inv->v[6] = ((mv3 * m->v[7]) - (mv4 * m->v[6])) * invdet;
+    inv->v[7] = (-(mv0 * m->v[7]) + (mv1 * mv6)) * invdet;
+    inv->v[8] = ((mv0 * mv4) - (mv1 * mv3)) * invdet;
 }
 
 static void
 matrix_inverse4 (mat4 * inv, mat4 * m)
 {
+    float mv0 = m->v[0];
+    float mv1 = m->v[1];
+    float mv2 = m->v[2];
+    float mv3 = m->v[3];
+    float mv4 = m->v[4];
+    float mv5 = m->v[5];
+    float mv6 = m->v[6];
+    float mv7 = m->v[7];
+    float mv8 = m->v[8];
+    float mv9 = m->v[9];
+    float mv10 = m->v[10];
+    float mv12 = m->v[12];
+    float mv13 = m->v[13];
     float invdet = 1.0 / sweet_matrix_det4 (m);
 
     inv->v[0] = m->v[5] * (m->v[10] * m->v[15] - m->v[11] * m->v[14]) -
@@ -658,61 +688,61 @@ matrix_inverse4 (mat4 * inv, mat4 * m)
               m->v[2] * (m->v[9] * m->v[15] - m->v[11] * m->v[13]) -
               m->v[3] * (m->v[9] * m->v[14] - m->v[10] * m->v[13]);
 
-    inv->v[2] = m->v[1] * (m->v[6] * m->v[15] - m->v[7] * m->v[14]) -
+    inv->v[2] = mv1 * (m->v[6] * m->v[15] - m->v[7] * m->v[14]) -
               m->v[2] * (m->v[5] * m->v[15] - m->v[7] * m->v[13]) +
               m->v[3] * (m->v[5] * m->v[14] - m->v[6] * m->v[13]);
 
-    inv->v[3] = -m->v[1] * (m->v[6] * m->v[11] - m->v[7] * m->v[10]) +
-              m->v[2] * (m->v[5] * m->v[11] - m->v[7] * m->v[9]) -
+    inv->v[3] = -mv1 * (m->v[6] * m->v[11] - m->v[7] * m->v[10]) +
+              mv2 * (m->v[5] * m->v[11] - m->v[7] * m->v[9]) -
               m->v[3] * (m->v[5] * m->v[10] - m->v[6] * m->v[9]);
 
     inv->v[4] = -m->v[4] * (m->v[10] * m->v[15] - m->v[11] * m->v[14]) +
               m->v[6] * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) -
               m->v[7] * (m->v[8] * m->v[14] - m->v[10] * m->v[12]);
 
-    inv->v[5] = m->v[0] * (m->v[10] * m->v[15] - m->v[11] * m->v[14]) -
-              m->v[2] * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) +
-              m->v[3] * (m->v[8] * m->v[14] - m->v[10] * m->v[12]);
+    inv->v[5] = mv0 * (m->v[10] * m->v[15] - m->v[11] * m->v[14]) -
+              mv2 * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) +
+              mv3 * (m->v[8] * m->v[14] - m->v[10] * m->v[12]);
 
-    inv->v[6] = -m->v[0] * (m->v[6] * m->v[15] - m->v[7] * m->v[14]) +
-              m->v[2] * (m->v[4] * m->v[15] - m->v[7] * m->v[12]) -
-              m->v[3] * (m->v[4] * m->v[14] - m->v[6] * m->v[12]);
+    inv->v[6] = -mv0 * (m->v[6] * m->v[15] - m->v[7] * m->v[14]) +
+              mv2 * (mv4 * m->v[15] - m->v[7] * m->v[12]) -
+              mv3 * (mv4 * m->v[14] - m->v[6] * m->v[12]);
 
-    inv->v[7] = m->v[0] * (m->v[6] * m->v[11] - m->v[7] * m->v[10]) -
-              m->v[2] * (m->v[4] * m->v[11] - m->v[7] * m->v[8]) +
-              m->v[3] * (m->v[4] * m->v[10] - m->v[6] * m->v[8]);
+    inv->v[7] = mv0 * (mv6 * m->v[11] - m->v[7] * m->v[10]) -
+              mv2 * (mv4 * m->v[11] - m->v[7] * m->v[8]) +
+              mv3 * (mv4 * m->v[10] - mv6 * m->v[8]);
 
-    inv->v[8] = m->v[4] * (m->v[9] * m->v[15] - m->v[11] * m->v[13]) -
-              m->v[5] * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) +
-              m->v[7] * (m->v[8] * m->v[13] - m->v[9] * m->v[12]);
+    inv->v[8] = mv4 * (m->v[9] * m->v[15] - m->v[11] * m->v[13]) -
+              mv5 * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) +
+              mv7 * (m->v[8] * m->v[13] - m->v[9] * m->v[12]);
 
-    inv->v[9] = -m->v[0] * (m->v[9] * m->v[15] - m->v[11] * m->v[13]) +
-              m->v[1] * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) -
-              m->v[3] * (m->v[8] * m->v[13] - m->v[9] * m->v[12]);
+    inv->v[9] = -mv0 * (m->v[9] * m->v[15] - m->v[11] * m->v[13]) +
+              mv1 * (m->v[8] * m->v[15] - m->v[11] * m->v[12]) -
+              mv3 * (m->v[8] * m->v[13] - m->v[9] * m->v[12]);
 
-    inv->v[10] = m->v[0] * (m->v[5] * m->v[15] - m->v[7] * m->v[13]) -
-               m->v[1] * (m->v[4] * m->v[15] - m->v[7] * m->v[12]) +
-               m->v[3] * (m->v[4] * m->v[13] - m->v[5] * m->v[12]);
+    inv->v[10] = mv0 * (mv5 * m->v[15] - mv7 * m->v[13]) -
+               mv1 * (mv4 * m->v[15] - mv7 * m->v[12]) +
+               mv3 * (mv4 * m->v[13] - mv5 * m->v[12]);
 
-    inv->v[11] = -m->v[0] * (m->v[5] * m->v[11] - m->v[7] * m->v[9]) +
-               m->v[1] * (m->v[4] * m->v[11] - m->v[7] * m->v[8]) -
-               m->v[3] * (m->v[4] * m->v[9] - m->v[5] * m->v[8]);
+    inv->v[11] = -mv0 * (mv5 * m->v[11] - mv7 * mv9) +
+               mv1 * (mv4 * m->v[11] - mv7 * mv8) -
+               mv3 * (mv4 * mv9 - mv5 * mv8);
 
-    inv->v[12] = -m->v[4] * (m->v[9] * m->v[14] - m->v[10] * m->v[13]) +
-               m->v[5] * (m->v[8] * m->v[14] - m->v[10] * m->v[12]) -
-               m->v[6] * (m->v[8] * m->v[13] - m->v[9] * m->v[12]);
+    inv->v[12] = -mv4 * (mv9 * m->v[14] - mv10 * m->v[13]) +
+               mv5 * (mv8 * m->v[14] - mv10 * m->v[12]) -
+               mv6 * (mv8 * m->v[13] - mv9 * m->v[12]);
 
-    inv->v[13] = m->v[0] * (m->v[9] * m->v[14] - m->v[10] * m->v[13]) -
-               m->v[1] * (m->v[8] * m->v[14] - m->v[10] * m->v[12]) +
-               m->v[2] * (m->v[8] * m->v[13] - m->v[9] * m->v[12]);
+    inv->v[13] = mv0 * (mv9 * m->v[14] - mv10 * m->v[13]) -
+               mv1 * (mv8 * m->v[14] - mv10 * mv12) +
+               mv2 * (mv8 * m->v[13] - mv9 * mv12);
 
-    inv->v[14] = -m->v[0] * (m->v[5] * m->v[14] - m->v[6] * m->v[13]) +
-               m->v[1] * (m->v[4] * m->v[14] - m->v[6] * m->v[12]) -
-               m->v[2] * (m->v[4] * m->v[13] - m->v[5] * m->v[12]);
+    inv->v[14] = -mv0 * (mv5 * m->v[14] - mv6 * mv13) +
+               mv1 * (mv4 * m->v[14] - mv6 * mv12) -
+               mv2 * (mv4 * mv13 - mv5 * mv12);
 
-    inv->v[15] = m->v[0] * (m->v[5] * m->v[10] - m->v[6] * m->v[9]) -
-               m->v[1] * (m->v[4] * m->v[10] - m->v[6] * m->v[8]) +
-               m->v[2] * (m->v[4] * m->v[9] - m->v[5] * m->v[8]);
+    inv->v[15] = mv0 * (mv5 * mv10 - mv6 * mv9) -
+               mv1 * (mv4 * mv10 - mv6 * mv8) +
+               mv2 * (mv4 * mv9 - mv5 * mv8);
 
     inv->v[ 0] *= invdet;
     inv->v[ 1] *= invdet;
@@ -748,7 +778,7 @@ sweet_matrix_inverse4 (mat4 * inv, mat4 * m)
     if (sweet_math_approx_equals (detp, 0, 0.0000001))
     {
         /* In this case the standard algorithm is used */
-        matrix_inverse4 (m, inv);
+        matrix_inverse4 (inv, m);
         return;
     }
 
@@ -807,50 +837,95 @@ sweet_matrix_inverse4 (mat4 * inv, mat4 * m)
 void
 sweet_matrix_product2 (mat2 * p, mat2 * a, mat2 * b)
 {
+    float av0 = a->v[0];
+    float av1 = a->v[1];
+    float bv0 = b->v[0];
+    float bv2 = b->v[2];
     p->v[0] = a->v[0] * b->v[0] + a->v[2] * b->v[1];
-    p->v[1] = a->v[1] * b->v[0] + a->v[3] * b->v[1];
-    p->v[2] = a->v[0] * b->v[2] + a->v[2] * b->v[3];
-    p->v[3] = a->v[1] * b->v[2] + a->v[3] * b->v[3];
+    p->v[1] = a->v[1] * bv0 + a->v[3] * b->v[1];
+    p->v[2] = av0 * b->v[2] + a->v[2] * b->v[3];
+    p->v[3] = av1 * bv2 + a->v[3] * b->v[3];
 }
 
 void
 sweet_matrix_product3 (mat3 * p, mat3 * a, mat3 * b)
 {
+    float bv0 = b->v[0];
+    float bv1 = b->v[1];
+    float bv3 = b->v[3];
+    float bv4 = b->v[4];
+    float bv6 = b->v[6];
+    float bv7 = b->v[7];
+    float av0 = a->v[0];
+    float av1 = a->v[1];
+    float av2 = a->v[2];
+    float av3 = a->v[3];
+    float av4 = a->v[4];
+    float av5 = a->v[5];
+
     p->v[0] = a->v[0] * b->v[0] + a->v[3] * b->v[1] + a->v[6] * b->v[2];
-    p->v[1] = a->v[1] * b->v[0] + a->v[4] * b->v[1] + a->v[7] * b->v[2];
-    p->v[2] = a->v[2] * b->v[0] + a->v[5] * b->v[1] + a->v[8] * b->v[2];
+    p->v[1] = a->v[1] * bv0 + a->v[4] * b->v[1] + a->v[7] * b->v[2];
+    p->v[2] = a->v[2] * bv0 + a->v[5] * bv1 + a->v[8] * b->v[2];
 
-    p->v[3] = a->v[0] * b->v[3] + a->v[3] * b->v[4] + a->v[6] * b->v[5];
-    p->v[4] = a->v[1] * b->v[3] + a->v[4] * b->v[4] + a->v[7] * b->v[5];
-    p->v[5] = a->v[2] * b->v[3] + a->v[5] * b->v[4] + a->v[8] * b->v[5];
+    p->v[3] = av0 * b->v[3] + a->v[3] * b->v[4] + a->v[6] * b->v[5];
+    p->v[4] = av1 * bv3 + a->v[4] * b->v[4] + a->v[7] * b->v[5];
+    p->v[5] = av2 * bv3 + a->v[5] * bv4 + a->v[8] * b->v[5];
 
-    p->v[6] = a->v[0] * b->v[6] + a->v[3] * b->v[7] + a->v[6] * b->v[8];
-    p->v[7] = a->v[1] * b->v[6] + a->v[4] * b->v[7] + a->v[7] * b->v[8];
-    p->v[8] = a->v[2] * b->v[6] + a->v[5] * b->v[7] + a->v[8] * b->v[8];
+    p->v[6] = av0 * b->v[6] + av3 * b->v[7] + a->v[6] * b->v[8];
+    p->v[7] = av1 * bv6 + av4 * b->v[7] + a->v[7] * b->v[8];
+    p->v[8] = av2 * bv6 + av5 * bv7 + a->v[8] * b->v[8];
 }
 
 void
 sweet_matrix_product4 (mat4 * p, mat4 * a, mat4 * b)
 {
+    float bv0 = b->v[0];
+    float bv1 = b->v[1];
+    float bv2 = b->v[2];
+    float bv3 = b->v[3];
+    float bv4 = b->v[4];
+    float bv5 = b->v[5];
+    float bv6 = b->v[6];
+    float bv7 = b->v[7];
+    float bv8 = b->v[8];
+    float bv9 = b->v[9];
+    float bv10 = b->v[10];
+    float bv11 = b->v[11];
+    float bv12 = b->v[12];
+    float bv13 = b->v[13];
+    float bv14 = b->v[14];
+    float av0 = a->v[0];
+    float av1 = a->v[1];
+    float av2 = a->v[2];
+    float av3 = a->v[3];
+    float av4 = a->v[4];
+    float av5 = a->v[5];
+    float av6 = a->v[6];
+    float av7 = a->v[7];
+    float av8 = a->v[8];
+    float av9 = a->v[9];
+    float av10 = a->v[10];
+    float av11 = a->v[11];
+
     p->v[0] = a->v[0] * b->v[0] + a->v[4] * b->v[1] + a->v[8] * b->v[2] + a->v[12] * b->v[3];
-    p->v[1] = a->v[1] * b->v[0] + a->v[5] * b->v[1] + a->v[9] * b->v[2] + a->v[13] * b->v[3];
-    p->v[2] = a->v[2] * b->v[0] + a->v[6] * b->v[1] + a->v[10] * b->v[2] + a->v[14] * b->v[3];
-    p->v[3] = a->v[3] * b->v[0] + a->v[7] * b->v[1] + a->v[11] * b->v[2] + a->v[15] * b->v[3];
+    p->v[1] = a->v[1] * bv0 + a->v[5] * b->v[1] + a->v[9] * b->v[2] + a->v[13] * b->v[3];
+    p->v[2] = a->v[2] * bv0 + a->v[6] * bv1 + a->v[10] * b->v[2] + a->v[14] * b->v[3];
+    p->v[3] = a->v[3] * bv0 + a->v[7] * bv1 + a->v[11] * bv2 + a->v[15] * b->v[3];
 
-    p->v[4] = a->v[0] * b->v[4] + a->v[4] * b->v[5] + a->v[8] * b->v[6] + a->v[12] * b->v[7];
-    p->v[5] = a->v[1] * b->v[4] + a->v[5] * b->v[5] + a->v[9] * b->v[6] + a->v[13] * b->v[7];
-    p->v[6] = a->v[2] * b->v[4] + a->v[6] * b->v[5] + a->v[10] * b->v[6] + a->v[14] * b->v[7];
-    p->v[7] = a->v[3] * b->v[4] + a->v[7] * b->v[5] + a->v[11] * b->v[6] + a->v[15] * b->v[7];
+    p->v[4] = av0 * b->v[4] + a->v[4] * b->v[5] + a->v[8] * b->v[6] + a->v[12] * b->v[7];
+    p->v[5] = av1 * b->v[4] + a->v[5] * b->v[5] + a->v[9] * b->v[6] + a->v[13] * b->v[7];
+    p->v[6] = av2 * bv4 + a->v[6] * bv5 + a->v[10] * b->v[6] + a->v[14] * b->v[7];
+    p->v[7] = av3 * bv4 + a->v[7] * bv5 + a->v[11] * bv6 + a->v[15] * b->v[7];
 
-    p->v[8] = a->v[0] * b->v[8] + a->v[4] * b->v[9] + a->v[8] * b->v[10] + a->v[12] * b->v[11];
-    p->v[9] = a->v[1] * b->v[8] + a->v[5] * b->v[9] + a->v[9] * b->v[10] + a->v[13] * b->v[11];
-    p->v[10] = a->v[2] * b->v[8] + a->v[6] * b->v[9] + a->v[10] * b->v[10] + a->v[14] * b->v[11];
-    p->v[11] = a->v[3] * b->v[8] + a->v[7] * b->v[9] + a->v[11] * b->v[10] + a->v[15] * b->v[11];
+    p->v[8] = av0 * b->v[8] + av4 * b->v[9] + a->v[8] * b->v[10] + a->v[12] * b->v[11];
+    p->v[9] = av1 * bv8 + av5 * b->v[9] + a->v[9] * b->v[10] + a->v[13] * b->v[11];
+    p->v[10] = av2 * bv8 + av6 * bv9 + a->v[10] * b->v[10] + a->v[14] * b->v[11];
+    p->v[11] = av3 * bv8 + av7 * bv9 + a->v[11] * bv10 + a->v[15] * b->v[11];
 
-    p->v[12] = a->v[0] * b->v[12] + a->v[4] * b->v[13] + a->v[8] * b->v[14] + a->v[12] * b->v[15];
-    p->v[13] = a->v[1] * b->v[12] + a->v[5] * b->v[13] + a->v[9] * b->v[14] + a->v[13] * b->v[15];
-    p->v[14] = a->v[2] * b->v[12] + a->v[6] * b->v[13] + a->v[10] * b->v[14] + a->v[14] * b->v[15];
-    p->v[15] = a->v[3] * b->v[12] + a->v[7] * b->v[13] + a->v[11] * b->v[14] + a->v[15] * b->v[15];
+    p->v[12] = av0 * b->v[12] + av4 * b->v[13] + av8 * b->v[14] + a->v[12] * b->v[15];
+    p->v[13] = av1 * bv12 + av5 * b->v[13] + av9 * b->v[14] + a->v[13] * b->v[15];
+    p->v[14] = av2 * bv12 + av6 * bv13 + av10 * b->v[14] + a->v[14] * b->v[15];
+    p->v[15] = av3 * bv12 + av7 * bv13 + av11 * bv14 + a->v[15] * b->v[15];
 }
 
 void
diff --git a/sweet_matrix.h b/sweet_matrix.h
index e79a9b0..110088e 100644
--- a/sweet_matrix.h
+++ b/sweet_matrix.h
@@ -1,5 +1,5 @@
 /*
- *   Sweet is a small library for basic math and small matrix operation.
+ *   Sweet is a small library for basic math and small matrix operations.
  *   Copyright 2014 Luc Girod.
  *
  * This library is free software: you can redistribute it and/or
@@ -45,8 +45,8 @@ mat4 sweet_matrix_identity4 ();
 /** @param zFar as float */
 /** @return Perspective matrix 4x4 */
 mat4 sweet_matrix_frustum (float left, float right,
-                          float bottom, float top,
-                          float zNear, float zFar);
+                           float bottom, float top,
+                           float zNear, float zFar);
 
 /** Perspective */
 /** @param fovy of the camera in degree as flaot */
@@ -65,8 +65,8 @@ mat4 sweet_matrix_perspective (double fovy, double aspect, double z_near, double
 /** @param zFar as float */
 /** @return Orthogonal matrix 4x4 */
 mat4 sweet_matrix_ortho (float left, float right,
-                        float bottom, float top,
-                        float zNear, float zFar);
+                         float bottom, float top,
+                         float zNear, float zFar);
 
 /** @return camera matrix */
 mat4 sweet_matrix_look_at (vec3 pos, vec3 dir, vec3 up);
diff --git a/sweet_matrix_stack.c b/sweet_matrix_stack.c
new file mode 100644
index 0000000..b2b1dbb
--- /dev/null
+++ b/sweet_matrix_stack.c
@@ -0,0 +1,199 @@
+/*
+ *   Sweet is a small library for basic math and small matrix operations.
+ *   Copyright 2014 Luc Girod.
+ *
+ * This library is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation, either
+ * version 3 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <stdlib.h>
+#include "sweet_matrix_stack.h"
+
+matrix_stack2
+sweet_matrix_stack2_new ()
+{
+    matrix_stack2 m;
+    m.position = 0;
+    m.matrix[0] = sweet_matrix_identity2 ();
+    return m; 
+}
+
+matrix_stack3
+sweet_matrix_stack3_new ()
+{
+    matrix_stack3 m;
+    m.position = 0;
+    m.matrix[0] = sweet_matrix_identity3 ();
+    return m; 
+}
+
+matrix_stack4
+sweet_matrix_stack4_new ()
+{
+    matrix_stack4 m;
+    m.position = 0;
+    m.matrix[0] = sweet_matrix_identity4 ();
+    return m; 
+}
+
+int
+sweet_matrix_stack2_push (matrix_stack2 * ms)
+{
+    if (ms->position >= SWEET_MATRIX_STACK_SIZE) { return 1; }
+
+    ms->matrix[ms->position + 1] = ms->matrix[ms->position];
+    ms->position++;
+    return 0;
+}
+
+int
+sweet_matrix_stack3_push (matrix_stack3 * ms)
+{
+    if (ms->position >= SWEET_MATRIX_STACK_SIZE) { return 1; }
+
+    ms->matrix[ms->position + 1] = ms->matrix[ms->position];
+    ms->position++;
+    return 0;
+}
+
+int
+sweet_matrix_stack4_push (matrix_stack4 * ms)
+{
+    if (ms->position >= SWEET_MATRIX_STACK_SIZE) { return 1; }
+
+    ms->matrix[ms->position + 1] = ms->matrix[ms->position];
+    ms->position = ms->position++;
+    return 0;
+}
+
+
+int
+sweet_matrix_stack2_pop (matrix_stack2 * ms)
+{
+    if (ms->position <= 0) { return 1; }
+
+    ms->position--;
+    return 0;
+}
+
+int
+sweet_matrix_stack3_pop (matrix_stack3 * ms)
+{
+    if (ms->position <= 0) { return 1; }
+
+    ms->position--;
+    return 0;
+}
+
+int
+sweet_matrix_stack4_pop (matrix_stack4 * ms)
+{
+    if (ms->position <= 0) { return 1; }
+
+    ms->position--;
+    return 0;
+}
+
+void
+sweet_matrix_stack2_mult (matrix_stack2 * ms, mat2 * m)
+{
+    sweet_matrix_product2 (ms->matrix + ms->position,
+                           ms->matrix + ms->position, m);
+}
+
+void
+sweet_matrix_stack3_mult (matrix_stack3 * ms, mat3 * m)
+{
+    sweet_matrix_product3 (ms->matrix + ms->position,
+                           ms->matrix + ms->position, m);
+}
+
+void
+sweet_matrix_stack4_mult (matrix_stack4 * ms, mat4 * m)
+{
+    sweet_matrix_product4 (ms->matrix + ms->position,
+                           ms->matrix + ms->position, m);
+}
+
+void
+sweet_matrix_stack2_identity (matrix_stack2 * ms)
+{
+    ms->matrix[ms->position] = sweet_matrix_identity2 ();
+}
+
+void
+sweet_matrix_stack3_identity (matrix_stack3 * ms)
+{
+    ms->matrix[ms->position] = sweet_matrix_identity3 ();
+}
+
+void
+sweet_matrix_stack4_identity (matrix_stack4 * ms)
+{
+    ms->matrix[ms->position] = sweet_matrix_identity4 ();
+}
+
+void
+sweet_matrix_stack2_set (matrix_stack2 * ms, mat2 * m)
+{
+    ms->matrix[ms->position] = *m;
+}
+
+void
+sweet_matrix_stack3_set (matrix_stack3 * ms, mat3 * m)
+{
+    ms->matrix[ms->position] = *m;
+}
+
+void
+sweet_matrix_stack4_set (matrix_stack4 * ms, mat4 * m)
+{
+    ms->matrix[ms->position] = *m;
+}
+
+mat2 *
+sweet_matrix_stack2_get_matrix_pointer (matrix_stack2 * ms)
+{
+    return (ms->matrix + ms->position);
+}
+
+mat3 *
+sweet_matrix_stack3_get_matrix_pointer (matrix_stack3 * ms)
+{
+    return (ms->matrix + ms->position);
+}
+
+mat4 *
+sweet_matrix_stack4_get_matrix_pointer (matrix_stack4 * ms)
+{
+    return (ms->matrix + ms->position);
+}
+    
+mat2
+sweet_matrix_stack2_get_matrix (matrix_stack2 * ms)
+{
+    return ms->matrix[ms->position];
+}
+
+mat3
+sweet_matrix_stack3_get_matrix (matrix_stack3 * ms)
+{
+    return ms->matrix[ms->position];
+}
+
+mat4
+sweet_matrix_stack4_get_matrix (matrix_stack4 * ms)
+{
+    return ms->matrix[ms->position];
+}
+
diff --git a/sweet_matrix_stack.h b/sweet_matrix_stack.h
new file mode 100644
index 0000000..b4ae0c3
--- /dev/null
+++ b/sweet_matrix_stack.h
@@ -0,0 +1,92 @@
+/*
+ *   Sweet is a small library for basic math and small matrix operations.
+ *   Copyright 2014 Luc Girod.
+ *
+ * This library is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation, either
+ * version 3 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SWEET_MATRIX_STACK
+#define SWEET_MATRIX_STACK
+
+#include "sweet_types.h"
+#include "sweet_matrix.h"
+
+#define SWEET_MATRIX_STACK_SIZE 128
+
+typedef struct matrix_stack2
+{
+    mat2 matrix[SWEET_MATRIX_STACK_SIZE];
+    unsigned int position;
+
+}matrix_stack2;
+
+typedef struct matrix_stack3
+{
+    mat3 matrix[SWEET_MATRIX_STACK_SIZE];
+    unsigned int position;
+
+}matrix_stack3;
+
+typedef struct matrix_stack4
+{
+    mat4 matrix[SWEET_MATRIX_STACK_SIZE];
+    unsigned int position;
+
+}matrix_stack4;
+
+
+matrix_stack2 sweet_matrix_stack2_new ();
+
+matrix_stack3 sweet_matrix_stack3_new ();
+
+matrix_stack4 sweet_matrix_stack4_new ();
+
+int sweet_matrix_stack2_push (matrix_stack2 * ms);
+
+int sweet_matrix_stack3_push (matrix_stack3 * ms);
+
+int sweet_matrix_stack4_push (matrix_stack4 * ms);
+
+int sweet_matrix_stack2_pop (matrix_stack2 * ms);
+
+int sweet_matrix_stack3_pop (matrix_stack3 * ms);
+
+int sweet_matrix_stack4_pop (matrix_stack4 * ms);
+
+void sweet_matrix_stack2_mult (matrix_stack2 * ms, mat2 * m);
+
+void sweet_matrix_stack3_mult (matrix_stack3 * ms, mat3 * m);
+
+void sweet_matrix_stack4_mult (matrix_stack4 * ms, mat4 * m);
+
+void sweet_matrix_stack2_set (matrix_stack2 * ms, mat2 * m);
+
+void sweet_matrix_stack3_set (matrix_stack3 * ms, mat3 * m);
+
+void sweet_matrix_stack4_set (matrix_stack4 * ms, mat4 * m);
+
+mat2 * sweet_matrix_stack2_get_matrix_pointer (matrix_stack2 * m);
+
+mat3 * sweet_matrix_stack3_get_matrix_pointer (matrix_stack3 * m);
+
+mat4 * sweet_matrix_stack4_get_matrix_pointer (matrix_stack4 * m);
+    
+mat2 sweet_matrix_stack2_get_matrix (matrix_stack2 * m);
+
+mat3 sweet_matrix_stack3_get_matrix (matrix_stack3 * m);
+
+mat4 sweet_matrix_stack4_get_matrix (matrix_stack4 * m);
+
+#endif
+
diff --git a/sweet_types.h b/sweet_types.h
index 1854150..a4d3c9f 100644
--- a/sweet_types.h
+++ b/sweet_types.h
@@ -1,5 +1,5 @@
 /*
- *   Sweet is a small library for basic math and small matrix operation.
+ *   Sweet is a small library for basic math and small matrix operations.
  *   Copyright 2014 Luc Girod.
  *
  * This library is free software: you can redistribute it and/or