From 6db7c98ab607db5329f7c94e3524af67fd80143a Mon Sep 17 00:00:00 2001
From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Sat, 29 Apr 2000 02:13:29 +0000
Subject: [PATCH] - Extended check (FIXME: may under certain conditions rout
 oom).

---
 check/check_matrices.cpp | 57 +++++++++++++++++++++++++++++++------
 check/checks.h           |  3 +-
 check/genex.cpp          | 61 +++++++++++++++++++++++++++++-----------
 3 files changed, 94 insertions(+), 27 deletions(-)

diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index 71893bb1..c1fa6257 100644
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -22,14 +22,14 @@
 
 #include "checks.h"
 
-// determinants of some sparse symbolic size x size matrices over
-// an integral domain.
+/* determinants of some sparse symbolic matrices with coefficients in
+ * an integral domain. */
 static unsigned integdom_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a");
     
-    for (int size=3; size<17; ++size) {
+    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:
@@ -37,7 +37,7 @@ static unsigned integdom_matrix_determinants(void)
         }
         for (int c=0; c<size; ++c) {
             // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant vanishes:
+            // to guarantee that the determinant is zero:
             A.set(size-1,c,A(0,c)-A(size-2,c));
         }
         if (!A.determinant().is_zero()) {
@@ -51,22 +51,60 @@ static unsigned integdom_matrix_determinants(void)
     return result;
 }
 
+/* determinants of some sparse 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));
+        }
+        if (!A.determinant().is_zero()) {
+            clog << "Determinant of " << size << "x" << size << " matrix "
+                 << endl << A << endl
+                 << "was not found to vanish!" << endl;
+            ++result;
+        }
+    }
+    
+    return result;
+}
+
+/* Some quite wild determinants with functions and stuff like that. */
+static unsigned wild_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
-    for (int size=3; size<13; ++size) {
+    for (int size=3; size<6; ++size) {
         matrix A(size,size);
         for (int r=0; r<size-1; ++r) {
             // populate one element in each row:
-            // FIXME: the line using sparse_tree() should be used:
-            // A.set(r,unsigned(rand()%size),sparse_tree(a, b, c, 3, true, true)/sparse_tree(a, b, c, 2, true, true));
-            A.set(r,unsigned(rand()%size),dense_univariate_poly(a,4)/dense_univariate_poly(a,2));
+            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 vanishes:
+            // to guarantee that the determinant is zero:
             A.set(size-1,c,A(0,c)-A(size-2,c));
         }
         if (!A.determinant().is_zero()) {
@@ -89,6 +127,7 @@ unsigned check_matrices(void)
     
     result += integdom_matrix_determinants();  cout << '.' << flush;
     result += rational_matrix_determinants();  cout << '.' << flush;
+    result += wild_matrix_determinants();  cout << '.' << flush;
     
     if (!result) {
         cout << " passed " << endl;
diff --git a/check/checks.h b/check/checks.h
index 4c79e2f3..d95ad1a2 100644
--- a/check/checks.h
+++ b/check/checks.h
@@ -36,7 +36,8 @@ using namespace GiNaC;
 const ex dense_univariate_poly(const symbol & x, unsigned degree);
 const ex dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree);
 const ex sparse_tree(const symbol & x, const symbol & y, const symbol & z,
-                     int level, bool trig = false, bool rational = true);
+                     int level,
+                     bool trig = false, bool rational = true, bool complex = false);
 
 // prototypes for all individual checks should be unsigned fcn():
 unsigned check_numeric();
diff --git a/check/genex.cpp b/check/genex.cpp
index dacb40eb..18c2ded4 100644
--- a/check/genex.cpp
+++ b/check/genex.cpp
@@ -58,11 +58,13 @@ dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree)
     return bipoly;
 }
 
+/* Chose a randum symbol or number from the argument list. */
 const ex
 random_symbol(const symbol & x,
               const symbol & y,
               const symbol & z,
-              bool rational = true)
+              bool rational = true,
+              bool complex = false)
 {
     ex e;
     switch (abs(rand()) % 4) {
@@ -76,13 +78,15 @@ random_symbol(const symbol & x,
             e = z;
             break;
         case 3: {
-            int c1 = rand() % 20 - 10;
-            int c2 = rand() % 20 - 10;
-            if (c1 == 0) c1 = 1;
-            if (c2 == 0) c2 = 1;
+            int c1;
+            do { c1 = rand()%20 - 10; } while (!c1);
+            int c2;
+            do { c2 = rand()%20 - 10; } while (!c2);
             if (!rational)
                 c2 = 1;
-            e = numeric(c1) / numeric(c2);
+            e = numeric(c1, c2);
+            if (complex && !(rand()%5))
+                e = e*I;
             break;
         }
     }
@@ -96,24 +100,33 @@ sparse_tree(const symbol & x,
             const symbol & z,
             int level,
             bool trig = false,    // true includes trigonomatric functions
-            bool rational = true) // false includes coefficients in Q
+            bool rational = true, // false excludes coefficients in Q
+            bool complex = false) // true includes complex numbers
 {
     if (level == 0)
-        return random_symbol(x,y,z,rational);
-    switch (abs(rand()) % 7) {
+        return random_symbol(x,y,z,rational,complex);
+    switch (abs(rand()) % 10) {
         case 0:
         case 1:
-            return add(sparse_tree(x,y,z,level-1, trig, rational),
-                       sparse_tree(x,y,z,level-1, trig, rational));
         case 2:
         case 3:
-            return mul(sparse_tree(x,y,z,level-1, trig, rational),
+            return add(sparse_tree(x,y,z,level-1, trig, rational),
                        sparse_tree(x,y,z,level-1, trig, rational));
         case 4:
         case 5:
-            return power(sparse_tree(x,y,z,level-1, trig, rational),
-                         abs(rand() % 4));
         case 6:
+            return mul(sparse_tree(x,y,z,level-1, trig, rational),
+                       sparse_tree(x,y,z,level-1, trig, rational));
+        case 7:
+        case 8: {
+            ex powbase;
+            do {
+                powbase = sparse_tree(x,y,z,level-1, trig, rational);
+            } while (powbase.is_zero());
+            return pow(powbase, abs(rand() % 4));
+            break;
+        }
+        case 9:
             if (trig) {
                 switch (abs(rand()) % 4) {
                     case 0:
@@ -122,10 +135,24 @@ sparse_tree(const symbol & x,
                         return cos(sparse_tree(x,y,z,level-1, trig, rational));
                     case 2:
                         return exp(sparse_tree(x,y,z,level-1, trig, rational));
-                    case 3:
-                        return log(sparse_tree(x,y,z,level-1, trig, rational));
+                    case 3: {
+                        ex logex;
+                        do {
+                            ex logarg;
+                            do {
+                                logarg = sparse_tree(x,y,z,level-1, trig, rational);
+                            } while (logarg.is_zero());
+                            // Keep the evaluator from accidentally plugging an
+                            // unwanted I in the tree:
+                            if (!complex && logarg.info(info_flags::negative))
+                                logarg = -logarg;
+                            logex = log(logarg);
+                        } while (logex.is_zero());
+                        return logex;
+                        break;
+                    }
                 }
             } else
-                return random_symbol(x,y,z,rational);
+                return random_symbol(x,y,z,rational,complex);
     }
 }
-- 
2.44.0