From 3ecdb295b5bdc8441f647174513bb670f9c939c9 Mon Sep 17 00:00:00 2001
From: Luc <luc@gluc.ch>
Date: Sun, 19 Aug 2018 04:56:27 +0200
Subject: [PATCH] fix in eigen vector

---
 sweet_matrix.c | 102 ++++++++++++++++++++++++++-----------------------
 1 file changed, 55 insertions(+), 47 deletions(-)

diff --git a/sweet_matrix.c b/sweet_matrix.c
index 2991fa5..174348a 100644
--- a/sweet_matrix.c
+++ b/sweet_matrix.c
@@ -1254,66 +1254,74 @@ sweet_matrix_solve3 (mat3 * m, vec3 * v)
     return 0;
 }
 
-static int is_identity (mat3 * m)
+static int is_diagonal (mat3 * m)
 {
-    if (!sweet_math_approx_equals (m->v[0], 1.0, EPSILON)) { return 0; }
     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[4], 1.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; }
-    if (!sweet_math_approx_equals (m->v[8], 1.0, EPSILON)) { return 0; }
     return 1;
 }
 
-static void get_eigen_vector (vec3 * v, mat3 * a, float h)
+static int col_n_zero (mat3 * m, int i)
 {
-    char free_value;
-    a->v[0] = a->v[0] - h;
-    a->v[4] = a->v[4] - h;
-    a->v[8] = a->v[8] - h;
+   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 int row_n_zero (mat3 * m, int i)
+{
+   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));
+}
+
+static void get_eigen_vector (vec3 * v, mat3 a, float h)
+{
+    a.v[0] = a.v[0] - h;
+    a.v[4] = a.v[4] - h;
+    a.v[8] = a.v[8] - h;
 
     v->x = 0;
     v->y = 0;
     v->z = 0;
-    free_value = 0;
 
-    sweet_matrix_gaussian_elimination3 (a, v);
+    sweet_matrix_gaussian_elimination3 (&a, v);
 
-    if (sweet_math_approx_equals (a->v[0], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[1], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[2], 0, EPSILON))
-    {
-        free_value = 1;
+    if (col_n_zero (&a, 0))
         v->x = 1;
-    }
-    if (sweet_math_approx_equals (a->v[3], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[4], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[5], 0, EPSILON))
-    {
-        free_value = 1;
-        v->y = 1;
-    }
-    if (sweet_math_approx_equals (a->v[6], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[7], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[8], 0, EPSILON))
-    {
-        free_value = 1;
+    else if (col_n_zero(&a, 1))
+        v->y = 1; 
+    else if (col_n_zero(&a, 2))
         v->z = 1;
-    }
-
-    if (!free_value &&
-        sweet_math_approx_equals (a->v[2], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[5], 0, EPSILON) &&
-        sweet_math_approx_equals (a->v[8], 0, EPSILON))
+    else
     {
-       v->z = 1;
-       v->y = - a->v[7];
-       v->x = - a->v[3] * v->y - a->v[6];
+        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]);
+        }
     }
     *v = sweet_vector_normalize3 (*v);
 }
@@ -1326,7 +1334,7 @@ sweet_matrix_eigen_vector3 (vec3 * vectors, mat3 * m)
     int n;
     float h1, h2, h3;
 
-    if (is_identity (m))
+    if (is_diagonal (m))
     {
         vectors[0] = sweet_vector_new3 (1,0,0);
         vectors[1] = sweet_vector_new3 (0,1,0);
@@ -1337,16 +1345,16 @@ sweet_matrix_eigen_vector3 (vec3 * vectors, mat3 * m)
     n = sweet_matrix_eigen_value3 (&h1, &h2, &h3, m);
 
     /* vector associated with h1 */
-    a = *m;
-    get_eigen_vector (&vectors[0], &a, h1);
+    get_eigen_vector (&vectors[0], *m, h1);
 
     /* vector associated with h2 */
-    a = *m;
-    get_eigen_vector (&vectors[1], &a, h2);
+    get_eigen_vector (&vectors[1], *m, h2);
 
     /* vector associated with h3 */
-    a = *m;
-    get_eigen_vector (&vectors[2], &a, h3);
+    if (n > 2)
+        get_eigen_vector (&vectors[2], *m, h3);
+    else
+        vectors[2] = sweet_vector_cross (vectors[0], vectors[1]);
     return n;
 }