From: Richard Kreckel 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 . --- 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 {