From 9df145c8bfa8ce9f2cbe6c05673481b6ca4c4c22 Mon Sep 17 00:00:00 2001
From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Wed, 26 Jul 2000 00:49:22 +0000
Subject: [PATCH] - Made determinant_algo (in flags.h) really work. - Fixed
 minor bugs in matrix.cpp. - Small Documantation updates.

---
 INSTALL                  |   1 +
 check/check_matrices.cpp | 114 ++++++++++++------
 check/exam_lsolve.cpp    |  37 +++++-
 config.guess             |  13 +-
 doc/tutorial/ginac.texi  |   2 +-
 ginac/flags.h            |   6 +-
 ginac/inifcns.cpp        | 110 ++++++++++++-----
 ginac/inifcns.h          |   7 +-
 ginac/inifcns_gamma.cpp  |  30 ++---
 ginac/matrix.cpp         | 250 ++++++++++++++++++++-------------------
 ginac/matrix.h           |  16 +--
 ginac/numeric.h          |   1 -
 12 files changed, 358 insertions(+), 229 deletions(-)

diff --git a/INSTALL b/INSTALL
index 3730e91c..10e0b58f 100644
--- a/INSTALL
+++ b/INSTALL
@@ -113,3 +113,4 @@ to compile, install and work properly:
  < 5.14.39    | `VERBOTEN' by license (please bite your favorite lawyer)
  < 5.14.39,40 | compiles but does not feel happy at all (inconsistent!)
    5.14.41    | tested on egcs 1.1.1, gcc 2.95.2: only minor weirdnesses
+   5.14.44    | G__cpp_ginaccint.C needs manual fixes, doesn't work well
diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index c1fa6257..3e702b8e 100644
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -31,15 +31,13 @@ static unsigned integdom_matrix_determinants(void)
     
     for (int size=3; size<20; ++size) {
         matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
+        // populate one element in each row:
+        for (int r=0; r<size-1; ++r)
             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
+        // set the last row to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int c=0; c<size; ++c)
             A.set(size-1,c,A(0,c)-A(size-2,c));
-        }
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
@@ -51,29 +49,30 @@ static unsigned integdom_matrix_determinants(void)
     return result;
 }
 
-/* determinants of some sparse symbolic matrices with multivariate rational
- * function coefficients. */
+/* determinants of some symbolic matrices with multivariate rational function
+ * coefficients. */
 static unsigned rational_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
-
+    
     for (int size=3; size<8; ++size) {
         matrix A(size,size);
         for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 4, false, false, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, false, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
+            // populate one or two elements in each row:
+            for (int ec=0; ec<2; ++ec) {
+                ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
+                ex denom;
+                do {
+                    denom = sparse_tree(a, b, c, rand()%2, false, false, false);
+                } while (denom.is_zero());
+                A.set(r,unsigned(rand()%size),numer/denom);
+            }
         }
+        // set the last row to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int co=0; co<size; ++co)
+            A.set(size-1,co,A(0,co)-A(size-2,co));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
@@ -85,28 +84,29 @@ static unsigned rational_matrix_determinants(void)
     return result;
 }
 
-/* Some quite wild determinants with functions and stuff like that. */
-static unsigned wild_matrix_determinants(void)
+/* Some quite funny determinants with functions and stuff like that inside. */
+static unsigned funny_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
-    for (int size=3; size<6; ++size) {
+    for (int size=3; size<7; ++size) {
         matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 3, true, true, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, true, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
+        for (int co=0; co<size-1; ++co) {
+            // populate one or two elements in each row:
+            for (int ec=0; ec<2; ++ec) {
+                ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
+                ex denom;
+                do {
+                    denom = sparse_tree(a, b, c, rand()%2, false, true, false);
+                } while (denom.is_zero());
+                A.set(unsigned(rand()%size),co,numer/denom);
+            }
         }
+        // set the last column to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int ro=0; ro<size; ++ro)
+            A.set(ro,size-1,A(ro,0)-A(ro,size-2));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
@@ -118,6 +118,41 @@ static unsigned wild_matrix_determinants(void)
     return result;
 }
 
+/* compare results from different determinant algorithms.*/
+static unsigned compare_matrix_determinants(void)
+{
+    unsigned result = 0;
+    symbol a("a");
+    
+    for (int size=2; size<6; ++size) {
+        matrix A(size,size);
+        for (int co=0; co<size; ++co) {
+            for (int ro=0; ro<size; ++ro) {
+                // populate some elements
+                ex elem = 0;
+                if (rand()%(size-1) == 0)
+                    elem = sparse_tree(a, a, a, rand()%3, false, true, false);
+                A.set(ro,co,elem);
+            }
+        }
+        ex det_gauss = A.determinant(determinant_algo::gauss);
+        ex det_laplace = A.determinant(determinant_algo::laplace);
+        ex det_bareiss = A.determinant(determinant_algo::bareiss);
+        if ((det_gauss-det_laplace).normal() != 0 ||
+            (det_bareiss-det_laplace).normal() != 0) {
+            clog << "Determinant of " << size << "x" << size << " matrix "
+                 << endl << A << endl
+                 << "is inconsistent between different algorithms:" << endl
+                 << "Gauss elimination:   " << det_gauss << endl
+                 << "Minor elimination:   " << det_laplace << endl
+                 << "Fraction-free elim.: " << det_bareiss << endl;
+            ++result;
+        }
+    }
+    
+    return result;
+}
+
 unsigned check_matrices(void)
 {
     unsigned result = 0;
@@ -127,7 +162,8 @@ unsigned check_matrices(void)
     
     result += integdom_matrix_determinants();  cout << '.' << flush;
     result += rational_matrix_determinants();  cout << '.' << flush;
-    result += wild_matrix_determinants();  cout << '.' << flush;
+    result += funny_matrix_determinants();  cout << '.' << flush;
+    result += compare_matrix_determinants();  cout << '.' << flush;
     
     if (!result) {
         cout << " passed " << endl;
diff --git a/check/exam_lsolve.cpp b/check/exam_lsolve.cpp
index bcbc25fe..3deddbc2 100644
--- a/check/exam_lsolve.cpp
+++ b/check/exam_lsolve.cpp
@@ -32,7 +32,7 @@ static unsigned exam_lsolve1(void)
     eq = (3*x+5 == numeric(8));
     aux = lsolve(eq, x);
     if (aux != 1) {
-        result++;
+        ++result;
         clog << "solution of 3*x+5==8 erroneously returned "
              << aux << endl;
     }
@@ -60,7 +60,7 @@ static unsigned exam_lsolve2a(void)
     // It should have returned [x==(3+b^2)/(a+b),y==(3-a*b)/(a+b)]
     if (!normal(sol_x - (3+pow(b,2))/(a+b)).is_zero() ||
         !normal(sol_y - (3-a*b)/(a+b)).is_zero()) {
-        result++;
+        ++result;
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }
@@ -88,7 +88,7 @@ static unsigned exam_lsolve2b(void)
     // It should have returned [x==43/17,y==-10/17]
     if (!(sol_x - numeric(43,17)).is_zero() ||
         !(sol_y - numeric(-10,17)).is_zero()) {
-        result++;
+        ++result;
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }
@@ -116,7 +116,35 @@ static unsigned exam_lsolve2c(void)
     // It should have returned [x==-3/2*I,y==-1/2]
     if (!(sol_x - numeric(-3,2)*I).is_zero() ||
         !(sol_y - numeric(-1,2)).is_zero()) {
-        result++;
+        ++result;
+        clog << "solution of the system " << eqns << " for " << vars
+             << " erroneously returned " << sol << endl;
+    }
+    
+    return result;
+}
+
+static unsigned exam_lsolve2S(void)
+{
+    // A degenerate example that went wrong in GiNaC 0.6.2.
+    unsigned result = 0;
+    symbol x("x"), y("y"), t("t");
+    lst eqns, vars;
+    ex sol;
+    
+    // Create the linear system [0*x+0*y==0,0*x+1*y==t]...
+    eqns.append(0*x+0*y==0).append(0*x+1*y==t);
+    // ...to be solved for [x,y]...
+    vars.append(x).append(y);
+    // ...and solve it:
+    sol = lsolve(eqns, vars);
+    ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
+    ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
+    
+    // It should have returned [x==x,y==t]
+    if (!(sol_x - x).is_zero() ||
+        !(sol_y - t).is_zero()) {
+        ++result;
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }
@@ -135,6 +163,7 @@ unsigned exam_lsolve(void)
     result += exam_lsolve2a();  cout << '.' << flush;
     result += exam_lsolve2b();  cout << '.' << flush;
     result += exam_lsolve2c();  cout << '.' << flush;
+    result += exam_lsolve2S();  cout << '.' << flush;
     
     if (!result) {
         cout << " passed " << endl;
diff --git a/config.guess b/config.guess
index 49949643..323f8c9c 100755
--- a/config.guess
+++ b/config.guess
@@ -3,7 +3,7 @@
 #   Copyright (C) 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000
 #   Free Software Foundation, Inc.
 
-version='2000-06-13'
+version='2000-07-24'
 
 # This file is free software; you can redistribute it and/or modify it
 # under the terms of the GNU General Public License as published by
@@ -1088,6 +1088,17 @@ EOF
     DS/*:UNIX_System_V:*:*)
 	echo ${UNAME_MACHINE}-${UNAME_SYSTEM}-${UNAME_RELEASE}
 	exit 0 ;;
+    *:Plan9:*:*)
+	# "uname -m" is not consistent, so use $cputype instead. 386
+	# is converted to i386 for consistency with other x86
+	# operating systems.
+	if test "$cputype" = "386"; then
+	    UNAME_MACHINE=i386
+	else
+	    UNAME_MACHINE="$cputype"
+	fi
+	echo ${UNAME_MACHINE}-unknown-plan9
+	exit 0 ;;
 esac
 
 #echo '(No uname command or uname output not recognized.)' 1>&2
diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 536bf897..bc330c7d 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -2075,7 +2075,7 @@ negative real axis where the points on the axis itself belong to the
 upper part (i.e. continuous with quadrant II).  The inverse
 trigonometric and hyperbolic functions are not defined for complex
 arguments by the C++ standard, however.  Here, we follow the conventions
-used by CLN, which in turn follow the carefully structured definitions
+used by CLN, which in turn follow the carefully designed definitions
 in the Common Lisp standard.  Hopefully, future revisions of the C++
 standard incorporate these functions in the complex domain in a manner
 compatible with Common Lisp.
diff --git a/ginac/flags.h b/ginac/flags.h
index 77e3ec01..41a890b1 100644
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -39,9 +39,11 @@ public:
          };
 };
 
-class determinant_options {
+class determinant_algo {
 public:
-    enum { laplace,
+    enum { automatic,
+           gauss,
+           laplace,
            bareiss
          };
 };
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index e0d1e91e..b560100f 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -45,21 +45,21 @@ namespace GiNaC {
 // absolute value
 //////////
 
-static ex abs_evalf(const ex & x)
+static ex abs_evalf(const ex & arg)
 {
     BEGIN_TYPECHECK
-        TYPECHECK(x,numeric)
-    END_TYPECHECK(abs(x))
+        TYPECHECK(arg,numeric)
+    END_TYPECHECK(abs(arg))
     
-    return abs(ex_to_numeric(x));
+    return abs(ex_to_numeric(arg));
 }
 
-static ex abs_eval(const ex & x)
+static ex abs_eval(const ex & arg)
 {
-    if (is_ex_exactly_of_type(x, numeric))
-        return abs(ex_to_numeric(x));
+    if (is_ex_exactly_of_type(arg, numeric))
+        return abs(ex_to_numeric(arg));
     else
-        return abs(x).hold();
+        return abs(arg).hold();
 }
 
 REGISTER_FUNCTION(abs, eval_func(abs_eval).
@@ -70,41 +70,41 @@ REGISTER_FUNCTION(abs, eval_func(abs_eval).
 // Complex sign
 //////////
 
-static ex csgn_evalf(const ex & x)
+static ex csgn_evalf(const ex & arg)
 {
     BEGIN_TYPECHECK
-        TYPECHECK(x,numeric)
-    END_TYPECHECK(csgn(x))
+        TYPECHECK(arg,numeric)
+    END_TYPECHECK(csgn(arg))
     
-    return csgn(ex_to_numeric(x));
+    return csgn(ex_to_numeric(arg));
 }
 
-static ex csgn_eval(const ex & x)
+static ex csgn_eval(const ex & arg)
 {
-    if (is_ex_exactly_of_type(x, numeric))
-        return csgn(ex_to_numeric(x));
+    if (is_ex_exactly_of_type(arg, numeric))
+        return csgn(ex_to_numeric(arg));
     
-    else if (is_ex_exactly_of_type(x, mul)) {
-        numeric oc = ex_to_numeric(x.op(x.nops()-1));
+    else if (is_ex_exactly_of_type(arg, mul)) {
+        numeric oc = ex_to_numeric(arg.op(arg.nops()-1));
         if (oc.is_real()) {
             if (oc > 0)
                 // csgn(42*x) -> csgn(x)
-                return csgn(x/oc).hold();
+                return csgn(arg/oc).hold();
             else
                 // csgn(-42*x) -> -csgn(x)
-                return -csgn(x/oc).hold();
+                return -csgn(arg/oc).hold();
         }
         if (oc.real().is_zero()) {
             if (oc.imag() > 0)
                 // csgn(42*I*x) -> csgn(I*x)
-                return csgn(I*x/oc).hold();
+                return csgn(I*arg/oc).hold();
             else
                 // csgn(-42*I*x) -> -csgn(I*x)
-                return -csgn(I*x/oc).hold();
+                return -csgn(I*arg/oc).hold();
         }
 	}
    
-    return csgn(x).hold();
+    return csgn(arg).hold();
 }
 
 static ex csgn_series(const ex & arg,
@@ -113,13 +113,10 @@ static ex csgn_series(const ex & arg,
                       unsigned options)
 {
     const ex arg_pt = arg.subs(rel);
-    if (arg_pt.info(info_flags::numeric)) {
-        if (ex_to_numeric(arg_pt).real().is_zero())
-            throw (std::domain_error("csgn_series(): on imaginary axis"));
-        epvector seq;
-        seq.push_back(expair(csgn(arg_pt), _ex0()));
-        return pseries(rel,seq);
-    }
+    if (arg_pt.info(info_flags::numeric) &&
+        ex_to_numeric(arg_pt).real().is_zero())
+        throw (std::domain_error("csgn_series(): on imaginary axis"));
+    
     epvector seq;
     seq.push_back(expair(csgn(arg_pt), _ex0()));
     return pseries(rel,seq);
@@ -129,6 +126,61 @@ REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
                         evalf_func(csgn_evalf).
                         series_func(csgn_series));
 
+
+//////////
+// Eta function: log(x*y) == log(x) + log(y) + eta(x,y).
+//////////
+
+static ex eta_evalf(const ex & x, const ex & y)
+{
+    BEGIN_TYPECHECK
+        TYPECHECK(x,numeric)
+        TYPECHECK(y,numeric)
+    END_TYPECHECK(eta(x,y))
+        
+    numeric xim = imag(ex_to_numeric(x));
+    numeric yim = imag(ex_to_numeric(y));
+    numeric xyim = imag(ex_to_numeric(x*y));
+    return evalf(I/4*Pi)*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1));
+}
+
+static ex eta_eval(const ex & x, const ex & y)
+{
+    if (is_ex_exactly_of_type(x, numeric) &&
+        is_ex_exactly_of_type(y, numeric)) {
+        // don't call eta_evalf here because it would call Pi.evalf()!
+        numeric xim = imag(ex_to_numeric(x));
+        numeric yim = imag(ex_to_numeric(y));
+        numeric xyim = imag(ex_to_numeric(x*y));
+        return (I/4)*Pi*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1));
+    }
+    
+    return eta(x,y).hold();
+}
+
+static ex eta_series(const ex & arg1,
+                     const ex & arg2,
+                     const relational & rel,
+                     int order,
+                     unsigned options)
+{
+    const ex arg1_pt = arg1.subs(rel);
+    const ex arg2_pt = arg2.subs(rel);
+    if (ex_to_numeric(arg1_pt).imag().is_zero() ||
+        ex_to_numeric(arg2_pt).imag().is_zero() ||
+        ex_to_numeric(arg1_pt*arg2_pt).imag().is_zero()) {
+        throw (std::domain_error("eta_series(): on discontinuity"));
+    }
+    epvector seq;
+    seq.push_back(expair(eta(arg1_pt,arg2_pt), _ex0()));
+    return pseries(rel,seq);
+}
+
+REGISTER_FUNCTION(eta, eval_func(eta_eval).
+                       evalf_func(eta_evalf).
+                       series_func(eta_series));
+
+
 //////////
 // dilogarithm
 //////////
diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 9781d935..95ed15a1 100644
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -36,6 +36,9 @@ DECLARE_FUNCTION_1P(abs)
 /** Complex sign. */
 DECLARE_FUNCTION_1P(csgn)
 
+/** Eta function: log(a*b) == log(a) + log(b) + eta(a, b). */
+DECLARE_FUNCTION_2P(eta)
+
 /** Sine. */
 DECLARE_FUNCTION_1P(sin)
 
@@ -87,7 +90,7 @@ DECLARE_FUNCTION_1P(Li2)
 /** Trilogarithm. */
 DECLARE_FUNCTION_1P(Li3)
 
-// overloading at work: we cannot use the macros
+// overloading at work: we cannot use the macros here
 /** Riemann's Zeta-function. */
 extern const unsigned function_index_zeta1;
 inline function zeta(const ex & p1) {
@@ -106,7 +109,7 @@ DECLARE_FUNCTION_1P(tgamma)
 /** Beta-function. */
 DECLARE_FUNCTION_2P(beta)
 
-// overloading at work: we cannot use the macros
+// overloading at work: we cannot use the macros here
 /** Psi-function (aka digamma-function). */
 extern const unsigned function_index_psi1;
 inline function psi(const ex & p1) {
diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index bf60c1e5..aa1cab4c 100644
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -63,11 +63,10 @@ static ex lgamma_eval(const ex & x)
         // trap integer arguments:
         if (x.info(info_flags::integer)) {
             // lgamma(n) -> log((n-1)!) for postitive n
-            if (x.info(info_flags::posint)) {
+            if (x.info(info_flags::posint))
                 return log(factorial(x.exadd(_ex_1())));
-            } else {
+            else
                 throw (pole_error("lgamma_eval(): logarithmic pole",0));
-            }
         }
         //  lgamma_evalf should be called here once it becomes available
     }
@@ -98,15 +97,16 @@ static ex lgamma_series(const ex & arg,
     // from which follows
     //   series(lgamma(x),x==-m,order) ==
     //   series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order);
-    // This, however, seems to fail utterly because you run into branch-cut
-    // problems.  Somebody ought to implement it some day using an asymptotic
-    // series for tgamma:
     const ex arg_pt = arg.subs(rel);
     if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
     // if we got here we have to care for a simple pole of tgamma(-m):
-    throw (std::overflow_error("lgamma_series: please implement my at the poles"));
-    return _ex0();  // not reached
+    numeric m = -ex_to_numeric(arg_pt);
+    ex recur;
+    for (numeric p; p<=m; ++p)
+        recur += log(arg+p);
+    cout << recur << endl;
+    return (lgamma(arg+m+_ex1())-recur).series(rel, order, options);
 }
 
 
@@ -202,7 +202,7 @@ static ex tgamma_series(const ex & arg,
     ex ser_denom = _ex1();
     for (numeric p; p<=m; ++p)
         ser_denom *= arg+p;
-    return (tgamma(arg+m+_ex1())/ser_denom).series(rel, order+1);
+    return (tgamma(arg+m+_ex1())/ser_denom).series(rel, order+1, options);
 }
 
 
@@ -299,21 +299,21 @@ static ex beta_series(const ex & arg1,
         throw do_taylor();  // caught by function::series()
     // trap the case where arg1 is on a pole:
     if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive))
-        arg1_ser = tgamma(arg1+*s).series(rel,order);
+        arg1_ser = tgamma(arg1+*s).series(rel, order, options);
     else
         arg1_ser = tgamma(arg1).series(rel,order);
     // trap the case where arg2 is on a pole:
     if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive))
-        arg2_ser = tgamma(arg2+*s).series(rel,order);
+        arg2_ser = tgamma(arg2+*s).series(rel, order, options);
     else
         arg2_ser = tgamma(arg2).series(rel,order);
     // trap the case where arg1+arg2 is on a pole:
     if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive))
-        arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel,order);
+        arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel, order, options);
     else
         arg1arg2_ser = tgamma(arg2+arg1).series(rel,order);
     // compose the result (expanding all the terms):
-    return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel,order).expand();
+    return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel, order, options).expand();
 }
 
 
@@ -410,7 +410,7 @@ static ex psi1_series(const ex & arg,
     ex recur;
     for (numeric p; p<=m; ++p)
         recur += power(arg+p,_ex_1());
-    return (psi(arg+m+_ex1())-recur).series(rel, order);
+    return (psi(arg+m+_ex1())-recur).series(rel, order, options);
 }
 
 const unsigned function_index_psi1 =
@@ -536,7 +536,7 @@ static ex psi2_series(const ex & n,
     for (numeric p; p<=m; ++p)
         recur += power(arg+p,-n+_ex_1());
     recur *= factorial(n)*power(_ex_1(),n);
-    return (psi(n, arg+m+_ex1())-recur).series(rel, order);
+    return (psi(n, arg+m+_ex1())-recur).series(rel, order, options);
 }
 
 const unsigned function_index_psi2 =
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index feae32fe..52b4de36 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -152,7 +152,7 @@ void matrix::archive(archive_node &n) const
     exvector::const_iterator i = m.begin(), iend = m.end();
     while (i != iend) {
         n.add_ex("m", *i);
-        i++;
+        ++i;
     }
 }
 
@@ -174,15 +174,13 @@ void matrix::print(std::ostream & os, unsigned upper_precedence) const
     os << "[[ ";
     for (unsigned r=0; r<row-1; ++r) {
         os << "[[";
-        for (unsigned c=0; c<col-1; ++c) {
+        for (unsigned c=0; c<col-1; ++c)
             os << m[r*col+c] << ",";
-        }
         os << m[col*(r+1)-1] << "]], ";
     }
     os << "[[";
-    for (unsigned c=0; c<col-1; ++c) {
+    for (unsigned c=0; c<col-1; ++c)
         os << m[(row-1)*col+c] << ",";
-    }
     os << m[row*col-1] << "]] ]]";
 }
 
@@ -192,15 +190,13 @@ void matrix::printraw(std::ostream & os) const
     os << "matrix(" << row << "," << col <<",";
     for (unsigned r=0; r<row-1; ++r) {
         os << "(";
-        for (unsigned c=0; c<col-1; ++c) {
+        for (unsigned c=0; c<col-1; ++c)
             os << m[r*col+c] << ",";
-        }
         os << m[col*(r-1)-1] << "),";
     }
     os << "(";
-    for (unsigned c=0; c<col-1; ++c) {
+    for (unsigned c=0; c<col-1; ++c)
         os << m[(row-1)*col+c] << ",";
-    }
     os << m[row*col-1] << "))";
 }
 
@@ -219,6 +215,9 @@ ex matrix::op(int i) const
 /** returns matrix entry at position (i/col, i%col). */
 ex & matrix::let_op(int i)
 {
+    GINAC_ASSERT(i>=0);
+    GINAC_ASSERT(i<nops());
+    
     return m[i];
 }
 
@@ -226,9 +225,9 @@ ex & matrix::let_op(int i)
 ex matrix::expand(unsigned options) const
 {
     exvector tmp(row*col);
-    for (unsigned i=0; i<row*col; ++i) {
-        tmp[i]=m[i].expand(options);
-    }
+    for (unsigned i=0; i<row*col; ++i)
+        tmp[i] = m[i].expand(options);
+    
     return matrix(row, col, tmp);
 }
 
@@ -242,9 +241,9 @@ bool matrix::has(const ex & other) const
     if (is_equal(*other.bp)) return true;
     
     // search all the elements
-    for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
+    for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r)
         if ((*r).has(other)) return true;
-    }
+    
     return false;
 }
 
@@ -263,12 +262,10 @@ ex matrix::eval(int level) const
     
     // eval() entry by entry
     exvector m2(row*col);
-    --level;    
-    for (unsigned r=0; r<row; ++r) {
-        for (unsigned c=0; c<col; ++c) {
+    --level;
+    for (unsigned r=0; r<row; ++r)
+        for (unsigned c=0; c<col; ++c)
             m2[r*col+c] = m[r*col+c].eval(level);
-        }
-    }
     
     return (new matrix(row, col, m2))->setflag(status_flags::dynallocated |
                                                status_flags::evaluated );
@@ -291,11 +288,10 @@ ex matrix::evalf(int level) const
     // evalf() entry by entry
     exvector m2(row*col);
     --level;
-    for (unsigned r=0; r<row; ++r) {
-        for (unsigned c=0; c<col; ++c) {
+    for (unsigned r=0; r<row; ++r)
+        for (unsigned c=0; c<col; ++c)
             m2[r*col+c] = m[r*col+c].evalf(level);
-        }
-    }
+    
     return matrix(row, col, m2);
 }
 
@@ -343,11 +339,9 @@ matrix matrix::add(const matrix & other) const
     exvector sum(this->m);
     exvector::iterator i;
     exvector::const_iterator ci;
-    for (i=sum.begin(), ci=other.m.begin();
-         i!=sum.end();
-         ++i, ++ci) {
+    for (i=sum.begin(), ci=other.m.begin(); i!=sum.end(); ++i, ++ci)
         (*i) += (*ci);
-    }
+    
     return matrix(row,col,sum);
 }
 
@@ -363,11 +357,9 @@ matrix matrix::sub(const matrix & other) const
     exvector dif(this->m);
     exvector::iterator i;
     exvector::const_iterator ci;
-    for (i=dif.begin(), ci=other.m.begin();
-         i!=dif.end();
-         ++i, ++ci) {
+    for (i=dif.begin(), ci=other.m.begin(); i!=dif.end(); ++i, ++ci)
         (*i) -= (*ci);
-    }
+    
     return matrix(row,col,dif);
 }
 
@@ -438,7 +430,7 @@ matrix matrix::transpose(void) const
 
 /** Determinant of square matrix.  This routine doesn't actually calculate the
  *  determinant, it only implements some heuristics about which algorithm to
- *  call.  If all the elements of the matrix are elements of an integral domain
+ *  run.  If all the elements of the matrix are elements of an integral domain
  *  the determinant is also in that integral domain and the result is expanded
  *  only.  If one or more elements are from a quotient field the determinant is
  *  usually also in that quotient field and the result is normalized before it
@@ -446,85 +438,105 @@ matrix matrix::transpose(void) const
  *  [[a/(a-b),1],[b/(a-b),1]] is returned as unity.  (In this respect, it
  *  behaves like MapleV and unlike Mathematica.)
  *
+ *  @param     algo allows to chose an algorithm
  *  @return    the determinant as a new expression
- *  @exception logic_error (matrix not square) */
-ex matrix::determinant(void) const
+ *  @exception logic_error (matrix not square)
+ *  @see       determinant_algo */
+ex matrix::determinant(unsigned algo) const
 {
     if (row!=col)
         throw (std::logic_error("matrix::determinant(): matrix not square"));
     GINAC_ASSERT(row*col==m.capacity());
     if (this->row==1)  // continuation would be pointless
         return m[0];
-    
-    // Gather some information about the matrix:
+        
+    // Gather some statistical information about this matrix:
     bool numeric_flag = true;
     bool normal_flag = false;
-    unsigned sparse_count = 0;  // count non-zero elements
+    unsigned sparse_count = 0;  // counts non-zero elements
     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
-        if (!(*r).is_zero())
+        lst srl;  // symbol replacement list
+        ex rtest = (*r).to_rational(srl);
+        if (!rtest.is_zero())
             ++sparse_count;
-        if (!(*r).info(info_flags::numeric))
+        if (!rtest.info(info_flags::numeric))
             numeric_flag = false;
-        if ((*r).info(info_flags::rational_function) &&
-            !(*r).info(info_flags::crational_polynomial))
+        if (!rtest.info(info_flags::crational_polynomial) &&
+             rtest.info(info_flags::rational_function))
             normal_flag = true;
     }
     
-    // Purely numeric matrix handled by Gauss elimination
-    if (numeric_flag) {
-        ex det = 1;
-        matrix tmp(*this);
-        int sign = tmp.gauss_elimination();
-        for (int d=0; d<row; ++d)
-            det *= tmp.m[d*col+d];
-        return (sign*det);
-    }
-    
-    // Does anybody know when a matrix is really sparse?
-    // Maybe <~row/2.2 nonzero elements average in a row?
-    if (5*sparse_count<=row*col) {
-        // copy *this:
-        matrix tmp(*this);
-        int sign;
-        sign = tmp.fraction_free_elimination(true);
-        if (normal_flag)
-            return (sign*tmp.m[row*col-1]).normal();
-        else
-            return (sign*tmp.m[row*col-1]).expand();
+    // Here is the heuristics in case this routine has to decide:
+    if (algo == determinant_algo::automatic) {
+        // Minor expansion is generally a good starting point:
+        algo = determinant_algo::laplace;
+        // Does anybody know when a matrix is really sparse?
+        // Maybe <~row/2.236 nonzero elements average in a row?
+        if (5*sparse_count<=row*col)
+            algo = determinant_algo::bareiss;
+        // Purely numeric matrix can be handled by Gauss elimination.
+        // This overrides any prior decisions.
+        if (numeric_flag)
+            algo = determinant_algo::gauss;
     }
     
-    // Now come the minor expansion schemes.  We always develop such that the
-    // smallest minors (i.e, the trivial 1x1 ones) are on the rightmost column.
-    // For this to be efficient it turns out that the emptiest columns (i.e.
-    // the ones with most zeros) should be the ones on the right hand side.
-    // Therefore we presort the columns of the matrix:
-    typedef std::pair<unsigned,unsigned> uintpair;  // # of zeros, column
-    std::vector<uintpair> c_zeros;  // number of zeros in column
-    for (unsigned c=0; c<col; ++c) {
-        unsigned acc = 0;
-        for (unsigned r=0; r<row; ++r)
-            if (m[r*col+c].is_zero())
-                ++acc;
-        c_zeros.push_back(uintpair(acc,c));
-    }
-    sort(c_zeros.begin(),c_zeros.end());
-    // unfortunately std::vector<uintpair> can't be used for permutation_sign.
-    std::vector<unsigned> pre_sort;
-    for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
-        pre_sort.push_back(i->second);
-    int sign = permutation_sign(pre_sort);
-    exvector result(row*col);  // represents sorted matrix
-    unsigned c = 0;
-    for (std::vector<unsigned>::iterator i=pre_sort.begin();
-         i!=pre_sort.end();
-         ++i,++c) {
-        for (unsigned r=0; r<row; ++r)
-            result[r*col+c] = m[r*col+(*i)];
+    switch(algo) {
+        case determinant_algo::gauss: {
+            ex det = 1;
+            matrix tmp(*this);
+            int sign = tmp.gauss_elimination();
+            for (unsigned d=0; d<row; ++d)
+                det *= tmp.m[d*col+d];
+            if (normal_flag)
+                return (sign*det).normal();
+            else
+                return (sign*det).expand();
+        }
+        case determinant_algo::bareiss: {
+            matrix tmp(*this);
+            int sign;
+            sign = tmp.fraction_free_elimination(true);
+            if (normal_flag)
+                return (sign*tmp.m[row*col-1]).normal();
+            else
+                return (sign*tmp.m[row*col-1]).expand();
+        }
+        case determinant_algo::laplace:
+        default: {
+            // This is the minor expansion scheme.  We always develop such
+            // that the smallest minors (i.e, the trivial 1x1 ones) are on the
+            // rightmost column.  For this to be efficient it turns out that
+            // the emptiest columns (i.e. the ones with most zeros) should be
+            // the ones on the right hand side.  Therefore we presort the
+            // columns of the matrix:
+            typedef std::pair<unsigned,unsigned> uintpair;
+            std::vector<uintpair> c_zeros;  // number of zeros in column
+            for (unsigned c=0; c<col; ++c) {
+                unsigned acc = 0;
+                for (unsigned r=0; r<row; ++r)
+                    if (m[r*col+c].is_zero())
+                        ++acc;
+                c_zeros.push_back(uintpair(acc,c));
+            }
+            sort(c_zeros.begin(),c_zeros.end());
+            std::vector<unsigned> pre_sort;
+            for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
+                pre_sort.push_back(i->second);
+            int sign = permutation_sign(pre_sort);
+            exvector result(row*col);  // represents sorted matrix
+            unsigned c = 0;
+            for (std::vector<unsigned>::iterator i=pre_sort.begin();
+                 i!=pre_sort.end();
+                 ++i,++c) {
+                for (unsigned r=0; r<row; ++r)
+                    result[r*col+c] = m[r*col+(*i)];
+            }
+            
+            if (normal_flag)
+                return sign*matrix(row,col,result).determinant_minor().normal();
+            return sign*matrix(row,col,result).determinant_minor();
+        }
     }
-    
-    if (normal_flag)
-        return sign*matrix(row,col,result).determinant_minor().normal();
-    return sign*matrix(row,col,result).determinant_minor();
 }
 
 
@@ -647,16 +659,6 @@ matrix matrix::inverse(void) const
     return tmp;
 }
 
-// superfluous helper function, to be removed:
-void matrix::swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
-{
-    ensure_if_modifiable();
-    
-    ex tmp = (*this)(r1,c1);
-    set(r1,c1,(*this)(r2,c2));
-    set(r2,c2,tmp);
-}
-
 
 /** 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'
@@ -669,7 +671,7 @@ void matrix::swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
 matrix matrix::fraction_free_elim(const matrix & vars,
                                   const matrix & rhs) const
 {
-    // FIXME: use implementation of matrix::fraction_free_elimination
+    // FIXME: use implementation of matrix::fraction_free_elimination instead!
     if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
         throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
     
@@ -687,29 +689,29 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned k=0; (k<n)&&(r<m); ++k) {
         // find a nonzero pivot
         unsigned p;
-        for (p=r; (p<m)&&(a(p,k).is_zero()); ++p) {}
+        for (p=r; (p<m)&&(a.m[p*a.cols()+k].is_zero()); ++p) {}
         // pivot is in row p
         if (p<m) {
             if (p!=r) {
                 // swap rows p and r
                 for (unsigned j=k; j<n; ++j)
-                    a.swap(p,j,r,j);
-                b.swap(p,0,r,0);
+                    a.m[p*a.cols()+j].swap(a.m[r*a.cols()+j]);
+                b.m[p*b.cols()].swap(b.m[r*b.cols()]);
                 // keep track of sign changes due to row exchange
                 sign *= -1;
             }
             for (unsigned i=r+1; i<m; ++i) {
                 for (unsigned j=k+1; j<n; ++j) {
-                    a.set(i,j,(a(r,k)*a(i,j)
-                              -a(r,j)*a(i,k))/divisor);
-                    a.set(i,j,a(i,j).normal());
+                    a.set(i,j,(a.m[r*a.cols()+k]*a.m[i*a.cols()+j]
+                              -a.m[r*a.cols()+j]*a.m[i*a.cols()+k])/divisor);
+                    a.set(i,j,a.m[i*a.cols()+j].normal());
                 }
-                b.set(i,0,(a(r,k)*b(i,0)
-                          -b(r,0)*a(i,k))/divisor);
-                b.set(i,0,b(i,0).normal());
-                a.set(i,k,0);
+                b.set(i,0,(a.m[r*a.cols()+k]*b.m[i*b.cols()]
+                          -b.m[r*b.cols()]*a.m[i*a.cols()+k])/divisor);
+                b.set(i,0,b.m[i*b.cols()].normal());
+                a.set(i,k,_ex0());
             }
-            divisor = a(r,k);
+            divisor = a.m[r*a.cols()+k];
             ++r;
         }
     }
@@ -720,8 +722,8 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned r=0; r<m; ++r) {
         int zero_in_this_row=0;
         for (unsigned c=0; c<n; ++c) {
-            if (a(r,c).is_zero())
-               zero_in_this_row++;
+            if (a.m[r*a.cols()+c].is_zero())
+               ++zero_in_this_row;
             else
                 break;
         }
@@ -739,26 +741,26 @@ matrix matrix::fraction_free_elim(const matrix & vars,
             first_non_zero++;
         if (first_non_zero>n) {
             // row consists only of zeroes, corresponding rhs must be 0 as well
-            if (!b(r,0).is_zero()) {
+            if (!b.m[r*b.cols()].is_zero()) {
                 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
             }
         } else {
             // assign solutions for vars between first_non_zero+1 and
             // last_assigned_sol-1: free parameters
             for (unsigned c=first_non_zero; c<last_assigned_sol-1; ++c)
-                sol.set(c,0,vars(c,0));
-            ex e = b(r,0);
+                sol.set(c,0,vars.m[c*vars.cols()]);
+            ex e = b.m[r*b.cols()];
             for (unsigned c=first_non_zero; c<n; ++c)
-                e -= a(r,c)*sol(c,0);
+                e -= a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
             sol.set(first_non_zero-1,0,
-                    (e/a(r,first_non_zero-1)).normal());
+                    (e/(a.m[r*a.cols()+(first_non_zero-1)])).normal());
             last_assigned_sol = first_non_zero;
         }
     }
     // assign solutions for vars between 1 and
     // last_assigned_sol-1: free parameters
     for (unsigned c=0; c<last_assigned_sol-1; ++c)
-        sol.set(c,0,vars(c,0));
+        sol.set(c,0,vars.m[c*vars.cols()]);
     
 #ifdef DO_GINAC_ASSERT
     // test solution with echelon matrix
@@ -1159,7 +1161,7 @@ int matrix::pivot(unsigned ro, bool symbolic)
     
     if (symbolic) {  // search first non-zero
         for (unsigned r=ro; r<row; ++r) {
-            if (!m[r*col+ro].is_zero()) {
+            if (!m[r*col+ro].expand().is_zero()) {
                 k = r;
                 break;
             }
diff --git a/ginac/matrix.h b/ginac/matrix.h
index cef3dc36..b737a38c 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -27,11 +27,6 @@
 #include "basic.h"
 #include "ex.h"
 
-namespace std {
-    // forward declaration, so <stdexcept> need not be included:
-    class range_error;
-}
-
 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
 #endif // ndef NO_NAMESPACE_GINAC
@@ -82,9 +77,9 @@ protected:
     
     // non-virtual functions in this class
 public:
-    unsigned rows(void) const        //! get number of rows.
+    unsigned rows(void) const        //! Get number of rows.
         { return row; }
-    unsigned cols(void) const        //! get number of columns.
+    unsigned cols(void) const        //! Get number of columns.
         { return col; }
     matrix add(const matrix & other) const;
     matrix sub(const matrix & other) const;
@@ -92,7 +87,7 @@ public:
     const ex & operator() (unsigned ro, unsigned co) const;
     matrix & set(unsigned ro, unsigned co, ex value);
     matrix transpose(void) const;
-    ex determinant(void) const;
+    ex determinant(unsigned options = determinant_algo::automatic) const;
     ex trace(void) const;
     ex charpoly(const symbol & lambda) const;
     matrix inverse(void) const;
@@ -105,7 +100,6 @@ protected:
     int division_free_elimination(void);
     int fraction_free_elimination(bool det = false);
     int pivot(unsigned ro, bool symbolic=true);
-    void swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2);
     
 // member variables
 protected:
@@ -147,8 +141,8 @@ inline unsigned cols(const matrix & m)
 inline matrix transpose(const matrix & m)
 { return m.transpose(); }
 
-inline ex determinant(const matrix & m)
-{ return m.determinant(); }
+inline ex determinant(const matrix & m, unsigned options = determinant_algo::automatic)
+{ return m.determinant(options); }
 
 inline ex trace(const matrix & m)
 { return m.trace(); }
diff --git a/ginac/numeric.h b/ginac/numeric.h
index d55e1d13..c5086f43 100644
--- a/ginac/numeric.h
+++ b/ginac/numeric.h
@@ -83,7 +83,6 @@ class numeric : public basic
     friend const numeric atanh(const numeric & x);
     friend const numeric Li2(const numeric & x);
     friend const numeric zeta(const numeric & x);
-    // friend const numeric bernoulli(const numeric & n);
     friend const numeric fibonacci(const numeric & n);
     friend numeric abs(const numeric & x);
     friend numeric mod(const numeric & a, const numeric & b);
-- 
2.45.0