diff --git a/Makefile b/Makefile
index a5a6140..fac5263 100644
--- a/Makefile
+++ b/Makefile
@@ -10,15 +10,9 @@ OBJS_WITH_SHORT = sweet_math_short.o sweet_matrix_short.o
 CFLAGS = -pedantic -fPIC -Wall
 LD = -shared
 
-ifeq ($(TARGET),short)
 all: $(PROG)
 $(PROG): $(OBJS) $(OBJS_WITH_SHORT)
 	$(CC) $(LD) $(OBJS) $(OBJS_WITH_SHORT) -o $(PROG)
-else
-all: $(PROG)
-$(PROG): $(OBJS)
-	$(CC) $(LD) $(OBJS) -o $(PROG)
-endif
 
 .c.o:
 	$(CC) $(CFLAGS) -c $*.c
diff --git a/sweet_math.c b/sweet_math.c
index a75487b..869b324 100644
--- a/sweet_math.c
+++ b/sweet_math.c
@@ -17,6 +17,7 @@
  */
 
 #include <stdio.h>
+#include <stdlib.h>
 #include <math.h>
 #include <stdarg.h>
 #include "sweet_math.h"
@@ -78,10 +79,10 @@ sweet_math_quadratic_polynomial (vec2 * r, float a, float b, float c)
 
 
 #define FRAC_1_ON_3   0.3333333333333333333333333333333333333333333333333333333333333333333
-#define FRAC_2_ON_3   0.6666666666666666666666666666666666666666666666666666666666666666666
-#define FRAC_1_ON_6   0.1666666666666666666666666666666666666666666666666666666666666666666
+#define FRAC_2_ON_3   0.6666666666666666666666666666666666666666666666666666666666666666667
+#define FRAC_1_ON_6   0.1666666666666666666666666666666666666666666666666666666666666666667
 #define FRAC_1_ON_9   0.1111111111111111111111111111111111111111111111111111111111111111111
-#define FRAC_2_ON_27  0.0740740740740740740740740740740740740740740740740740740740740740740
+#define FRAC_2_ON_27  0.0740740740740740740740740740740740740740740740740740740740740740741
 #define FRAC_1_ON_27  0.0370370370370370370370370370370370370370370370370370370370370370370
 
 int
@@ -117,6 +118,8 @@ sweet_math_cubic_polynomial (vec3 * roots, float t, float a, float b, float c)
     if (D < 0)
     {
         roots->x = (r + s) - a_over_3;
+        roots->y = -r - a_over_3;
+        roots->z = roots->y;
         return 1;
     }
     else if (sweet_math_approx_equals (D, 0, EPSILON))
@@ -486,6 +489,33 @@ sweet_vector_rescale4 (vec4 v, float size)
     return v;
 }
 
+vec2
+sweet_vector_inverse2 (vec2 v)
+{
+    v.x = - v.x;
+    v.y = - v.y;
+    return v;
+}
+
+vec3
+sweet_vector_inverse3 (vec3 v)
+{
+    v.x = - v.x;
+    v.y = - v.y;
+    v.z = - v.z;
+    return v;
+}
+
+vec4
+sweet_vector_inverse4 (vec4 v)
+{
+    v.x = - v.x;
+    v.y = - v.y;
+    v.z = - v.z;
+    v.w = - v.w;
+    return v;
+}
+
 vec2
 sweet_vector_add2 (vec2 v1, vec2 v2)
 {
@@ -598,6 +628,24 @@ sweet_vector_angle3 (vec3 a, vec3 b)
     return acos (dot);
 }
 
+vec3 
+sweet_vector_normalized_random ()
+{
+    vec3 v;
+    v.x = (float)rand() / (float)RAND_MAX;
+    v.y = (float)rand() / (float)RAND_MAX;
+    v.z = (float)rand() / (float)RAND_MAX;
+    return sweet_vector_normalize3 (v);
+}
+
+vec3
+sweet_vector_random_ortho(vec3 v)
+{
+    vec3 u = sweet_vector_normalized_random ();
+    v = sweet_vector_cross(u, v);
+    return sweet_vector_normalize3(v);
+}
+
 vec2
 sweet_vector_product2 (vec2 v1, vec2 v2)
 {
diff --git a/sweet_math.h b/sweet_math.h
index 1679de9..a3c76a0 100644
--- a/sweet_math.h
+++ b/sweet_math.h
@@ -185,6 +185,13 @@ vec2 sweet_vector_rescale2 (vec2 vector, float size);
 vec3 sweet_vector_rescale3 (vec3 vector, float size);
 vec4 sweet_vector_rescale4 (vec4 vector, float size);
 
+/** Inverse */
+/** @param Vector as vector*/
+/** @return - vector */
+vec2 sweet_vector_inverse2 (vec2 vector);
+vec3 sweet_vector_inverse3 (vec3 vector);
+vec4 sweet_vector_inverse4 (vec4 vector);
+
 /** Add */
 /** @param Vector1 as vector */
 /** @param Vector2 as vector */
@@ -216,6 +223,10 @@ vec4 sweet_vector_sub4 (vec4 vector1, vec4 vector2);
 float sweet_vector_angle2 (vec2 vector1, vec2 vector2);
 float sweet_vector_angle3 (vec3 vector1, vec3 vector2);
 
+/** Random */
+vec3 sweet_vector_normalized_random ();
+vec3 sweet_vector_random_ortho(vec3 v);
+
 /** Mult */
 /** @param Vector1 as vector */
 /** @param Vector2 as vector */
diff --git a/sweet_math_short.h b/sweet_math_short.h
index 745dfd7..e2dddfb 100644
--- a/sweet_math_short.h
+++ b/sweet_math_short.h
@@ -81,6 +81,9 @@ vec4 vec4_new_v4 (vec4 v);
 #define vec2_rescale sweet_vector_rescale2 
 #define vec3_rescale sweet_vector_rescale3 
 #define vec4_rescale sweet_vector_rescale4 
+#define vec2_inverse sweet_vector_inverse2 
+#define vec3_inverse sweet_vector_inverse3 
+#define vec4_inverse sweet_vector_inverse4 
 #define vec2_add sweet_vector_add2 
 #define vec3_add sweet_vector_add3 
 #define vec4_add sweet_vector_add4 
diff --git a/sweet_matrix.c b/sweet_matrix.c
index 174348a..4e3b761 100644
--- a/sweet_matrix.c
+++ b/sweet_matrix.c
@@ -16,8 +16,6 @@
  * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-
-
 #include <string.h>
 #include <math.h>
 #include "sweet_matrix.h"
@@ -86,7 +84,6 @@ sweet_matrix_null4 (void)
     return m;
 }
 
-
 mat2
 sweet_matrix_identity2 (void)
 {
@@ -149,13 +146,12 @@ sweet_matrix_identity4 (void)
     return m;
 }
 
-
 mat2 
 sweet_matrix_ortho_basis2 (vec2 u, vec2 v)
 {
     mat2 m;
     m.v[0] = u.x;
-    m.v[3] = v.x;
+    m.v[1] = v.x;
 
     m.v[2] = u.y;
     m.v[3] = v.y;
@@ -182,28 +178,55 @@ sweet_matrix_ortho_basis3 (vec3 u, vec3 v, vec3 w)
     return m;
 }
 
+mat4
+sweet_matrix_ortho_basis3h (vec3 p, vec3 u, vec3 v, vec3 w)
+{
+    mat4 m;
+    m.v[0] = u.x;
+    m.v[1] = v.x;
+    m.v[2] = w.x;
+    m.v[3] = 0;
+
+    m.v[4] = u.y;
+    m.v[5] = v.y;
+    m.v[6] = w.y;
+    m.v[7] = 0;
+
+    m.v[8] = u.z;
+    m.v[9] = v.z;
+    m.v[10] = w.z;
+    m.v[11] = 0;
+
+    m.v[12] = -p.x * u.x - p.y * u.y - p.z * u.z;
+    m.v[13] = -p.x * v.x - p.y * v.y - p.z * v.z;
+    m.v[14] = -p.x * w.x - p.y * w.y - p.z * w.z;
+    m.v[15] = 1;
+
+    return m;
+}
+
 mat4 
 sweet_matrix_ortho_basis4 (vec4 u, vec4 v, vec4 w, vec4 t)
 {
     mat4 m;
     m.v[0] = u.x;
-    m.v[1] = v.y;
-    m.v[2] = w.z;
-    m.v[3] = t.w;
+    m.v[1] = v.x;
+    m.v[2] = w.x;
+    m.v[3] = t.x;
 
-    m.v[4] = u.x;
+    m.v[4] = u.y;
     m.v[5] = v.y;
-    m.v[6] = w.z;
-    m.v[7] = t.w;
+    m.v[6] = w.y;
+    m.v[7] = t.y;
 
-    m.v[8]  = u.x;
-    m.v[9]  = v.y;
+    m.v[8]  = u.z;
+    m.v[9]  = v.z;
     m.v[10] = w.z;
-    m.v[11] = t.w;
+    m.v[11] = t.z;
 
-    m.v[12] = u.x;
-    m.v[13] = v.y;
-    m.v[14] = w.z;
+    m.v[12] = u.w;
+    m.v[13] = v.w;
+    m.v[14] = w.w;
     m.v[15] = t.w;
 
     return m;
@@ -331,6 +354,14 @@ sweet_matrix_look_at (vec3 pos, vec3 dir, vec3 up)
     return m;
 }
 
+mat4
+sweet_matrix_look_at_ortho (vec3 pos, vec3 dir, vec3 up, vec3 right)
+{
+    mat4 m;
+    dir = sweet_vector_scale3 (dir, -1);
+    return sweet_matrix_ortho_basis3h (pos, right, up, dir);
+}
+
 mat3
 sweet_matrix_quat_rotation3 (quaternion q)
 {
@@ -495,13 +526,13 @@ sweet_matrix_translation2h (vec2 v)
     m.v[1] = 0.0;
     m.v[2] = 0.0;
 
-    m.v[4] = 0.0;
-    m.v[5] = 1.0;
-    m.v[6] = 0.0;
+    m.v[3] = 0.0;
+    m.v[4] = 1.0;
+    m.v[5] = 0.0;
 
-    m.v[8] = v.x;
-    m.v[9] = v.y;
-    m.v[10] = 1.0;
+    m.v[6] = v.x;
+    m.v[7] = v.y;
+    m.v[8] = 1.0;
 
     return m;
 }
@@ -987,6 +1018,62 @@ sweet_matrix_product4 (mat4 * p, mat4 * a, mat4 * b)
     p->v[15] = av3 * bv12 + av7 * bv13 + av11 * bv14 + a->v[15] * b->v[15];
 }
 
+void
+sweet_matrix_add_diag3 (mat3 * out, mat3 * in, float x, float y, float z)
+{
+    out->v[0] = in->v[0] + x;
+    out->v[1] = in->v[1];
+    out->v[2] = in->v[2];
+    out->v[3] = in->v[3];
+    out->v[4] = in->v[4] + y;
+    out->v[5] = in->v[5];
+    out->v[6] = in->v[6];
+    out->v[6] = in->v[7];
+    out->v[8] = in->v[8] + z;
+}
+
+void
+sweet_matrix_substract_diag3 (mat3 * out, mat3 * in, float x, float y, float z)
+{
+    out->v[0] = in->v[0] - x;
+    out->v[1] = in->v[1];
+    out->v[2] = in->v[2];
+    out->v[3] = in->v[3];
+    out->v[4] = in->v[4] - y;
+    out->v[5] = in->v[5];
+    out->v[6] = in->v[6];
+    out->v[6] = in->v[7];
+    out->v[8] = in->v[8] - z;
+}
+
+void
+sweet_matrix_add3 (mat3 * out, mat3 * m1, mat3 * m2)
+{
+    out->v[0] = m1->v[0] + m2->v[0];
+    out->v[1] = m1->v[1] + m2->v[1];
+    out->v[2] = m1->v[2] + m2->v[2];
+    out->v[3] = m1->v[3] + m2->v[3];
+    out->v[4] = m1->v[4] + m2->v[4];
+    out->v[5] = m1->v[5] + m2->v[5];
+    out->v[6] = m1->v[6] + m2->v[6];
+    out->v[7] = m1->v[7] + m2->v[7];
+    out->v[8] = m1->v[8] + m2->v[8];
+}
+
+void
+sweet_matrix_substract3 (mat3 * out, mat3 * m1, mat3 * m2)
+{
+    out->v[0] = m1->v[0] - m2->v[0];
+    out->v[1] = m1->v[1] - m2->v[1];
+    out->v[2] = m1->v[2] - m2->v[2];
+    out->v[3] = m1->v[3] - m2->v[3];
+    out->v[4] = m1->v[4] - m2->v[4];
+    out->v[5] = m1->v[5] - m2->v[5];
+    out->v[6] = m1->v[6] - m2->v[6];
+    out->v[7] = m1->v[7] - m2->v[7];
+    out->v[8] = m1->v[8] - m2->v[8];
+}
+
 void
 sweet_matrix_sub3 (mat2 * sub, mat3 * m)
 {
@@ -1089,18 +1176,31 @@ sweet_matrix_eigen_value3 (float * h1, float * h2, float * h3, mat3 * m)
     float b;
     float c;
     int n;
-    vec3 v;
+    vec3 v = sweet_vector_new3 (0,0,0);
 
     a = m->v[0] + m->v[4] + m->v[8];
 
     b = m->v[1] * m->v[3] + m->v[2] * m->v[6] + m->v[7] * m->v[5] -
         m->v[0] * m->v[4] - m->v[0] * m->v[8] - m->v[4] * m->v[8];
 
-    c = m->v[0] * m->v[4] * m->v[8] - m->v[0] * m->v[7] * m->v[5] -
-        m->v[1] * m->v[3] * m->v[8] + m->v[1] * m->v[6] * m->v[5] +
-        m->v[2] * m->v[3] * m->v[7] - m->v[2] * m->v[6] * m->v[4];
+    c = m->v[0] * (m->v[4] * m->v[8] - m->v[5] * m->v[7]) -
+        m->v[1] * (m->v[3] * m->v[8] - m->v[5] * m->v[6]) +
+        m->v[2] * (m->v[3] * m->v[7] - m->v[4] * m->v[6]);
 
     n = sweet_math_cubic_polynomial (&v, -1, a, b, c);
+    if (v.y > v.x)
+    {
+        float t = v.x; v.x = v.y; v.y = t;
+    }
+    if (v.z > v.y)
+    {
+        float t = v.y; v.y = v.z; v.z = t;
+    } 
+    if (v.y > v.x)
+    {
+        float t = v.x; v.x = v.y; v.y = t;
+    }
+
     *h1 = v.x;
     *h2 = v.y;
     *h3 = v.z;
@@ -1205,7 +1305,7 @@ static void substract_row (mat3 * a, float * s, int l1, int l2, int i)
 
 }
 
-#define EPSILON 0.00005
+#define EPSILON 1e-05f 
 
 int
 sweet_matrix_gaussian_elimination3 (mat3 * a, vec3 * v)
@@ -1254,107 +1354,178 @@ sweet_matrix_solve3 (mat3 * m, vec3 * v)
     return 0;
 }
 
-static int is_diagonal (mat3 * m)
+static void build_Q0 (mat3 * q, mat3 * m)
 {
-    if (!sweet_math_approx_equals (m->v[1], 0.0, EPSILON)) { return 0; }
-    if (!sweet_math_approx_equals (m->v[2], 0.0, EPSILON)) { return 0; }
-    if (!sweet_math_approx_equals (m->v[3], 0.0, EPSILON)) { return 0; }
-    if (!sweet_math_approx_equals (m->v[5], 0.0, EPSILON)) { return 0; }
-    if (!sweet_math_approx_equals (m->v[6], 0.0, EPSILON)) { return 0; }
-    if (!sweet_math_approx_equals (m->v[7], 0.0, EPSILON)) { return 0; }
-    return 1;
+    float a = sqrt (m->v[0] * m->v[0] + m->v[1] * m->v[1] + m->v[2] * m->v[2]);
+    if (m->v[0] > 0) { a = - a; }
+
+    vec3 v = sweet_vector_new3 (m->v[0] - a, m->v[1], m->v[2]);
+    float norm = sweet_vector_norm3 (v);
+    if (norm > 0)
+        v = sweet_vector_scale3 (v, 1.0 / norm);
+
+    q->v[0] = 1 - 2 * v.x * v.x;
+    q->v[1] =   - 2 * v.y * v.x;
+    q->v[2] =   - 2 * v.z * v.x;
+
+    q->v[3] =   - 2 * v.x * v.y;
+    q->v[4] = 1 - 2 * v.y * v.y;
+    q->v[5] =   - 2 * v.z * v.y;
+
+    q->v[6] =   - 2 * v.x * v.z;
+    q->v[7] =   - 2 * v.y * v.z;
+    q->v[8] = 1 - 2 * v.z * v.z;
 }
 
-static int col_n_zero (mat3 * m, int i)
-{
-   i = i * 3;
-   return (sweet_math_approx_equals (m->v[i], 0, EPSILON) &&
-           sweet_math_approx_equals (m->v[i+1], 0, EPSILON) &&
-           sweet_math_approx_equals (m->v[i+2], 0, EPSILON));
+static void build_Q1 (mat3 * q, mat3 * m)
+{ 
+    float a = sqrt (m->v[4] * m->v[4] + m->v[5] * m->v[5]);
+    if (m->v[4] > 0) { a = -a; }
+
+    vec2 v = sweet_vector_new2 (m->v[4] - a, m->v[5]);
+    float norm = sweet_vector_norm2 (v);
+    if (norm > 0)
+        v = sweet_vector_scale2 (v, 1.0 / norm);
+
+    q->v[0] = 1;
+    q->v[1] = 0;
+    q->v[2] = 0;
+    q->v[3] = 0;
+    q->v[4] = 1 - 2 * v.x * v.x;
+    q->v[5] =   - 2 * v.y * v.x;
+    q->v[6] = 0;
+    q->v[7] =   - 2 * v.x * v.y;
+    q->v[8] = 1 - 2 * v.y * v.y;
 }
 
-static int row_n_zero (mat3 * m, int i)
+void sweet_matrix_QR3 (mat3 * Q, mat3 * R, mat3 * m)
 {
-   return (sweet_math_approx_equals (m->v[i], 0, EPSILON) &&
-           sweet_math_approx_equals (m->v[i+3], 0, EPSILON) &&
-           sweet_math_approx_equals (m->v[i+6], 0, EPSILON));
+    mat3 tmp;
+    mat3 Q0;
+    mat3 Q1;
+
+    build_Q0 (&Q0, m);
+    sweet_matrix_product3 (&tmp, &Q0, m);
+    build_Q1 (&Q1, &tmp);
+    sweet_matrix_product3 (Q, &Q0, &Q1);
+
+    sweet_matrix_transpose3 (&tmp, Q);
+    sweet_matrix_product3 (R, &tmp, m);
 }
 
-static void get_eigen_vector (vec3 * v, mat3 a, float h)
+void sweet_matrix_hess3 (mat3 * H, mat3 * Q, mat3 * m)
 {
-    a.v[0] = a.v[0] - h;
-    a.v[4] = a.v[4] - h;
-    a.v[8] = a.v[8] - h;
+    mat3 tmp;
 
-    v->x = 0;
-    v->y = 0;
-    v->z = 0;
+    float a = sqrt (m->v[1] * m->v[1] + m->v[2] * m->v[2]);
+    if (m->v[1] > 0)
+        a = - a;
 
-    sweet_matrix_gaussian_elimination3 (&a, v);
+    vec2 v = sweet_vector_new2 (m->v[1] - a, m->v[2]);
+    float norm = sweet_vector_norm2 (v);
+    if (norm > 0)
+        v = sweet_vector_scale2 (v, 1.0/norm);
 
-    if (col_n_zero (&a, 0))
-        v->x = 1;
-    else if (col_n_zero(&a, 1))
-        v->y = 1; 
-    else if (col_n_zero(&a, 2))
-        v->z = 1;
-    else
+    Q->v[0] = 1;
+    Q->v[1] = 0;
+    Q->v[2] = 0;
+    Q->v[3] = 0;
+    Q->v[4] = 1 - 2 * v.x * v.x;
+    Q->v[5] =   - 2 * v.x * v.y;
+    Q->v[6] = 0;
+    Q->v[7] =   - 2 * v.x * v.y;
+    Q->v[8] = 1 - 2 * v.y * v.y;
+
+    sweet_matrix_product3 (&tmp, Q, m);
+    sweet_matrix_product3 (H, &tmp, Q);
+}
+
+static void shift (mat3 * m, float s)
+{
+    m->v[0] = m->v[0] + s;
+    m->v[4] = m->v[4] + s;
+    m->v[8] = m->v[8] + s;
+}
+
+
+void print_matrix (char * title, mat3 * m)
+{
+    printf ("%s:\n[%.10f      %.10f      %.10f]\n[%.10f      %.10f      %.10f]\n[%.10f      %.10f      %.10f]\n",
+            title,
+            m->v[0], m->v[3], m->v[6],
+            m->v[1], m->v[4], m->v[7],
+            m->v[2], m->v[5], m->v[8]);
+}
+
+static void sort3 (float * val, vec3 * vec)
+{
+    if (val[1] > val[0])
     {
-        if (sweet_math_approx_equals (a.v[2], 0, EPSILON) &&
-            sweet_math_approx_equals (a.v[5], 0, EPSILON))
-        {
-            if (sweet_math_approx_equals (a.v[8], 0, EPSILON))
-            {
-                v->z = 1;
-                v->y = - a.v[7];
-                v->x = (-v->y * a.v[3] - a.v[6]);
-            }
-            else
-            {
-                v->z = 0;
-                v->y = 1;
-                v->x = (-a.v[3]);
-            }
-        } 
-        else if (row_n_zero (&a, 1))
-        {
-            v->y = 1;
-            v->z = 0;
-            v->x = (-v->y * a.v[3]);
-        }
+        float t = val[0];
+        vec3 v = vec[0];
+        val[0] = val[1];
+        val[1] = t;
+        vec[0] = vec[1];
+        vec[1] = v;
     }
-    *v = sweet_vector_normalize3 (*v);
+    if (val[2] > val[1])
+    {
+        float t = val[1];
+        vec3 v = vec[1];
+        val[1] = val[2];
+        val[2] = t;
+        vec[1] = vec[2];
+        vec[2] = v;    
+    }
+    /*
+    if (val[1] > val[0])
+    { 
+        float t = val[0];
+        vec3 v = vec[0];
+        val[0] = val[1];
+        val[1] = t;
+        vec[0] = vec[1];
+        vec[1] = v;
+    }*/
 }
 
-int
-sweet_matrix_eigen_vector3 (vec3 * vectors, mat3 * m)
+void sweet_matrix_eigen_vectors3 (vec3 * eigen_vectors, float * eigen_values, mat3 * A, int it)
 {
-    /* (A - hI)V = 0 */
-    mat3 a;
-    int n;
-    float h1, h2, h3;
+    int i;
+    mat3 H, S, Q, W, R;
+    sweet_matrix_hess3 (&H, &W, A);
 
-    if (is_diagonal (m))
+    float s = H.v[8];
+    shift (&H, -H.v[8]);
+
+    sweet_matrix_QR3 (&Q, &R, &H);
+    sweet_matrix_product3 (&H, &R, &Q);
+    S = Q;
+    for (i = 0; i < it; i++)
     {
-        vectors[0] = sweet_vector_new3 (1,0,0);
-        vectors[1] = sweet_vector_new3 (0,1,0);
-        vectors[2] = sweet_vector_new3 (0,0,1);
-        return 3;
+        sweet_matrix_QR3 (&Q, &R, &H);
+        sweet_matrix_product3 (&H, &R, &Q);
+        sweet_matrix_product3 (&S, &S, &Q);
     }
 
-    n = sweet_matrix_eigen_value3 (&h1, &h2, &h3, m);
+    shift (&H, s);
+    sweet_matrix_product3 (&S, &W, &S);
 
-    /* vector associated with h1 */
-    get_eigen_vector (&vectors[0], *m, h1);
+    eigen_values[0] = H.v[0];
+    eigen_values[1] = H.v[4];
+    eigen_values[2] = H.v[8];
+    eigen_vectors[0] = sweet_vector_new3 (S.v[0], S.v[1], S.v[2]);
+    eigen_vectors[1] = sweet_vector_new3 (S.v[3], S.v[4], S.v[5]);
+    eigen_vectors[2] = sweet_vector_new3 (S.v[6], S.v[7], S.v[8]);
+    sort3 (eigen_values, eigen_vectors);
 
-    /* vector associated with h2 */
-    get_eigen_vector (&vectors[1], *m, h2);
+/*
+    printf ("%f %f %f\n", eigen_vectors[0].x, eigen_vectors[0].y, eigen_vectors[0].z);
+    printf ("%f %f %f\n", eigen_vectors[1].x, eigen_vectors[1].y, eigen_vectors[1].z);
+    printf ("%f %f %f\n", eigen_vectors[2].x, eigen_vectors[2].y, eigen_vectors[2].z);
 
-    /* vector associated with h3 */
-    if (n > 2)
-        get_eigen_vector (&vectors[2], *m, h3);
-    else
-        vectors[2] = sweet_vector_cross (vectors[0], vectors[1]);
-    return n;
+    printf ("%f %f %f\n", eigen_values[0], eigen_values[1], eigen_values[2]);
+    */
 }
 
+
diff --git a/sweet_matrix.h b/sweet_matrix.h
index 79d4165..2e69a8e 100644
--- a/sweet_matrix.h
+++ b/sweet_matrix.h
@@ -41,6 +41,7 @@ mat4 sweet_matrix_identity4 (void);
 /** @return transform matrix to the new basis */
 mat2 sweet_matrix_ortho_basis2 (vec2 u, vec2 v);
 mat3 sweet_matrix_ortho_basis3 (vec3 u, vec3 v, vec3 w);
+mat4 sweet_matrix_ortho_basis3h (vec3 p, vec3 u, vec3 v, vec3 w);
 mat4 sweet_matrix_ortho_basis4 (vec4 u, vec4 v, vec4 w, vec4 t);
 
 /** Frustum */
@@ -78,6 +79,8 @@ mat4 sweet_matrix_ortho (float left, float right,
 /** @return camera matrix */
 mat4 sweet_matrix_look_at (vec3 pos, vec3 dir, vec3 up);
 
+mat4 sweet_matrix_look_at_ortho (vec3 pos, vec3 dir, vec3 up, vec3 right);
+
 /** Rotation from quaterion */
 /** @param quaternion representing rotation */
 /** @return 3d Rotation matrix */
@@ -170,7 +173,8 @@ int sweet_matrix_gaussian_elimination3 (mat3 * m, vec3 * v);
 int sweet_matrix_solve3 (mat3 * a, vec3 * v);
 
 /** Eigen Vector */
-int sweet_matrix_eigen_vector3 (vec3 * vectors, mat3 * m);
+void sweet_matrix_eigen_vectors3 (vec3 * eigen_vectors, float * eigen_values, mat3 * A, int it);
+
 
 #endif
 
diff --git a/sweet_matrix_short.h b/sweet_matrix_short.h
index 9bbf6ac..ee2126f 100644
--- a/sweet_matrix_short.h
+++ b/sweet_matrix_short.h
@@ -51,11 +51,13 @@ mat4 mat4_new_4v (vec4 e1, vec4 e2, vec4 e3, vec4 e4);
 #define identity4 sweet_matrix_identity4 
 #define ortho_basis2 sweet_matrix_ortho_basis2
 #define ortho_basis3 sweet_matrix_ortho_basis3
+#define ortho_basis3h sweet_matrix_ortho_basis3h
 #define ortho_basis4 sweet_matrix_ortho_basis4
 #define matrix_frustum sweet_matrix_frustum 
 #define matrix_perspective sweet_matrix_perspective 
 #define matrix_ortho sweet_matrix_ortho 
 #define matrix_look_at sweet_matrix_look_at 
+#define matrix_look_at_ortho sweet_matrix_look_at_ortho
 #define matrix_quat_rotation3 sweet_matrix_quat_rotation3 
 #define matrix_quat_rotation3h sweet_matrix_quat_rotation3h 
 #define rotation2 sweet_matrix_rotation2 
@@ -91,7 +93,7 @@ mat4 mat4_new_4v (vec4 e1, vec4 e2, vec4 e3, vec4 e4);
 #define mat3_eigen_value sweet_matrix_eigen_value3 
 #define mat3_gaussian_elimination sweet_matrix_gaussian_elimination3 
 #define mat3_solve sweet_matrix_solve3 
-#define mat3_eigen_vector sweet_matrix_eigen_vector3 
+#define mat3_eigen_vector sweet_matrix_eigen_vectors3 
 
 #endif