From: Richard Kreckel <kreckel@ginac.de>
Date: Sun, 28 Jan 2018 23:07:17 +0000 (+0100)
Subject: Make matrix::solve() work with non-normal zeros.
X-Git-Tag: release_1-7-3~6
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=d0ff428fb5b7bb565a0aea69e05e5705d840c16d;hp=afb0ccaa49a0cca001d854594e09125a58434123

Make matrix::solve() work with non-normal zeros.

Normalize elements of augmented matrix before checking if they are zero.
This ensures that we don't divide by a non-normal zero. Also added a test.

This bug was reported by Vitaly Magerya <vmagerya@gmail.com>.
---

diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp
index 586c52f5..7770e6bf 100644
--- a/check/exam_matrices.cpp
+++ b/check/exam_matrices.cpp
@@ -1,3 +1,4 @@
+
 /** @file exam_matrices.cpp
  *
  *  Here we examine manipulations on GiNaC's symbolic matrices. */
@@ -275,6 +276,35 @@ static unsigned matrix_rank()
 	return result;	
 }
 
+unsigned matrix_solve_nonnormal()
+{
+	symbol a("a"), b("b"), c("c"), x("x");
+	// This matrix has a non-normal zero element!
+	matrix mx {{1,0,0},
+	           {0,1/(x+1)-(x-1)/(x*x-1),1},
+	           {0,0,0}};
+	matrix zero {{0}, {0}, {0}};
+	matrix vars {{a}, {b}, {c}};
+	try {
+		matrix sol_gauss   = mx.solve(vars, zero, solve_algo::gauss);
+		matrix sol_divfree = mx.solve(vars, zero, solve_algo::divfree);
+		matrix sol_bareiss = mx.solve(vars, zero, solve_algo::bareiss);
+		if (sol_gauss != sol_divfree  ||  sol_gauss != sol_bareiss) {
+			clog << "different solutions while solving "
+			     << mx << " * " << vars << " == " << zero << endl
+			     << "gauss:   " << sol_gauss << endl
+			     << "divfree: " << sol_divfree << endl
+			     << "bareiss: " << sol_bareiss << endl;
+			return 1;
+		}
+	} catch (const exception & e) {
+		clog << "exception thrown while solving "
+		     << mx << " * " << vars << " == " << zero << endl;
+		return 1;
+	}
+	return 0;
+}
+
 static unsigned matrix_misc()
 {
 	unsigned result = 0;
@@ -337,6 +367,7 @@ unsigned exam_matrices()
 	result += matrix_solve2();  cout << '.' << flush;
 	result += matrix_evalm();  cout << "." << flush;
 	result += matrix_rank();  cout << "." << flush;
+	result += matrix_solve_nonnormal();  cout << "." << flush;
 	result += matrix_misc();  cout << '.' << flush;
 	
 	return result;
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 80ca3e45..bf7f8347 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1054,11 +1054,11 @@ matrix matrix::solve(const matrix & vars,
 		unsigned last_assigned_sol = n+1;
 		for (int r=m-1; r>=0; --r) {
 			unsigned fnz = 1;    // first non-zero in row
-			while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].is_zero()))
+			while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].normal().is_zero()))
 				++fnz;
 			if (fnz>n) {
 				// row consists only of zeros, corresponding rhs must be 0, too
-				if (!aug.m[r*(n+p)+n+co].is_zero()) {
+				if (!aug.m[r*(n+p)+n+co].normal().is_zero()) {
 					throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
 				}
 			} else {