[gnumeric] EIGEN: Detect underflow of pivot element.

[Morten Welinder <mortenw src gnome org>]

commit 194e8b76466ec03058aca443072717939b65eac1
Author: Morten Welinder <terra gnome org>
Date:   Mon Jan 7 21:06:19 2013 -0500

    EIGEN: Detect underflow of pivot element.

 ChangeLog      |    5 +++++
 src/mathfunc.c |   24 +++++++++++++++---------
 2 files changed, 20 insertions(+), 9 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 6c1d080..723fb81 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-01-07  Morten Welinder  <terra gnome org>
+
+	* src/mathfunc.c (gnm_matrix_eigen): Detect underflow of the
+	pivot.  Use gnm_hypot where possible.
+
 2013-01-04  Morten Welinder  <terra gnome org>
 
 	* src/gui-clipboard.c (gnm_x_claim_clipboard): Take a GdkDisplay
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 5dc59ef..4df86aa 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -6528,11 +6528,13 @@ static void
 gnm_matrix_eigen_update (guint k, gnm_float t, gnm_float *eigenvalues, gboolean *changed, guint *state)
 {
 	gnm_float y = eigenvalues[k];
+	gboolean unchanged;
 	eigenvalues[k] += t;
-	if (changed[k] && y == eigenvalues[k]) {
+	unchanged = (y == eigenvalues[k]);
+	if (changed[k] && unchanged) {
 		changed[k] = FALSE;
 		(*state)--;
-	} else if ((!changed[k]) && (y != eigenvalues[k])) {
+	} else if (!changed[k] && !unchanged) {
 		changed[k] = TRUE;
 		(*state)++;
 	}
@@ -6577,19 +6579,23 @@ gnm_matrix_eigen (gnm_float **matrix, gnm_float **eigenvectors, gnm_float *eigen
 			g_print ("gnm_matrix_eigen exceeded iterations\n");
 			return FALSE;
 		}
-		for (k = 1; k < (usize-1); k++)
+		for (k = 1; k < usize - 1; k++)
 			if (gnm_abs (matrix[k][ind[k]]) > gnm_abs (matrix[m][ind[m]]))
 				m = k;
 		l = ind[m];
 		pivot = matrix[m][l];
 		/* pivot is (m,l) */
+		if (pivot == 0) {
+			g_printerr ("gnm_matrix_eigen underflow in pivot.\n");
+			break;
+		}
 
-		y = (eigenvalues[l] - eigenvalues[m])/2;
-		t = gnm_abs (y) + gnm_sqrt (pivot*pivot+y*y);
-		s = gnm_sqrt (pivot*pivot+t*t);
-		c = t/s;
-		s = pivot/s;
-		t = pivot * pivot /t;
+		y = (eigenvalues[l] - eigenvalues[m]) / 2;
+		t = gnm_abs (y) + gnm_hypot (pivot, y);
+		s = gnm_hypot (pivot, t);
+		c = t / s;
+		s = pivot / s;
+		t = pivot * pivot / t;
 		if (y < 0) {
 			s = -s;
 			t = -t;