[GiNaC-list] [patch] Fix a failure in matrix::division_free_elimination()

[Vitaly Magerya]

Hi, folks. I've run into another case of matrix::solve() failing.
This time it seems like the division-free elimination doesn't
handle un-normal zeros properly: currently it only does expand()
(instead of normal()) on the results of it's calculations, but
then uses matrix::pivot() to find the pivot, and that routine
also uses expand().is_zero() for zero testing. This breaks when
fractions are present in the initial matrix. Not every time, but
at least in one example.

Note that this problem also affects matrix::solve() with 'algo'
set to solve_algo::automatic, which is how I encountered it.

I'm attaching the example added to the tests, and also the fix
that makes matrix::division_free_elimination() use normal() in
stead of expand(). An alternative would be to use normal() in
the matrix::pivot() routine, but I decided against it because:
1) that way the results of the normal() call would be discarded;
2) that could slow down e.g. matrix::gauss_elimination() which
   also uses matrix::pivot(), but normalizes expressions itself.

Of course, using normal() makes division_free_elimination only
truly "division-free" if the input matrix didn't have fractions.
In other cases, GCDs will be computed.
-------------- next part --------------
From 710b1b41b477463ac8668134c970a0f9d90b2a45 Mon Sep 17 00:00:00 2001
From: Vitaly Magerya <magv at tx97.net>
Date: Tue, 9 Oct 2018 13:32:29 +0200
Subject: [PATCH 1/2] Add a test case for matrix::solve that breaks for
 solve_algo::divfree

---
 check/exam_matrices.cpp | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp
index 7770e6bf..3ef7c9d5 100644
--- a/check/exam_matrices.cpp
+++ b/check/exam_matrices.cpp
@@ -215,6 +215,34 @@ static unsigned matrix_solve2()
 	return result;
 }
 
+static unsigned matrix_solve3()
+{
+	unsigned result = 0;
+	symbol x("x");
+	symbol t1("t1"), t2("t2"), t3("t3");
+	matrix A = {
+        {3+6*x, 6*(x+x*x)/(2+3*x), 0},
+        {-(2+7*x+6*x*x)/x, -2-2*x, 0},
+        {-2*(2+3*x)/(1+2*x), -6*x/(1+2*x), 1+4*x}
+    };
+	matrix B = {{0}, {0}, {0}};
+	matrix X = {{t1}, {t2}, {t3}};
+	matrix cmp = {{1, 0},
+	              {3, 2},
+	              {3, 4}};
+    for (auto algo : vector<int>({
+        solve_algo::gauss, solve_algo::divfree, solve_algo::bareiss, solve_algo::markowitz
+    })) {
+        matrix sol(A.solve(X, B, algo));
+        if (!normal((A*sol - B).evalm()).is_zero_matrix()) {
+            clog << "Solving " << A << " * " << X << " == " << B << " with algo=" << algo << endl
+                 << "erroneously returned " << sol << endl;
+            result += 1;
+        }
+    }
+	return result;
+}
+
 static unsigned matrix_evalm()
 {
 	unsigned result = 0;
@@ -365,6 +393,7 @@ unsigned exam_matrices()
 	result += matrix_invert2();  cout << '.' << flush;
 	result += matrix_invert3();  cout << '.' << flush;
 	result += matrix_solve2();  cout << '.' << flush;
+	result += matrix_solve3();  cout << '.' << flush;
 	result += matrix_evalm();  cout << "." << flush;
 	result += matrix_rank();  cout << "." << flush;
 	result += matrix_solve_nonnormal();  cout << "." << flush;
-- 
2.17.1

-------------- next part --------------
From 107f42727b2de9040a8b50ee2a2bf5d7c2408ac9 Mon Sep 17 00:00:00 2001
From: Vitaly Magerya <magv at tx97.net>
Date: Tue, 9 Oct 2018 14:19:59 +0200
Subject: [PATCH 2/2] Fix matrix::division_free_elimination for unnormal zeros

This fix makes the division_free_elimination function not really
division-free, but at least it returns the correct result now.
---
 ginac/matrix.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index f5c99a37..dc2105ae 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1470,7 +1470,7 @@ int matrix::division_free_elimination(const bool det)
 				sign = -sign;
 			for (unsigned r2=r0+1; r2<m; ++r2) {
 				for (unsigned c=c0+1; c<n; ++c)
-					this->m[r2*n+c] = (this->m[r0*n+c0]*this->m[r2*n+c] - this->m[r2*n+c0]*this->m[r0*n+c]).expand();
+					this->m[r2*n+c] = (this->m[r0*n+c0]*this->m[r2*n+c] - this->m[r2*n+c0]*this->m[r0*n+c]).normal();
 				// fill up left hand side with zeros
 				for (unsigned c=r0; c<=c0; ++c)
 					this->m[r2*n+c] = _ex0;
-- 
2.17.1