From: Vitaly Magerya <vmagerya@gmail.com>
Date: Fri, 12 Oct 2018 18:45:27 +0000 (+0200)
Subject: Handle un-normal zeros properly in the division-free elimination.
X-Git-Tag: release_1-7-5~4
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=96be7a17790700998e2d071e00ad5c1ffb4a2d3e;ds=sidebyside

Handle un-normal zeros properly in the division-free elimination.

Call .normal() instead of just .expand(), such that matrix::pivot() finds
the pivot.

Note that this problem also affects matrix::solve() with 'algo'
set to solve_algo::automatic.

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.
---

diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp
index 7770e6bf..1f402ba9 100644
--- a/check/exam_matrices.cpp
+++ b/check/exam_matrices.cpp
@@ -215,6 +215,31 @@ 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}};
+	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 +390,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;
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;