- Completely restructured the checks in subdir check/.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Tue, 29 Feb 2000 17:28:13 +0000 (17:28 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Tue, 29 Feb 2000 17:28:13 +0000 (17:28 +0000)
- Optimized Laplace expansion of determinants to not compute anything more
  than once.
- Moved int permutation_sign(vector<T> s) from matrix.cpp to utils.h.
- A handful of other fixes assiciated with the above.

51 files changed:
INSTALL
Makefile.in
NEWS
check/Makefile.am
check/Makefile.in
check/check_inifcns.cpp [moved from check/inifcns_consist.cpp with 93% similarity]
check/check_lsolve.cpp [new file with mode: 0644]
check/check_matrices.cpp [new file with mode: 0644]
check/check_numeric.cpp [new file with mode: 0644]
check/checks.cpp [new file with mode: 0644]
check/checks.h [moved from check/check.h with 57% similarity]
check/checks.ref [new file with mode: 0644]
check/exam_differentiation.cpp [moved from check/differentiation.cpp with 90% similarity]
check/exam_lsolve.cpp [moved from check/linear_solve.cpp with 81% similarity]
check/exam_matrices.cpp [moved from check/matrix_checks.cpp with 86% similarity]
check/exam_misc.cpp [moved from check/expand_subs.cpp with 60% similarity]
check/exam_noncommut.cpp [moved from check/lortensor_check.cpp with 87% similarity]
check/exam_normalization.cpp [moved from check/normalization.cpp with 69% similarity]
check/exam_numeric.cpp [moved from check/numeric_consist.cpp with 87% similarity]
check/exam_paranoia.cpp [moved from check/paranoia_check.cpp with 86% similarity]
check/exam_polygcd.cpp [moved from check/poly_gcd.cpp with 87% similarity]
check/exam_powerlaws.cpp [moved from check/powerlaws.cpp with 93% similarity]
check/exam_pseries.cpp [moved from check/series_expansion.cpp with 90% similarity]
check/exams.cpp [new file with mode: 0644]
check/exams.h [new file with mode: 0644]
check/exams.ref [new file with mode: 0644]
check/genex.cpp [new file with mode: 0644]
check/numeric_output.cpp [deleted file]
check/result.ref [deleted file]
check/run_checks
check/run_exams [new file with mode: 0755]
check/run_times [new file with mode: 0755]
check/time_dennyfliegner.cpp [new file with mode: 0644]
check/time_gammaseries.cpp [new file with mode: 0644]
check/timer.cpp [moved from check/fcntimer.cpp with 63% similarity]
check/times.cpp [moved from check/main.cpp with 56% similarity]
check/times.h [new file with mode: 0644]
check/times.ref [new file with mode: 0644]
cint/Makefile.in
config.h.in
doc/Makefile.in
doc/reference/Makefile.in
doc/tutorial/Makefile.in
ginac/Makefile.in
ginac/inifcns.cpp
ginac/matrix.cpp
ginac/matrix.h
ginac/remember.h
ginac/utils.h
ginsh/Makefile.in
tools/Makefile.in

diff --git a/INSTALL b/INSTALL
index aee07da..a7ef9be 100644 (file)
--- a/INSTALL
+++ b/INSTALL
@@ -108,5 +108,5 @@ Good luck!
 
   Known to work with:  |  Known not to work with:
 -----------------------+----------------------------
-  Cint 5.14.25         |  Cint 5.14.24
-  Cint 5.14.26         |  Cint 5.14.29
+  Cint 5.14.31         |  Cint before 5.14.29
+
index 7e2ca91..1ab8391 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
diff --git a/NEWS b/NEWS
index fb41acc..601e575 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,13 @@
 This file records noteworthy changes.
 
+0.5.4 (?? 2000)
+* Algorithms in class matrix (determinant and solve) were replaced by
+  less brain-dead ones and should now have much better performance.
+* Checks were reorganized and split up into three parts:
+  a) exams (small tests with predefined input)
+  b) checks (lenghty consistency checks)
+  c) timings (for crude benchmarking)
+
 0.5.3 (23 February 2000)
 * A more flexible scheme for registering functions was implemented,
   allowing for remembering, too.
index 87976d8..96ec408 100644 (file)
@@ -1,12 +1,17 @@
 ## Process this file with automake to produce Makefile.in
-TESTS = run_checks
-check_PROGRAMS = check_ginac
-check_ginac_SOURCES = paranoia_check.cpp numeric_output.cpp \
-  numeric_consist.cpp powerlaws.cpp expand_subs.cpp inifcns_consist.cpp \
-  differentiation.cpp poly_gcd.cpp normalization.cpp matrix_checks.cpp \
-  linear_solve.cpp series_expansion.cpp lortensor_check.cpp \
-  fcntimer.cpp main.cpp check.h
-check_ginac_LDADD = ../ginac/libginac.la
+TESTS = run_exams run_checks run_times
+check_PROGRAMS = exams checks times
+exams_SOURCES = exam_paranoia.cpp exam_numeric.cpp exam_powerlaws.cpp \
+  exam_differentiation.cpp exam_polygcd.cpp exam_normalization.cpp \
+  exam_pseries.cpp exam_matrices.cpp exam_lsolve.cpp exam_noncommut.cpp \
+  exam_misc.cpp exams.cpp exams.h
+exams_LDADD = ../ginac/libginac.la
+checks_SOURCES =  check_numeric.cpp check_inifcns.cpp check_matrices.cpp \
+  check_lsolve.cpp genex.cpp checks.cpp checks.h
+checks_LDADD = ../ginac/libginac.la
+times_SOURCES = time_dennyfliegner.cpp time_gammaseries.cpp timer.cpp \
+  times.cpp times.h
+times_LDADD = ../ginac/libginac.la
 INCLUDES = -I$(srcdir)/../ginac
-CLEANFILES = result.out
-EXTRA_DIST = result.ref run_checks
+CLEANFILES = exams.out checks.out times.out
+EXTRA_DIST = exams.ref checks.ref times.ref run_exams run_checks run_times
index 0a17c44..0410e06 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
@@ -102,14 +103,20 @@ TUTORIAL_TARGETS = @TUTORIAL_TARGETS@
 VERSION = @VERSION@
 YACC = @YACC@
 
-TESTS = run_checks
-check_PROGRAMS = check_ginac
-check_ginac_SOURCES = paranoia_check.cpp numeric_output.cpp   numeric_consist.cpp powerlaws.cpp expand_subs.cpp inifcns_consist.cpp   differentiation.cpp poly_gcd.cpp normalization.cpp matrix_checks.cpp   linear_solve.cpp series_expansion.cpp lortensor_check.cpp   fcntimer.cpp main.cpp check.h
+TESTS = run_exams run_checks run_times
+check_PROGRAMS = exams checks times
+exams_SOURCES = exam_paranoia.cpp exam_numeric.cpp exam_powerlaws.cpp   exam_differentiation.cpp exam_polygcd.cpp exam_normalization.cpp   exam_pseries.cpp exam_matrices.cpp exam_lsolve.cpp exam_noncommut.cpp   exam_misc.cpp exams.cpp exams.h
 
-check_ginac_LDADD = ../ginac/libginac.la
+exams_LDADD = ../ginac/libginac.la
+checks_SOURCES = check_numeric.cpp check_inifcns.cpp check_matrices.cpp   check_lsolve.cpp genex.cpp checks.cpp checks.h
+
+checks_LDADD = ../ginac/libginac.la
+times_SOURCES = time_dennyfliegner.cpp time_gammaseries.cpp timer.cpp   times.cpp times.h
+
+times_LDADD = ../ginac/libginac.la
 INCLUDES = -I$(srcdir)/../ginac
-CLEANFILES = result.out
-EXTRA_DIST = result.ref run_checks
+CLEANFILES = exams.out checks.out times.out
+EXTRA_DIST = exams.ref checks.ref times.ref run_exams run_checks run_times
 mkinstalldirs = $(SHELL) $(top_srcdir)/mkinstalldirs
 CONFIG_HEADER = ../config.h
 CONFIG_CLEAN_FILES = 
@@ -118,12 +125,19 @@ DEFS = @DEFS@ -I. -I$(srcdir) -I..
 CPPFLAGS = @CPPFLAGS@
 LDFLAGS = @LDFLAGS@
 LIBS = @LIBS@
-check_ginac_OBJECTS =  paranoia_check.o numeric_output.o \
-numeric_consist.o powerlaws.o expand_subs.o inifcns_consist.o \
-differentiation.o poly_gcd.o normalization.o matrix_checks.o \
-linear_solve.o series_expansion.o lortensor_check.o fcntimer.o main.o
-check_ginac_DEPENDENCIES =  ../ginac/libginac.la
-check_ginac_LDFLAGS = 
+exams_OBJECTS =  exam_paranoia.o exam_numeric.o exam_powerlaws.o \
+exam_differentiation.o exam_polygcd.o exam_normalization.o \
+exam_pseries.o exam_matrices.o exam_lsolve.o exam_noncommut.o \
+exam_misc.o exams.o
+exams_DEPENDENCIES =  ../ginac/libginac.la
+exams_LDFLAGS = 
+checks_OBJECTS =  check_numeric.o check_inifcns.o check_matrices.o \
+check_lsolve.o genex.o checks.o
+checks_DEPENDENCIES =  ../ginac/libginac.la
+checks_LDFLAGS = 
+times_OBJECTS =  time_dennyfliegner.o time_gammaseries.o timer.o times.o
+times_DEPENDENCIES =  ../ginac/libginac.la
+times_LDFLAGS = 
 CXXFLAGS = @CXXFLAGS@
 CXXCOMPILE = $(CXX) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS)
 LTCXXCOMPILE = $(LIBTOOL) --mode=compile $(CXX) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS)
@@ -141,14 +155,16 @@ DISTFILES = $(DIST_COMMON) $(SOURCES) $(HEADERS) $(TEXINFOS) $(EXTRA_DIST)
 
 TAR = tar
 GZIP_ENV = --best
-DEP_FILES =  .deps/differentiation.P .deps/expand_subs.P \
-.deps/fcntimer.P .deps/inifcns_consist.P .deps/linear_solve.P \
-.deps/lortensor_check.P .deps/main.P .deps/matrix_checks.P \
-.deps/normalization.P .deps/numeric_consist.P .deps/numeric_output.P \
-.deps/paranoia_check.P .deps/poly_gcd.P .deps/powerlaws.P \
-.deps/series_expansion.P
-SOURCES = $(check_ginac_SOURCES)
-OBJECTS = $(check_ginac_OBJECTS)
+DEP_FILES =  .deps/check_inifcns.P .deps/check_lsolve.P \
+.deps/check_matrices.P .deps/check_numeric.P .deps/checks.P \
+.deps/exam_differentiation.P .deps/exam_lsolve.P .deps/exam_matrices.P \
+.deps/exam_misc.P .deps/exam_noncommut.P .deps/exam_normalization.P \
+.deps/exam_numeric.P .deps/exam_paranoia.P .deps/exam_polygcd.P \
+.deps/exam_powerlaws.P .deps/exam_pseries.P .deps/exams.P .deps/genex.P \
+.deps/time_dennyfliegner.P .deps/time_gammaseries.P .deps/timer.P \
+.deps/times.P
+SOURCES = $(exams_SOURCES) $(checks_SOURCES) $(times_SOURCES)
+OBJECTS = $(exams_OBJECTS) $(checks_OBJECTS) $(times_OBJECTS)
 
 all: all-redirect
 .SUFFIXES:
@@ -202,9 +218,17 @@ distclean-libtool:
 
 maintainer-clean-libtool:
 
-check_ginac: $(check_ginac_OBJECTS) $(check_ginac_DEPENDENCIES)
-       @rm -f check_ginac
-       $(CXXLINK) $(check_ginac_LDFLAGS) $(check_ginac_OBJECTS) $(check_ginac_LDADD) $(LIBS)
+exams: $(exams_OBJECTS) $(exams_DEPENDENCIES)
+       @rm -f exams
+       $(CXXLINK) $(exams_LDFLAGS) $(exams_OBJECTS) $(exams_LDADD) $(LIBS)
+
+checks: $(checks_OBJECTS) $(checks_DEPENDENCIES)
+       @rm -f checks
+       $(CXXLINK) $(checks_LDFLAGS) $(checks_OBJECTS) $(checks_LDADD) $(LIBS)
+
+times: $(times_OBJECTS) $(times_DEPENDENCIES)
+       @rm -f times
+       $(CXXLINK) $(times_LDFLAGS) $(times_OBJECTS) $(times_LDADD) $(LIBS)
 .cpp.o:
        $(CXXCOMPILE) -c $<
 .cpp.lo:
similarity index 93%
rename from check/inifcns_consist.cpp
rename to check/check_inifcns.cpp
index c30ddef..fb7d999 100644 (file)
@@ -1,4 +1,4 @@
-/** @file inifcns_consist.cpp
+/** @file check_inifcns.cpp
  *
  *  This test routine applies assorted tests on initially known higher level
  *  functions. */
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "checks.h"
 
 /* Some tests on the sine trigonometric function. */
 static unsigned inifcns_consist_sin(void)
 {
     unsigned result = 0;
-    bool errorflag;
+    bool errorflag = false;
     
     // sin(n*Pi) == 0?
     errorflag = false;
@@ -68,7 +64,7 @@ static unsigned inifcns_consist_sin(void)
     errorflag = false;
     ex argument;
     numeric epsilon(double(1e-8));
-    for (int n=-240; n<=240; ++n) {
+    for (int n=-340; n<=340; ++n) {
         argument = n*Pi/60;
         if (abs(sin(evalf(argument))-evalf(sin(argument)))>epsilon) {
             clog << "sin(" << argument << ") returns "
@@ -123,7 +119,7 @@ static unsigned inifcns_consist_cos(void)
     errorflag = false;
     ex argument;
     numeric epsilon(double(1e-8));
-    for (int n=-240; n<=240; ++n) {
+    for (int n=-340; n<=340; ++n) {
         argument = n*Pi/60;
         if (abs(cos(evalf(argument))-evalf(cos(argument)))>epsilon) {
             clog << "cos(" << argument << ") returns "
@@ -151,7 +147,7 @@ static unsigned inifcns_consist_tan(void)
     errorflag = false;
     ex argument;
     numeric epsilon(double(1e-8));
-    for (int n=-240; n<=240; ++n) {
+    for (int n=-340; n<=340; ++n) {
         if (!(n%30) && (n%60))  // skip poles
             ++n;
         argument = n*Pi/60;
@@ -332,26 +328,26 @@ static unsigned inifcns_consist_zeta(void)
     return result;
 }
 
-unsigned inifcns_consist(void)
+unsigned check_inifcns(void)
 {
     unsigned result = 0;
 
-    cout << "checking consistency of symbolic functions..." << flush;
+    cout << "checking consistency of symbolic functions" << flush;
     clog << "---------consistency of symbolic functions:" << endl;
     
-    result += inifcns_consist_sin();
-    result += inifcns_consist_cos();
-    result += inifcns_consist_tan();
-    result += inifcns_consist_trans();
-    result += inifcns_consist_gamma();
-    result += inifcns_consist_psi();
-    result += inifcns_consist_zeta();
+    result += inifcns_consist_sin();  cout << '.' << flush;
+    result += inifcns_consist_cos();  cout << '.' << flush;
+    result += inifcns_consist_tan();  cout << '.' << flush;
+    result += inifcns_consist_trans();  cout << '.' << flush;
+    result += inifcns_consist_gamma();  cout << '.' << flush;
+    result += inifcns_consist_psi();  cout << '.' << flush;
+    result += inifcns_consist_zeta();  cout << '.' << flush;
 
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     
     return result;
diff --git a/check/check_lsolve.cpp b/check/check_lsolve.cpp
new file mode 100644 (file)
index 0000000..4d1f237
--- /dev/null
@@ -0,0 +1,128 @@
+/** @file check_lsolve.cpp
+ *
+ *  These test routines do some simple checks on solving linear systems of
+ *  symbolic equations. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "checks.h"
+
+static unsigned lsolve1(int size)
+{
+    // A dense size x size matrix in dense univariate random polynomials
+    // of order 4.
+    unsigned result = 0;
+    symbol a("a");
+    ex sol;
+    
+    // Create two dense linear matrices A and B where all entries are random
+    // univariate polynomials 
+    matrix A(size,size), B(size,2), X(size,2);
+    for (int ro=0; ro<size; ++ro) {
+        for (int co=0; co<size; ++co)
+            A.set(ro,co,dense_univariate_poly(a, 5));
+        for (int co=0; co<2; ++co)
+            B.set(ro,co,dense_univariate_poly(a, 5));
+    }
+    if (A.determinant().is_zero())
+        clog << "lsolve1: singular system!" << endl;
+    
+    // Solve the system A*X==B:
+    X = A.old_solve(B);
+    
+    // check the result:
+    bool errorflag = false;
+    matrix Aux(size,2);
+    Aux = A.mul(X).sub(B);
+    for (int ro=0; ro<size && !errorflag; ++ro)
+        for (int co=0; co<2; ++co)
+            if (!(Aux(ro,co)).normal().is_zero())
+                errorflag = true;
+    if (errorflag) {
+        clog << "Our solve method claims that A*X==B, with matrices" << endl
+             << "A == " << A << endl
+             << "X == " << X << endl
+             << "B == " << B << endl;
+        ++result;
+    }
+    return result;
+}
+
+static unsigned lsolve2(int size)
+{
+    // A dense size x size matrix in dense bivariate random polynomials
+    // of order 2.
+    unsigned result = 0;
+    symbol a("a"), b("b");
+    ex sol;
+    
+    // Create two dense linear matrices A and B where all entries are dense random
+    // bivariate polynomials:
+    matrix A(size,size), B(size,2), X(size,2);
+    for (int ro=0; ro<size; ++ro) {
+        for (int co=0; co<size; ++co)
+            A.set(ro,co,dense_bivariate_poly(a,b,2));
+        for (int co=0; co<2; ++co)
+            B.set(ro,co,dense_bivariate_poly(a,b,2));
+    }
+    if (A.determinant().is_zero())
+        clog << "lsolve2: singular system!" << endl;
+    
+    // Solve the system A*X==B:
+    X = A.old_solve(B);
+    
+    // check the result:
+    bool errorflag = false;
+    matrix Aux(size,2);
+    Aux = A.mul(X).sub(B);
+    for (int ro=0; ro<size && !errorflag; ++ro)
+        for (int co=0; co<2; ++co)
+            if (!(Aux(ro,co)).normal().is_zero())
+                errorflag = true;
+    if (errorflag) {
+        clog << "Our solve method claims that A*X==B, with matrices" << endl
+             << "A == " << A << endl
+             << "X == " << X << endl
+             << "B == " << B << endl;
+        ++result;
+    }
+    return result;
+}
+
+unsigned check_lsolve(void)
+{
+    unsigned result = 0;
+    
+    cout << "checking linear solve" << flush;
+    clog << "---------linear solve:" << endl;
+    
+    //result += lsolve1(2);  cout << '.' << flush;
+    //result += lsolve1(3);  cout << '.' << flush;
+    //result += lsolve2(2);  cout << '.' << flush;
+    //result += lsolve2(3);  cout << '.' << flush;
+    
+    if (!result) {
+        cout << " passed " << endl;
+        clog << "(no output)" << endl;
+    } else {
+        cout << " failed " << endl;
+    }
+    
+    return result;
+}
diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
new file mode 100644 (file)
index 0000000..ec7efef
--- /dev/null
@@ -0,0 +1,70 @@
+/** @file check_matrices.cpp
+ *
+ *  Here we test manipulations on GiNaC's symbolic matrices. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "checks.h"
+
+// determinants of some sparse symbolic size x size matrices
+static unsigned matrix_determinants(void)
+{
+    unsigned result = 0;
+    symbol a("a");
+
+    for (int size=3; size<16; ++size) {
+        matrix A(size,size);
+        for (int c=0; c<size; ++c) {
+            for (int r=0;r<size-1; ++r)
+                // populate 10 percent of the entries, the rest remains 0:
+                if (!(int)(10.0*rand()/(RAND_MAX+1.0)))
+                    A.set(r,c,dense_univariate_poly(a,5));
+            // set the last line to a linear combination of two other lines
+            // to guarantee that the determinant vanishes:
+            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;
+}
+
+unsigned check_matrices(void)
+{
+    unsigned result = 0;
+    
+    cout << "checking symbolic matrix manipulations" << flush;
+    clog << "---------symbolic matrix manipulations:" << endl;
+    
+    result += matrix_determinants();  cout << '.' << flush;
+    
+    if (!result) {
+        cout << " passed " << endl;
+        clog << "(no output)" << endl;
+    } else {
+        cout << " failed " << endl;
+    }
+    
+    return result;
+}
diff --git a/check/check_numeric.cpp b/check/check_numeric.cpp
new file mode 100644 (file)
index 0000000..811ba52
--- /dev/null
@@ -0,0 +1,120 @@
+/** @file check_numeric.cpp
+ *
+ *  These exams creates some numbers and check the result of several boolean
+ *  tests on these numbers like is_integer() etc... */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "checks.h"
+
+#ifndef NO_NAMESPACE_GINAC
+using namespace GiNaC;
+#endif // ndef NO_NAMESPACE_GINAC
+
+/* Simple and maybe somewhat pointless consistency tests of assorted tests and
+ * conversions. */
+static unsigned check_numeric1(void)
+{
+    unsigned result = 0;
+    bool errorflag = false;
+    int re_q, im_q;
+    
+    // Check some numerator and denominator calculations:
+    for (int i=0; i<200; ++i) {
+        do { re_q = rand(); } while (re_q == 0);
+        do { im_q = rand(); } while (im_q == 0);
+        numeric r(rand()-RAND_MAX/2, re_q);
+        numeric i(rand()-RAND_MAX/2, im_q);
+        numeric z = r + I*i;
+        numeric p = numer(z);
+        numeric q = denom(z);
+        numeric res = p/q;
+        if (res != z) {
+            clog << z << " erroneously transformed into " 
+                 << p << "/" << q << " by numer() and denom()" << endl;
+            errorflag = true;
+        }
+    }
+    if (errorflag)
+        ++result;
+    
+    return result;
+}
+
+static unsigned check_numeric2(void)
+{
+    unsigned result = 0;
+    bool errorflag = false;
+    int i_num, i_den;
+    
+    // Check non-nested radicals (n/d)^(m/n) in ex wrapper class:
+    for (int i=0; i<200; ++i) {  // FIXME: run to ~200
+        for (int j=2; j<13; ++j) {
+            // construct an exponent 1/j...
+            numeric nm(1,j);
+            nm += numeric(int(20.0*rand()/(RAND_MAX+1.0))-10);
+            // ...a numerator...
+            do { i_num = rand(); } while (i_num == 0);
+            numeric num(i_num);
+            // ...and a denominator.
+            do { i_den = (rand())/100; } while (i_den == 0);
+            numeric den(i_den);
+            // construct the radicals:
+            ex radical = pow(ex(num)/ex(den),ex(nm));
+            numeric floating = pow(num/den,nm);
+            // test the result:
+            if (is_ex_of_type(radical,numeric)) {
+                clog << "(" << num << "/" << den << ")^(" << nm
+                     << ") should have been a product, instead it's "
+                     << radical << endl;
+                errorflag = true;
+            }
+            numeric ratio = ex_to_numeric(evalf(radical))/floating;
+            if (ratio>1.0001 && ratio<0.9999) {
+                clog << "(" << num << "/" << den << ")^(" << nm
+                     << ") erroneously evaluated to " << radical;
+                errorflag = true;
+            }
+        }
+    }
+    if (errorflag)
+        ++result;
+    
+    return result;
+}
+
+unsigned check_numeric(void)
+{
+    unsigned result = 0;
+    
+    cout << "checking consistency of numeric types" << flush;
+    clog << "---------consistency of numeric types:" << endl;
+    
+    result += check_numeric1();  cout << '.' << flush;
+    result += check_numeric2();  cout << '.' << flush;
+    
+    if (!result) {
+        cout << " passed " << endl;
+        clog << "(no output)" << endl;
+    } else {
+        cout << " failed " << endl;
+    }
+    
+    return result;
+}
diff --git a/check/checks.cpp b/check/checks.cpp
new file mode 100644 (file)
index 0000000..7a979f0
--- /dev/null
@@ -0,0 +1,79 @@
+/** @file checks.cpp
+ *
+ *  Main program that calls the individual tests. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include <stdexcept>
+#include <iostream>
+#include <time.h>
+
+#include "checks.h"
+
+int main()
+{
+    unsigned result = 0;
+    
+    srand((unsigned)time(NULL));
+    
+    try {
+        for (int i=0; i<1; ++i)
+            result += check_numeric();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        for (int i=0; i<1; ++i)
+            result += check_inifcns();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        for (int i=0; i<1; ++i)
+            result += check_matrices();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        for (int i=0; i<1; ++i)
+            result += check_lsolve();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    if (result) {
+        cout << "Error: something went wrong. ";
+        if (result == 1) {
+            cout << "(one failure)" << endl;
+        } else {
+            cout << "(" << result << " individual failures)" << endl;
+        }
+        cout << "please check check.out against check.ref for more details."
+             << endl << "happy debugging!" << endl;
+    }
+    
+    return result;
+}
similarity index 57%
rename from check/check.h
rename to check/checks.h
index 7bdfe2c..02fb4eb 100644 (file)
@@ -1,4 +1,4 @@
-/** @file check.h
+/** @file checks.h
  *
  *  Prototypes for all individual checks. */
 
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#ifndef CHECK_H
-#define CHECK_H
-
-// fcntimer is defined in timer.cpp and used for timing check functions only:
-unsigned fcntimer(unsigned fcn());
-
-// prototypes for all individual checks must be unsigned fcn() in order to be
-// able to use fcntimer() as a wrapper:
-unsigned paranoia_check();
-unsigned numeric_output();
-unsigned numeric_consist();
-unsigned powerlaws();
-unsigned expand_subs();
-unsigned inifcns_consist();
-unsigned differentiation();
-unsigned poly_gcd();
-unsigned normalization();
-unsigned matrix_checks();
-unsigned linear_solve();
-unsigned series_expansion();
-unsigned lortensor_check();
-
-#endif // ndef CHECK_H
+#ifndef CHECKS_H
+#define CHECKS_H
+
+// For rand() and friends:
+#include <stdlib.h>
+
+#include "ginac.h"
+
+#ifndef NO_NAMESPACE_GINAC
+using namespace GiNaC;
+#endif // ndef NO_NAMESPACE_GINAC
+
+// prototypes for the expression generating functions in:
+const ex dense_univariate_poly(const symbol & x, unsigned degree);
+const ex dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree);
+
+// prototypes for all individual checks should be unsigned fcn():
+unsigned check_numeric();
+unsigned check_inifcns();
+unsigned check_matrices();
+unsigned check_lsolve();
+
+#endif // ndef CHECKS_H
diff --git a/check/checks.ref b/check/checks.ref
new file mode 100644 (file)
index 0000000..a109f51
--- /dev/null
@@ -0,0 +1,8 @@
+---------consistency of numeric types:
+(no output)
+---------consistency of symbolic functions:
+(no output)
+---------symbolic matrix manipulations:
+(no output)
+---------linear solve:
+(no output)
similarity index 90%
rename from check/differentiation.cpp
rename to check/exam_differentiation.cpp
index 307bd59..0b1e7f1 100644 (file)
@@ -1,4 +1,4 @@
-/** @file differentiation.cpp
+/** @file exam_differentiation.cpp
  *
  *  Tests for symbolic differentiation, including various functions. */
 
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 static unsigned check_diff(const ex &e, const symbol &x,
                            const ex &d, unsigned nth=1)
@@ -59,7 +55,7 @@ static unsigned check_diff(const ex &e, const symbol &x,
 }
 
 // Simple (expanded) polynomials
-static unsigned differentiation1(void)
+static unsigned exam_differentiation1(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y");
@@ -111,7 +107,7 @@ static unsigned differentiation1(void)
 }
 
 // Trigonometric functions
-static unsigned differentiation2(void)
+static unsigned exam_differentiation2(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y"), a("a"), b("b");
@@ -160,7 +156,7 @@ static unsigned differentiation2(void)
 }
     
 // exp function
-static unsigned differentiation3(void)
+static unsigned exam_differentiation3(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y"), a("a"), b("b");
@@ -188,7 +184,7 @@ static unsigned differentiation3(void)
 }
 
 // log functions
-static unsigned differentiation4(void)
+static unsigned exam_differentiation4(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y"), a("a"), b("b");
@@ -218,7 +214,7 @@ static unsigned differentiation4(void)
 }
 
 // Functions with two variables
-static unsigned differentiation5(void)
+static unsigned exam_differentiation5(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y"), a("a"), b("b");
@@ -253,7 +249,7 @@ static unsigned differentiation5(void)
 }
 
 // Series
-static unsigned differentiation6(void)
+static unsigned exam_differentiation6(void)
 {
     symbol x("x");
     ex e, d, ed;
@@ -272,25 +268,25 @@ static unsigned differentiation6(void)
     return 0;
 }
 
-unsigned differentiation(void)
+unsigned exam_differentiation(void)
 {
     unsigned result = 0;
     
-    cout << "checking symbolic differentiation..." << flush;
-    clog << "---------symbolic differentiation:" << endl;
+    cout << "examining symbolic differentiation" << flush;
+    clog << "----------symbolic differentiation:" << endl;
     
-    result += differentiation1();
-    result += differentiation2();
-    result += differentiation3();
-    result += differentiation4();
-    result += differentiation5();
-    result += differentiation6();
+    result += exam_differentiation1();  cout << '.' << flush;
+    result += exam_differentiation2();  cout << '.' << flush;
+    result += exam_differentiation3();  cout << '.' << flush;
+    result += exam_differentiation4();  cout << '.' << flush;
+    result += exam_differentiation5();  cout << '.' << flush;
+    result += exam_differentiation6();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     return result;
 }
similarity index 81%
rename from check/linear_solve.cpp
rename to check/exam_lsolve.cpp
index d006beb..dfa71da 100644 (file)
@@ -1,7 +1,6 @@
-/** @file linear_solve.cpp
+/** @file exam_lsolve.cpp
  *
- * These test routines do some simple checks on solving linear systems of
- * symbolic equations. */
+ *  These exams test solving small linear systems of symbolic equations. */
 
 /*
  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
+#include "exams.h"
 
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
-
-static unsigned lsolve1(void)
+static unsigned exam_lsolve1(void)
 {
     // A trivial example.
     unsigned result = 0;
@@ -45,9 +40,9 @@ static unsigned lsolve1(void)
     return result;
 }
 
-static unsigned lsolve2a(void)
+static unsigned exam_lsolve2a(void)
 {
-    // An example from the Maple online-help.
+    // An example from the Maple online help.
     unsigned result = 0;
     symbol a("a"), b("b"), x("x"), y("y");
     lst eqns, vars;
@@ -73,9 +68,9 @@ static unsigned lsolve2a(void)
     return result;
 }
 
-static unsigned lsolve2b(void)
+static unsigned exam_lsolve2b(void)
 {
-    // A boring example from Mathematica's online-help.
+    // A boring example from Mathematica's online help.
     unsigned result = 0;
     symbol x("x"), y("y");
     lst eqns, vars;
@@ -101,9 +96,9 @@ static unsigned lsolve2b(void)
     return result;
 }
 
-static unsigned lsolve2c(void)
+static unsigned exam_lsolve2c(void)
 {
-    // An example from the Maple online-help.
+    // An example from the Maple online help.
     unsigned result = 0;
     symbol x("x"), y("y");
     lst eqns, vars;
@@ -129,23 +124,23 @@ static unsigned lsolve2c(void)
     return result;
 }
 
-unsigned linear_solve(void)
+unsigned exam_lsolve(void)
 {
     unsigned result = 0;
     
-    cout << "checking linear solve..." << flush;
-    clog << "---------linear solve:" << endl;
+    cout << "examining linear solve" << flush;
+    clog << "----------linear solve:" << endl;
     
-    result += lsolve1();
-    result += lsolve2a();
-    result += lsolve2b();
-    result += lsolve2c();
+    result += exam_lsolve1();  cout << '.' << flush;
+    result += exam_lsolve2a();  cout << '.' << flush;
+    result += exam_lsolve2b();  cout << '.' << flush;
+    result += exam_lsolve2c();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     
     return result;
similarity index 86%
rename from check/matrix_checks.cpp
rename to check/exam_matrices.cpp
index 9e907ec..81fed94 100644 (file)
@@ -1,6 +1,6 @@
-/** @file matrix_checks.cpp
+/** @file exam_matrices.cpp
  *
- *  Here we test manipulations on GiNaC's symbolic matrices. */
+ *  Here we examine manipulations on GiNaC's symbolic matrices. */
 
 /*
  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
  */
 
 #include <stdexcept>
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 static unsigned matrix_determinants(void)
 {
     unsigned result = 0;
     ex det;
-    matrix m1(1,1), m2(2,2), m3(3,3);
+    matrix m1(1,1), m2(2,2), m3(3,3), m4(4,4);
     symbol a("a"), b("b"), c("c");
     symbol d("d"), e("e"), f("f");
     symbol g("g"), h("h"), i("i");
@@ -87,6 +83,15 @@ static unsigned matrix_determinants(void)
         ++result;
     }
 
+    // check sparse symbolic 4x4 matrix determinant
+    m4.set(0,1,a).set(1,0,b).set(3,2,c).set(2,3,d);
+    det = m4.determinant();
+    if (det != a*b*c*d) {
+        clog << "determinant of 4x4 matrix " << m4
+             << " erroneously returned " << det << endl;
+        ++result;
+    }
+    
     // check characteristic polynomial
     m3.set(0,0,a).set(0,1,-2).set(0,2,2);
     m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
@@ -204,8 +209,7 @@ static unsigned matrix_misc(void)
     bool caught=false;
     try {
         m5 = inverse(m4);
-    }
-    catch (std::runtime_error err) {
+    } catch (std::runtime_error err) {
         caught=true;
     }
     if (!caught) {
@@ -217,24 +221,24 @@ static unsigned matrix_misc(void)
     return result;
 }
 
-unsigned matrix_checks(void)
+unsigned exam_matrices(void)
 {
     unsigned result = 0;
     
-    cout << "checking symbolic matrix manipulations..." << flush;
-    clog << "---------symbolic matrix manipulations:" << endl;
-    
-    result += matrix_determinants();
-    result += matrix_invert1();
-    result += matrix_invert2();
-    result += matrix_invert3();
-    result += matrix_misc();
+    cout << "examining symbolic matrix manipulations" << flush;
+    clog << "----------symbolic matrix manipulations:" << endl;
+
+    result += matrix_determinants();  cout << '.' << flush;
+    result += matrix_invert1();  cout << '.' << flush;
+    result += matrix_invert2();  cout << '.' << flush;
+    result += matrix_invert3();  cout << '.' << flush;
+    result += matrix_misc();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     
     return result;
similarity index 60%
rename from check/expand_subs.cpp
rename to check/exam_misc.cpp
index 419cea7..86bd12a 100644 (file)
@@ -1,17 +1,6 @@
-/** @file expand_subs.cpp
+/** @file exam_misc.cpp
  *
- *  The first test routine implements Denny Fliegner's quick consistency check:
- *     e = (a0 + a1 + a2 + a3 + ...)^2
- *     expand e
- *     substitute a0 by (-a2 - a3 - ...) in e
- *     expand e
- *  after which e should be just a1^2.
- *  In addition, a simpler modification is tested in the second test:
- *     e = (a0 + a1)^200
- *     expand e
- *     substitute a0 by -a1 in e
- *  after which e should return 0 (without expanding). */
-
+ */
 /*
  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
  *
  */
 
 
-#include "ginac.h"
+#include "exams.h"
 
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
-
-#define VECSIZE 100
-
-static unsigned expand_subs1(void)
+#define VECSIZE 30
+static unsigned exam_expand_subs(void)
 {
+    unsigned result = 0;
     symbol a1("a1");
     symbol a[VECSIZE];
     ex e, aux;
-
+    
     a[1] = a1;
     for (unsigned i=0; i<VECSIZE; ++i) {
         e = e + a[i];
     }
-
+    
     // prepare aux so it will swallow anything but a1^2:
     aux = -e + a[0] + a[1];
     e = expand(subs(expand(pow(e, 2)), a[0] == aux));
-
+    
     if (e != pow(a1,2)) {
         clog << "Denny Fliegner's quick consistency check erroneously returned "
              << e << "." << endl;
-        return 1;
+        ++result;
     }
-    return 0;
+    
+    return result;
 }
 
-static unsigned expand_subs2(void)
+/*  A simple modification of Denny Fliegner's three step consistency test:
+ *  1)  e = (a0 + a1)^200
+ *  2)  expand e
+ *  3)  substitute a0 by -a1 in e
+ *  after which e should return 0 (without expanding). */
+static unsigned exam_expand_subs2(void)
 {
+    unsigned result = 0;
     symbol a("a"), b("b");
     ex e, f;
-
-    // Here the final expand() should be superflous. For no particular reason
-    // at all, we don't use the wrapper-functions but the methods instead:
+    
     e = pow(a+b,200).expand();
     f = e.subs(a == -b);
 
     if (f != 0) {
         clog << "e = pow(a+b,200).expand(); f = e.subs(a == -b); erroneously returned "
              << f << " instead of simplifying to 0." << endl;
-        return 1;
+        ++result;
     }
-    return 0;
+    
+    return result;
 }
 
-unsigned expand_subs(void)
+unsigned exam_misc(void)
 {
     unsigned result = 0;
-
-    cout << "checking commutative expansion and substitution..." << flush;
-    clog << "---------commutative expansion and substitution:" << endl;
     
-    result += expand_subs1();
-    result += expand_subs2();
+    cout << "examining miscellaneous other things" << flush;
+    clog << "----------miscellaneous other things:" << endl;
+    
+    result += exam_expand_subs();  cout << '.' << flush;
+    result += exam_expand_subs2();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
-
+    
     return result;
 }
similarity index 87%
rename from check/lortensor_check.cpp
rename to check/exam_noncommut.cpp
index c153ee4..634fad7 100644 (file)
@@ -1,4 +1,4 @@
-/** @file lortensor_check.cpp
+/** @file exam_noncommut.cpp
  *
  *  Here we test manipulations on GiNaC's lortensors. */
 
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include <stdexcept>
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 static unsigned lortensor_check1(void)
 {
@@ -90,21 +85,21 @@ static unsigned lortensor_check2(void)
     return result;
 }
 
-unsigned lortensor_check(void)
+unsigned exam_noncommut(void)
 {
     unsigned result = 0;
     
-    cout << "checking manipulations of lortensor objects..." << flush;
-    clog << "---------manipulations of lortensor objects:" << endl;
+    cout << "examining behaviour of noncommutative objects" << flush;
+    clog << "----------behaviour of noncommutative objects:" << endl;
     
-    result += lortensor_check1();
-    result += lortensor_check2();
+    result += lortensor_check1();  cout << '.' << flush;
+    result += lortensor_check2();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     
     return result;    
similarity index 69%
rename from check/normalization.cpp
rename to check/exam_normalization.cpp
index 814c45b..98e1d7c 100644 (file)
@@ -1,4 +1,4 @@
-/** @file normalization.cpp
+/** @file exam_normalization.cpp
  *
  *  Rational function normalization test suite. */
 
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 static symbol w("w"), x("x"), y("y"), z("z");
 
@@ -39,7 +35,7 @@ static unsigned check_normal(const ex &e, const ex &d)
     return 0;
 }
 
-static unsigned normal1(void)
+static unsigned exam_normal1(void)
 {
     unsigned result = 0;
     ex e, d;
@@ -62,20 +58,28 @@ static unsigned normal1(void)
     e = pow(x, -1) + x/(x+1);
     d = (pow(x, 2)+x+1)/(x*(x+1));
     result += check_normal(e, d);
-    
-    // Fraction cancellation
-       e = numeric(1)/2 * z * (2*x + 2*y);
-       d = z * (x + y);
-       result += check_normal(e, d);
 
-       e = numeric(1)/6 * z * (3*x + 3*y) * (2*x + 2*w);
-       d = z * (x + y) * (x + w);
-       result += check_normal(e, d);
-
-       e = (3*x + 3*y) * (w/3 + z/3);
-       d = (x + y) * (w + z);
-       result += check_normal(e, d);
+    return result;
+}
 
+static unsigned exam_normal2(void)
+{
+    unsigned result = 0;
+    ex e, d;
+    
+    // Fraction cancellation
+    e = numeric(1)/2 * z * (2*x + 2*y);
+    d = z * (x + y);
+    result += check_normal(e, d);
+    
+    e = numeric(1)/6 * z * (3*x + 3*y) * (2*x + 2*w);
+    d = z * (x + y) * (x + w);
+    result += check_normal(e, d);
+    
+    e = (3*x + 3*y) * (w/3 + z/3);
+    d = (x + y) * (w + z);
+    result += check_normal(e, d);
+    
     e = (pow(x, 2) - pow(y, 2)) / pow(x-y, 3);
     d = (x + y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
     result += check_normal(e, d);
@@ -83,16 +87,24 @@ static unsigned normal1(void)
     e = (pow(x, -1) + x) / (pow(x , 2) * 2 + 2);
     d = pow(x * 2, -1);
     result += check_normal(e, d);
+    
+    // Fraction cancellation with rational coefficients
+    e = (pow(x, 2) - pow(y, 2)) / pow(x/2 - y/2, 3);
+    d = (8 * x + 8 * y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
+    result += check_normal(e, d);
+    
+    // Fraction cancellation with rational coefficients
+    e = z/5 * (x/7 + y/10) / (x/14 + y/20);
+    d = 2*z/5;
+    result += check_normal(e, d);
+    
+    return result;
+}
 
-       // Fraction cancellation with rational coefficients
-       e = (pow(x, 2) - pow(y, 2)) / pow(x/2 - y/2, 3);
-       d = (8 * x + 8 * y) / (pow(x, 2) + pow(y, 2) - x * y * 2);
-       result += check_normal(e, d);
-
-       // Fraction cancellation with rational coefficients
-       e = z/5 * (x/7 + y/10) / (x/14 + y/20);
-       d = 2*z/5;
-       result += check_normal(e, d);
+static unsigned exam_normal3(void)
+{
+    unsigned result = 0;
+    ex e, d;
     
     // Distribution of powers
     e = pow(x/y, 2);
@@ -109,6 +121,14 @@ static unsigned normal1(void)
     d = pow((pow(x, 2) + 1) / x, numeric(1)/2);
     result += check_normal(e, d);
     
+    return result;
+}
+
+static unsigned exam_normal4(void)
+{
+    unsigned result = 0;
+    ex e, d;
+    
     // Replacement of functions with temporary symbols and fraction cancellation
     e = pow(sin(x), 2) - pow(cos(x), 2);
     e /= sin(x) + cos(x);
@@ -137,20 +157,24 @@ static unsigned normal1(void)
     return result;
 }
 
-unsigned normalization(void)
+unsigned exam_normalization(void)
 {
     unsigned result = 0;
     
-    cout << "checking rational function normalization..." << flush;
-    clog << "---------rational function normalization:" << endl;
+    cout << "examining rational function normalization" << flush;
+    clog << "----------rational function normalization:" << endl;
     
-    result += normal1();
+    result += exam_normal1();  cout << '.' << flush;
+    result += exam_normal2();  cout << '.' << flush;
+    result += exam_normal3();  cout << '.' << flush;
+    result += exam_normal4();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
+    
     return result;
 }
similarity index 87%
rename from check/numeric_consist.cpp
rename to check/exam_numeric.cpp
index 5839da1..d3c1264 100644 (file)
@@ -1,7 +1,7 @@
-/** @file numeric_consist.cpp
+/** @file exam_numeric.cpp
  *
- *  This test routine creates some numbers and check the result of several
- *  boolean tests on these numbers like is_integer() etc... */
+ *  These exams creates some numbers and check the result of several boolean
+ *  tests on these numbers like is_integer() etc... */
 
 /*
  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include <stdlib.h>
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 /* Simple and maybe somewhat pointless consistency tests of assorted tests and
  * conversions. */
-static unsigned numeric_consist1(void)
+static unsigned exam_numeric1(void)
 {
     unsigned result = 0;
     numeric test_int1(42);
@@ -107,23 +102,6 @@ static unsigned numeric_consist1(void)
         ++result;
     }
     
-    // Check some numerator and denominator calculations:
-    for (int i=0; i<10; ++i) {
-        int re_q, im_q;
-        do { re_q = rand(); } while (re_q == 0);
-        do { im_q = rand(); } while (im_q == 0);
-        numeric r(rand()-RAND_MAX/2, re_q);
-        numeric i(rand()-RAND_MAX/2, im_q);
-        numeric z = r + I*i;
-        numeric p = numer(z);
-        numeric q = denom(z);
-        numeric res = p/q;
-        if (res != z) {
-            clog << z << " erroneously transformed into " 
-                 << p << "/" << q << " by numer() and denom()" << endl;
-            ++result;
-        }
-    }    
     return result;
 }
 
@@ -132,7 +110,7 @@ static unsigned numeric_consist1(void)
  * Implementing a workaround sadly introduced another bug on May 28th 1999
  * that was fixed on May 31st.  The workaround turned out to be stupid and
  * the original bug in CLN was finally killed on September 2nd. */
-static unsigned numeric_consist2(void)
+static unsigned exam_numeric2(void)
 {
     unsigned result = 0;
     
@@ -171,7 +149,7 @@ static unsigned numeric_consist2(void)
 
 /* Assorted tests to ensure some crucial functions behave exactly as specified
  * in the documentation. */
-static unsigned numeric_consist3(void)
+static unsigned exam_numeric3(void)
 {
     unsigned result = 0;
     numeric calc_rem, calc_quo;
@@ -289,7 +267,7 @@ static unsigned numeric_consist3(void)
 
 /* Now we perform some less trivial checks about several functions which should
  * return exact numbers if possible. */
-static unsigned numeric_consist4(void)
+static unsigned exam_numeric4(void)
 {
     unsigned result = 0;
     bool passed;
@@ -323,23 +301,23 @@ static unsigned numeric_consist4(void)
     return result;
 }
 
-unsigned numeric_consist(void)
+unsigned exam_numeric(void)
 {
     unsigned result = 0;
-
-    cout << "checking consistency of numeric types..." << flush;
-    clog << "---------consistency of numeric types:" << endl;
     
-    result += numeric_consist1();
-    result += numeric_consist2();
-    result += numeric_consist3();
-    result += numeric_consist4();
-
+    cout << "examining consistency of numeric types" << flush;
+    clog << "----------consistency of numeric types:" << endl;
+    
+    result += exam_numeric1();  cout << '.' << flush;
+    result += exam_numeric2();  cout << '.' << flush;
+    result += exam_numeric3();  cout << '.' << flush;
+    result += exam_numeric4();  cout << '.' << flush;
+    
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     
     return result;
similarity index 86%
rename from check/paranoia_check.cpp
rename to check/exam_paranoia.cpp
index 125fb73..7ad2cad 100644 (file)
@@ -1,4 +1,4 @@
-/** @file paranoia_check.cpp
+/** @file exam_paranoia.cpp
  *
  *  This set of tests checks for some of GiNaC's oopses which showed up during
  *  development.  Things were evaluated wrongly and so.  Such a sick behaviour
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 // The very first pair of historic problems had its roots in power.cpp and was
 // finally resolved on April 27th 1999. (Fixing the first on April 23rd
 // actually introduced the second.)
-static unsigned paranoia_check1(void)
+static unsigned exam_paranoia1(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y"), z("z");
@@ -63,7 +59,7 @@ static unsigned paranoia_check1(void)
 // And here the second oops which showed up until May 17th 1999.  It had to do
 // with lexicographic canonicalization and thus showed up only if the variables
 // had the names as given here:
-static unsigned paranoia_check2(void)
+static unsigned exam_paranoia2(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y"), z("z");
@@ -99,7 +95,7 @@ static unsigned paranoia_check2(void)
 // The third bug was introduced on May 18th 1999, discovered on May 19 and
 // fixed that same day.  It worked when x was substituted by 1 but not with
 // other numbers:
-static unsigned paranoia_check3(void)
+static unsigned exam_paranoia3(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y");
@@ -128,7 +124,7 @@ static unsigned paranoia_check3(void)
 }
 
 // The fourth bug was also discovered on May 19th 1999 and fixed immediately:
-static unsigned paranoia_check4(void)
+static unsigned exam_paranoia4(void)
 {
     unsigned result = 0;
     symbol x("x");
@@ -153,7 +149,7 @@ static unsigned paranoia_check4(void)
 }
 
 // The fifth oops was discovered on May 20th 1999 and fixed a day later:
-static unsigned paranoia_check5(void)
+static unsigned exam_paranoia5(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y");
@@ -172,7 +168,7 @@ static unsigned paranoia_check5(void)
 }
 
 // This one was discovered on Jun 1st 1999 and fixed the same day:
-static unsigned paranoia_check6(void)
+static unsigned exam_paranoia6(void)
 {
     unsigned result = 0;
     symbol x("x");
@@ -190,7 +186,7 @@ static unsigned paranoia_check6(void)
 
 // This one was introduced on June 1st 1999 by some aggressive manual
 // optimization. Discovered and fixed on June 2nd.
-static unsigned paranoia_check7(void)
+static unsigned exam_paranoia7(void)
 {
     unsigned result = 0;
     symbol x("x"), y("y");
@@ -210,7 +206,7 @@ static unsigned paranoia_check7(void)
 // This one was a result of the rewrite of mul::max_coefficient when we
 // introduced the overall_coefficient field in expairseq objects on Oct 1st
 // 1999. Fixed on Oct 4th.
-static unsigned paranoia_check8(void)
+static unsigned exam_paranoia8(void)
 {
     unsigned result = 0;
     symbol x("x");
@@ -236,7 +232,7 @@ static unsigned paranoia_check8(void)
 // Z[X]. multiply_lcm() forgot to multiply the x-linear term with the LCM of
 // the coefficient's denominators (2 in this case).  Introduced on Jan 25th
 // 2000 and fixed on Jan 31th.
-static unsigned paranoia_check9(void)
+static unsigned exam_paranoia9(void)
 {
     unsigned result = 0;
     symbol x("x");
@@ -255,7 +251,7 @@ static unsigned paranoia_check9(void)
 // and on Feb 13th 2000 I found out that things like 2^(3/2) throw an exception
 // "power::eval(): pow(0,0) is undefined" instead of simplifying to 2*2^(1/2).
 // It was fixed that same day.
-static unsigned paranoia_check10(void)
+static unsigned exam_paranoia10(void)
 {
     unsigned result = 0;
     
@@ -280,7 +276,7 @@ static unsigned paranoia_check10(void)
 // add::normal() forgot to multiply the denominator of the overall_coeff of
 // its expanded and normalized children with the denominator of the expanded
 // child (did you get this? Well, never mind...). Fixed on Feb 21th 2000.
-static unsigned paranoia_check11(void)
+static unsigned exam_paranoia11(void)
 {
     unsigned result = 0;
        symbol x("x");
@@ -296,31 +292,31 @@ static unsigned paranoia_check11(void)
     return result;
 }
 
-unsigned paranoia_check(void)
+unsigned exam_paranoia(void)
 {
     unsigned result = 0;
-
-    cout << "checking several ex-bugs just out of pure paranoia..." << flush;
-    clog << "---------several ex-bugs just out of pure paranoia:" << endl;
-
-    result += paranoia_check1();
-    result += paranoia_check2();
-    result += paranoia_check3();
-    result += paranoia_check4();
-    result += paranoia_check5();
-    result += paranoia_check6();
-    result += paranoia_check7();
-    result += paranoia_check8();
-    result += paranoia_check9();
-    result += paranoia_check10();
-    result += paranoia_check11();
-
+    
+    cout << "examining several historic failures just out of paranoia" << flush;
+    clog << "----------several historic failures:" << endl;
+    
+    result += exam_paranoia1();  cout << '.' << flush;
+    result += exam_paranoia2();  cout << '.' << flush;
+    result += exam_paranoia3();  cout << '.' << flush;
+    result += exam_paranoia4();  cout << '.' << flush;
+    result += exam_paranoia5();  cout << '.' << flush;
+    result += exam_paranoia6();  cout << '.' << flush;
+    result += exam_paranoia7();  cout << '.' << flush;
+    result += exam_paranoia8();  cout << '.' << flush;
+    result += exam_paranoia9();  cout << '.' << flush;
+    result += exam_paranoia10();  cout << '.' << flush;
+    result += exam_paranoia11();  cout << '.' << flush;
+    
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
-
+    
     return result;
 }
similarity index 87%
rename from check/poly_gcd.cpp
rename to check/exam_polygcd.cpp
index 0e6780c..34aea62 100644 (file)
@@ -1,4 +1,4 @@
-/** @file poly_gcd.cpp
+/** @file exam_polygcd.cpp
  *
  *  Some test with polynomial GCD calculations. See also the checks for
  *  rational function normalization in normalization.cpp. */
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
+#include "exams.h"
 
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
-
-const int MAX_VARIABLES = 5;
+const int MAX_VARIABLES = 3;
 
 static symbol x("x"), z("z");
 static symbol y[MAX_VARIABLES];
@@ -215,7 +211,7 @@ static unsigned poly_gcd7(void)
        ex p = x - y * z + 1;
        ex q = x - y + z * 3;
 
-       for (int j=1; j<=3; j++) {
+       for (int j=1; j<=MAX_VARIABLES; j++) {
                for (int k=j+1; k<=4; k++) {
                        ex d = pow(p, j) * pow(q, j);
                        ex f = pow(p, j) * pow(q, k);
@@ -230,28 +226,29 @@ static unsigned poly_gcd7(void)
        return 0;
 }
 
-unsigned poly_gcd(void)
+unsigned exam_polygcd(void)
 {
     unsigned result = 0;
-
-       cout << "checking polynomial GCD computation..." << flush;
-       clog << "---------polynomial GCD computation:" << endl;
-
-       result += poly_gcd1();
-       result += poly_gcd2();
-       result += poly_gcd3();
-//     result += poly_gcd3p(); // takes extremely long (PRS "worst" case)
-       result += poly_gcd4();
-       result += poly_gcd5();
-       result += poly_gcd5p();
-       result += poly_gcd6();
-       result += poly_gcd7();
-
+    
+       cout << "examining polynomial GCD computation" << flush;
+       clog << "----------polynomial GCD computation:" << endl;
+    
+       result += poly_gcd1();  cout << '.' << flush;
+       result += poly_gcd2();  cout << '.' << flush;
+       result += poly_gcd3();  cout << '.' << flush;
+    result += poly_gcd3p();     cout << '.' << flush; // takes extremely long (PRS "worst" case)
+       result += poly_gcd4();  cout << '.' << flush;
+       result += poly_gcd5();  cout << '.' << flush;
+       result += poly_gcd5p();  cout << '.' << flush;
+       result += poly_gcd6();  cout << '.' << flush;
+       result += poly_gcd7();  cout << '.' << flush;
+    
        if (!result) {
-               cout << " passed ";
+               cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-               cout << " failed ";
+               cout << " failed " << endl;
     }
+    
        return result;
 }
similarity index 93%
rename from check/powerlaws.cpp
rename to check/exam_powerlaws.cpp
index e885369..74e95de 100644 (file)
@@ -1,4 +1,4 @@
-/** @file powerlaws.cpp
+/** @file exam_powerlaws.cpp
  *
  *  Tests for power laws.  You shouldn't try to draw much inspiration from
  *  this code, it is a sanity check rather deeply rooted in GiNaC's classes. */
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
+#include "exams.h"
 
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
-
-static unsigned powerlaws1(void)
+static unsigned exam_powerlaws1(void)
 {
     // (x^a)^b = x^(a*b)
     
@@ -103,7 +99,7 @@ static unsigned powerlaws1(void)
     return 0;
 }
 
-static unsigned powerlaws2(void)
+static unsigned exam_powerlaws2(void)
 {
     // (a*x)^b = a^b * x^b
     
@@ -201,7 +197,7 @@ static unsigned powerlaws2(void)
     return 0;
 }
 
-static unsigned powerlaws3(void)
+static unsigned exam_powerlaws3(void)
 {
     // numeric evaluation
 
@@ -240,7 +236,7 @@ static unsigned powerlaws3(void)
     return 0;
 }
 
-static unsigned powerlaws4(void)
+static unsigned exam_powerlaws4(void)
 {
     // test for mul::eval()
 
@@ -264,23 +260,24 @@ static unsigned powerlaws4(void)
     return 0;
 }
 
-unsigned powerlaws(void)
+unsigned exam_powerlaws(void)
 {
     unsigned result = 0;
     
-    cout << "checking power laws..." << flush;
-    clog << "---------power laws:" << endl;
+    cout << "examining power laws" << flush;
+    clog << "----------power laws:" << endl;
     
-    result += powerlaws1();
-    result += powerlaws2();
-    result += powerlaws3();
-    result += powerlaws4();
+    result += exam_powerlaws1();  cout << '.' << flush;
+    result += exam_powerlaws2();  cout << '.' << flush;
+    result += exam_powerlaws3();  cout << '.' << flush;
+    result += exam_powerlaws4();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
+    
     return result;
 }
similarity index 90%
rename from check/series_expansion.cpp
rename to check/exam_pseries.cpp
index 4cd1697..efcc0fc 100644 (file)
@@ -1,4 +1,4 @@
-/** @file series_expansion.cpp
+/** @file exam_pseries.cpp
  *
  *  Series expansion test (Laurent and Taylor series). */
 
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
+#include "exams.h"
 
 static symbol x("x");
 
@@ -43,7 +39,7 @@ static unsigned check_series(const ex &e, const ex &point, const ex &d, int orde
 }
 
 // Series expansion
-static unsigned series1(void)
+static unsigned exam_series1(void)
 {
     unsigned result = 0;
     ex e, d;
@@ -106,7 +102,7 @@ static unsigned series1(void)
 }
 
 // Series addition
-static unsigned series2(void)
+static unsigned exam_series2(void)
 {
     unsigned result = 0;
     ex e, d;
@@ -119,7 +115,7 @@ static unsigned series2(void)
 }
 
 // Series multiplication
-static unsigned series3(void)
+static unsigned exam_series3(void)
 {
     unsigned result = 0;
     ex e, d;
@@ -132,7 +128,7 @@ static unsigned series3(void)
 }
 
 // Order term handling
-static unsigned series4(void)
+static unsigned exam_series4(void)
 {
     unsigned result = 0;
     ex e, d;
@@ -152,7 +148,7 @@ static unsigned series4(void)
 }
 
 // Series of special functions
-static unsigned series5(void)
+static unsigned exam_series5(void)
 {
     unsigned result = 0;
     ex e, d;
@@ -199,24 +195,24 @@ static unsigned series5(void)
     return result;
 }
 
-unsigned series_expansion(void)
+unsigned exam_pseries(void)
 {
     unsigned result = 0;
     
-    cout << "checking series expansion..." << flush;
-    clog << "---------series expansion:" << endl;
+    cout << "examining series expansion" << flush;
+    clog << "----------series expansion:" << endl;
     
-    result += series1();
-    result += series2();
-    result += series3();
-    result += series4();
-    result += series5();
+    result += exam_series1();  cout << '.' << flush;
+    result += exam_series2();  cout << '.' << flush;
+    result += exam_series3();  cout << '.' << flush;
+    result += exam_series4();  cout << '.' << flush;
+    result += exam_series5();  cout << '.' << flush;
     
     if (!result) {
-        cout << " passed ";
+        cout << " passed " << endl;
         clog << "(no output)" << endl;
     } else {
-        cout << " failed ";
+        cout << " failed " << endl;
     }
     return result;
 }
diff --git a/check/exams.cpp b/check/exams.cpp
new file mode 100644 (file)
index 0000000..31b5cf7
--- /dev/null
@@ -0,0 +1,121 @@
+/** @file exams.cpp
+ *
+ *  Main program that calls all individual exams. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include <stdexcept>
+#include <iostream>
+
+#include "exams.h"
+
+int main()
+{
+    unsigned result = 0;
+    
+    try {
+        result += exam_paranoia();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_numeric();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_powerlaws();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_differentiation();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_polygcd();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_normalization();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_pseries();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_matrices();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_lsolve();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_noncommut();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        result += exam_misc();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    if (result) {
+        cout << "Error: something went wrong. ";
+        if (result == 1) {
+            cout << "(one failure)" << endl;
+        } else {
+            cout << "(" << result << " individual failures)" << endl;
+        }
+        cout << "please check exam.out against exam.ref for more details."
+             << endl << "happy debugging!" << endl;
+    }
+    
+    return result;
+}
diff --git a/check/exams.h b/check/exams.h
new file mode 100644 (file)
index 0000000..3582d53
--- /dev/null
@@ -0,0 +1,45 @@
+/** @file exams.h
+ *
+ *  Prototypes for all individual exams. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef EXAMS_H
+#define EXAMS_H
+
+#include "ginac.h"
+
+#ifndef NO_NAMESPACE_GINAC
+using namespace GiNaC;
+#endif // ndef NO_NAMESPACE_GINAC
+
+// prototypes for all individual checks should be unsigned fcn():
+unsigned exam_paranoia();
+unsigned exam_numeric();
+unsigned exam_powerlaws();
+unsigned exam_differentiation();
+unsigned exam_polygcd();
+unsigned exam_normalization();
+unsigned exam_pseries();
+unsigned exam_matrices();
+unsigned exam_lsolve();
+unsigned exam_noncommut();
+unsigned exam_misc();
+
+#endif // ndef EXAMS_H
diff --git a/check/exams.ref b/check/exams.ref
new file mode 100644 (file)
index 0000000..e15e241
--- /dev/null
@@ -0,0 +1,22 @@
+----------several historic failures:
+(no output)
+----------consistency of numeric types:
+(no output)
+----------power laws:
+(no output)
+----------symbolic differentiation:
+(no output)
+----------polynomial GCD computation:
+(no output)
+----------rational function normalization:
+(no output)
+----------series expansion:
+(no output)
+----------symbolic matrix manipulations:
+(no output)
+----------linear solve:
+(no output)
+----------behaviour of noncommutative objects:
+(no output)
+----------miscellaneous other things:
+(no output)
diff --git a/check/genex.cpp b/check/genex.cpp
new file mode 100644 (file)
index 0000000..77ea18f
--- /dev/null
@@ -0,0 +1,58 @@
+/** @file genex.cpp
+ *
+ *  Provides some routines for generating expressions that are later used as input
+ *  in the consistency checks. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+// For rand() and friends:
+#include <stdlib.h>
+
+#include "ginac.h"
+
+#ifndef NO_NAMESPACE_GINAC
+using namespace GiNaC;
+#endif // ndef NO_NAMESPACE_GINAC
+
+/* Create a dense univariate random polynomial in x.
+ * (of the form 9 - 22*a - 17*a^2 + 14*a^3 + 7*a^4 + 7a^5 if degree==5) */
+const ex
+dense_univariate_poly(const symbol & x, unsigned degree)
+{
+    ex unipoly;
+    
+    for (unsigned i=0; i<=degree; ++i)
+        unipoly += numeric((rand()-RAND_MAX/2))*pow(x,i);
+    
+    return unipoly;
+}
+
+/* Create a dense bivariate random polynomial in x1 and x2.
+ * (of the form 9 + 52*x1 - 27*x1^2 + 84*x2 + 7*x2^2 - 12*x1*x2 if degree ==2) */
+const ex
+dense_bivariate_poly(const symbol & x1, const symbol & x2, unsigned degree)
+{
+    ex bipoly;
+    
+    for (unsigned i1=0; i1<=degree; ++i1)
+        for (unsigned i2=0; i2<=degree-i1; ++i2)
+            bipoly += numeric((rand()-RAND_MAX/2))*pow(x1,i1)*pow(x2,i2);
+    
+    return bipoly;
+}
diff --git a/check/numeric_output.cpp b/check/numeric_output.cpp
deleted file mode 100644 (file)
index 6b3a259..0000000
+++ /dev/null
@@ -1,61 +0,0 @@
-/** @file numeric_output.cpp
- *
- *  Test output of numeric types. */
-
-/*
- *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
- *
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation; either version 2 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
- */
-
-#include "ginac.h"
-
-#ifndef NO_NAMESPACE_GINAC
-using namespace GiNaC;
-#endif // ndef NO_NAMESPACE_GINAC
-
-unsigned numeric_output(void)
-{
-    unsigned result = 0;
-    
-    cout << "checking output of numeric types..." << flush;
-    clog << "---------output of numeric types:" << endl;
-    
-    unsigned long Digits_before = Digits;
-    Digits = 222;
-    clog << "Using " << Digits << " digits" << endl;
-    clog << Pi << " evalfs to: " << Pi.evalf() << endl;
-    clog << Catalan << " evalfs to: " << Catalan.evalf() << endl;
-    clog << EulerGamma << " evalfs to: " << EulerGamma.evalf() << endl;
-    clog << "Complex integers: ";
-    clog << "{(0,0)=" << ex(0 + 0*I) << "} ";
-    clog << "{(1,0)=" << ex(1 + 0*I) << "} ";
-    clog << "{(1,1)=" << ex(1 + 1*I) << "} ";
-    clog << "{(0,1)=" << ex(0 + 1*I) << "} ";
-    clog << "{(-1,1)=" << ex(-1 + 1*I) << "} ";
-    clog << "{(-1,0)=" << ex(-1 + 0*I) << "} ";
-    clog << "{(-1,-1)=" << ex(-1 - 1*I) << "} ";
-    clog << "{(0,-1)=" << ex(0 - 1*I) << "} ";
-    clog << "{(1,-1)=" << ex(1 - 1*I) << "} " << endl;
-    Digits = Digits_before;
-    
-    if (! result) {
-        cout << " passed ";
-    } else {
-        cout << " failed ";
-    }
-    
-    return result;
-}
diff --git a/check/result.ref b/check/result.ref
deleted file mode 100644 (file)
index b4fa587..0000000
+++ /dev/null
@@ -1,30 +0,0 @@
----------several ex-bugs just out of pure paranoia:
-(no output)
----------output of numeric types:
-Using 222 digits
-Pi evalfs to: 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648235L0
-Catalan evalfs to: 0.915965594177219015054603514932384110774149374281672134266498119621763019776254769479356512926115106248574422619196199579035898803325859059431594737481158406995332028773319460519038727478164087865909024706484152163000228727640942388L0
-EulerGamma evalfs to: 0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504014486542836224173997644923536253500333742937337737673942792595258247094916008735203948165670853233151777L0
-Complex integers: {(0,0)=0} {(1,0)=1} {(1,1)=1+I} {(0,1)=I} {(-1,1)=-1+I} {(-1,0)=-1} {(-1,-1)=-1-I} {(0,-1)=-I} {(1,-1)=1-I} 
----------consistency of numeric types:
-(no output)
----------power laws:
-(no output)
----------commutative expansion and substitution:
-(no output)
----------consistency of symbolic functions:
-(no output)
----------symbolic differentiation:
-(no output)
----------polynomial GCD computation:
-(no output)
----------rational function normalization:
-(no output)
----------symbolic matrix manipulations:
-(no output)
----------linear solve:
-(no output)
----------series expansion:
-(no output)
----------manipulations of lortensor objects:
-(no output)
index 8d3adcb..92ba1a2 100755 (executable)
@@ -1,5 +1,4 @@
 #! /bin/sh
-echo "Running checks..."
-./check_ginac 2>result.out
-echo "Comparing output..."
-cmp ${srcdir}/result.ref result.out
+echo "GiNaC will now run through some rather costly consistency checks:"
+./checks 2>checks.out
+cmp ${srcdir}/checks.ref checks.out
diff --git a/check/run_exams b/check/run_exams
new file mode 100755 (executable)
index 0000000..8c8efa5
--- /dev/null
@@ -0,0 +1,4 @@
+#! /bin/sh
+echo "GiNaC will now take an exam with specific input (like a pupils' exam):"
+./exams 2>exams.out
+cmp ${srcdir}/exams.ref exams.out
diff --git a/check/run_times b/check/run_times
new file mode 100755 (executable)
index 0000000..3687f1e
--- /dev/null
@@ -0,0 +1,4 @@
+#! /bin/sh
+echo "GiNaC will now run through some basic timings:"
+./times 2>times.out
+cmp ${srcdir}/times.ref times.out
diff --git a/check/time_dennyfliegner.cpp b/check/time_dennyfliegner.cpp
new file mode 100644 (file)
index 0000000..23c7a95
--- /dev/null
@@ -0,0 +1,96 @@
+/** @file time_dennyfliegner.cpp
+ *
+ *  The first test routine implements Denny Fliegner's quick consistency check:
+ *  1)  e = (a0 + a1 + a2 + a3 + ...)^2, in expanded form
+ *  2)  substitute a0 by (-a2 - a3 - ...) in e
+ *  3)  expand e
+ *  after which e should be just a1^2. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "times.h"
+
+#define VECSIZE 200
+
+static unsigned expand_subs(unsigned size)
+{
+    unsigned result = 0;
+    symbol a1("a1");
+    symbol a[VECSIZE];
+    ex e, aux;
+    
+    a[1] = a1;
+    for (unsigned i=0; i<size; ++i) {
+        e = e + a[i];
+    }
+    
+    // prepare aux so it will swallow anything but a1^2:
+    aux = -e + a[0] + a[1];
+    e = expand(subs(expand(pow(e, 2)), a[0] == aux));
+    
+    if (e != pow(a1,2)) {
+        clog << "Denny Fliegner's quick consistency check erroneously returned "
+             << e << "." << endl;
+        ++result;
+    }
+    
+    return result;
+}
+
+unsigned time_dennyfliegner(void)
+{
+    unsigned result = 0;
+    
+    cout << "timing commutative expansion and substitution" << flush;
+    clog << "-------commutative expansion and substitution:" << endl;
+    
+    vector<unsigned> sizes;
+    vector<double> times;
+    timer rolex;
+    
+    sizes.push_back(40);
+    sizes.push_back(60);
+    sizes.push_back(100);
+    sizes.push_back(150);
+    
+    for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
+        rolex.start();
+        result += expand_subs(*i);  cout << '.' << flush;
+        times.push_back(rolex.read());
+    }
+    
+    if (!result) {
+        cout << " passed ";
+        clog << "(no output)" << endl;
+    } else {
+        cout << " failed ";
+    }
+    // print the report:
+    cout << endl << "    size:  ";
+    for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
+        cout << '\t' << (*i);
+    }
+    cout << endl << "    time/s:";
+    for (vector<double>::iterator i=times.begin(); i!=times.end(); ++i) {
+        cout << '\t' << (*i);
+    }
+    cout << endl;
+    
+    return result;
+}
diff --git a/check/time_gammaseries.cpp b/check/time_gammaseries.cpp
new file mode 100644 (file)
index 0000000..afc9f53
--- /dev/null
@@ -0,0 +1,86 @@
+/** @file time_gammaseries.cpp
+ *
+ *  Some timings on series expansion of the gamma function around a pole. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "times.h"
+
+unsigned gammaseries(unsigned order)
+{
+    unsigned result = 0;
+    symbol x;
+    
+    ex myseries = series(gamma(x),x,0,order);
+    // compute the last coefficient numerically:
+    ex last_coeff = myseries.coeff(x,order-1).evalf();
+    // compute a bound for that coefficient using a variation of the leading
+    // term in Stirling's formula:
+    ex bound = evalf(exp(ex(-.57721566490153286*(order-1)))/(order-1));
+    if (evalf(abs((last_coeff-pow(-1,order))/bound)) > numeric(1)) {
+        clog << "The " << order-1
+             << "th order coefficient in the power series expansion of gamma(0) was erroneously found to be "
+             << last_coeff << ", violating a simple estimate." << endl;
+        ++result;
+    }
+    
+    return result;
+}
+
+unsigned time_gammaseries(void)
+{
+    unsigned result = 0;
+    
+    cout << "timing Laurent series expansion of gamma function" << flush;
+    clog << "-------Laurent series expansion of gamma function:" << endl;
+    
+    vector<unsigned> sizes;
+    vector<double> times;
+    timer omega;
+    
+    sizes.push_back(10);
+    sizes.push_back(15);
+    sizes.push_back(20);
+    sizes.push_back(25);
+    
+    for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
+        omega.start();
+        result += gammaseries(*i);  cout << '.' << flush;
+        times.push_back(omega.read());
+    }
+    
+    if (!result) {
+        cout << " passed ";
+        clog << "(no output)" << endl;
+    } else {
+        cout << " failed ";
+    }
+    // print the report:
+    cout << endl << "    order: ";
+    for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
+        cout << '\t' << (*i);
+    }
+    cout << endl << "    time/s:";
+    for (vector<double>::iterator i=times.begin(); i!=times.end(); ++i) {
+        cout << '\t' << (*i);
+    }
+    cout << endl;
+    
+    return result;
+}
similarity index 63%
rename from check/fcntimer.cpp
rename to check/timer.cpp
index b040c7e..c4c21a3 100644 (file)
@@ -1,6 +1,6 @@
-/** @file fcntimer.cpp
+/** @file timer.cpp
  *
- *  Function execution timer. */
+ *  A simple stop watch class. */
 
 /*
  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
+#include "times.h"
 
-#include <stdio.h>
-#include <sys/resource.h>
+timer::timer(void) : on(false)
+{
+    getrusage(RUSAGE_SELF, &used1);
+    getrusage(RUSAGE_SELF, &used2);
+}
 
+void timer::start(void)
+{
+    on = true;
+    getrusage(RUSAGE_SELF, &used1);
+    getrusage(RUSAGE_SELF, &used2);
+}
 
-// fcntimer() is a little wrapper around GiNaC's automated checks.  All those
-// functions are passed void and return unsigned.  fcntimer() accepts one such
-// function fcn(), returns its result and as a side-effect prints to stdout how
-// much CPU time was consumed by fcn's execution in the fashion "(0.07s)\n".
-unsigned fcntimer(unsigned fcn())
+void timer::stop(void)
 {
-    unsigned fcnresult;
-    struct rusage used1, used2;
-    double elapsed;
+    on = false;
+    getrusage(RUSAGE_SELF, &used2);
+}
 
-    // time the execution of the function:
+void timer::reset(void)
+{
     getrusage(RUSAGE_SELF, &used1);
-    fcnresult = fcn();
     getrusage(RUSAGE_SELF, &used2);
+}
 
-    // add elapsed user and system time in microseconds:
+double timer::read(void)
+{
+    double elapsed;
+    if (this->running())
+        getrusage(RUSAGE_SELF, &used2);
     elapsed = ((used2.ru_utime.tv_sec - used1.ru_utime.tv_sec) +
                (used2.ru_stime.tv_sec - used1.ru_stime.tv_sec) +
                (used2.ru_utime.tv_usec - used1.ru_utime.tv_usec) / 1e6 +
                (used2.ru_stime.tv_usec - used1.ru_stime.tv_usec) / 1e6);
+    // round to 10ms for safety:
+    return 0.01*int(elapsed*100+0.5);
+}
 
-    printf("(%.2fs)\n", elapsed);
-
-    return fcnresult;
+bool timer::running(void)
+{
+    return on;
 }
similarity index 56%
rename from check/main.cpp
rename to check/times.cpp
index d3fbea9..d2b40e8 100644 (file)
@@ -1,6 +1,6 @@
-/** @file main.cpp
+/** @file times.cpp
  *
- *  Main program that calls all individual tests. */
+ *  Main program that calls the individual timings. */
 
 /*
  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
 #include <stdexcept>
 #include <iostream>
 
-#include "check.h"
+#include "times.h"
 
 int main()
 {
     unsigned result = 0;
     
     try {
-        for (int i=0; i<1; ++i) {
-            result += fcntimer(paranoia_check);
-            result += fcntimer(numeric_output);
-            result += fcntimer(numeric_consist);
-            result += fcntimer(powerlaws);
-            result += fcntimer(expand_subs);
-            result += fcntimer(inifcns_consist);
-            result += fcntimer(differentiation);
-            result += fcntimer(poly_gcd);
-            result += fcntimer(normalization);
-            result += fcntimer(matrix_checks);
-            result += fcntimer(linear_solve);
-            result += fcntimer(series_expansion);
-            result += fcntimer(lortensor_check);
-        }
+        for (int i=0; i<1; ++i)
+            result += time_dennyfliegner();
     } catch (const exception &e) {
-        cout << "error: caught an exception: " << e.what() << endl;
-        result++;
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
+    }
+    
+    try {
+        for (int i=0; i<1; ++i)
+            result += time_gammaseries();
+    } catch (const exception &e) {
+        cout << "Error: caught exception " << e.what() << endl;
+        ++result;
     }
     
     if (result) {
-        cout << "error: something went wrong. ";
+        cout << "Error: something went wrong. ";
         if (result == 1) {
             cout << "(one failure)" << endl;
         } else {
             cout << "(" << result << " individual failures)" << endl;
         }
-        cout << "please check result.out against result.ref for more details."
+        cout << "please check times.out against times.ref for more details."
              << endl << "happy debugging!" << endl;
     }
-
+    
     return result;
 }
diff --git a/check/times.h b/check/times.h
new file mode 100644 (file)
index 0000000..bb54399
--- /dev/null
@@ -0,0 +1,53 @@
+/** @file times.h
+ *
+ *  Prototypes for all individual timings. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef CHECKS_H
+#define CHECKS_H
+
+#include <sys/resource.h>
+#include <stdlib.h>
+#include <vector>
+
+#include "ginac.h"
+
+#ifndef NO_NAMESPACE_GINAC
+using namespace GiNaC;
+#endif // ndef NO_NAMESPACE_GINAC
+
+class timer {
+public:
+    timer();
+    void start(void);
+    void stop(void);
+    void reset(void);
+    double read(void);
+    bool running(void);
+private:
+    bool on;
+    struct rusage used1, used2;
+};
+
+// prototypes for all individual timings should be unsigned fcn():
+unsigned time_dennyfliegner();
+unsigned time_gammaseries();
+
+#endif // ndef CHECKS_H
diff --git a/check/times.ref b/check/times.ref
new file mode 100644 (file)
index 0000000..4aef1c9
--- /dev/null
@@ -0,0 +1,4 @@
+-------commutative expansion and substitution:
+(no output)
+-------Laurent series expansion of gamma function:
+(no output)
index 43dc2d8..a494a33 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
index f4d7fb1..8709ddc 100644 (file)
@@ -18,6 +18,9 @@
 /* The number of bytes in a long long.  */
 #undef SIZEOF_LONG_LONG
 
+/* Define if you have the tgetent function.  */
+#undef HAVE_TGETENT
+
 /* Define if you have the <algorithm> header file.  */
 #undef HAVE_ALGORITHM
 
index e153dbe..f7d7934 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
index 3ed3518..021b10e 100644 (file)
@@ -87,6 +87,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
index 0964c1a..ddb7149 100644 (file)
@@ -86,6 +86,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
index 0408108..e4eca61 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
index bc401f4..55da24d 100644 (file)
@@ -212,10 +212,10 @@ ex lsolve(const ex &eqns, const ex &symbols)
     matrix vars(symbols.nops(),1);
     
     for (unsigned r=0; r<eqns.nops(); r++) {
-        ex eq=eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
-        ex linpart=eq;
+        ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
+        ex linpart = eq;
         for (unsigned c=0; c<symbols.nops(); c++) {
-            ex co=eq.coeff(ex_to_symbol(symbols.op(c)),1);
+            ex co = eq.coeff(ex_to_symbol(symbols.op(c)),1);
             linpart -= co*symbols.op(c);
             sys.set(r,c,co);
         }
@@ -226,18 +226,16 @@ ex lsolve(const ex &eqns, const ex &symbols)
     // test if system is linear and fill vars matrix
     for (unsigned i=0; i<symbols.nops(); i++) {
         vars.set(i,0,symbols.op(i));
-        if (sys.has(symbols.op(i))) {
+        if (sys.has(symbols.op(i)))
             throw(std::logic_error("lsolve: system is not linear"));
-        }
-        if (rhs.has(symbols.op(i))) {
+        if (rhs.has(symbols.op(i)))
             throw(std::logic_error("lsolve: system is not linear"));
-        }
     }
     
     //matrix solution=sys.solve(rhs);
     matrix solution;
     try {
-        solution=sys.fraction_free_elim(vars,rhs);
+        solution = sys.fraction_free_elim(vars,rhs);
     } catch (const runtime_error & e) {
         // probably singular matrix (or other error)
         // return empty solution list
index e15feb6..8fbe862 100644 (file)
  */
 
 #include <algorithm>
+#include <map>
 #include <stdexcept>
 
 #include "matrix.h"
 #include "archive.h"
 #include "utils.h"
 #include "debugmsg.h"
+#include "numeric.h"
 
 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
@@ -301,23 +303,21 @@ ex matrix::evalf(int level) const
 int matrix::compare_same_type(const basic & other) const
 {
     GINAC_ASSERT(is_exactly_of_type(other, matrix));
-    const matrix & o=static_cast<matrix &>(const_cast<basic &>(other));
+    const matrix & o = static_cast<matrix &>(const_cast<basic &>(other));
     
     // compare number of rows
-    if (row != o.rows()) {
+    if (row != o.rows())
         return row < o.rows() ? -1 : 1;
-    }
     
     // compare number of columns
-    if (col != o.cols()) {
+    if (col != o.cols())
         return col < o.cols() ? -1 : 1;
-    }
     
     // equal number of rows and columns, compare individual elements
     int cmpval;
     for (unsigned r=0; r<row; ++r) {
         for (unsigned c=0; c<col; ++c) {
-            cmpval=((*this)(r,c)).compare(o(r,c));
+            cmpval = ((*this)(r,c)).compare(o(r,c));
             if (cmpval!=0) return cmpval;
         }
     }
@@ -415,7 +415,7 @@ matrix & matrix::set(unsigned ro, unsigned co, ex value)
     }
     
     ensure_if_modifiable();
-    m[ro*col+co]=value;
+    m[ro*col+co] = value;
     return *this;
 }
 
@@ -425,135 +425,18 @@ matrix matrix::transpose(void) const
 {
     exvector trans(col*row);
     
-    for (unsigned r=0; r<col; ++r) {
-        for (unsigned c=0; c<row; ++c) {
+    for (unsigned r=0; r<col; ++r)
+        for (unsigned c=0; c<row; ++c)
             trans[r*row+c] = m[c*col+r];
-        }
-    }
-    return matrix(col,row,trans);
-}
-
-/* Determiant of purely numeric matrix, using pivoting. This routine is only
- * called internally by matrix::determinant(). */
-ex determinant_numeric(const matrix & M)
-{
-    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
-    matrix tmp(M);
-    ex det=_ex1();
-    ex piv;
-    
-    for (unsigned r1=0; r1<M.rows(); ++r1) {
-        int indx = tmp.pivot(r1);
-        if (indx == -1) {
-            return _ex0();
-        }
-        if (indx != 0) {
-            det *= _ex_1();
-        }
-        det = det * tmp.m[r1*M.cols()+r1];
-        for (unsigned r2=r1+1; r2<M.rows(); ++r2) {
-            piv = tmp.m[r2*M.cols()+r1] / tmp.m[r1*M.cols()+r1];
-            for (unsigned c=r1+1; c<M.cols(); c++) {
-                tmp.m[r2*M.cols()+c] -= piv * tmp.m[r1*M.cols()+c];
-            }
-        }
-    }
-    return det;
-}
-
-// Compute the sign of a permutation of a vector of things, used internally
-// by determinant_symbolic_perm() where it is instantiated for int.
-template <typename T>
-int permutation_sign(vector<T> s)
-{
-    if (s.size() < 2)
-        return 0;
-    int sigma=1;
-    for (typename vector<T>::iterator i=s.begin(); i!=s.end()-1; ++i) {
-        for (typename vector<T>::iterator j=i+1; j!=s.end(); ++j) {
-            if (*i == *j)
-                return 0;
-            if (*i > *j) {
-                iter_swap(i,j);
-                sigma = -sigma;
-            }
-        }
-    }
-    return sigma;
-}
-
-/** Determinant built by application of the full permutation group. This
- *  routine is only called internally by matrix::determinant(). */
-ex determinant_symbolic_perm(const matrix & M)
-{
-    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
-    
-    if (M.rows()==1) {  // speed things up
-        return M(0,0);
-    }
-    
-    ex det;
-    ex term;
-    vector<unsigned> sigma(M.cols());
-    for (unsigned i=0; i<M.cols(); ++i) sigma[i]=i;
-    
-    do {
-        term = M(sigma[0],0);
-        for (unsigned i=1; i<M.cols(); ++i) term *= M(sigma[i],i);
-        det += permutation_sign(sigma)*term;
-    } while (next_permutation(sigma.begin(), sigma.end()));
-    
-    return det;
-}
-
-/** Recursive determiant for small matrices having at least one symbolic entry.
- *  This algorithm is also known as Laplace-expansion. This routine is only
- *  called internally by matrix::determinant(). */
-ex determinant_symbolic_minor(const matrix & M)
-{
-    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
-    
-    if (M.rows()==1) {  // end of recursion
-        return M(0,0);
-    }
-    if (M.rows()==2) {  // speed things up
-        return (M(0,0)*M(1,1)-
-                M(1,0)*M(0,1));
-    }
-    if (M.rows()==3) {  // speed things up even a little more
-        return ((M(2,1)*M(0,2)-M(2,2)*M(0,1))*M(1,0)+
-                (M(1,2)*M(0,1)-M(1,1)*M(0,2))*M(2,0)+
-                (M(2,2)*M(1,1)-M(2,1)*M(1,2))*M(0,0));
-    }
     
-    ex det;
-    matrix minorM(M.rows()-1,M.cols()-1);
-    for (unsigned r1=0; r1<M.rows(); ++r1) {
-        // assemble the minor matrix
-        for (unsigned r=0; r<minorM.rows(); ++r) {
-            for (unsigned c=0; c<minorM.cols(); ++c) {
-                if (r<r1) {
-                    minorM.set(r,c,M(r,c+1));
-                } else {
-                    minorM.set(r,c,M(r+1,c+1));
-                }
-            }
-        }
-        // recurse down
-        if (r1%2) {
-            det -= M(r1,0) * determinant_symbolic_minor(minorM);
-        } else {
-            det += M(r1,0) * determinant_symbolic_minor(minorM);
-        }
-    }
-    return det;
+    return matrix(col,row,trans);
 }
 
 /*  Leverrier algorithm for large matrices having at least one symbolic entry.
  *  This routine is only called internally by matrix::determinant(). The
- *  algorithm is deemed bad for symbolic matrices since it returns expressions
- *  that are very hard to canonicalize. */
-/*ex determinant_symbolic_leverrier(const matrix & M)
+ *  algorithm is very bad for symbolic matrices since it returns expressions
+ *  that are quite hard to expand. */
+/*ex matrix::determinant_symbolic_leverrier(const matrix & M)
  *{
  *    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
  *    
@@ -595,15 +478,15 @@ ex matrix::determinant(bool normalized) const
     // check, if there are non-numeric entries in the matrix:
     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
         if (!(*r).info(info_flags::numeric)) {
-            if (normalized) {
-                return determinant_symbolic_minor(*this).normal();
-            } else {
-                return determinant_symbolic_perm(*this);
-            }
+            if (normalized)
+                // return determinant_symbolic_minor().normal();
+                return determinant_symbolic_minor().normal();
+            else
+                return determinant_symbolic_perm();
         }
     }
     // if it turns out that all elements are numeric
-    return determinant_numeric(*this);
+    return determinant_numeric();
 }
 
 /** Trace of a matrix.
@@ -617,9 +500,9 @@ ex matrix::trace(void) const
     }
     
     ex tr;
-    for (unsigned r=0; r<col; ++r) {
+    for (unsigned r=0; r<col; ++r)
         tr += m[r*col+r];
-    }
+
     return tr;
 }
 
@@ -638,9 +521,9 @@ ex matrix::charpoly(const ex & lambda) const
     }
     
     matrix M(*this);
-    for (unsigned r=0; r<col; ++r) {
+    for (unsigned r=0; r<col; ++r)
         M.m[r*col+r] -= lambda;
-    }
+    
     return (M.determinant());
 }
 
@@ -657,9 +540,9 @@ matrix matrix::inverse(void) const
     
     matrix tmp(row,col);
     // set tmp to the unit matrix
-    for (unsigned i=0; i<col; ++i) {
+    for (unsigned i=0; i<col; ++i)
         tmp.m[i*col+i] = _ex1();
-    }
+
     // create a copy of this matrix
     matrix cpy(*this);
     for (unsigned r1=0; r1<row; ++r1) {
@@ -690,27 +573,30 @@ matrix matrix::inverse(void) const
     return tmp;
 }
 
+// superfluous helper function
 void matrix::ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
 {
     ensure_if_modifiable();
     
-    ex tmp=ffe_get(r1,c1);
+    ex tmp = ffe_get(r1,c1);
     ffe_set(r1,c1,ffe_get(r2,c2));
     ffe_set(r2,c2,tmp);
 }
 
+// superfluous helper function
 void matrix::ffe_set(unsigned r, unsigned c, ex e)
 {
     set(r-1,c-1,e);
 }
 
+// superfluous helper function
 ex matrix::ffe_get(unsigned r, unsigned c) const
 {
     return operator()(r-1,c-1);
 }
 
 /** Solve a set of equations for an m x n matrix by fraction-free Gaussian
- *  elimination. Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
+ *  elimination.  Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
  *  by Keith O. Geddes et al.
  *
  *  @param vars n x p matrix
@@ -720,25 +606,19 @@ ex matrix::ffe_get(unsigned r, unsigned c) const
 matrix matrix::fraction_free_elim(const matrix & vars,
                                   const matrix & rhs) const
 {
-    if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col)) {
-        throw (std::logic_error("matrix::solve(): incompatible matrices"));
-    }
+    // FIXME: implement a Sasaki-Murao scheme which avoids division at all!
+    if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
+        throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
     
-    matrix a(*this); // make a copy of the matrix
-    matrix b(rhs);     // make a copy of the rhs vector
-
-    /*
-    cout << "before" << endl;
-    cout << "a=" << a << endl;
-    cout << "b=" << b << endl;
-    */
+    matrix a(*this);  // make a copy of the matrix
+    matrix b(rhs);    // make a copy of the rhs vector
     
     // given an m x n matrix a, reduce it to upper echelon form
-    unsigned m=a.row;
-    unsigned n=a.col;
-    int sign=1;
-    ex divisor=1;
-    unsigned r=1;
+    unsigned m = a.row;
+    unsigned n = a.col;
+    int sign = 1;
+    ex divisor = 1;
+    unsigned r = 1;
     
     // eliminate below row r, with pivot in column k
     for (unsigned k=1; (k<=n)&&(r<=m); ++k) {
@@ -749,12 +629,11 @@ matrix matrix::fraction_free_elim(const matrix & vars,
         if (p<=m) {
             if (p!=r) {
                 // switch rows p and r
-                for (unsigned j=k; j<=n; ++j) {
+                for (unsigned j=k; j<=n; ++j)
                     a.ffe_swap(p,j,r,j);
-                }
                 b.ffe_swap(p,1,r,1);
                 // keep track of sign changes due to row exchange
-                sign=-sign;
+                sign = -sign;
             }
             for (unsigned i=r+1; i<=m; ++i) {
                 for (unsigned j=k+1; j<=n; ++j) {
@@ -767,12 +646,12 @@ matrix matrix::fraction_free_elim(const matrix & vars,
                 b.ffe_set(i,1,b.ffe_get(i,1).normal() /*.normal() */ );
                 a.ffe_set(i,k,0);
             }
-            divisor=a.ffe_get(r,k);
+            divisor = a.ffe_get(r,k);
             r++;
         }
     }
     // optionally compute the determinant for square or augmented matrices
-    // if (r==m+1) { det=sign*divisor; } else { det=0; }
+    // if (r==m+1) { det = sign*divisor; } else { det = 0; }
     
     /*
     for (unsigned r=1; r<=m; ++r) {
@@ -785,21 +664,20 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     
 #ifdef DO_GINAC_ASSERT
     // test if we really have an upper echelon matrix
-    int zero_in_last_row=-1;
+    int zero_in_last_row = -1;
     for (unsigned r=1; r<=m; ++r) {
         int zero_in_this_row=0;
         for (unsigned c=1; c<=n; ++c) {
-            if (a.ffe_get(r,c).is_equal(_ex0())) {
+            if (a.ffe_get(r,c).is_equal(_ex0()))
                zero_in_this_row++;
-            } else {
+            else
                 break;
-            }
         }
         GINAC_ASSERT((zero_in_this_row>zero_in_last_row)||(zero_in_this_row=n));
-        zero_in_last_row=zero_in_this_row;
+        zero_in_last_row = zero_in_this_row;
     }
 #endif // def DO_GINAC_ASSERT
-
+    
     /*
     cout << "after" << endl;
     cout << "a=" << a << endl;
@@ -808,12 +686,11 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     
     // assemble solution
     matrix sol(n,1);
-    unsigned last_assigned_sol=n+1;
+    unsigned last_assigned_sol = n+1;
     for (unsigned r=m; r>0; --r) {
-        unsigned first_non_zero=1;
-        while ((first_non_zero<=n)&&(a.ffe_get(r,first_non_zero).is_zero())) {
+        unsigned first_non_zero = 1;
+        while ((first_non_zero<=n)&&(a.ffe_get(r,first_non_zero).is_zero()))
             first_non_zero++;
-        }
         if (first_non_zero>n) {
             // row consists only of zeroes, corresponding rhs must be 0 as well
             if (!b.ffe_get(r,1).is_zero()) {
@@ -825,36 +702,26 @@ matrix matrix::fraction_free_elim(const matrix & vars,
             for (unsigned c=first_non_zero+1; c<=last_assigned_sol-1; ++c) {
                 sol.ffe_set(c,1,vars.ffe_get(c,1));
             }
-            ex e=b.ffe_get(r,1);
+            ex e = b.ffe_get(r,1);
             for (unsigned c=first_non_zero+1; c<=n; ++c) {
                 e=e-a.ffe_get(r,c)*sol.ffe_get(c,1);
             }
             sol.ffe_set(first_non_zero,1,
                         (e/a.ffe_get(r,first_non_zero)).normal());
-            last_assigned_sol=first_non_zero;
+            last_assigned_sol = first_non_zero;
         }
     }
     // assign solutions for vars between 1 and
     // last_assigned_sol-1: free parameters
-    for (unsigned c=1; c<=last_assigned_sol-1; ++c) {
+    for (unsigned c=1; c<=last_assigned_sol-1; ++c)
         sol.ffe_set(c,1,vars.ffe_get(c,1));
-    }
-
-    /*
-    for (unsigned c=1; c<=n; ++c) {
-        cout << vars.ffe_get(c,1) << "->" << sol.ffe_get(c,1) << endl;
-    }
-    */
-
-    // cout << "sol=" << sol << endl;
     
 #ifdef DO_GINAC_ASSERT
     // test solution with echelon matrix
     for (unsigned r=1; r<=m; ++r) {
-        ex e=0;
-        for (unsigned c=1; c<=n; ++c) {
-            e=e+a.ffe_get(r,c)*sol.ffe_get(c,1);
-        }
+        ex e = 0;
+        for (unsigned c=1; c<=n; ++c)
+            e = e+a.ffe_get(r,c)*sol.ffe_get(c,1);
         if (!(e-b.ffe_get(r,1)).normal().is_zero()) {
             cout << "e=" << e;
             cout << "b.ffe_get(" << r<<",1)=" << b.ffe_get(r,1) << endl;
@@ -862,25 +729,24 @@ matrix matrix::fraction_free_elim(const matrix & vars,
         }
         GINAC_ASSERT((e-b.ffe_get(r,1)).normal().is_zero());
     }
-
+    
     // test solution with original matrix
     for (unsigned r=1; r<=m; ++r) {
-        ex e=0;
-        for (unsigned c=1; c<=n; ++c) {
-            e=e+ffe_get(r,c)*sol.ffe_get(c,1);
-        }
+        ex e = 0;
+        for (unsigned c=1; c<=n; ++c)
+            e = e+ffe_get(r,c)*sol.ffe_get(c,1);
         try {
-        if (!(e-rhs.ffe_get(r,1)).normal().is_zero()) {
-            cout << "e=" << e << endl;
-            e.printtree(cout);
-            ex en=e.normal();
-            cout << "e.normal()=" << en << endl;
-            en.printtree(cout);
-            cout << "rhs.ffe_get(" << r<<",1)=" << rhs.ffe_get(r,1) << endl;
-            cout << "diff=" << (e-rhs.ffe_get(r,1)).normal() << endl;
-        }
+            if (!(e-rhs.ffe_get(r,1)).normal().is_zero()) {
+                cout << "e=" << e << endl;
+                e.printtree(cout);
+                ex en = e.normal();
+                cout << "e.normal()=" << en << endl;
+                en.printtree(cout);
+                cout << "rhs.ffe_get(" << r<<",1)=" << rhs.ffe_get(r,1) << endl;
+                cout << "diff=" << (e-rhs.ffe_get(r,1)).normal() << endl;
+            }
         } catch (...) {
-            ex xxx=e-rhs.ffe_get(r,1);
+            ex xxx = e - rhs.ffe_get(r,1);
             cerr << "xxx=" << xxx << endl << endl;
         }
         GINAC_ASSERT((e-rhs.ffe_get(r,1)).normal().is_zero());
@@ -888,76 +754,281 @@ matrix matrix::fraction_free_elim(const matrix & vars,
 #endif // def DO_GINAC_ASSERT
     
     return sol;
-}   
+}
+
+/** Solve a set of equations for an m x n matrix.
+ *
+ *  @param vars n x p matrix
+ *  @param rhs m x p matrix
+ *  @exception logic_error (incompatible matrices)
+ *  @exception runtime_error (singular matrix) */
+matrix matrix::solve(const matrix & vars,
+                     const matrix & rhs) const
+{
+    if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
+        throw (std::logic_error("matrix::solve(): incompatible matrices"));
     
-/** Solve simultaneous set of equations. */
-matrix matrix::solve(const matrix & v) const
+    throw (std::runtime_error("FIXME: need implementation."));
+}
+
+/** Old and obsolete interface: */
+matrix matrix::old_solve(const matrix & v) const
 {
-    if (!(row == col && col == v.row)) {
+    if ((v.row != col) || (col != v.row))
         throw (std::logic_error("matrix::solve(): incompatible matrices"));
-    }
     
-    // build the extended matrix of *this with v attached to the right
+    // build the augmented matrix of *this with v attached to the right
     matrix tmp(row,col+v.col);
     for (unsigned r=0; r<row; ++r) {
-        for (unsigned c=0; c<col; ++c) {
-            tmp.m[r*tmp.col+c] = m[r*col+c];
-        }
-        for (unsigned c=0; c<v.col; ++c) {
+        for (unsigned c=0; c<col; ++c)
+            tmp.m[r*tmp.col+c] = this->m[r*col+c];
+        for (unsigned c=0; c<v.col; ++c)
             tmp.m[r*tmp.col+c+col] = v.m[r*v.col+c];
+    }
+    // cout << "augmented: " << tmp << endl;
+    tmp.gauss_elimination();
+    // cout << "degaussed: " << tmp << endl;
+    // assemble the solution matrix
+    exvector sol(v.row*v.col);
+    for (unsigned c=0; c<v.col; ++c) {
+        for (unsigned r=row; r>0; --r) {
+            for (unsigned i=r; i<col; ++i)
+                sol[(r-1)*v.col+c] -= tmp.m[(r-1)*tmp.col+i]*sol[i*v.col+c];
+            sol[(r-1)*v.col+c] += tmp.m[(r-1)*tmp.col+col+c];
+            sol[(r-1)*v.col+c] = (sol[(r-1)*v.col+c]/tmp.m[(r-1)*tmp.col+(r-1)]).normal();
         }
     }
+    return matrix(v.row, v.col, sol);
+}
+
+// protected
+
+/** Determinant of purely numeric matrix, using pivoting.
+ *
+ *  @see       matrix::determinant() */
+ex matrix::determinant_numeric(void) const
+{
+    matrix tmp(*this);
+    ex det = _ex1();
+    ex piv;
+    
     for (unsigned r1=0; r1<row; ++r1) {
         int indx = tmp.pivot(r1);
-        if (indx == -1) {
-            throw (std::runtime_error("matrix::solve(): singular matrix"));
-        }
-        for (unsigned c=r1; c<tmp.col; ++c) {
-            tmp.m[r1*tmp.col+c] /= tmp.m[r1*tmp.col+r1];
-        }
+        if (indx == -1)
+            return _ex0();
+        if (indx != 0)
+            det *= _ex_1();
+        det = det * tmp.m[r1*col+r1];
         for (unsigned r2=r1+1; r2<row; ++r2) {
-            for (unsigned c=r1; c<tmp.col; ++c) {
-                tmp.m[r2*tmp.col+c]
-                    -= tmp.m[r2*tmp.col+r1] * tmp.m[r1*tmp.col+c];
+            piv = tmp.m[r2*col+r1] / tmp.m[r1*col+r1];
+            for (unsigned c=r1+1; c<col; c++) {
+                tmp.m[r2*col+c] -= piv * tmp.m[r1*col+c];
             }
         }
     }
+    return det;
+}
+
+/** Recursive determinant for small matrices having at least one symbolic
+ *  entry.  The basic algorithm, known as Laplace-expansion, is enhanced by
+ *  some bookkeeping to avoid calculation of the same submatrices ("minors")
+ *  more than once.  According to W.M.Gentleman and S.C.Johnson this algorithm
+ *  is better than elimination schemes for sparse multivariate polynomials and
+ *  also for dense univariate polynomials once the dimesion becomes larger
+ *  than 7.
+ *
+ *  @see matrix::determinant() */
+ex matrix::determinant_symbolic_minor(void) const
+{
+    // for small matrices the algorithm does not make sense:
+    if (this->row==1)
+        return m[0];
+    if (this->row==2)
+        return (m[0]*m[3]-m[2]*m[1]);
+    if (this->row==3)
+        return ((m[4]*m[8]-m[5]*m[7])*m[0]-
+                (m[1]*m[8]-m[2]*m[7])*m[3]+
+                (m[1]*m[5]-m[4]*m[2])*m[6]);
     
-    // assemble the solution matrix
-    exvector sol(v.row*v.col);
-    for (unsigned c=0; c<v.col; ++c) {
-        for (unsigned r=col-1; r>=0; --r) {
-            sol[r*v.col+c] = tmp[r*tmp.col+c];
-            for (unsigned i=r+1; i<col; ++i) {
-                sol[r*v.col+c]
-                    -= tmp[r*tmp.col+i] * sol[i*v.col+c];
+    // This algorithm can best be understood by looking at a naive
+    // implementation of Laplace-expansion, like this one:
+    // ex det;
+    // matrix minorM(this->row-1,this->col-1);
+    // for (unsigned r1=0; r1<this->row; ++r1) {
+    //     // shortcut if element(r1,0) vanishes
+    //     if (m[r1*col].is_zero())
+    //         continue;
+    //     // assemble the minor matrix
+    //     for (unsigned r=0; r<minorM.rows(); ++r) {
+    //         for (unsigned c=0; c<minorM.cols(); ++c) {
+    //             if (r<r1)
+    //                 minorM.set(r,c,m[r*col+c+1]);
+    //             else
+    //                 minorM.set(r,c,m[(r+1)*col+c+1]);
+    //         }
+    //     }
+    //     // recurse down and care for sign:
+    //     if (r1%2)
+    //         det -= m[r1*col] * minorM.determinant_symbolic_minor();
+    //     else
+    //         det += m[r1*col] * minorM.determinant_symbolic_minor();
+    // }
+    // return det;
+    // What happens is that while proceeding down many of the minors are
+    // computed more than once.  In particular, there are binomial(n,k)
+    // kxk minors and each one is computed factorial(n-k) times.  Therefore
+    // it is reasonable to store the results of the minors.  We proceed from
+    // right to left.  At each column c we only need to retrieve the minors
+    // calculated in step c-1.  We therefore only have to store at most 
+    // 2*binomial(n,n/2) minors.
+    
+    // we store our subminors in maps, keys being the rows they arise from
+    typedef map<vector<unsigned>,class ex> Rmap;
+    typedef map<vector<unsigned>,class ex>::value_type Rmap_value;
+    Rmap A, B;
+    ex det;
+    vector<unsigned> Pkey;    // Unique flipper counter for partitioning into minors
+    Pkey.reserve(this->col);
+    vector<unsigned> Mkey;    // key for minor determinant (a subpartition of Pkey)
+    Mkey.reserve(this->col-1);
+    // initialize A with last column:
+    for (unsigned r=0; r<this->col; ++r) {
+        Pkey.erase(Pkey.begin(),Pkey.end());
+        Pkey.push_back(r);
+        A.insert(Rmap_value(Pkey,m[this->col*r+this->col-1]));
+    }
+    // proceed from right to left through matrix
+    for (int c=this->col-2; c>=0; --c) {
+        Pkey.erase(Pkey.begin(),Pkey.end());  // don't change capacity
+        Mkey.erase(Mkey.begin(),Mkey.end());
+        for (unsigned i=0; i<this->col-c; ++i)
+            Pkey.push_back(i);
+        unsigned fc = 0;  // controls logic for our strange flipper counter
+        do {
+            A.insert(Rmap_value(Pkey,_ex0()));
+            det = _ex0();
+            for (unsigned r=0; r<this->col-c; ++r) {
+                // maybe there is nothing to do?
+                if (m[Pkey[r]*this->col+c].is_zero())
+                    continue;
+                // create the sorted key for all possible minors
+                Mkey.erase(Mkey.begin(),Mkey.end());
+                for (unsigned i=0; i<this->col-c; ++i)
+                    if (i!=r)
+                        Mkey.push_back(Pkey[i]);
+                // Fetch the minors and compute the new determinant
+                if (r%2)
+                    det -= m[Pkey[r]*this->col+c]*A[Mkey];
+                else
+                    det += m[Pkey[r]*this->col+c]*A[Mkey];
             }
-        }
+            // Store the new determinant at its place in B:
+            B.insert(Rmap_value(Pkey,det));
+            // increment our strange flipper counter
+            for (fc=this->col-c; fc>0; --fc) {
+                ++Pkey[fc-1];
+                if (Pkey[fc-1]<fc+c)
+                    break;
+            }
+            if (fc<this->col-c)
+                for (unsigned j=fc; j<this->col-c; ++j)
+                    Pkey[j] = Pkey[j-1]+1;
+        } while(fc);
+        // change the role of A and B:
+        A = B;
+        B.clear();
     }
-    return matrix(v.row, v.col, sol);
+    
+    return det;
 }
 
-// protected
+/** Determinant built by application of the full permutation group.  This
+ *  routine is only called internally by matrix::determinant(). */
+ex matrix::determinant_symbolic_perm(void) const
+{
+    if (rows()==1)  // speed things up
+        return m[0];
+    
+    ex det;
+    ex term;
+    vector<unsigned> sigma(col);
+    for (unsigned i=0; i<col; ++i)
+        sigma[i]=i;
+    
+    do {
+        term = (*this)(sigma[0],0);
+        for (unsigned i=1; i<col; ++i)
+            term *= (*this)(sigma[i],i);
+        det += permutation_sign(sigma)*term;
+    } while (next_permutation(sigma.begin(), sigma.end()));
+    
+    return det;
+}
+
+/** Perform the steps of an ordinary Gaussian elimination to bring the matrix
+ *  into an upper echelon form.
+ *
+ *  @return sign is 1 if an even number of rows was swapped, -1 if an odd
+ *  number of rows was swapped and 0 if the matrix is singular. */
+int matrix::gauss_elimination(void)
+{
+    int sign = 1;
+    ensure_if_modifiable();
+    for (unsigned r1=0; r1<row-1; ++r1) {
+        int indx = pivot(r1);
+        if (indx == -1)
+            return 0;  // Note: leaves *this in a messy state.
+        if (indx > 0)
+            sign = -sign;
+        for (unsigned r2=r1+1; r2<row; ++r2) {
+            for (unsigned c=r1+1; c<col; ++c)
+                this->m[r2*col+c] -= this->m[r2*col+r1]*this->m[r1*col+c]/this->m[r1*col+r1];
+            for (unsigned c=0; c<=r1; ++c)
+                this->m[r2*col+c] = _ex0();
+        }
+    }
+    return sign;
+}
 
 /** Partial pivoting method.
- *  Usual pivoting returns the index to the element with the largest absolute
- *  value and swaps the current row with the one where the element was found.
- *  Here it does the same with the first non-zero element. (This works fine,
- *  but may be far from optimal for numerics.) */
-int matrix::pivot(unsigned ro)
+ *  Usual pivoting (symbolic==false) returns the index to the element with the
+ *  largest absolute value in column ro and swaps the current row with the one
+ *  where the element was found.  With (symbolic==true) it does the same thing
+ *  with the first non-zero element.
+ *
+ *  @param ro is the row to be inspected
+ *  @param symbolic signal if we want the first non-zero element to be pivoted
+ *  (true) or the one with the largest absolute value (false).
+ *  @return 0 if no interchange occured, -1 if all are zero (usually signaling
+ *  a degeneracy) and positive integer k means that rows ro and k were swapped.
+ */
+int matrix::pivot(unsigned ro, bool symbolic)
 {
-    unsigned k=ro;
+    unsigned k = ro;
     
-    for (unsigned r=ro; r<row; ++r) {
-        if (!m[r*col+ro].is_zero()) {
-            k = r;
-            break;
+    if (symbolic) {  // search first non-zero
+        for (unsigned r=ro; r<row; ++r) {
+            if (!m[r*col+ro].is_zero()) {
+                k = r;
+                break;
+            }
+        }
+    } else {  // search largest
+        numeric tmp(0);
+        numeric maxn(-1);
+        for (unsigned r=ro; r<row; ++r) {
+            GINAC_ASSERT(is_ex_of_type(m[r*col+ro],numeric));
+            if ((tmp = abs(ex_to_numeric(m[r*col+ro]))) > maxn &&
+                !tmp.is_zero()) {
+                maxn = tmp;
+                k = r;
+            }
         }
     }
-    if (m[k*col+ro].is_zero()) {
+    if (m[k*col+ro].is_zero())
         return -1;
-    }
     if (k!=ro) {  // swap rows
+        ensure_if_modifiable();
         for (unsigned c=0; c<col; ++c) {
             m[k*col+c].swap(m[ro*col+c]);
         }
index 38e9164..3a6f0e5 100644 (file)
@@ -37,10 +37,7 @@ class matrix : public basic
     GINAC_DECLARE_REGISTERED_CLASS(matrix, basic)
 
 // friends
-    friend ex determinant_numeric(const matrix & m);
-    friend ex determinant_symbolic_perm(const matrix & m);
-    friend ex determinant_symbolic_minor(const matrix & m);
-
+// (none)
 // member functions
 
     // default constructor, destructor, copy constructor, assignment operator
@@ -95,9 +92,16 @@ public:
     ex charpoly(const ex & lambda) const;
     matrix inverse(void) const;
     matrix fraction_free_elim(const matrix & vars, const matrix & v) const;
-    matrix solve(const matrix & v) const;
+    matrix solve(const matrix & vars, const matrix & rhs) const;
+    matrix old_solve(const matrix & v) const;  // FIXME: may be removed
 protected:
-    int pivot(unsigned ro);
+    ex determinant_numeric(void) const;
+    ex determinant_symbolic_minor(void) const;
+    ex determinant_symbolic_perm(void) const;
+    int gauss_elimination(void);
+    int fraction_free_elimination(void);
+    int pivot(unsigned ro, bool symbolic=true);
+private:  // FIXME: these should be obsoleted
     void ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2);
     void ffe_set(unsigned r, unsigned c, ex e);
     ex ffe_get(unsigned r, unsigned c) const;
index 07e8ec9..403e3ed 100644 (file)
@@ -21,6 +21,9 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
+#ifndef __GINAC_REMEMBER_H__
+#define __GINAC_REMEMBER_H__
+
 #include <vector>
 #include <list>
 
@@ -65,7 +68,7 @@ protected:
 };
 
 /** The remember table is organized like an n-fold associative cache
-    in a microprocessor. The table has a width of 's' (which is rounded
+    in a microprocessor.  The table has a width of 's' (which is rounded
     to table_size, some power of 2 near 's', internally) and a depth of 'as'
     (unless you choose that entries are never discarded). The place where
     an entry is stored depends on the hashvalue of the parameters of the
@@ -97,3 +100,5 @@ protected:
 #ifndef NO_NAMESPACE_GINAC
 } // namespace GiNaC
 #endif // ndef NO_NAMESPACE_GINAC
+
+#endif // ndef __GINAC_REMEMBER_H__
index fa6f54a..7018bc9 100644 (file)
@@ -138,6 +138,26 @@ OutputIterator mymerge3(InputIterator1 first1, InputIterator1 last1,
     }
 }
 
+// Compute the sign of a permutation of a vector of things.
+template <typename T>
+int permutation_sign(vector<T> s)
+{
+    if (s.size() < 2)
+        return 0;
+    int sigma = 1;
+    for (typename vector<T>::iterator i=s.begin(); i!=s.end()-1; ++i) {
+        for (typename vector<T>::iterator j=i+1; j!=s.end(); ++j) {
+            if (*i == *j)
+                return 0;
+            if (*i > *j) {
+                iter_swap(i,j);
+                sigma = -sigma;
+            }
+        }
+    }
+    return sigma;
+}
+
 // Collection of `construct on first use' wrappers for safely avoiding
 // internal object replication without running into the `static
 // initialization order fiasco'.  This chest of numbers helps speed up
index 52d9402..5bd33c8 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@
index 13f92f3..0f42a2c 100644 (file)
@@ -84,6 +84,7 @@ GINSH_LIBS = @GINSH_LIBS@
 LATEX = @LATEX@
 LEX = @LEX@
 LIBGINACCINT = @LIBGINACCINT@
+LIBTERMCAP = @LIBTERMCAP@
 LIBTOOL = @LIBTOOL@
 LN_S = @LN_S@
 LT_AGE = @LT_AGE@