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 586c52f57673a0bc20a33e89ba152f01697e31f7..7770e6bf829beaaa3101b965bb57571d9e07242b 100644 (file)
--- 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. */
 /** @file exam_matrices.cpp
  *
  *  Here we examine manipulations on GiNaC's symbolic matrices. */


@@ -275,6 +276,35 @@ static unsigned matrix_rank()
        return result;  
 }
 
        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;
 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_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;
        result += matrix_misc();  cout << '.' << flush;
        
        return result;

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 80ca3e459ff42a4002920fd68466054e0398c42a..bf7f8347c9a1f7b31645595d408ce602cf77081d 100644 (file)
--- 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
                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
                                ++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 {
                                        throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
                                }
                        } else {