X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=check%2Fcheck_matrices.cpp;h=3e702b8e963c174bdf37a896b91f89d7f5003d18;hp=c1fa6257599938bc57c015e363a6822bf876a9c2;hb=9df145c8bfa8ce9f2cbe6c05673481b6ca4c4c22;hpb=fe9dbfb9947b24149b3ce7dd9285f27ab286cbd7

diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index c1fa6257..3e702b8e 100644
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -31,15 +31,13 @@ static unsigned integdom_matrix_determinants(void)
     
     for (int size=3; size<20; ++size) {
         matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
+        // populate one element in each row:
+        for (int r=0; r<size-1; ++r)
             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
+        // set the last row to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int c=0; c<size; ++c)
             A.set(size-1,c,A(0,c)-A(size-2,c));
-        }
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
@@ -51,29 +49,30 @@ static unsigned integdom_matrix_determinants(void)
     return result;
 }
 
-/* determinants of some sparse symbolic matrices with multivariate rational
- * function coefficients. */
+/* determinants of some symbolic matrices with multivariate rational function
+ * coefficients. */
 static unsigned rational_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
-
+    
     for (int size=3; size<8; ++size) {
         matrix A(size,size);
         for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 4, false, false, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, false, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
+            // populate one or two elements in each row:
+            for (int ec=0; ec<2; ++ec) {
+                ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
+                ex denom;
+                do {
+                    denom = sparse_tree(a, b, c, rand()%2, false, false, false);
+                } while (denom.is_zero());
+                A.set(r,unsigned(rand()%size),numer/denom);
+            }
         }
+        // set the last row to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int co=0; co<size; ++co)
+            A.set(size-1,co,A(0,co)-A(size-2,co));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
@@ -85,28 +84,29 @@ static unsigned rational_matrix_determinants(void)
     return result;
 }
 
-/* Some quite wild determinants with functions and stuff like that. */
-static unsigned wild_matrix_determinants(void)
+/* Some quite funny determinants with functions and stuff like that inside. */
+static unsigned funny_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
-    for (int size=3; size<6; ++size) {
+    for (int size=3; size<7; ++size) {
         matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 3, true, true, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, true, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
+        for (int co=0; co<size-1; ++co) {
+            // populate one or two elements in each row:
+            for (int ec=0; ec<2; ++ec) {
+                ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
+                ex denom;
+                do {
+                    denom = sparse_tree(a, b, c, rand()%2, false, true, false);
+                } while (denom.is_zero());
+                A.set(unsigned(rand()%size),co,numer/denom);
+            }
         }
+        // set the last column to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int ro=0; ro<size; ++ro)
+            A.set(ro,size-1,A(ro,0)-A(ro,size-2));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
@@ -118,6 +118,41 @@ static unsigned wild_matrix_determinants(void)
     return result;
 }
 
+/* compare results from different determinant algorithms.*/
+static unsigned compare_matrix_determinants(void)
+{
+    unsigned result = 0;
+    symbol a("a");
+    
+    for (int size=2; size<6; ++size) {
+        matrix A(size,size);
+        for (int co=0; co<size; ++co) {
+            for (int ro=0; ro<size; ++ro) {
+                // populate some elements
+                ex elem = 0;
+                if (rand()%(size-1) == 0)
+                    elem = sparse_tree(a, a, a, rand()%3, false, true, false);
+                A.set(ro,co,elem);
+            }
+        }
+        ex det_gauss = A.determinant(determinant_algo::gauss);
+        ex det_laplace = A.determinant(determinant_algo::laplace);
+        ex det_bareiss = A.determinant(determinant_algo::bareiss);
+        if ((det_gauss-det_laplace).normal() != 0 ||
+            (det_bareiss-det_laplace).normal() != 0) {
+            clog << "Determinant of " << size << "x" << size << " matrix "
+                 << endl << A << endl
+                 << "is inconsistent between different algorithms:" << endl
+                 << "Gauss elimination:   " << det_gauss << endl
+                 << "Minor elimination:   " << det_laplace << endl
+                 << "Fraction-free elim.: " << det_bareiss << endl;
+            ++result;
+        }
+    }
+    
+    return result;
+}
+
 unsigned check_matrices(void)
 {
     unsigned result = 0;
@@ -127,7 +162,8 @@ unsigned check_matrices(void)
     
     result += integdom_matrix_determinants();  cout << '.' << flush;
     result += rational_matrix_determinants();  cout << '.' << flush;
-    result += wild_matrix_determinants();  cout << '.' << flush;
+    result += funny_matrix_determinants();  cout << '.' << flush;
+    result += compare_matrix_determinants();  cout << '.' << flush;
     
     if (!result) {
         cout << " passed " << endl;