[GiNaC-list] [patch] Fix a failure in matrix::division_free_elimination()
Vitaly Magerya
vmagerya at gmail.com
Tue Oct 9 14:44:20 CEST 2018
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
More information about the GiNaC-list
mailing list