From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Thu, 10 Aug 2000 17:14:39 +0000 (+0000)
Subject: - Complete revamp of methods in class matrix.  Some redundant (and poor)
X-Git-Tag: release_0-6-4~1
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=db5765dc91202851739e196ba11bfccb0b3fe7bc;hp=55d35dcf72dc411c8265628fcad2bd67d320a8c9

- Complete revamp of methods in class matrix.  Some redundant (and poor)
  implementations of elimination schemes were thrown out.  The code is now
  highly orthogonal, more flexible and much more efficient.  The function
  lsolve() now falls back to matrix::solve(), which in turn calls different
  elimination strategies.  matrix::fraction_free_elim() was thrown out for
  good.  matrix::pivot() has now two integer arguments (ro and co) to allow
  pivoting when eliminating to echelon form.
- Some new checks to guard against regressions in the new matrix code.
- Some minor bug-fixes in documentation.
- Bumped up version from 0.6.3 to 0.6.4.
- ...and the usual changes triggered by other changes.
---

diff --git a/INSTALL b/INSTALL
index 10e0b58f..584801e3 100644
--- a/INSTALL
+++ b/INSTALL
@@ -56,13 +56,13 @@ The time the "make" step takes depends heavily on optimization levels.
 Large amounts of memory (>128MB) will be required by the compiler,
 also depending on optimization.  To give you a rough idea of what you
 have to expect the following table may be helpful.  It was measured on
-a P-III/500MHz with "enough" memory:
+an Athlon/800MHz with "enough" memory:
 
 step:         | GCC optimization level:  | comment:
               |  -O0   |  -O1   |  -O2   |
 --------------+--------+--------+--------+-------------------------------
-make          |  ~4m   |  ~6m   |  ~10m  | building shared and static lib
-make check    |  ~3m   |  ~3m   |  ~7m   | mostly due to compilation
+make          |   ~4m  |   ~5m  |   ~6m  | building shared and static lib
+make check    |  ~20m  |  ~11m  |  ~12m  | largely due to compilation
 
 
 COMMON PROBLEMS
@@ -112,5 +112,6 @@ 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
+   5.14.41    | tested on egcs 1.1.1, gcc 2.95.2: some weirdnesses
+   5.14.44    | does not compile: G__cpp_ginaccint.C needs manual fixes
+   5.14.45    | compiles fine, still some weirdnesses during execution
diff --git a/NEWS b/NEWS
index 7ee51a54..4a119ee6 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,12 @@
 This file records noteworthy changes.
 
+0.6.4 (10 August 2000)
+* Complete revamp of methods in class matrix.  Some redundant (and poor)
+  implementations of elimination schemes were thrown out.  The code is now
+  highly orthogonal, more flexible and much more efficient.
+* Some long standing and quite nasty bugs were discovered and fixed in the
+  following functions: add::normal(), heur_gcd(), sr_gcd() and Order_eval().
+
 0.6.3 (25 July 2000)
 * Derivatives are now assembled in a slightly different manner (i.e. they
   might 'look' different on first sight).  Under certain circumstances this
diff --git a/check/check_lsolve.cpp b/check/check_lsolve.cpp
index 4d1f2373..1873bb2d 100644
--- a/check/check_lsolve.cpp
+++ b/check/check_lsolve.cpp
@@ -23,85 +23,134 @@
 
 #include "checks.h"
 
-static unsigned lsolve1(int size)
+static unsigned check_matrix_solve(unsigned m, unsigned n, unsigned p,
+                                   unsigned degree)
 {
-    // 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));
+    const symbol a("a");
+    matrix A(m,n);
+    matrix B(m,p);
+    // set the first min(m,n) rows of A and B
+    for (unsigned ro=0; (ro<m)&&(ro<n); ++ro) {
+        for (unsigned co=0; co<n; ++co)
+            A.set(ro,co,dense_univariate_poly(a,degree));
+        for (unsigned co=0; co<p; ++co)
+            B.set(ro,co,dense_univariate_poly(a,degree));
     }
-    if (A.determinant().is_zero())
-        clog << "lsolve1: singular system!" << endl;
-    
+    // repeat excessive rows of A and B to avoid excessive construction of
+    // overdetermined linear systems
+    for (unsigned ro=n; ro<m; ++ro) {
+        for (unsigned co=0; co<n; ++co)
+            A.set(ro,co,A(ro-1,co));
+        for (unsigned co=0; co<p; ++co)
+            B.set(ro,co,B(ro-1,co));
+    }
+    // create a vector of n*p symbols all named "xrc" where r and c are ints
+    vector<symbol> x;
+    matrix X(n,p);
+    for (unsigned i=0; i<n; ++i) {
+        for (unsigned j=0; j<p; ++j) {
+            char buf[4];
+            ostrstream(buf,sizeof(buf)) << i << j << ends;
+            x.push_back(symbol(string("x")+buf));
+            X.set(i,j,x[p*i+j]);
+        }
+    }
+    matrix sol(n,p);
     // Solve the system A*X==B:
-    X = A.old_solve(B);
+    try {
+        sol = A.solve(X, B);
+    } catch (const exception & err) {  // catch runtime_error
+        // Presumably, the coefficient matrix A was degenerate
+        string errwhat = err.what();
+        if (errwhat == "matrix::solve(): inconsistent linear system")
+            return 0;
+        else
+            clog << "caught exception: " << errwhat << endl;
+        throw;
+    }
     
-    // check the result:
+    // check the result with our original matrix:
     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())
+    for (unsigned ro=0; ro<m; ++ro) {
+        for (unsigned pco=0; pco<p; ++pco) {
+            ex e = 0;
+            for (unsigned co=0; co<n; ++co)
+            e += A(ro,co)*sol(co,pco);
+            if (!(e-B(ro,pco)).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
+             << "X == " << sol << endl
              << "B == " << B << endl;
-        ++result;
+        return 1;
     }
-    return result;
+    
+    return 0;
 }
 
-static unsigned lsolve2(int size)
+static unsigned check_inifcns_lsolve(unsigned n)
 {
-    // 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));
+    for (int repetition=0; repetition<100; ++repetition) {
+        // create two size n vectors of symbols, one for the coefficients
+        // a[0],..,a[n], one for indeterminates x[0]..x[n]:
+        vector<symbol> a;
+        vector<symbol> x;
+        for (unsigned i=0; i<n; ++i) {
+            char buf[3];
+            ostrstream(buf,sizeof(buf)) << i << ends;
+            a.push_back(symbol(string("a")+buf));
+            x.push_back(symbol(string("x")+buf));
+        }
+        lst eqns;  // equation list
+        lst vars;  // variable list
+        ex sol;    // solution
+        // Create a random linear system...
+        for (unsigned i=0; i<n; ++i) {
+            ex lhs = rand()%201-100;
+            ex rhs = rand()%201-100;
+            for (unsigned j=0; j<n; ++j) {
+                // ...with small coefficients to give degeneracy a chance...
+                lhs += a[j]*(rand()%21-10);
+                rhs += x[j]*(rand()%21-10);
+            }
+            eqns.append(lhs==rhs);
+            vars.append(x[i]);
+        }
+        // ...solve it...
+        sol = lsolve(eqns, vars);
+        
+        // ...and check the solution:
+        if (sol.nops() == 0) {
+            // no solution was found
+            // is the coefficient matrix really, really, really degenerate?
+            matrix coeffmat(n,n);
+            for (unsigned ro=0; ro<n; ++ro)
+                for (unsigned co=0; co<n; ++co)
+                    coeffmat.set(ro,co,eqns.op(co).rhs().coeff(a[co],1));
+            if (!coeffmat.determinant().is_zero()) {
+                ++result;
+                clog << "solution of the system " << eqns << " for " << vars
+                     << " was not found" << endl;
+            }
+        } else {
+            // insert the solution into rhs of out equations
+            bool errorflag = false;
+            for (unsigned i=0; i<n; ++i)
+                if (eqns.op(i).rhs().subs(sol) != eqns.op(i).lhs())
+                    errorflag = true;
+            if (errorflag) {
+                ++result;
+                clog << "solution of the system " << eqns << " for " << vars
+                     << " erroneously returned " << sol << endl;
+            }
+        }
     }
-    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;
 }
 
@@ -112,11 +161,34 @@ unsigned check_lsolve(void)
     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;
+    // solve some numeric linear systems
+    for (unsigned n=1; n<12; ++n)
+        result += check_matrix_solve(n, n, 1, 0);
+    cout << '.' << flush;
+    // solve some underdetermined numeric systems
+    for (unsigned n=1; n<12; ++n)
+        result += check_matrix_solve(n+1, n, 1, 0);
+    cout << '.' << flush;
+    // solve some overdetermined numeric systems
+    for (unsigned n=1; n<12; ++n)
+        result += check_matrix_solve(n, n+1, 1, 0);
+    cout << '.' << flush;
+    // solve some multiple numeric systems
+    for (unsigned n=1; n<12; ++n)
+        result += check_matrix_solve(n, n, n/3+1, 0);
+    cout << '.' << flush;
+    // solve some symbolic linear systems
+    for (unsigned n=1; n<7; ++n)
+        result += check_matrix_solve(n, n, 1, 2);
+    cout << '.' << flush;
     
+    // check lsolve, the wrapper function around matrix::solve()
+    result += check_inifcns_lsolve(2);  cout << '.' << flush;
+    result += check_inifcns_lsolve(3);  cout << '.' << flush;
+    result += check_inifcns_lsolve(4);  cout << '.' << flush;
+    result += check_inifcns_lsolve(5);  cout << '.' << flush;
+    result += check_inifcns_lsolve(6);  cout << '.' << flush;
+        
     if (!result) {
         cout << " passed " << endl;
         clog << "(no output)" << endl;
diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index 3e702b8e..62b74ba2 100644
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -29,14 +29,14 @@ static unsigned integdom_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a");
     
-    for (int size=3; size<20; ++size) {
+    for (unsigned size=3; size<20; ++size) {
         matrix A(size,size);
         // populate one element in each row:
-        for (int r=0; r<size-1; ++r)
+        for (unsigned r=0; r<size-1; ++r)
             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
         // 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)
+        for (unsigned 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 "
@@ -56,11 +56,11 @@ static unsigned rational_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
-    for (int size=3; size<8; ++size) {
+    for (unsigned size=3; size<8; ++size) {
         matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
+        for (unsigned r=0; r<size-1; ++r) {
             // populate one or two elements in each row:
-            for (int ec=0; ec<2; ++ec) {
+            for (unsigned ec=0; ec<2; ++ec) {
                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
                 ex denom;
                 do {
@@ -71,7 +71,7 @@ static unsigned rational_matrix_determinants(void)
         }
         // 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)
+        for (unsigned 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 "
@@ -90,11 +90,11 @@ static unsigned funny_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
-    for (int size=3; size<7; ++size) {
+    for (unsigned size=3; size<7; ++size) {
         matrix A(size,size);
-        for (int co=0; co<size-1; ++co) {
+        for (unsigned co=0; co<size-1; ++co) {
             // populate one or two elements in each row:
-            for (int ec=0; ec<2; ++ec) {
+            for (unsigned ec=0; ec<2; ++ec) {
                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
                 ex denom;
                 do {
@@ -103,9 +103,9 @@ static unsigned funny_matrix_determinants(void)
                 A.set(unsigned(rand()%size),co,numer/denom);
             }
         }
-        // set the last column to a linear combination of two other lines
+        // set the last column to a linear combination of two other columns
         // to guarantee that the determinant is zero:
-        for (int ro=0; ro<size; ++ro)
+        for (unsigned 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 "
@@ -124,27 +124,30 @@ static unsigned compare_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a");
     
-    for (int size=2; size<6; ++size) {
+    for (unsigned size=2; size<7; ++size) {
         matrix A(size,size);
-        for (int co=0; co<size; ++co) {
-            for (int ro=0; ro<size; ++ro) {
+        for (unsigned co=0; co<size; ++co) {
+            for (unsigned ro=0; ro<size; ++ro) {
                 // populate some elements
                 ex elem = 0;
-                if (rand()%(size-1) == 0)
+                if (rand()%(size/2) == 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_divfree = A.determinant(determinant_algo::divfree);
         ex det_bareiss = A.determinant(determinant_algo::bareiss);
         if ((det_gauss-det_laplace).normal() != 0 ||
-            (det_bareiss-det_laplace).normal() != 0) {
+            (det_bareiss-det_laplace).normal() != 0 ||
+            (det_divfree-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
+                 << "Division-free elim.: " << det_divfree << endl
                  << "Fraction-free elim.: " << det_bareiss << endl;
             ++result;
         }
@@ -153,6 +156,43 @@ static unsigned compare_matrix_determinants(void)
     return result;
 }
 
+static unsigned symbolic_matrix_inverse(void)
+{
+    unsigned result = 0;
+    symbol a("a"), b("b"), c("c");
+    
+    for (unsigned size=2; size<5; ++size) {
+        matrix A(size,size);
+        do {
+            for (unsigned co=0; co<size; ++co) {
+                for (unsigned ro=0; ro<size; ++ro) {
+                    // populate some elements
+                    ex elem = 0;
+                    if (rand()%(size/2) == 0)
+                        elem = sparse_tree(a, b, c, rand()%2, false, true, false);
+                    A.set(ro,co,elem);
+                }
+            }
+        } while (A.determinant() == 0);
+        matrix B = A.inverse();
+        matrix C = A.mul(B);
+        bool ok = true;
+        for (unsigned ro=0; ro<size; ++ro)
+            for (unsigned co=0; co<size; ++co)
+                if (C(ro,co).normal() != (ro==co?1:0))
+                    ok = false;
+        if (!ok) {
+            clog << "Inverse of " << size << "x" << size << " matrix "
+                 << endl << A << endl
+                 << "erroneously returned: "
+                 << endl << B << endl;
+            ++result;
+        }
+    }
+    
+    return result;
+}
+
 unsigned check_matrices(void)
 {
     unsigned result = 0;
@@ -164,6 +204,7 @@ unsigned check_matrices(void)
     result += rational_matrix_determinants();  cout << '.' << flush;
     result += funny_matrix_determinants();  cout << '.' << flush;
     result += compare_matrix_determinants();  cout << '.' << flush;
+    result += symbolic_matrix_inverse();  cout << '.' << flush;
     
     if (!result) {
         cout << " passed " << endl;
diff --git a/check/checks.cpp b/check/checks.cpp
index 7a979f08..b7827381 100644
--- a/check/checks.cpp
+++ b/check/checks.cpp
@@ -71,7 +71,7 @@ int main()
         } else {
             cout << "(" << result << " individual failures)" << endl;
         }
-        cout << "please check check.out against check.ref for more details."
+        cout << "please check checks.out against check.ref for more details."
              << endl << "happy debugging!" << endl;
     }
     
diff --git a/check/checks.h b/check/checks.h
index d95ad1a2..5827f223 100644
--- a/check/checks.h
+++ b/check/checks.h
@@ -27,6 +27,7 @@
 #include <stdlib.h>
 
 #include "ginac.h"
+using namespace std;
 
 #ifndef NO_NAMESPACE_GINAC
 using namespace GiNaC;
diff --git a/check/exam_differentiation.cpp b/check/exam_differentiation.cpp
index 085ff451..1a5aabb1 100644
--- a/check/exam_differentiation.cpp
+++ b/check/exam_differentiation.cpp
@@ -67,40 +67,19 @@ static unsigned exam_differentiation1(void)
     e = (e1 * e2).expand();
     
     // d e / dx:
-    d = 121 - 55*pow(x,-2) - 66*pow(x,-3) - 30*pow(x,-3)*pow(y,-2)
-        - 42*pow(x,-3)*pow(y,-1) - 78*pow(x,-3)*y
-        - 102*pow(x,-3)*pow(y,2) - 25*pow(x,-2) * pow(y,-2)
-        - 35*pow(x,-2)*pow(y,-1) - 65*pow(x,-2)*y
-        - 85*pow(x,-2)*pow(y,2) + 77*pow(y,-1) + 143*y + 187*pow(y,2)
-        + 130*x*pow(y,-2) + 182*pow(y,-1)*x + 338*x*y + 442*x*pow(y,2)
-        + 55*pow(y,-2) + 286*x;
+    d = ex("121-55/x^2-66/x^3-30/x^3/y^2-42/x^3/y-78/x^3*y-102/x^3*y^2-25/x^2/y^2-35/x^2/y-65/x^2*y-85/x^2*y^2+77/y+143*y+187*y^2+130*x/y^2+182/y*x+338*x*y+442*x*y^2+55/y^2+286*x",lst(x,y));
     result += check_diff(e, x, d);
     
     // d e / dy:
-    d = 91 - 30*pow(x,-2)*pow(y,-3) - 21*pow(x,-2)*pow(y,-2)
-        + 39*pow(x,-2) + 102*pow(x,-2)*y - 50*pow(x,-1)*pow(y,-3)
-        - 35*pow(x,-1)*pow(y,-2) + 65*pow(x,-1) + 170*pow(x,-1)*y
-        - 77*pow(y,-2)*x + 143*x + 374*x*y - 130*pow(y,-3)*pow(x,2)
-        - 91*pow(y,-2)*pow(x,2) + 169*pow(x,2) + 442*pow(x,2)*y
-        - 110*pow(y,-3)*x - 70*pow(y,-3) + 238*y - 49*pow(y,-2);
+    d = ex("91-30/x^2/y^3-21/x^2/y^2+39/x^2+102/x^2*y-50/x/y^3-35/x/y^2+65/x+170/x*y-77*x/y^2+143*x+374*x*y-130/y^3*x^2-91/y^2*x^2+169*x^2+442*x^2*y-110/y^3*x-70/y^3+238*y-49/y^2",lst(x,y));
     result += check_diff(e, y, d);
     
     // d^2 e / dx^2:
-    d = 286 + 90*pow(x,-4)*pow(y,-2) + 126*pow(x,-4)*pow(y,-1)
-        + 234*pow(x,-4)*y + 306*pow(x,-4)*pow(y,2)
-        + 50*pow(x,-3)*pow(y,-2) + 70*pow(x,-3)*pow(y,-1)
-        + 130*pow(x,-3)*y + 170*pow(x,-3)*pow(y,2)
-        + 130*pow(y,-2) + 182*pow(y,-1) + 338*y + 442*pow(y,2)
-        + 198*pow(x,-4) + 110*pow(x,-3);
+    d = ex("286+90/x^4/y^2+126/x^4/y+234/x^4*y+306/x^4*y^2+50/x^3/y^2+70/x^3/y+130/x^3*y+170/x^3*y^2+130/y^2+182/y+338*y+442*y^2+198/x^4+110/x^3",lst(x,y));
     result += check_diff(e, x, d, 2);
     
     // d^2 e / dy^2:
-    d = 238 + 90*pow(x,-2)*pow(y,-4) + 42*pow(x,-2)*pow(y,-3)
-        + 102*pow(x,-2) + 150*pow(x,-1)*pow(y,-4)
-        + 70*pow(x,-1)*pow(y,-3) + 170*pow(x,-1) + 330*x*pow(y,-4)
-        + 154*x*pow(y,-3) + 374*x + 390*pow(x,2)*pow(y,-4)
-        + 182*pow(x,2)*pow(y,-3) + 442*pow(x,2) + 210*pow(y,-4)
-        + 98*pow(y,-3);
+    d = ex("238+90/x^2/y^4+42/x^2/y^3+102/x^2+150/x/y^4+70/x/y^3+170/x+330*x/y^4+154*x/y^3+374*x+390*x^2/y^4+182*x^2/y^3+442*x^2+210/y^4+98/y^3",lst(x,y));
     result += check_diff(e, y, d, 2);
     
     return result;
diff --git a/check/exam_lsolve.cpp b/check/exam_lsolve.cpp
index 3deddbc2..04eb438b 100644
--- a/check/exam_lsolve.cpp
+++ b/check/exam_lsolve.cpp
@@ -86,8 +86,8 @@ static unsigned exam_lsolve2b(void)
     ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
     
     // 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()) {
+    if ((sol_x != numeric(43,17)) ||
+        (sol_y != numeric(-10,17))) {
         ++result;
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
@@ -114,8 +114,8 @@ static unsigned exam_lsolve2c(void)
     ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
     
     // 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()) {
+    if ((sol_x != numeric(-3,2)*I) ||
+        (sol_y != numeric(-1,2))) {
         ++result;
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
@@ -142,8 +142,39 @@ static unsigned exam_lsolve2S(void)
     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()) {
+    if ((sol_x != x) ||
+        (sol_y != t)) {
+        ++result;
+        clog << "solution of the system " << eqns << " for " << vars
+             << " erroneously returned " << sol << endl;
+    }
+    
+    return result;
+}
+
+static unsigned exam_lsolve3S(void)
+{
+    // A degenerate example that went wrong while trying to improve elimination
+    unsigned result = 0;
+    symbol b("b"), c("c");
+    symbol x("x"), y("y"), z("z");
+    lst eqns, vars;
+    ex sol;
+    
+    // Create the linear system [y+z==b,-y+z==c] with one additional row...
+    eqns.append(ex(0)==ex(0)).append(b==z+y).append(c==z-y);
+    // ...to be solved for [x,y,z]...
+    vars.append(x).append(y).append(z);
+    // ...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)
+    ex sol_z = sol.op(2).rhs();  // rhs of solution for third variable (z)
+    
+    // It should have returned [x==x,y==t,]
+    if ((sol_x != x) ||
+        (sol_y != (b-c)/2) ||
+        (sol_z != (b+c)/2)) {
         ++result;
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
@@ -164,6 +195,7 @@ unsigned exam_lsolve(void)
     result += exam_lsolve2b();  cout << '.' << flush;
     result += exam_lsolve2c();  cout << '.' << flush;
     result += exam_lsolve2S();  cout << '.' << flush;
+    result += exam_lsolve3S();  cout << '.' << flush;
     
     if (!result) {
         cout << " passed " << endl;
diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp
index db61e3e1..9f17bb3d 100644
--- a/check/exam_matrices.cpp
+++ b/check/exam_matrices.cpp
@@ -180,6 +180,42 @@ static unsigned matrix_invert3(void)
     return result;
 }
 
+static unsigned matrix_solve2(void)
+{
+    // check the solution of the multiple system A*X = B:
+    //     [ 1  2 -1 ] [ x0 y0 ]   [ 4 0 ]
+    //     [ 1  4 -2 ]*[ x1 y1 ] = [ 7 0 ]
+    //     [ a -2  2 ] [ x2 y2 ]   [ a 4 ]
+    unsigned result = 0;
+    symbol a("a");
+    symbol x0("x0"), x1("x1"), x2("x2");
+    symbol y0("y0"), y1("y1"), y2("y2");
+    matrix A(3,3);
+    A.set(0,0,1).set(0,1,2).set(0,2,-1);
+    A.set(1,0,1).set(1,1,4).set(1,2,-2);
+    A.set(2,0,a).set(2,1,-2).set(2,2,2);
+    matrix B(3,2);
+    B.set(0,0,4).set(1,0,7).set(2,0,a);
+    B.set(0,1,0).set(1,1,0).set(2,1,4);
+    matrix X(3,2);
+    X.set(0,0,x0).set(1,0,x1).set(2,0,x2);
+    X.set(0,1,y0).set(1,1,y1).set(2,1,y2);
+    matrix cmp(3,2);
+    cmp.set(0,0,1).set(1,0,3).set(2,0,3);
+    cmp.set(0,1,0).set(1,1,2).set(2,1,4);
+    matrix sol(A.solve(X, B));
+    for (unsigned ro=0; ro<3; ++ro)
+        for (unsigned co=0; co<2; ++co)
+            if (cmp(ro,co) != sol(ro,co))
+                result = 1;
+    if (result) {
+        clog << "Solving " << A << " * " << X << " == " << B << endl
+             << "erroneously returned " << sol << endl;
+    }
+    
+    return result;
+}
+
 static unsigned matrix_misc(void)
 {
     unsigned result = 0;
@@ -242,6 +278,7 @@ unsigned exam_matrices(void)
     result += matrix_invert1();  cout << '.' << flush;
     result += matrix_invert2();  cout << '.' << flush;
     result += matrix_invert3();  cout << '.' << flush;
+    result += matrix_solve2();  cout << '.' << flush;
     result += matrix_misc();  cout << '.' << flush;
     
     if (!result) {
diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index 6fe0f097..2c6fe6ae 100644
--- a/check/exam_paranoia.cpp
+++ b/check/exam_paranoia.cpp
@@ -318,15 +318,21 @@ static unsigned exam_paranoia13(void)
 {
     unsigned result = 0;
 	symbol a("a"), b("b"), c("c");
-
+    
 	ex e = (b*a-c*a)/(4-a);
-	ex f = e.normal();
 	ex d = (c*a-b*a)/(a-4);
-
-	if (!(f - d).expand().is_zero()) {
-		clog << "normal(" << e << ") returns " << f << " instead of " << d << endl;
-		++result;
-	}
+    
+    try {
+        ex f = e.normal();    
+        if (!(f - d).expand().is_zero()) {
+            clog << "normal(" << e << ") returns " << f
+                 << " instead of " << d << endl;
+            ++result;
+        }
+    } catch (const exception &err) {
+        clog << "normal(" << e << ") throws " << err.what() << endl;
+        ++result;
+    }
 	return result;
 }
 
diff --git a/check/exams.h b/check/exams.h
index cd5f4b9b..416564ae 100644
--- a/check/exams.h
+++ b/check/exams.h
@@ -24,6 +24,7 @@
 #define EXAMS_H
 
 #include "ginac.h"
+using namespace std;
 
 #ifndef NO_NAMESPACE_GINAC
 using namespace GiNaC;
diff --git a/check/time_dennyfliegner.cpp b/check/time_dennyfliegner.cpp
index 0d787853..db9ff80c 100644
--- a/check/time_dennyfliegner.cpp
+++ b/check/time_dennyfliegner.cpp
@@ -31,20 +31,23 @@
 static unsigned expand_subs(unsigned size)
 {
     unsigned result = 0;
-    symbol a1("a1");
-    symbol a[VECSIZE];
+    // create a vector of size symbols named "a0", "a1", ...
+    vector<symbol> a;
+    for (unsigned i=0; i<size; ++i) {
+        char buf[4];
+        ostrstream(buf,sizeof(buf)) << i << ends;
+        a.push_back(symbol(string("a")+buf));
+    }
     ex e, aux;
     
-    a[1] = a1;
-    for (unsigned i=0; i<size; ++i) {
+    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));
+    e = pow(e,2).expand().subs(a[0]==aux).expand();
     
-    if (e != pow(a1,2)) {
+    if (e != pow(a[1],2)) {
         clog << "Denny Fliegner's quick consistency check erroneously returned "
              << e << "." << endl;
         ++result;
diff --git a/check/times.h b/check/times.h
index 74d3474e..90aac753 100644
--- a/check/times.h
+++ b/check/times.h
@@ -27,6 +27,7 @@
 #include <stdlib.h>
 #include <strstream>
 #include <vector>
+using namespace std;
 
 #include "ginac.h"
 
diff --git a/configure b/configure
index 1b903beb..d5c55789 100755
--- a/configure
+++ b/configure
@@ -545,9 +545,9 @@ fi
 
 GINACLIB_MAJOR_VERSION=0
 GINACLIB_MINOR_VERSION=6
-GINACLIB_MICRO_VERSION=3
+GINACLIB_MICRO_VERSION=4
 GINACLIB_INTERFACE_AGE=0
-GINACLIB_BINARY_AGE=0
+GINACLIB_BINARY_AGE=1
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 
diff --git a/configure.in b/configure.in
index 1f17c6a9..8d0ecaf3 100644
--- a/configure.in
+++ b/configure.in
@@ -19,9 +19,9 @@ dnl (don't we all *love* M4?)...
 
 GINACLIB_MAJOR_VERSION=0
 GINACLIB_MINOR_VERSION=6
-GINACLIB_MICRO_VERSION=3
+GINACLIB_MICRO_VERSION=4
 GINACLIB_INTERFACE_AGE=0
-GINACLIB_BINARY_AGE=0
+GINACLIB_BINARY_AGE=1
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 AC_SUBST(GINACLIB_MAJOR_VERSION)
diff --git a/doc/reference/DoxyfileHTML b/doc/reference/DoxyfileHTML
index 678525d1..aa44d49d 100644
--- a/doc/reference/DoxyfileHTML
+++ b/doc/reference/DoxyfileHTML
@@ -1,4 +1,4 @@
-# Doxyfile 1.0.0
+# Doxyfile 1.2.0
 
 # This file describes the settings to be used by doxygen for a project
 #
@@ -11,161 +11,221 @@
 # General configuration options
 #---------------------------------------------------------------------------
 
-# The PROJECT_NAME tag is a single word (or a sequence of word surrounded
+# The PROJECT_NAME tag is a single word (or a sequence of words surrounded 
 # by quotes) that should identify the project. 
 
-PROJECT_NAME         = GiNaC
+PROJECT_NAME          = GiNaC
 
-# The PROJECT_NUMBER tag can be used to enter a project or revision number.
+# The PROJECT_NUMBER tag can be used to enter a project or revision number. 
 # This could be handy for archiving the generated documentation or 
 # if some version control system is used.
 
-PROJECT_NUMBER       =
+PROJECT_NUMBER        = 
 
 # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) 
 # base path where the generated documentation will be put. 
 # If a relative path is entered, it will be relative to the location 
 # where doxygen was started. If left blank the current directory will be used.
 
-OUTPUT_DIRECTORY     = .
+OUTPUT_DIRECTORY      = .
 
-# The OUTPUT_LANGUAGE tag is used to specify the language in which all
-# documentation generated by doxygen is written. Doxygen will use this
-# information to generate all constant output in the proper language.
+# The OUTPUT_LANGUAGE tag is used to specify the language in which all 
+# documentation generated by doxygen is written. Doxygen will use this 
+# information to generate all constant output in the proper language. 
 # The default language is English, other supported languages are: 
-# Dutch, French, Italian, Czech, Swedish, German and Japanese
+# Dutch, French, Italian, Czech, Swedish, German, Finnish, Japanese, 
+# Spanish, Russian, Croatian and Polish.
 
-OUTPUT_LANGUAGE      = English
+OUTPUT_LANGUAGE       = English
 
-# The QUIET tag can be used to turn on/off the messages that are generated
-# by doxygen. Possible values are YES and NO. If left blank NO is used.
+# The DISABLE_INDEX tag can be used to turn on/off the condensed index at 
+# top of each HTML page. The value NO (the default) enables the index and 
+# the value YES disables it. 
 
-QUIET                = NO
+DISABLE_INDEX         = NO
 
-# The WARNINGS tag can be used to turn on/off the warning messages that are
-# generated by doxygen. Possible values are YES and NO. If left blank
-# NO is used.
+# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in 
+# documentation are documented, even if no documentation was available. 
+# Private class members and static file members will be hidden unless 
+# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES 
 
-WARNINGS             = YES
+EXTRACT_ALL           = YES
 
-# The DISABLE_INDEX tag can be used to turn on/off the condensed index at
-# top of each HTML page. The value NO (the default) enables the index and
-# the value YES disables it.
+# If the EXTRACT_PRIVATE tag is set to YES all private members of a class 
+# will be included in the documentation. 
 
-DISABLE_INDEX        = NO
+EXTRACT_PRIVATE       = YES
 
-# If the EXTRACT_ALL tag is set to YES all classes and functions will be
-# included in the documentation, even if no documentation was available.
+# If the EXTRACT_STATIC tag is set to YES all static members of a file 
+# will be included in the documentation. 
 
-EXTRACT_ALL          = YES
+EXTRACT_STATIC        = YES
 
-# If the EXTRACT_PRIVATE tag is set to YES all private members of a class
-# will be included in the documentation.
+# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all 
+# undocumented members of documented classes, files or namespaces. 
+# If set to NO (the default) these members will be included in the 
+# various overviews, but no documentation section is generated. 
+# This option has no effect if EXTRACT_ALL is enabled. 
 
-EXTRACT_PRIVATE      = YES
+HIDE_UNDOC_MEMBERS    = NO
 
-# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all
-# undocumented members inside documented classes or files.
+# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all 
+# undocumented classes that are normally visible in the class hierarchy. 
+# If set to NO (the default) these class will be included in the various 
+# overviews. This option has no effect if EXTRACT_ALL is enabled. 
 
-HIDE_UNDOC_MEMBERS   = NO
+HIDE_UNDOC_CLASSES    = NO
 
-# If the HIDE_UNDOC_CLASSESS tag is set to YES, Doxygen will hide all
-# undocumented classes.
-
-HIDE_UNDOC_CLASSES   = NO
-
-# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will
+# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will 
 # include brief member descriptions after the members that are listed in 
-# the file and class documentation (similar to JavaDoc).
-# Set to NO to disable this.
+# the file and class documentation (similar to JavaDoc). 
+# Set to NO to disable this. 
 
-BRIEF_MEMBER_DESC    = YES
+BRIEF_MEMBER_DESC     = YES
 
-# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend
-# the brief description of a member or function before the detailed description.
+# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend 
+# the brief description of a member or function before the detailed description. 
 # Note: if both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the 
-# brief descriptions will be completely suppressed.
+# brief descriptions will be completely suppressed. 
 
-REPEAT_BRIEF         = YES
+REPEAT_BRIEF          = YES
 
-# If the ALWAYS_DETAILS_SEC and REPEAT_BRIEF tags are both set to YES then
-# Doxygen will generate a detailed section even if there is only a brief
-# description.
+# If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then 
+# Doxygen will generate a detailed section even if there is only a brief 
+# description. 
 
-ALWAYS_DETAILED_SEC  = NO
+ALWAYS_DETAILED_SEC   = NO
 
-# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full
-# path before files name in the file list and in the header files. If set
-# to NO the shortest path that makes the file name unique will be used.
+# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full 
+# path before files name in the file list and in the header files. If set 
+# to NO the shortest path that makes the file name unique will be used. 
 
-FULL_PATH_NAMES      = NO
+FULL_PATH_NAMES       = NO
 
-# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag
-# can be used to strip a user defined part of the path. Stripping is
-# only done if one of the specified strings matches the left-hand part of
-# the path.
+# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag 
+# can be used to strip a user defined part of the path. Stripping is 
+# only done if one of the specified strings matches the left-hand part of 
+# the path. It is allowed to use relative paths in the argument list.
 
-STRIP_FROM_PATH      =
+STRIP_FROM_PATH       = 
 
-# The INTERNAL_DOCS tag determines if documentation
+# The INTERNAL_DOCS tag determines if documentation 
 # that is typed after a \internal command is included. If the tag is set 
-# to NO (the default) then the documentation will be excluded.
-# Set it to YES to include the internal documentation.
+# to NO (the default) then the documentation will be excluded. 
+# Set it to YES to include the internal documentation. 
+
+INTERNAL_DOCS         = NO
+
+# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will 
+# generate a class diagram (in Html and LaTeX) for classes with base or 
+# super classes. Setting the tag to NO turns the diagrams off. 
+
+CLASS_DIAGRAMS        = YES
+
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will 
+# be generated. Documented entities will be cross-referenced with these sources. 
+
+SOURCE_BROWSER        = YES
+
+# Setting the INLINE_SOURCES tag to YES will include the body 
+# of functions and classes directly in the documentation. 
+
+INLINE_SOURCES        = NO
+
+# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct 
+# doxygen to hide any special comment blocks from generated source code 
+# fragments. Normal C and C++ comments will always remain visible. 
 
-INTERNAL_DOCS        = NO
+STRIP_CODE_COMMENTS   = YES
 
-# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will
-# generate a class diagram (in Html and LaTeX) for classes with base or
-# super classes. Setting the tag to NO turns the diagrams off.
+# If the CASE_SENSE_NAMES tag is set to NO (the default) then Doxygen 
+# will only generate file names in lower case letters. If set to 
+# YES upper case letters are also allowed. This is useful if you have 
+# classes or files whose names only differ in case and if your file system 
+# supports case sensitive file names. 
 
-CLASS_DIAGRAMS       = YES
+CASE_SENSE_NAMES      = YES
 
-# If the SOURCE_BROWSER tag is set to YES then a list of source files will
-# be generated. Documented entities will be cross-referenced with these sources.
+# If the HIDE_SCOPE_NAMES tag is set to NO (the default) then Doxygen 
+# will show members with their full class and namespace scopes in the 
+# documentation. If set to YES the scope will be hidden. 
 
-SOURCE_BROWSER       = YES
+HIDE_SCOPE_NAMES      = NO
 
-# Setting the INLINE_SOURCES tag to YES will include the body
-# of functions and classes directly in the documentation.
+# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen 
+# will generate a verbatim copy of the header file for each class for 
+# which an include is specified. Set to NO to disable this. 
 
-INLINE_SOURCES       = NO
+VERBATIM_HEADERS      = NO
 
-# If the CASE_SENSE_NAMES tag is set to NO (the default) then Doxygen
-# will only generate file names in lower case letters. If set to
-# YES upper case letters are also allowed. This is useful if you have
-# classes or files whose names only differ in case and if your file system
-# supports case sensitive file names.
+# If the SHOW_INCLUDE_FILES tag is set to YES (the default) then Doxygen 
+# will put list of the files that are included by a file in the documentation 
+# of that file. 
 
-CASE_SENSE_NAMES     = YES
+SHOW_INCLUDE_FILES    = YES
 
-# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen
-# will generate a verbatim copy of the header file for each class for
-# which an include is specified. Set to NO to disable this.
+# If the JAVADOC_AUTOBRIEF tag is set to YES (the default) then Doxygen 
+# will interpret the first line (until the first dot) of a JavaDoc-style 
+# comment as the brief description. If set to NO, the Javadoc-style will 
+# behave just like the Qt-style comments. 
 
-VERBATIM_HEADERS     = NO
+JAVADOC_AUTOBRIEF     = YES
 
-# If the JAVADOC_AUTOBRIEF tag is set to YES (the default) then Doxygen
-# will interpret the first line (until the first dot) of a JavaDoc-style
-# comment as the brief description. If set to NO, the Javadoc-style will
-# behave just like the Qt-style comments.
+# If the INHERIT_DOCS tag is set to YES (the default) then an undocumented 
+# member inherits the documentation from any documented member that it 
+# reimplements. 
 
-JAVADOC_AUTOBRIEF    = YES
+INHERIT_DOCS          = YES
 
-# if the INHERIT_DOCS tag is set to YES (the default) then an undocumented
-# member inherits the documentation from any documented member that it
-# reimplements.
+# If the INLINE_INFO tag is set to YES (the default) then a tag [inline] 
+# is inserted in the documentation for inline members. 
 
-INHERIT_DOCS         = YES
+INLINE_INFO           = YES
 
-# if the INLINE_INFO tag is set to YES (the default) then a tag [inline]
-# is inserted in the documentation for inline members.
+# If the SORT_MEMBER_DOCS tag is set to YES (the default) then doxygen 
+# will sort the (detailed) documentation of file and class members 
+# alphabetically by member name. If set to NO the members will appear in 
+# declaration order. 
 
-INLINE_INFO          = YES
+SORT_MEMBER_DOCS      = NO
 
-# the TAB_SIZE tag can be used to set the number of spaces in a tab.
-# Doxygen uses this value to replace tabs by spaces in code fragments.
+# The TAB_SIZE tag can be used to set the number of spaces in a tab. 
+# Doxygen uses this value to replace tabs by spaces in code fragments. 
 
-TAB_SIZE             = 4
+TAB_SIZE              = 4
+
+# The ENABLE_SECTIONS tag can be used to enable conditional 
+# documentation sections, marked by \if sectionname ... \endif. 
+
+ENABLED_SECTIONS      = 
+
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+
+# The QUIET tag can be used to turn on/off the messages that are generated 
+# by doxygen. Possible values are YES and NO. If left blank NO is used. 
+
+QUIET                 = NO
+
+# The WARNINGS tag can be used to turn on/off the warning messages that are 
+# generated by doxygen. Possible values are YES and NO. If left blank 
+# NO is used. 
+
+WARNINGS              = YES
+
+# If WARN_IF_UNDOCUMENTED is set to YES, then doxygen will generate warnings 
+# for undocumented members. If EXTRACT_ALL is set to YES then this flag will 
+# automatically be disabled. 
+
+WARN_IF_UNDOCUMENTED  = YES
+
+# The WARN_FORMAT tag determines the format of the warning messages that 
+# doxygen can produce. The string should contain the $file, $line, and $text 
+# tags, which will be replaced by the file and line number from which the 
+# warning originated and the warning text. 
+
+WARN_FORMAT           = "$file:$line: $text"
 
 #---------------------------------------------------------------------------
 # configuration options related to the input files
@@ -174,293 +234,428 @@ TAB_SIZE             = 4
 # The INPUT tag can be used to specify the files and/or directories that contain 
 # documented source files. You may enter file names like "myfile.cpp" or 
 # directories like "/usr/src/myproject". Separate the files or directories 
-# with spaces.
+# with spaces. 
 
-INPUT                = ../../ginac
+INPUT                 = ../../ginac
 
 # If the value of the INPUT tag contains directories, you can use the 
 # FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp 
 # and *.h) to filter out the source-files in the directories. If left 
-# blank all files are included.
+# blank all files are included. 
 
-FILE_PATTERNS        = *.cpp *.h
+FILE_PATTERNS         = *.cpp *.h
 
-# The RECURSIVE tag can be used to turn specify whether or not subdirectories
-# should be searched for input files as well. Possible values are YES and NO.
-# If left blank NO is used.
+# The RECURSIVE tag can be used to turn specify whether or not subdirectories 
+# should be searched for input files as well. Possible values are YES and NO. 
+# If left blank NO is used. 
 
-RECURSIVE            = NO
+RECURSIVE             = NO
 
-# The EXCLUDE tag can be used to specify files and/or directories that should
+# The EXCLUDE tag can be used to specify files and/or directories that should 
 # excluded from the INPUT source files. This way you can easily exclude a 
-# subdirectory from a directory tree whose root is specified with the INPUT tag.
+# subdirectory from a directory tree whose root is specified with the INPUT tag. 
 
-EXCLUDE              =
+EXCLUDE               = 
 
-# If the value of the INPUT tag contains directories, you can use the
-# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude
-# certain files from those directories.
+# If the value of the INPUT tag contains directories, you can use the 
+# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude 
+# certain files from those directories. 
 
-EXCLUDE_PATTERNS     =
+EXCLUDE_PATTERNS      = 
 
 # The EXAMPLE_PATH tag can be used to specify one or more files or 
 # directories that contain example code fragments that are included (see 
-# the \include command).
+# the \include command). 
 
-EXAMPLE_PATH         =
+EXAMPLE_PATH          = 
 
-# If the value of the EXAMPLE_PATH tag contains directories, you can use the
+# If the value of the EXAMPLE_PATH tag contains directories, you can use the 
 # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp 
 # and *.h) to filter out the source-files in the directories. If left 
-# blank all files are included.
+# blank all files are included. 
 
-EXAMPLE_PATTERNS     =
+EXAMPLE_PATTERNS      = 
 
 # The IMAGE_PATH tag can be used to specify one or more files or 
 # directories that contain image that are included in the documentation (see 
-# the \image command).
+# the \image command). 
 
-IMAGE_PATH           =
+IMAGE_PATH            = 
 
-# The INPUT_FILTER tag can be used to specify a program that doxygen should
+# The INPUT_FILTER tag can be used to specify a program that doxygen should 
 # invoke to filter for each input file. Doxygen will invoke the filter program 
-# by executing (via popen()) the command <filter> <input-file>, where <filter>
-# is the value of the INPUT_FILTER tag, and <input-file> is the name of an
-# input file. Doxygen will then use the output that the filter program writes
-# to standard output.
+# by executing (via popen()) the command <filter> <input-file>, where <filter> 
+# is the value of the INPUT_FILTER tag, and <input-file> is the name of an 
+# input file. Doxygen will then use the output that the filter program writes 
+# to standard output. 
 
-INPUT_FILTER         =
+INPUT_FILTER          = 
+
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+
+# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index 
+# of all compounds will be generated. Enable this if the project 
+# contains a lot of classes, structs, unions or interfaces. 
+
+ALPHABETICAL_INDEX    = YES
+
+# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then 
+# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns 
+# in which this list will be split (can be a number in the range [1..20]) 
+
+COLS_IN_ALPHA_INDEX   = 5
+
+# In case all classes in a project start with a common prefix, all 
+# classes will be put under the same header in the alphabetical index. 
+# The IGNORE_PREFIX tag can be used to specify one or more prefixes that 
+# should be ignored while generating the index headers. 
+
+IGNORE_PREFIX         = 
 
 #---------------------------------------------------------------------------
 # configuration options related to the HTML output
 #---------------------------------------------------------------------------
 
-# If the GENERATE_HTML tag is set to YES (the default) Doxygen will
-# generate HTML output
+# If the GENERATE_HTML tag is set to YES (the default) Doxygen will 
+# generate HTML output. 
 
-GENERATE_HTML        = YES
+GENERATE_HTML         = YES
 
-# The HTML_OUTPUT tag is used to specify where the HTML docs will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be
-# put in front of it. If left blank `html' will be used as the default path.
+# The HTML_OUTPUT tag is used to specify where the HTML docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `html' will be used as the default path. 
 
-HTML_OUTPUT          = .
+HTML_OUTPUT           = .
 
 # The HTML_HEADER tag can be used to specify a personal HTML header for 
 # each generated HTML page. If it is left blank doxygen will generate a 
 # standard header.
 
-HTML_HEADER          =
+HTML_HEADER           = 
 
 # The HTML_FOOTER tag can be used to specify a personal HTML footer for 
 # each generated HTML page. If it is left blank doxygen will generate a 
 # standard footer.
 
-HTML_FOOTER          = Doxyfooter
+HTML_FOOTER           = Doxyfooter
 
-# The HTML_STYLESHEET tag can be used to specify a user defined cascading
+# The HTML_STYLESHEET tag can be used to specify a user defined cascading 
 # style sheet that is used by each HTML page. It can be used to 
-# fine-tune the look of the HTML output. If the tag is left blank doxygen
-# will generate a default style sheet
-
-HTML_STYLESHEET      =
-
-# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes,
-# files or namespaces will be aligned in HTML using tables. If set to
-# NO a bullet list will be used.
+# fine-tune the look of the HTML output. If the tag is left blank doxygen 
+# will generate a default style sheet 
 
-HTML_ALIGN_MEMBERS   = YES
+HTML_STYLESHEET       = 
 
-# If the GENERATE_HTMLHELP tag is set to YES, additional index files
-# will be generated that can be used as input for tools like the
-# Microsoft HTML help workshop to generate a compressed HTML help file (.chm)
-# of the generated HTML documentation.
+# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes, 
+# files or namespaces will be aligned in HTML using tables. If set to 
+# NO a bullet list will be used. 
 
-GENERATE_HTMLHELP    = NO
+HTML_ALIGN_MEMBERS    = YES
 
-# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index
-# of all compounds will be generated. Enable this if the project
-# contains a lot of classes, structs, unions or interfaces.
+# If the GENERATE_HTMLHELP tag is set to YES, additional index files 
+# will be generated that can be used as input for tools like the 
+# Microsoft HTML help workshop to generate a compressed HTML help file (.chm) 
+# of the generated HTML documentation. 
 
-ALPHABETICAL_INDEX   = YES
-
-# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then
-# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns
-# in which this list will be split (can be a number in the range [1..20])
-
-COLS_IN_ALPHA_INDEX  = 5
+GENERATE_HTMLHELP     = NO
 
 #---------------------------------------------------------------------------
 # configuration options related to the LaTeX output
 #---------------------------------------------------------------------------
 
-# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
-# generate Latex output.
+# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will 
+# generate Latex output. 
 
-GENERATE_LATEX       = NO
+GENERATE_LATEX        = NO
 
-# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be
-# put in front of it. If left blank `latex' will be used as the default path.
+# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `latex' will be used as the default path. 
 
-LATEX_OUTPUT         = latex
+LATEX_OUTPUT          = latex
 
-# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact
-# LaTeX documents. This may be useful for small projects and may help to
-# save some trees in general.
+# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact 
+# LaTeX documents. This may be useful for small projects and may help to 
+# save some trees in general. 
 
-COMPACT_LATEX        = NO
+COMPACT_LATEX         = NO
 
-# The PAPER_TYPE tag can be used to set the paper type that is used
+# The PAPER_TYPE tag can be used to set the paper type that is used 
 # by the printer. Possible values are: a4, a4wide, letter, legal and 
-# executive. If left blank a4wide will be used.
+# executive. If left blank a4wide will be used. 
 
-PAPER_TYPE           = a4wide
+PAPER_TYPE            = a4wide
 
-# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
-# packages that should be included in the LaTeX output.
+# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX 
+# packages that should be included in the LaTeX output. 
 
-EXTRA_PACKAGES       =
+EXTRA_PACKAGES        = 
 
 # The LATEX_HEADER tag can be used to specify a personal LaTeX header for 
-# the generated latex document. The header should contain everything until
+# the generated latex document. The header should contain everything until 
 # the first chapter. If it is left blank doxygen will generate a 
-# standard header. Notice: only use this tag if you know what you are doing!
+# standard header. Notice: only use this tag if you know what you are doing! 
+
+LATEX_HEADER          = 
+
+# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated 
+# is prepared for conversion to pdf (using ps2pdf). The pdf file will 
+# contain links (just like the HTML output) instead of page references 
+# This makes the output suitable for online browsing using a pdf viewer. 
 
-LATEX_HEADER         =
+PDF_HYPERLINKS        = NO
 
-# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated
-# is prepared for conversion to pdf (using ps2pdf). The pdf file will
-# contain links (just like the HTML output) instead of page references
-# This makes the output suitable for online browsing using a pdf viewer.
+# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \\batchmode. 
+# command to the generated LaTeX files. This will instruct LaTeX to keep 
+# running if errors occur, instead of asking the user for help. 
+# This option is also used when generating formulas in HTML. 
 
-PDF_HYPERLINKS       = NO
+LATEX_BATCHMODE       = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_RTF tag is set to YES Doxygen will generate RTF output 
+# For now this is experimental and is disabled by default. The RTF output 
+# is optimised for Word 97 and may not look too pretty with other readers 
+# or editors.
+
+GENERATE_RTF          = NO
+
+# The RTF_OUTPUT tag is used to specify where the RTF docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `rtf' will be used as the default path. 
+
+RTF_OUTPUT            = rtf
+
+# If the COMPACT_RTF tag is set to YES Doxygen generates more compact 
+# RTF documents. This may be useful for small projects and may help to 
+# save some trees in general. 
+
+COMPACT_RTF           = NO
+
+# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated 
+# will contain hyperlink fields. The RTF file will 
+# contain links (just like the HTML output) instead of page references. 
+# This makes the output suitable for online browsing using a WORD or other. 
+# programs which support those fields. 
+# Note: wordpad (write) and others do not support links. 
+
+RTF_HYPERLINKS        = NO
+
+# Load stylesheet definitions from file. Syntax is similar to doxygen's 
+# config file, i.e. a series of assigments. You only have to provide 
+# replacements, missing definitions are set to their default value. 
+
+RTF_STYLESHEET_FILE   = 
 
 #---------------------------------------------------------------------------
 # configuration options related to the man page output
 #---------------------------------------------------------------------------
 
-# If the GENERATE_MAN tag is set to YES (the default) Doxygen will
-# generate man pages
+# If the GENERATE_MAN tag is set to YES (the default) Doxygen will 
+# generate man pages 
 
-GENERATE_MAN         = NO
+GENERATE_MAN          = NO
 
-# The MAN_OUTPUT tag is used to specify where the man pages will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be
-# put in front of it. If left blank `man' will be used as the default path.
+# The MAN_OUTPUT tag is used to specify where the man pages will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `man' will be used as the default path. 
 
-MAN_OUTPUT           =
+MAN_OUTPUT            = man
 
-# The MAN_EXTENSION tag determines the extension that is added to
-# the generated man pages (default is the subroutine's section .3)
+# The MAN_EXTENSION tag determines the extension that is added to 
+# the generated man pages (default is the subroutine's section .3) 
 
-MAN_EXTENSION        = .3
+MAN_EXTENSION         = .3
 
 #---------------------------------------------------------------------------
-# Configuration options related to the preprocessor 
+# Configuration options related to the preprocessor   
 #---------------------------------------------------------------------------
 
-# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will
-# evaluate all C-preprocessor directives found in the sources and include
-# files.
+# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will 
+# evaluate all C-preprocessor directives found in the sources and include 
+# files. 
 
-ENABLE_PREPROCESSING = YES
+ENABLE_PREPROCESSING  = YES
 
-# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro
+# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro 
 # names in the source code. If set to NO (the default) only conditional 
-# compilation will be performed.
+# compilation will be performed. Macro expansion can be done in a controlled 
+# way by setting EXPAND_ONLY_PREDEF to YES. 
+
+MACRO_EXPANSION       = NO
+
+# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES 
+# then the macro expansion is limited to the macros specified with the 
+# PREDEFINED and EXPAND_AS_PREDEFINED tags. 
 
-MACRO_EXPANSION      = NO
+EXPAND_ONLY_PREDEF    = NO
 
-# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files
-# in the INCLUDE_PATH (see below) will be search if a #include is found.
+# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files 
+# in the INCLUDE_PATH (see below) will be search if a #include is found. 
 
-SEARCH_INCLUDES      = YES
+SEARCH_INCLUDES       = YES
 
-# The INCLUDE_PATH tag can be used to specify one or more directories that
-# contain include files that are not input files but should be processed by
-# the preprocessor.
+# The INCLUDE_PATH tag can be used to specify one or more directories that 
+# contain include files that are not input files but should be processed by 
+# the preprocessor. 
 
-INCLUDE_PATH         =
+INCLUDE_PATH          = 
 
-# The PREDEFINED tag can be used to specify one or more macro names that
-# are defined before the preprocessor is started (similar to the -D option of
-# gcc). The argument of the tag is a list of macros of the form: name
+# You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard 
+# patterns (like *.h and *.hpp) to filter out the header-files in the 
+# directories. If left blank, the patterns specified with FILE_PATTERNS will 
+# be used. 
+
+INCLUDE_FILE_PATTERNS = *.h
+
+# The PREDEFINED tag can be used to specify one or more macro names that 
+# are defined before the preprocessor is started (similar to the -D option of 
+# gcc). The argument of the tag is a list of macros of the form: name 
 # or name=definition (no spaces). If the definition and the = are 
-# omitted =1 is assumed.
+# omitted =1 is assumed. 
 
-PREDEFINED           = NO_NAMESPACE_GINAC
+PREDEFINED            = NO_NAMESPACE_GINAC
 
-# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES
-# then the macro expansion is limited to the macros specified with the
-# PREDEFINED tag.
+# If the MACRO_EXPANSION and EXPAND_PREDEF_ONLY tags are set to YES then 
+# this tag can be used to specify a list of macro names that should be expanded. 
+# The macro definition that is found in the sources will be used. 
+# Use the PREDEFINED tag if you want to use a different macro definition. 
 
-EXPAND_ONLY_PREDEF   = NO
+EXPAND_AS_DEFINED     = 
 
 #---------------------------------------------------------------------------
-# Configuration options related to external references 
+# Configuration::addtions related to external references   
 #---------------------------------------------------------------------------
 
 # The TAGFILES tag can be used to specify one or more tagfiles. 
 
-TAGFILES             =
+TAGFILES              = 
+
+# When a file name is specified after GENERATE_TAGFILE, doxygen will create 
+# a tag file that is based on the input files it reads. 
+
+GENERATE_TAGFILE      = 
+
+# If the ALLEXTERNALS tag is set to YES all external classes will be listed 
+# in the class index. If set to NO only the inherited external classes 
+# will be listed. 
+
+ALLEXTERNALS          = NO
+
+# The PERL_PATH should be the absolute path and name of the perl script 
+# interpreter (i.e. the result of `which perl'). 
+
+PERL_PATH             = /usr/bin/perl
+
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool   
+#---------------------------------------------------------------------------
+
+# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is 
+# available from the path. This tool is part of Graphviz, a graph visualization 
+# toolkit from AT&T and Lucent Bell Labs. The other options in this section 
+# have no effect if this option is set to NO (the default) 
+
+HAVE_DOT              = YES
+
+# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for each documented class showing the direct and 
+# indirect inheritance relations. Setting this tag to YES will force the 
+# the CLASS_DIAGRAMS tag to NO.
+
+CLASS_GRAPH           = YES
+
+# If the COLLABORATION_GRAPH and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for each documented class showing the direct and 
+# indirect implementation dependencies (inheritance, containment, and 
+# class references variables) of the class with other documented classes. 
+
+COLLABORATION_GRAPH   = YES
+
+# If the ENABLE_PREPROCESSING, INCLUDE_GRAPH, and HAVE_DOT tags are set to 
+# YES then doxygen will generate a graph for each documented file showing 
+# the direct and indirect include dependencies of the file with other 
+# documented files. 
+
+INCLUDE_GRAPH         = YES
+
+# If the ENABLE_PREPROCESSING, INCLUDED_BY_GRAPH, and HAVE_DOT tags are set to 
+# YES then doxygen will generate a graph for each documented header file showing 
+# the documented files that directly or indirectly include this file 
+
+INCLUDED_BY_GRAPH     = YES
+
+# If the GRAPHICAL_HIERARCHY and HAVE_DOT tags are set to YES then doxygen 
+# will graphical hierarchy of all classes instead of a textual one. 
+
+GRAPHICAL_HIERARCHY   = YES
 
-# When a file name is specified after GENERATE_TAGFILE, doxygen will create
-# a tag file that is based on the input files it reads.
+# The tag DOT_PATH can be used to specify the path where the dot tool can be 
+# found. If left blank, it is assumed the dot tool can be found on the path. 
 
-GENERATE_TAGFILE     =
+DOT_PATH              = 
 
-# If the ALLEXTERNALS tag is set to YES all external classes will be listed
-# in the class index. If set to NO only the inherited external classes
-# will be listed.
+# The MAX_DOT_GRAPH_WIDTH tag can be used to set the maximum allowed width 
+# (in pixels) of the graphs generated by dot. If a graph becomes larger than 
+# this value, doxygen will try to truncate the graph, so that it fits within 
+# the specified constraint. Beware that most browsers cannot cope with very 
+# large images. 
 
-ALLEXTERNALS         = NO
+MAX_DOT_GRAPH_WIDTH   = 1024
 
-# The PERL_PATH should be the absolute path and name of the perl script
-# interpreter (i.e. the result of `which perl').
+# The MAX_DOT_GRAPH_HEIGHT tag can be used to set the maximum allows height 
+# (in pixels) of the graphs generated by dot. If a graph becomes larger than 
+# this value, doxygen will try to truncate the graph, so that it fits within 
+# the specified constraint. Beware that most browsers cannot cope with very 
+# large images. 
 
-PERL_PATH            = /usr/bin/perl
+MAX_DOT_GRAPH_HEIGHT  = 1024
 
 #---------------------------------------------------------------------------
-# Configuration options related to the search engine 
+# Configuration::addtions related to the search engine   
 #---------------------------------------------------------------------------
 
 # The SEARCHENGINE tag specifies whether or not a search engine should be 
-# used. If set to NO the values of all tags below this one will be ignored.
+# used. If set to NO the values of all tags below this one will be ignored. 
 
-SEARCHENGINE         = NO
+SEARCHENGINE          = NO
 
-# The CGI_NAME tag should be the name of the CGI script that
-# starts the search engine (doxysearch) with the correct parameters.
-# A script with this name will be generated by doxygen.
+# The CGI_NAME tag should be the name of the CGI script that 
+# starts the search engine (doxysearch) with the correct parameters. 
+# A script with this name will be generated by doxygen. 
 
-CGI_NAME             = search.cgi
+CGI_NAME              = search.cgi
 
-# The CGI_URL tag should be the absolute URL to the directory where the
+# The CGI_URL tag should be the absolute URL to the directory where the 
 # cgi binaries are located. See the documentation of your http daemon for 
-# details.
+# details. 
 
-CGI_URL              =
+CGI_URL               = 
 
-# The DOC_URL tag should be the absolute URL to the directory where the
+# The DOC_URL tag should be the absolute URL to the directory where the 
 # documentation is located. If left blank the absolute path to the 
-# documentation, with file:// prepended to it, will be used.
+# documentation, with file:// prepended to it, will be used. 
 
-DOC_URL              =
+DOC_URL               = 
 
-# The DOC_ABSPATH tag should be the absolute path to the directory where the
-# documentation is located. If left blank the directory on the local machine
-# will be used.
+# The DOC_ABSPATH tag should be the absolute path to the directory where the 
+# documentation is located. If left blank the directory on the local machine 
+# will be used. 
 
-DOC_ABSPATH          =
+DOC_ABSPATH           = 
 
-# The BIN_ABSPATH tag must point to the directory where the doxysearch binary
-# is installed.
+# The BIN_ABSPATH tag must point to the directory where the doxysearch binary 
+# is installed. 
 
-BIN_ABSPATH          = /usr/local/bin/
+BIN_ABSPATH           = /usr/local/bin/
 
 # The EXT_DOC_PATHS tag can be used to specify one or more paths to 
-# documentation generated for other projects. This allows doxysearch to search
-# the documentation for these projects as well.
+# documentation generated for other projects. This allows doxysearch to search 
+# the documentation for these projects as well. 
 
-EXT_DOC_PATHS        =
+EXT_DOC_PATHS         = 
diff --git a/doc/reference/DoxyfileTEX b/doc/reference/DoxyfileTEX
index 9a97db23..2aa3ef42 100644
--- a/doc/reference/DoxyfileTEX
+++ b/doc/reference/DoxyfileTEX
@@ -1,4 +1,4 @@
-# Doxyfile 1.0.0
+# Doxyfile 1.2.0
 
 # This file describes the settings to be used by doxygen for a project
 #
@@ -11,161 +11,221 @@
 # General configuration options
 #---------------------------------------------------------------------------
 
-# The PROJECT_NAME tag is a single word (or a sequence of word surrounded
+# The PROJECT_NAME tag is a single word (or a sequence of words surrounded 
 # by quotes) that should identify the project. 
 
-PROJECT_NAME         = GiNaC
+PROJECT_NAME          = GiNaC
 
-# The PROJECT_NUMBER tag can be used to enter a project or revision number.
+# The PROJECT_NUMBER tag can be used to enter a project or revision number. 
 # This could be handy for archiving the generated documentation or 
 # if some version control system is used.
 
-PROJECT_NUMBER       =
+PROJECT_NUMBER        = 
 
 # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) 
 # base path where the generated documentation will be put. 
 # If a relative path is entered, it will be relative to the location 
 # where doxygen was started. If left blank the current directory will be used.
 
-OUTPUT_DIRECTORY     = .
+OUTPUT_DIRECTORY      = .
 
-# The OUTPUT_LANGUAGE tag is used to specify the language in which all
-# documentation generated by doxygen is written. Doxygen will use this
-# information to generate all constant output in the proper language.
+# The OUTPUT_LANGUAGE tag is used to specify the language in which all 
+# documentation generated by doxygen is written. Doxygen will use this 
+# information to generate all constant output in the proper language. 
 # The default language is English, other supported languages are: 
-# Dutch, French, Italian, Czech, Swedish, German and Japanese
+# Dutch, French, Italian, Czech, Swedish, German, Finnish, Japanese, 
+# Spanish, Russian, Croatian and Polish.
 
-OUTPUT_LANGUAGE      = English
+OUTPUT_LANGUAGE       = English
 
-# The QUIET tag can be used to turn on/off the messages that are generated
-# by doxygen. Possible values are YES and NO. If left blank NO is used.
+# The DISABLE_INDEX tag can be used to turn on/off the condensed index at 
+# top of each HTML page. The value NO (the default) enables the index and 
+# the value YES disables it. 
 
-QUIET                = NO
+DISABLE_INDEX         = NO
 
-# The WARNINGS tag can be used to turn on/off the warning messages that are
-# generated by doxygen. Possible values are YES and NO. If left blank
-# NO is used.
+# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in 
+# documentation are documented, even if no documentation was available. 
+# Private class members and static file members will be hidden unless 
+# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES 
 
-WARNINGS             = YES
+EXTRACT_ALL           = 
 
-# The DISABLE_INDEX tag can be used to turn on/off the condensed index at
-# top of each HTML page. The value NO (the default) enables the index and
-# the value YES disables it.
+# If the EXTRACT_PRIVATE tag is set to YES all private members of a class 
+# will be included in the documentation. 
 
-DISABLE_INDEX        = NO
+EXTRACT_PRIVATE       = YES
 
-# If the EXTRACT_ALL tag is set to YES all classes and functions will be
-# included in the documentation, even if no documentation was available.
+# If the EXTRACT_STATIC tag is set to YES all static members of a file 
+# will be included in the documentation. 
 
-EXTRACT_ALL          = NO
+EXTRACT_STATIC        = YES
 
-# If the EXTRACT_PRIVATE tag is set to YES all private members of a class
-# will be included in the documentation.
+# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all 
+# undocumented members of documented classes, files or namespaces. 
+# If set to NO (the default) these members will be included in the 
+# various overviews, but no documentation section is generated. 
+# This option has no effect if EXTRACT_ALL is enabled. 
 
-EXTRACT_PRIVATE      = YES
+HIDE_UNDOC_MEMBERS    = NO
 
-# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all
-# undocumented members inside documented classes or files.
+# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all 
+# undocumented classes that are normally visible in the class hierarchy. 
+# If set to NO (the default) these class will be included in the various 
+# overviews. This option has no effect if EXTRACT_ALL is enabled. 
 
-HIDE_UNDOC_MEMBERS   = NO
+HIDE_UNDOC_CLASSES    = NO
 
-# If the HIDE_UNDOC_CLASSESS tag is set to YES, Doxygen will hide all
-# undocumented classes.
-
-HIDE_UNDOC_CLASSES   = NO
-
-# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will
+# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will 
 # include brief member descriptions after the members that are listed in 
-# the file and class documentation (similar to JavaDoc).
-# Set to NO to disable this.
+# the file and class documentation (similar to JavaDoc). 
+# Set to NO to disable this. 
 
-BRIEF_MEMBER_DESC    = YES
+BRIEF_MEMBER_DESC     = YES
 
-# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend
-# the brief description of a member or function before the detailed description.
+# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend 
+# the brief description of a member or function before the detailed description. 
 # Note: if both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the 
-# brief descriptions will be completely suppressed.
+# brief descriptions will be completely suppressed. 
 
-REPEAT_BRIEF         = YES
+REPEAT_BRIEF          = YES
 
-# If the ALWAYS_DETAILS_SEC and REPEAT_BRIEF tags are both set to YES then
-# Doxygen will generate a detailed section even if there is only a brief
-# description.
+# If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then 
+# Doxygen will generate a detailed section even if there is only a brief 
+# description. 
 
-ALWAYS_DETAILED_SEC  = NO
+ALWAYS_DETAILED_SEC   = NO
 
-# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full
-# path before files name in the file list and in the header files. If set
-# to NO the shortest path that makes the file name unique will be used.
+# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full 
+# path before files name in the file list and in the header files. If set 
+# to NO the shortest path that makes the file name unique will be used. 
 
-FULL_PATH_NAMES      = NO
+FULL_PATH_NAMES       = NO
 
-# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag
-# can be used to strip a user defined part of the path. Stripping is
-# only done if one of the specified strings matches the left-hand part of
-# the path.
+# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag 
+# can be used to strip a user defined part of the path. Stripping is 
+# only done if one of the specified strings matches the left-hand part of 
+# the path. It is allowed to use relative paths in the argument list.
 
-STRIP_FROM_PATH      =
+STRIP_FROM_PATH       = 
 
-# The INTERNAL_DOCS tag determines if documentation
+# The INTERNAL_DOCS tag determines if documentation 
 # that is typed after a \internal command is included. If the tag is set 
-# to NO (the default) then the documentation will be excluded.
-# Set it to YES to include the internal documentation.
+# to NO (the default) then the documentation will be excluded. 
+# Set it to YES to include the internal documentation. 
+
+INTERNAL_DOCS         = NO
+
+# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will 
+# generate a class diagram (in Html and LaTeX) for classes with base or 
+# super classes. Setting the tag to NO turns the diagrams off. 
+
+CLASS_DIAGRAMS        = YES
+
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will 
+# be generated. Documented entities will be cross-referenced with these sources. 
+
+SOURCE_BROWSER        = 
+
+# Setting the INLINE_SOURCES tag to YES will include the body 
+# of functions and classes directly in the documentation. 
+
+INLINE_SOURCES        = NO
+
+# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct 
+# doxygen to hide any special comment blocks from generated source code 
+# fragments. Normal C and C++ comments will always remain visible. 
 
-INTERNAL_DOCS        = NO
+STRIP_CODE_COMMENTS   = YES
 
-# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will
-# generate a class diagram (in Html and LaTeX) for classes with base or
-# super classes. Setting the tag to NO turns the diagrams off.
+# If the CASE_SENSE_NAMES tag is set to NO (the default) then Doxygen 
+# will only generate file names in lower case letters. If set to 
+# YES upper case letters are also allowed. This is useful if you have 
+# classes or files whose names only differ in case and if your file system 
+# supports case sensitive file names. 
 
-CLASS_DIAGRAMS       = YES
+CASE_SENSE_NAMES      = YES
 
-# If the SOURCE_BROWSER tag is set to YES then a list of source files will
-# be generated. Documented entities will be cross-referenced with these sources.
+# If the HIDE_SCOPE_NAMES tag is set to NO (the default) then Doxygen 
+# will show members with their full class and namespace scopes in the 
+# documentation. If set to YES the scope will be hidden. 
 
-SOURCE_BROWSER       = NO
+HIDE_SCOPE_NAMES      = NO
 
-# Setting the INLINE_SOURCES tag to YES will include the body
-# of functions and classes directly in the documentation.
+# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen 
+# will generate a verbatim copy of the header file for each class for 
+# which an include is specified. Set to NO to disable this. 
 
-INLINE_SOURCES       = NO
+VERBATIM_HEADERS      = NO
 
-# If the CASE_SENSE_NAMES tag is set to NO (the default) then Doxygen
-# will only generate file names in lower case letters. If set to
-# YES upper case letters are also allowed. This is useful if you have
-# classes or files whose names only differ in case and if your file system
-# supports case sensitive file names.
+# If the SHOW_INCLUDE_FILES tag is set to YES (the default) then Doxygen 
+# will put list of the files that are included by a file in the documentation 
+# of that file. 
 
-CASE_SENSE_NAMES     = YES
+SHOW_INCLUDE_FILES    = YES
 
-# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen
-# will generate a verbatim copy of the header file for each class for
-# which an include is specified. Set to NO to disable this.
+# If the JAVADOC_AUTOBRIEF tag is set to YES (the default) then Doxygen 
+# will interpret the first line (until the first dot) of a JavaDoc-style 
+# comment as the brief description. If set to NO, the Javadoc-style will 
+# behave just like the Qt-style comments. 
 
-VERBATIM_HEADERS     = NO
+JAVADOC_AUTOBRIEF     = YES
 
-# If the JAVADOC_AUTOBRIEF tag is set to YES (the default) then Doxygen
-# will interpret the first line (until the first dot) of a JavaDoc-style
-# comment as the brief description. If set to NO, the Javadoc-style will
-# behave just like the Qt-style comments.
+# If the INHERIT_DOCS tag is set to YES (the default) then an undocumented 
+# member inherits the documentation from any documented member that it 
+# reimplements. 
 
-JAVADOC_AUTOBRIEF    = YES
+INHERIT_DOCS          = YES
 
-# if the INHERIT_DOCS tag is set to YES (the default) then an undocumented
-# member inherits the documentation from any documented member that it
-# reimplements.
+# If the INLINE_INFO tag is set to YES (the default) then a tag [inline] 
+# is inserted in the documentation for inline members. 
 
-INHERIT_DOCS         = YES
+INLINE_INFO           = YES
 
-# if the INLINE_INFO tag is set to YES (the default) then a tag [inline]
-# is inserted in the documentation for inline members.
+# If the SORT_MEMBER_DOCS tag is set to YES (the default) then doxygen 
+# will sort the (detailed) documentation of file and class members 
+# alphabetically by member name. If set to NO the members will appear in 
+# declaration order. 
 
-INLINE_INFO          = YES
+SORT_MEMBER_DOCS      = NO
 
-# the TAB_SIZE tag can be used to set the number of spaces in a tab.
-# Doxygen uses this value to replace tabs by spaces in code fragments.
+# The TAB_SIZE tag can be used to set the number of spaces in a tab. 
+# Doxygen uses this value to replace tabs by spaces in code fragments. 
 
-TAB_SIZE             = 4
+TAB_SIZE              = 4
+
+# The ENABLE_SECTIONS tag can be used to enable conditional 
+# documentation sections, marked by \if sectionname ... \endif. 
+
+ENABLED_SECTIONS      = 
+
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+
+# The QUIET tag can be used to turn on/off the messages that are generated 
+# by doxygen. Possible values are YES and NO. If left blank NO is used. 
+
+QUIET                 = NO
+
+# The WARNINGS tag can be used to turn on/off the warning messages that are 
+# generated by doxygen. Possible values are YES and NO. If left blank 
+# NO is used. 
+
+WARNINGS              = YES
+
+# If WARN_IF_UNDOCUMENTED is set to YES, then doxygen will generate warnings 
+# for undocumented members. If EXTRACT_ALL is set to YES then this flag will 
+# automatically be disabled. 
+
+WARN_IF_UNDOCUMENTED  = YES
+
+# The WARN_FORMAT tag determines the format of the warning messages that 
+# doxygen can produce. The string should contain the $file, $line, and $text 
+# tags, which will be replaced by the file and line number from which the 
+# warning originated and the warning text. 
+
+WARN_FORMAT           = "$file:$line: $text"
 
 #---------------------------------------------------------------------------
 # configuration options related to the input files
@@ -174,293 +234,428 @@ TAB_SIZE             = 4
 # The INPUT tag can be used to specify the files and/or directories that contain 
 # documented source files. You may enter file names like "myfile.cpp" or 
 # directories like "/usr/src/myproject". Separate the files or directories 
-# with spaces.
+# with spaces. 
 
-INPUT                = ../../ginac
+INPUT                 = ../../ginac
 
 # If the value of the INPUT tag contains directories, you can use the 
 # FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp 
 # and *.h) to filter out the source-files in the directories. If left 
-# blank all files are included.
+# blank all files are included. 
 
-FILE_PATTERNS        = *.cpp *.h
+FILE_PATTERNS         = *.cpp *.h
 
-# The RECURSIVE tag can be used to turn specify whether or not subdirectories
-# should be searched for input files as well. Possible values are YES and NO.
-# If left blank NO is used.
+# The RECURSIVE tag can be used to turn specify whether or not subdirectories 
+# should be searched for input files as well. Possible values are YES and NO. 
+# If left blank NO is used. 
 
-RECURSIVE            = NO
+RECURSIVE             = NO
 
-# The EXCLUDE tag can be used to specify files and/or directories that should
+# The EXCLUDE tag can be used to specify files and/or directories that should 
 # excluded from the INPUT source files. This way you can easily exclude a 
-# subdirectory from a directory tree whose root is specified with the INPUT tag.
+# subdirectory from a directory tree whose root is specified with the INPUT tag. 
 
-EXCLUDE              =
+EXCLUDE               = 
 
-# If the value of the INPUT tag contains directories, you can use the
-# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude
-# certain files from those directories.
+# If the value of the INPUT tag contains directories, you can use the 
+# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude 
+# certain files from those directories. 
 
-EXCLUDE_PATTERNS     =
+EXCLUDE_PATTERNS      = 
 
 # The EXAMPLE_PATH tag can be used to specify one or more files or 
 # directories that contain example code fragments that are included (see 
-# the \include command).
+# the \include command). 
 
-EXAMPLE_PATH         =
+EXAMPLE_PATH          = 
 
-# If the value of the EXAMPLE_PATH tag contains directories, you can use the
+# If the value of the EXAMPLE_PATH tag contains directories, you can use the 
 # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp 
 # and *.h) to filter out the source-files in the directories. If left 
-# blank all files are included.
+# blank all files are included. 
 
-EXAMPLE_PATTERNS     =
+EXAMPLE_PATTERNS      = 
 
 # The IMAGE_PATH tag can be used to specify one or more files or 
 # directories that contain image that are included in the documentation (see 
-# the \image command).
+# the \image command). 
 
-IMAGE_PATH           =
+IMAGE_PATH            = 
 
-# The INPUT_FILTER tag can be used to specify a program that doxygen should
+# The INPUT_FILTER tag can be used to specify a program that doxygen should 
 # invoke to filter for each input file. Doxygen will invoke the filter program 
-# by executing (via popen()) the command <filter> <input-file>, where <filter>
-# is the value of the INPUT_FILTER tag, and <input-file> is the name of an
-# input file. Doxygen will then use the output that the filter program writes
-# to standard output.
+# by executing (via popen()) the command <filter> <input-file>, where <filter> 
+# is the value of the INPUT_FILTER tag, and <input-file> is the name of an 
+# input file. Doxygen will then use the output that the filter program writes 
+# to standard output. 
 
-INPUT_FILTER         =
+INPUT_FILTER          = 
+
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+
+# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index 
+# of all compounds will be generated. Enable this if the project 
+# contains a lot of classes, structs, unions or interfaces. 
+
+ALPHABETICAL_INDEX    = YES
+
+# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then 
+# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns 
+# in which this list will be split (can be a number in the range [1..20]) 
+
+COLS_IN_ALPHA_INDEX   = 5
+
+# In case all classes in a project start with a common prefix, all 
+# classes will be put under the same header in the alphabetical index. 
+# The IGNORE_PREFIX tag can be used to specify one or more prefixes that 
+# should be ignored while generating the index headers. 
+
+IGNORE_PREFIX         = 
 
 #---------------------------------------------------------------------------
 # configuration options related to the HTML output
 #---------------------------------------------------------------------------
 
-# If the GENERATE_HTML tag is set to YES (the default) Doxygen will
-# generate HTML output
+# If the GENERATE_HTML tag is set to YES (the default) Doxygen will 
+# generate HTML output. 
 
-GENERATE_HTML        = NO
+GENERATE_HTML         = NO
 
-# The HTML_OUTPUT tag is used to specify where the HTML docs will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be
-# put in front of it. If left blank `html' will be used as the default path.
+# The HTML_OUTPUT tag is used to specify where the HTML docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `html' will be used as the default path. 
 
-HTML_OUTPUT          = .
+HTML_OUTPUT           = .
 
 # The HTML_HEADER tag can be used to specify a personal HTML header for 
 # each generated HTML page. If it is left blank doxygen will generate a 
 # standard header.
 
-HTML_HEADER          =
+HTML_HEADER           = 
 
 # The HTML_FOOTER tag can be used to specify a personal HTML footer for 
 # each generated HTML page. If it is left blank doxygen will generate a 
 # standard footer.
 
-HTML_FOOTER          = Doxyfooter
+HTML_FOOTER           = Doxyfooter
 
-# The HTML_STYLESHEET tag can be used to specify a user defined cascading
+# The HTML_STYLESHEET tag can be used to specify a user defined cascading 
 # style sheet that is used by each HTML page. It can be used to 
-# fine-tune the look of the HTML output. If the tag is left blank doxygen
-# will generate a default style sheet
-
-HTML_STYLESHEET      =
-
-# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes,
-# files or namespaces will be aligned in HTML using tables. If set to
-# NO a bullet list will be used.
+# fine-tune the look of the HTML output. If the tag is left blank doxygen 
+# will generate a default style sheet 
 
-HTML_ALIGN_MEMBERS   = YES
+HTML_STYLESHEET       = 
 
-# If the GENERATE_HTMLHELP tag is set to YES, additional index files
-# will be generated that can be used as input for tools like the
-# Microsoft HTML help workshop to generate a compressed HTML help file (.chm)
-# of the generated HTML documentation.
+# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes, 
+# files or namespaces will be aligned in HTML using tables. If set to 
+# NO a bullet list will be used. 
 
-GENERATE_HTMLHELP    = NO
+HTML_ALIGN_MEMBERS    = YES
 
-# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index
-# of all compounds will be generated. Enable this if the project
-# contains a lot of classes, structs, unions or interfaces.
+# If the GENERATE_HTMLHELP tag is set to YES, additional index files 
+# will be generated that can be used as input for tools like the 
+# Microsoft HTML help workshop to generate a compressed HTML help file (.chm) 
+# of the generated HTML documentation. 
 
-ALPHABETICAL_INDEX   = YES
-
-# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then
-# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns
-# in which this list will be split (can be a number in the range [1..20])
-
-COLS_IN_ALPHA_INDEX  = 5
+GENERATE_HTMLHELP     = NO
 
 #---------------------------------------------------------------------------
 # configuration options related to the LaTeX output
 #---------------------------------------------------------------------------
 
-# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
-# generate Latex output.
+# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will 
+# generate Latex output. 
 
-GENERATE_LATEX       = YES
+GENERATE_LATEX        = YES
 
-# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be
-# put in front of it. If left blank `latex' will be used as the default path.
+# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `latex' will be used as the default path. 
 
-LATEX_OUTPUT         = latex
+LATEX_OUTPUT          = latex
 
-# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact
-# LaTeX documents. This may be useful for small projects and may help to
-# save some trees in general.
+# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact 
+# LaTeX documents. This may be useful for small projects and may help to 
+# save some trees in general. 
 
-COMPACT_LATEX        = NO
+COMPACT_LATEX         = NO
 
-# The PAPER_TYPE tag can be used to set the paper type that is used
+# The PAPER_TYPE tag can be used to set the paper type that is used 
 # by the printer. Possible values are: a4, a4wide, letter, legal and 
-# executive. If left blank a4wide will be used.
+# executive. If left blank a4wide will be used. 
 
-PAPER_TYPE           = a4wide
+PAPER_TYPE            = a4wide
 
-# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
-# packages that should be included in the LaTeX output.
+# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX 
+# packages that should be included in the LaTeX output. 
 
-EXTRA_PACKAGES       =
+EXTRA_PACKAGES        = 
 
 # The LATEX_HEADER tag can be used to specify a personal LaTeX header for 
-# the generated latex document. The header should contain everything until
+# the generated latex document. The header should contain everything until 
 # the first chapter. If it is left blank doxygen will generate a 
-# standard header. Notice: only use this tag if you know what you are doing!
+# standard header. Notice: only use this tag if you know what you are doing! 
+
+LATEX_HEADER          = 
+
+# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated 
+# is prepared for conversion to pdf (using ps2pdf). The pdf file will 
+# contain links (just like the HTML output) instead of page references 
+# This makes the output suitable for online browsing using a pdf viewer. 
 
-LATEX_HEADER         =
+PDF_HYPERLINKS        = NO
 
-# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated
-# is prepared for conversion to pdf (using ps2pdf). The pdf file will
-# contain links (just like the HTML output) instead of page references
-# This makes the output suitable for online browsing using a pdf viewer.
+# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \\batchmode. 
+# command to the generated LaTeX files. This will instruct LaTeX to keep 
+# running if errors occur, instead of asking the user for help. 
+# This option is also used when generating formulas in HTML. 
 
-PDF_HYPERLINKS       = NO
+LATEX_BATCHMODE       = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_RTF tag is set to YES Doxygen will generate RTF output 
+# For now this is experimental and is disabled by default. The RTF output 
+# is optimised for Word 97 and may not look too pretty with other readers 
+# or editors.
+
+GENERATE_RTF          = NO
+
+# The RTF_OUTPUT tag is used to specify where the RTF docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `rtf' will be used as the default path. 
+
+RTF_OUTPUT            = rtf
+
+# If the COMPACT_RTF tag is set to YES Doxygen generates more compact 
+# RTF documents. This may be useful for small projects and may help to 
+# save some trees in general. 
+
+COMPACT_RTF           = NO
+
+# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated 
+# will contain hyperlink fields. The RTF file will 
+# contain links (just like the HTML output) instead of page references. 
+# This makes the output suitable for online browsing using a WORD or other. 
+# programs which support those fields. 
+# Note: wordpad (write) and others do not support links. 
+
+RTF_HYPERLINKS        = NO
+
+# Load stylesheet definitions from file. Syntax is similar to doxygen's 
+# config file, i.e. a series of assigments. You only have to provide 
+# replacements, missing definitions are set to their default value. 
+
+RTF_STYLESHEET_FILE   = 
 
 #---------------------------------------------------------------------------
 # configuration options related to the man page output
 #---------------------------------------------------------------------------
 
-# If the GENERATE_MAN tag is set to YES (the default) Doxygen will
-# generate man pages
+# If the GENERATE_MAN tag is set to YES (the default) Doxygen will 
+# generate man pages 
 
-GENERATE_MAN         = NO
+GENERATE_MAN          = NO
 
-# The MAN_OUTPUT tag is used to specify where the man pages will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be
-# put in front of it. If left blank `man' will be used as the default path.
+# The MAN_OUTPUT tag is used to specify where the man pages will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `man' will be used as the default path. 
 
-MAN_OUTPUT           =
+MAN_OUTPUT            = man
 
-# The MAN_EXTENSION tag determines the extension that is added to
-# the generated man pages (default is the subroutine's section .3)
+# The MAN_EXTENSION tag determines the extension that is added to 
+# the generated man pages (default is the subroutine's section .3) 
 
-MAN_EXTENSION        = .3
+MAN_EXTENSION         = .3
 
 #---------------------------------------------------------------------------
-# Configuration options related to the preprocessor 
+# Configuration options related to the preprocessor   
 #---------------------------------------------------------------------------
 
-# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will
-# evaluate all C-preprocessor directives found in the sources and include
-# files.
+# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will 
+# evaluate all C-preprocessor directives found in the sources and include 
+# files. 
 
-ENABLE_PREPROCESSING = YES
+ENABLE_PREPROCESSING  = YES
 
-# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro
+# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro 
 # names in the source code. If set to NO (the default) only conditional 
-# compilation will be performed.
+# compilation will be performed. Macro expansion can be done in a controlled 
+# way by setting EXPAND_ONLY_PREDEF to YES. 
+
+MACRO_EXPANSION       = NO
+
+# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES 
+# then the macro expansion is limited to the macros specified with the 
+# PREDEFINED and EXPAND_AS_PREDEFINED tags. 
 
-MACRO_EXPANSION      = NO
+EXPAND_ONLY_PREDEF    = NO
 
-# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files
-# in the INCLUDE_PATH (see below) will be search if a #include is found.
+# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files 
+# in the INCLUDE_PATH (see below) will be search if a #include is found. 
 
-SEARCH_INCLUDES      = YES
+SEARCH_INCLUDES       = YES
 
-# The INCLUDE_PATH tag can be used to specify one or more directories that
-# contain include files that are not input files but should be processed by
-# the preprocessor.
+# The INCLUDE_PATH tag can be used to specify one or more directories that 
+# contain include files that are not input files but should be processed by 
+# the preprocessor. 
 
-INCLUDE_PATH         =
+INCLUDE_PATH          = 
 
-# The PREDEFINED tag can be used to specify one or more macro names that
-# are defined before the preprocessor is started (similar to the -D option of
-# gcc). The argument of the tag is a list of macros of the form: name
+# You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard 
+# patterns (like *.h and *.hpp) to filter out the header-files in the 
+# directories. If left blank, the patterns specified with FILE_PATTERNS will 
+# be used. 
+
+INCLUDE_FILE_PATTERNS = *.h
+
+# The PREDEFINED tag can be used to specify one or more macro names that 
+# are defined before the preprocessor is started (similar to the -D option of 
+# gcc). The argument of the tag is a list of macros of the form: name 
 # or name=definition (no spaces). If the definition and the = are 
-# omitted =1 is assumed.
+# omitted =1 is assumed. 
 
-PREDEFINED           =
+PREDEFINED            = 
 
-# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES
-# then the macro expansion is limited to the macros specified with the
-# PREDEFINED tag.
+# If the MACRO_EXPANSION and EXPAND_PREDEF_ONLY tags are set to YES then 
+# this tag can be used to specify a list of macro names that should be expanded. 
+# The macro definition that is found in the sources will be used. 
+# Use the PREDEFINED tag if you want to use a different macro definition. 
 
-EXPAND_ONLY_PREDEF   = NO
+EXPAND_AS_DEFINED     = 
 
 #---------------------------------------------------------------------------
-# Configuration options related to external references 
+# Configuration::addtions related to external references   
 #---------------------------------------------------------------------------
 
 # The TAGFILES tag can be used to specify one or more tagfiles. 
 
-TAGFILES             =
+TAGFILES              = 
+
+# When a file name is specified after GENERATE_TAGFILE, doxygen will create 
+# a tag file that is based on the input files it reads. 
+
+GENERATE_TAGFILE      = 
+
+# If the ALLEXTERNALS tag is set to YES all external classes will be listed 
+# in the class index. If set to NO only the inherited external classes 
+# will be listed. 
+
+ALLEXTERNALS          = NO
+
+# The PERL_PATH should be the absolute path and name of the perl script 
+# interpreter (i.e. the result of `which perl'). 
+
+PERL_PATH             = /usr/bin/perl
+
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool   
+#---------------------------------------------------------------------------
+
+# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is 
+# available from the path. This tool is part of Graphviz, a graph visualization 
+# toolkit from AT&T and Lucent Bell Labs. The other options in this section 
+# have no effect if this option is set to NO (the default) 
+
+HAVE_DOT              = YES
+
+# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for each documented class showing the direct and 
+# indirect inheritance relations. Setting this tag to YES will force the 
+# the CLASS_DIAGRAMS tag to NO.
+
+CLASS_GRAPH           = YES
+
+# If the COLLABORATION_GRAPH and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for each documented class showing the direct and 
+# indirect implementation dependencies (inheritance, containment, and 
+# class references variables) of the class with other documented classes. 
+
+COLLABORATION_GRAPH   = YES
+
+# If the ENABLE_PREPROCESSING, INCLUDE_GRAPH, and HAVE_DOT tags are set to 
+# YES then doxygen will generate a graph for each documented file showing 
+# the direct and indirect include dependencies of the file with other 
+# documented files. 
+
+INCLUDE_GRAPH         = YES
+
+# If the ENABLE_PREPROCESSING, INCLUDED_BY_GRAPH, and HAVE_DOT tags are set to 
+# YES then doxygen will generate a graph for each documented header file showing 
+# the documented files that directly or indirectly include this file 
+
+INCLUDED_BY_GRAPH     = YES
+
+# If the GRAPHICAL_HIERARCHY and HAVE_DOT tags are set to YES then doxygen 
+# will graphical hierarchy of all classes instead of a textual one. 
+
+GRAPHICAL_HIERARCHY   = YES
 
-# When a file name is specified after GENERATE_TAGFILE, doxygen will create
-# a tag file that is based on the input files it reads.
+# The tag DOT_PATH can be used to specify the path where the dot tool can be 
+# found. If left blank, it is assumed the dot tool can be found on the path. 
 
-GENERATE_TAGFILE     =
+DOT_PATH              = 
 
-# If the ALLEXTERNALS tag is set to YES all external classes will be listed
-# in the class index. If set to NO only the inherited external classes
-# will be listed.
+# The MAX_DOT_GRAPH_WIDTH tag can be used to set the maximum allowed width 
+# (in pixels) of the graphs generated by dot. If a graph becomes larger than 
+# this value, doxygen will try to truncate the graph, so that it fits within 
+# the specified constraint. Beware that most browsers cannot cope with very 
+# large images. 
 
-ALLEXTERNALS         = NO
+MAX_DOT_GRAPH_WIDTH   = 1024
 
-# The PERL_PATH should be the absolute path and name of the perl script
-# interpreter (i.e. the result of `which perl').
+# The MAX_DOT_GRAPH_HEIGHT tag can be used to set the maximum allows height 
+# (in pixels) of the graphs generated by dot. If a graph becomes larger than 
+# this value, doxygen will try to truncate the graph, so that it fits within 
+# the specified constraint. Beware that most browsers cannot cope with very 
+# large images. 
 
-PERL_PATH            = /usr/bin/perl
+MAX_DOT_GRAPH_HEIGHT  = 1024
 
 #---------------------------------------------------------------------------
-# Configuration options related to the search engine 
+# Configuration::addtions related to the search engine   
 #---------------------------------------------------------------------------
 
 # The SEARCHENGINE tag specifies whether or not a search engine should be 
-# used. If set to NO the values of all tags below this one will be ignored.
+# used. If set to NO the values of all tags below this one will be ignored. 
 
-SEARCHENGINE         = NO
+SEARCHENGINE          = NO
 
-# The CGI_NAME tag should be the name of the CGI script that
-# starts the search engine (doxysearch) with the correct parameters.
-# A script with this name will be generated by doxygen.
+# The CGI_NAME tag should be the name of the CGI script that 
+# starts the search engine (doxysearch) with the correct parameters. 
+# A script with this name will be generated by doxygen. 
 
-CGI_NAME             = search.cgi
+CGI_NAME              = search.cgi
 
-# The CGI_URL tag should be the absolute URL to the directory where the
+# The CGI_URL tag should be the absolute URL to the directory where the 
 # cgi binaries are located. See the documentation of your http daemon for 
-# details.
+# details. 
 
-CGI_URL              =
+CGI_URL               = 
 
-# The DOC_URL tag should be the absolute URL to the directory where the
+# The DOC_URL tag should be the absolute URL to the directory where the 
 # documentation is located. If left blank the absolute path to the 
-# documentation, with file:// prepended to it, will be used.
+# documentation, with file:// prepended to it, will be used. 
 
-DOC_URL              =
+DOC_URL               = 
 
-# The DOC_ABSPATH tag should be the absolute path to the directory where the
-# documentation is located. If left blank the directory on the local machine
-# will be used.
+# The DOC_ABSPATH tag should be the absolute path to the directory where the 
+# documentation is located. If left blank the directory on the local machine 
+# will be used. 
 
-DOC_ABSPATH          =
+DOC_ABSPATH           = 
 
-# The BIN_ABSPATH tag must point to the directory where the doxysearch binary
-# is installed.
+# The BIN_ABSPATH tag must point to the directory where the doxysearch binary 
+# is installed. 
 
-BIN_ABSPATH          = /usr/local/bin/
+BIN_ABSPATH           = /usr/local/bin/
 
 # The EXT_DOC_PATHS tag can be used to specify one or more paths to 
-# documentation generated for other projects. This allows doxysearch to search
-# the documentation for these projects as well.
+# documentation generated for other projects. This allows doxysearch to search 
+# the documentation for these projects as well. 
 
-EXT_DOC_PATHS        =
+EXT_DOC_PATHS         = 
diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 28615659..97033860 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -2411,7 +2411,7 @@ enough to know how to differentiate.  But if the function you want to
 implement does have a pole somewhere in the complex plane, you need to
 write another method for Laurent expansion around that point.
 
-Now that all the ingrediences for @code{cos} have been set up, we need
+Now that all the ingredients for @code{cos} have been set up, we need
 to tell the system about it.  This is done by a macro and we are not
 going to descibe how it expands, please consult your preprocessor if you
 are curious:
diff --git a/doc/tutorial/stamp-vti b/doc/tutorial/stamp-vti
index 5458d424..74c77b59 100644
--- a/doc/tutorial/stamp-vti
+++ b/doc/tutorial/stamp-vti
@@ -1,3 +1,3 @@
-@set UPDATED 25 July 2000
-@set EDITION 0.6.3
-@set VERSION 0.6.3
+@set UPDATED 10 August 2000
+@set EDITION 0.6.4
+@set VERSION 0.6.4
diff --git a/doc/tutorial/version.texi b/doc/tutorial/version.texi
index 5458d424..74c77b59 100644
--- a/doc/tutorial/version.texi
+++ b/doc/tutorial/version.texi
@@ -1,3 +1,3 @@
-@set UPDATED 25 July 2000
-@set EDITION 0.6.3
-@set VERSION 0.6.3
+@set UPDATED 10 August 2000
+@set EDITION 0.6.4
+@set VERSION 0.6.4
diff --git a/ginac/archive.cpp b/ginac/archive.cpp
index 73df064a..d249f297 100644
--- a/ginac/archive.cpp
+++ b/ginac/archive.cpp
@@ -580,8 +580,8 @@ void archive_node::printraw(std::ostream &os) const
 	}
 }
 
-/** Create a dummy archive.  The intention is to fill archive_node's default ctor,
- *  which is currently a Cint-requirement. */
+/** Create a dummy archive.  The intention is to fill archive_node's default
+ *  ctor, which is currently a Cint-requirement. */
 archive* archive_node::dummy_ar_creator(void)
 {
     static archive* some_ar = new archive;
diff --git a/ginac/flags.h b/ginac/flags.h
index 41a890b1..24edd02a 100644
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -43,11 +43,21 @@ class determinant_algo {
 public:
     enum { automatic,
            gauss,
+           divfree,
            laplace,
            bareiss
          };
 };
 
+class solve_algo {
+public:
+    enum { automatic,
+           gauss,
+           divfree,
+           bareiss
+         };
+};
+
 class status_flags {
 public:
     enum { dynallocated              = 0x0001,
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 855c7317..677d92f6 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -433,7 +433,7 @@ ex lsolve(const ex &eqns, const ex &symbols)
     // solve a system of linear equations
     if (eqns.info(info_flags::relation_equal)) {
         if (!symbols.info(info_flags::symbol))
-            throw(std::invalid_argument("lsolve: 2nd argument must be a symbol"));
+            throw(std::invalid_argument("lsolve(): 2nd argument must be a symbol"));
         ex sol=lsolve(lst(eqns),lst(symbols));
         
         GINAC_ASSERT(sol.nops()==1);
@@ -444,19 +444,19 @@ ex lsolve(const ex &eqns, const ex &symbols)
     
     // syntax checks
     if (!eqns.info(info_flags::list)) {
-        throw(std::invalid_argument("lsolve: 1st argument must be a list"));
+        throw(std::invalid_argument("lsolve(): 1st argument must be a list"));
     }
     for (unsigned i=0; i<eqns.nops(); i++) {
         if (!eqns.op(i).info(info_flags::relation_equal)) {
-            throw(std::invalid_argument("lsolve: 1st argument must be a list of equations"));
+            throw(std::invalid_argument("lsolve(): 1st argument must be a list of equations"));
         }
     }
     if (!symbols.info(info_flags::list)) {
-        throw(std::invalid_argument("lsolve: 2nd argument must be a list"));
+        throw(std::invalid_argument("lsolve(): 2nd argument must be a list"));
     }
     for (unsigned i=0; i<symbols.nops(); i++) {
         if (!symbols.op(i).info(info_flags::symbol)) {
-            throw(std::invalid_argument("lsolve: 2nd argument must be a list of symbols"));
+            throw(std::invalid_argument("lsolve(): 2nd argument must be a list of symbols"));
         }
     }
     
@@ -473,7 +473,7 @@ ex lsolve(const ex &eqns, const ex &symbols)
             linpart -= co*symbols.op(c);
             sys.set(r,c,co);
         }
-        linpart=linpart.expand();
+        linpart = linpart.expand();
         rhs.set(r,0,-linpart);
     }
     
@@ -486,32 +486,21 @@ ex lsolve(const ex &eqns, const ex &symbols)
             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.solve(vars,rhs);
     } catch (const runtime_error & e) {
-        // probably singular matrix (or other error)
-        // return empty solution list
-        // cerr << e.what() << endl;
+        // Probably singular matrix or otherwise overdetermined system:
+        // It is consistent to return an empty list
         return lst();
-    }
-    
-    // return a list of equations
-    if (solution.cols()!=1) {
-        throw(std::runtime_error("lsolve: strange number of columns returned from matrix::solve"));
-    }
-    if (solution.rows()!=symbols.nops()) {
-        cout << "symbols.nops()=" << symbols.nops() << endl;
-        cout << "solution.rows()=" << solution.rows() << endl;
-        throw(std::runtime_error("lsolve: strange number of rows returned from matrix::solve"));
-    }
+    }    
+    GINAC_ASSERT(solution.cols()==1);
+    GINAC_ASSERT(solution.rows()==symbols.nops());
     
-    // return list of the form lst(var1==sol1,var2==sol2,...)
+    // return list of equations of the form lst(var1==sol1,var2==sol2,...)
     lst sollist;
-    for (unsigned i=0; i<symbols.nops(); i++) {
+    for (unsigned i=0; i<symbols.nops(); i++)
         sollist.append(symbols.op(i)==solution(i,0));
-    }
     
     return sollist;
 }
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index bf490f04..6240ce20 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -108,7 +108,7 @@ matrix::matrix(unsigned r, unsigned c)
     m.resize(r*c, _ex0());
 }
 
- // protected
+// protected
 
 /** Ctor from representation, for internal use only. */
 matrix::matrix(unsigned r, unsigned c, const exvector & m2)
@@ -369,17 +369,17 @@ matrix matrix::sub(const matrix & other) const
  *  @exception logic_error (incompatible matrices) */
 matrix matrix::mul(const matrix & other) const
 {
-    if (col != other.row)
+    if (this->cols() != other.rows())
         throw (std::logic_error("matrix::mul(): incompatible matrices"));
     
-    exvector prod(row*other.col);
+    exvector prod(this->rows()*other.cols());
     
-    for (unsigned r1=0; r1<rows(); ++r1) {
-        for (unsigned c=0; c<cols(); ++c) {
+    for (unsigned r1=0; r1<this->rows(); ++r1) {
+        for (unsigned c=0; c<this->cols(); ++c) {
             if (m[r1*col+c].is_zero())
                 continue;
-            for (unsigned r2=0; r2<other.col; ++r2)
-                prod[r1*other.col+r2] += m[r1*col+c] * other.m[c*other.col+r2];
+            for (unsigned r2=0; r2<other.cols(); ++r2)
+                prod[r1*other.col+r2] += (m[r1*col+c] * other.m[c*other.col+r2]).expand();
         }
     }
     return matrix(row, other.col, prod);
@@ -418,13 +418,13 @@ matrix & matrix::set(unsigned ro, unsigned co, ex value)
  *  represents the transposed. */
 matrix matrix::transpose(void) const
 {
-    exvector trans(col*row);
+    exvector trans(this->cols()*this->rows());
     
-    for (unsigned r=0; r<col; ++r)
-        for (unsigned c=0; c<row; ++c)
-            trans[r*row+c] = m[c*col+r];
+    for (unsigned r=0; r<this->cols(); ++r)
+        for (unsigned c=0; c<this->rows(); ++c)
+            trans[r*this->rows()+c] = m[c*this->cols()+r];
     
-    return matrix(col,row,trans);
+    return matrix(this->cols(),this->rows(),trans);
 }
 
 
@@ -447,9 +447,7 @@ 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 statistical information about this matrix:
     bool numeric_flag = true;
     bool normal_flag = false;
@@ -468,11 +466,11 @@ ex matrix::determinant(unsigned algo) const
     
     // Here is the heuristics in case this routine has to decide:
     if (algo == determinant_algo::automatic) {
-        // Minor expansion is generally a good starting point:
+        // Minor expansion is generally a good guess:
         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)
+        if (row>3 && 5*sparse_count<=row*col)
             algo = determinant_algo::bareiss;
         // Purely numeric matrix can be handled by Gauss elimination.
         // This overrides any prior decisions.
@@ -480,17 +478,27 @@ ex matrix::determinant(unsigned algo) const
             algo = determinant_algo::gauss;
     }
     
+    // Trap the trivial case here, since some algorithms don't like it
+    if (this->row==1) {
+        // for consistency with non-trivial determinants...
+        if (normal_flag)
+            return m[0].normal();
+        else
+            return m[0].expand();
+    }
+    
+    // Compute the determinant
     switch(algo) {
         case determinant_algo::gauss: {
             ex det = 1;
             matrix tmp(*this);
-            int sign = tmp.gauss_elimination();
+            int sign = tmp.gauss_elimination(true);
             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();
+                return (sign*det).normal().expand();
         }
         case determinant_algo::bareiss: {
             matrix tmp(*this);
@@ -501,6 +509,19 @@ ex matrix::determinant(unsigned algo) const
             else
                 return (sign*tmp.m[row*col-1]).expand();
         }
+        case determinant_algo::divfree: {
+            matrix tmp(*this);
+            int sign;
+            sign = tmp.division_free_elimination(true);
+            if (sign==0)
+                return _ex0();
+            ex det = tmp.m[row*col-1];
+            // factor out accumulated bogus slag
+            for (unsigned d=0; d<row-2; ++d)
+                for (unsigned j=0; j<row-d-2; ++j)
+                    det = (det/tmp.m[d*col+d]).normal();
+            return (sign*det);
+        }
         case determinant_algo::laplace:
         default: {
             // This is the minor expansion scheme.  We always develop such
@@ -533,8 +554,9 @@ ex matrix::determinant(unsigned algo) const
             }
             
             if (normal_flag)
-                return sign*matrix(row,col,result).determinant_minor().normal();
-            return sign*matrix(row,col,result).determinant_minor();
+                return (sign*matrix(row,col,result).determinant_minor()).normal();
+            else
+                return sign*matrix(row,col,result).determinant_minor();
         }
     }
 }
@@ -550,7 +572,6 @@ ex matrix::trace(void) const
 {
     if (row != col)
         throw (std::logic_error("matrix::trace(): matrix not square"));
-    GINAC_ASSERT(row*col==m.capacity());
     
     ex tr;
     for (unsigned r=0; r<col; ++r)
@@ -625,15 +646,19 @@ matrix matrix::inverse(void) const
     if (row != col)
         throw (std::logic_error("matrix::inverse(): matrix not square"));
     
+    // NOTE: the Gauss-Jordan elimination used here can in principle be
+    // replaced this by two clever calls to gauss_elimination() and some to
+    // transpose().  Wouldn't be more efficient (maybe less?), just more
+    // orthogonal.
     matrix tmp(row,col);
     // set tmp to the unit matrix
     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) {
-        int indx = cpy.pivot(r1);
+        int indx = cpy.pivot(r1, r1);
         if (indx == -1) {
             throw (std::runtime_error("matrix::inverse(): singular matrix"));
         }
@@ -651,202 +676,121 @@ matrix matrix::inverse(void) const
                 ex a2 = cpy.m[r2*col+r1];
                 for (unsigned c=0; c<col; ++c) {
                     cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
+                    if (!cpy.m[r2*col+c].info(info_flags::numeric))
+                        cpy.m[r2*col+c] = cpy.m[r2*col+c].normal();
                     tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
+                    if (!tmp.m[r2*col+c].info(info_flags::numeric))
+                        tmp.m[r2*col+c] = tmp.m[r2*col+c].normal();
                 }
             }
         }
     }
+    
     return 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'
- *  by Keith O. Geddes et al.
+/** Solve a linear system consisting of a m x n matrix and a m x p right hand
+ *  side by applying an elimination scheme to the augmented matrix.
  *
- *  @param vars n x p matrix
+ *  @param vars n x p matrix, all elements must be symbols 
  *  @param rhs m x p matrix
+ *  @return n x p solution matrix
  *  @exception logic_error (incompatible matrices)
- *  @exception runtime_error (singular matrix) */
-matrix matrix::fraction_free_elim(const matrix & vars,
-                                  const matrix & rhs) const
+ *  @exception invalid_argument (1st argument must be matrix of symbols)
+ *  @exception runtime_error (inconsistent linear system)
+ *  @see       solve_algo */
+matrix matrix::solve(const matrix & vars,
+                     const matrix & rhs,
+                     unsigned algo) const
 {
-    // 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"));
-    
-    matrix a(*this);  // make a copy of the matrix
-    matrix b(rhs);    // make a copy of the rhs vector
+    const unsigned m = this->rows();
+    const unsigned n = this->cols();
+    const unsigned p = rhs.cols();
     
-    // 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 = 0;
+    // syntax checks    
+    if ((rhs.rows() != m) || (vars.rows() != n) || (vars.col != p))
+        throw (std::logic_error("matrix::solve(): incompatible matrices"));
+    for (unsigned ro=0; ro<n; ++ro)
+        for (unsigned co=0; co<p; ++co)
+            if (!vars(ro,co).info(info_flags::symbol))
+                throw (std::invalid_argument("matrix::solve(): 1st argument must be matrix of symbols"));
     
-    // eliminate below row r, with pivot in column k
-    for (unsigned k=0; (k<n)&&(r<m); ++k) {
-        // find a nonzero pivot
-        unsigned 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.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.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.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.m[r*a.cols()+k];
-            ++r;
-        }
+    // build the augmented matrix of *this with rhs attached to the right
+    matrix aug(m,n+p);
+    for (unsigned r=0; r<m; ++r) {
+        for (unsigned c=0; c<n; ++c)
+            aug.m[r*(n+p)+c] = this->m[r*n+c];
+        for (unsigned c=0; c<p; ++c)
+            aug.m[r*(n+p)+c+n] = rhs.m[r*p+c];
     }
     
-#ifdef DO_GINAC_ASSERT
-    // test if we really have an upper echelon matrix
-    int zero_in_last_row = -1;
-    for (unsigned r=0; r<m; ++r) {
-        int zero_in_this_row=0;
-        for (unsigned c=0; c<n; ++c) {
-            if (a.m[r*a.cols()+c].is_zero())
-               ++zero_in_this_row;
-            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;
+    // Gather some statistical information about the augmented matrix:
+    bool numeric_flag = true;
+    for (exvector::const_iterator r=aug.m.begin(); r!=aug.m.end(); ++r) {
+        if (!(*r).info(info_flags::numeric))
+            numeric_flag = false;
     }
-#endif // def DO_GINAC_ASSERT
     
-    // assemble solution
-    matrix sol(n,1);
-    unsigned last_assigned_sol = n+1;
-    for (int r=m-1; r>=0; --r) {
-        unsigned first_non_zero = 1;
-        while ((first_non_zero<=n)&&(a(r,first_non_zero-1).is_zero()))
-            first_non_zero++;
-        if (first_non_zero>n) {
-            // row consists only of zeroes, corresponding rhs must be 0 as well
-            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.m[c*vars.cols()]);
-            ex e = b.m[r*b.cols()];
-            for (unsigned c=first_non_zero; c<n; ++c)
-                e -= a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
-            sol.set(first_non_zero-1,0,
-                    (e/(a.m[r*a.cols()+(first_non_zero-1)])).normal());
-            last_assigned_sol = first_non_zero;
-        }
+    // Here is the heuristics in case this routine has to decide:
+    if (algo == solve_algo::automatic) {
+        // Bareiss (fraction-free) elimination is generally a good guess:
+        algo = solve_algo::bareiss;
+        // For m<3, Bareiss elimination is equivalent to division free
+        // elimination but has more logistic overhead
+        if (m<3)
+            algo = solve_algo::divfree;
+        // This overrides any prior decisions.
+        if (numeric_flag)
+            algo = solve_algo::gauss;
     }
-    // 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.m[c*vars.cols()]);
     
-#ifdef DO_GINAC_ASSERT
-    // test solution with echelon matrix
-    for (unsigned r=0; r<m; ++r) {
-        ex e = 0;
-        for (unsigned c=0; c<n; ++c)
-            e += a(r,c)*sol(c,0);
-        if (!(e-b(r,0)).normal().is_zero()) {
-            cout << "e=" << e;
-            cout << "b(" << r <<",0)=" << b(r,0) << endl;
-            cout << "diff=" << (e-b(r,0)).normal() << endl;
-        }
-        GINAC_ASSERT((e-b(r,0)).normal().is_zero());
+    // Eliminate the augmented matrix:
+    switch(algo) {
+        case solve_algo::gauss:
+            aug.gauss_elimination();
+        case solve_algo::divfree:
+            aug.division_free_elimination();
+        case solve_algo::bareiss:
+        default:
+            aug.fraction_free_elimination();
     }
     
-    // test solution with original matrix
-    for (unsigned r=0; r<m; ++r) {
-        ex e = 0;
-        for (unsigned c=0; c<n; ++c)
-            e += this->m[r*cols()+c]*sol(c,0);
-        try {
-            if (!(e-rhs(r,0)).normal().is_zero()) {
-                cout << "e==" << e << endl;
-                e.printtree(cout);
-                ex en = e.normal();
-                cout << "e.normal()=" << en << endl;
-                en.printtree(cout);
-                cout << "rhs(" << r <<",0)=" << rhs(r,0) << endl;
-                cout << "diff=" << (e-rhs(r,0)).normal() << endl;
+    // assemble the solution matrix:
+    matrix sol(n,p);
+    for (unsigned co=0; co<p; ++co) {
+        unsigned last_assigned_sol = n+1;
+        for (int r=m-1; r>=0; --r) {
+            unsigned fnz = 1;    // first non-zero in row
+            while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].is_zero()))
+                ++fnz;
+            if (fnz>n) {
+                // row consists only of zeros, corresponding rhs must be 0, too
+                if (!aug.m[r*(n+p)+n+co].is_zero()) {
+                    throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
+                }
+            } else {
+                // assign solutions for vars between fnz+1 and
+                // last_assigned_sol-1: free parameters
+                for (unsigned c=fnz; c<last_assigned_sol-1; ++c)
+                    sol.set(c,co,vars.m[c*p+co]);
+                ex e = aug.m[r*(n+p)+n+co];
+                for (unsigned c=fnz; c<n; ++c)
+                    e -= aug.m[r*(n+p)+c]*sol.m[c*p+co];
+                sol.set(fnz-1,co,
+                        (e/(aug.m[r*(n+p)+(fnz-1)])).normal());
+                last_assigned_sol = fnz;
             }
-        } catch (...) {
-            ex xxx = e - rhs(r,0);
-            cerr << "xxx=" << xxx << endl << endl;
         }
-        GINAC_ASSERT((e-rhs(r,0)).normal().is_zero());
+        // assign solutions for vars between 1 and
+        // last_assigned_sol-1: free parameters
+        for (unsigned ro=0; ro<last_assigned_sol-1; ++ro)
+            sol.set(ro,co,vars(ro,co));
     }
-#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"));
-    
-    throw (std::runtime_error("FIXME: need implementation."));
-}
-
-/** Old and obsolete interface: */
-matrix matrix::old_solve(const matrix & v) const
-{
-    if ((v.row != col) || (col != v.row))
-        throw (std::logic_error("matrix::solve(): incompatible matrices"));
-    
-    // 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] = 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
 
@@ -863,11 +807,12 @@ matrix matrix::old_solve(const matrix & v) const
 ex matrix::determinant_minor(void) const
 {
     // for small matrices the algorithm does not make any sense:
-    if (this->row==1)
-        return m[0];
-    if (this->row==2)
+    const unsigned n = this->cols();
+    if (n==1)
+        return m[0].expand();
+    if (n==2)
         return (m[0]*m[3]-m[2]*m[1]).expand();
-    if (this->row==3)
+    if (n==3)
         return (m[0]*m[4]*m[8]-m[0]*m[5]*m[7]-
                 m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
                 m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
@@ -875,8 +820,8 @@ ex matrix::determinant_minor(void) const
     // 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) {
+    // matrix minorM(this->rows()-1,this->cols()-1);
+    // for (unsigned r1=0; r1<this->rows(); ++r1) {
     //     // shortcut if element(r1,0) vanishes
     //     if (m[r1*col].is_zero())
     //         continue;
@@ -906,10 +851,10 @@ ex matrix::determinant_minor(void) const
     
     // Unique flipper counter for partitioning into minors
     std::vector<unsigned> Pkey;
-    Pkey.reserve(this->col);
+    Pkey.reserve(n);
     // key for minor determinant (a subpartition of Pkey)
     std::vector<unsigned> Mkey;
-    Mkey.reserve(this->col-1);
+    Mkey.reserve(n-1);
     // we store our subminors in maps, keys being the rows they arise from
     typedef std::map<std::vector<unsigned>,class ex> Rmap;
     typedef std::map<std::vector<unsigned>,class ex>::value_type Rmap_value;
@@ -917,34 +862,34 @@ ex matrix::determinant_minor(void) const
     Rmap B;
     ex det;
     // initialize A with last column:
-    for (unsigned r=0; r<this->col; ++r) {
+    for (unsigned r=0; r<n; ++r) {
         Pkey.erase(Pkey.begin(),Pkey.end());
         Pkey.push_back(r);
-        A.insert(Rmap_value(Pkey,m[this->col*r+this->col-1]));
+        A.insert(Rmap_value(Pkey,m[n*(r+1)-1]));
     }
     // proceed from right to left through matrix
-    for (int c=this->col-2; c>=0; --c) {
+    for (int c=n-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)
+        for (unsigned i=0; i<n-c; ++i)
             Pkey.push_back(i);
         unsigned fc = 0;  // controls logic for our strange flipper counter
         do {
             det = _ex0();
-            for (unsigned r=0; r<this->col-c; ++r) {
+            for (unsigned r=0; r<n-c; ++r) {
                 // maybe there is nothing to do?
-                if (m[Pkey[r]*this->col+c].is_zero())
+                if (m[Pkey[r]*n+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)
+                for (unsigned i=0; i<n-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];
+                    det -= m[Pkey[r]*n+c]*A[Mkey];
                 else
-                    det += m[Pkey[r]*this->col+c]*A[Mkey];
+                    det += m[Pkey[r]*n+c]*A[Mkey];
             }
             // prevent build-up of deep nesting of expressions saves time:
             det = det.expand();
@@ -952,13 +897,13 @@ ex matrix::determinant_minor(void) const
             if (!det.is_zero())
                 B.insert(Rmap_value(Pkey,det));
             // increment our strange flipper counter
-            for (fc=this->col-c; fc>0; --fc) {
+            for (fc=n-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)
+            if (fc<n-c)
+                for (unsigned j=fc; j<n-c; ++j)
                     Pkey[j] = Pkey[j-1]+1;
         } while(fc);
         // next column, so change the role of A and B:
@@ -970,28 +915,51 @@ ex matrix::determinant_minor(void) const
 }
 
 
-/** Perform the steps of an ordinary Gaussian elimination to bring the matrix
- *  into an upper echelon form.
+/** Perform the steps of an ordinary Gaussian elimination to bring the m x n
+ *  matrix into an upper echelon form.  The algorithm is ok for matrices
+ *  with numeric coefficients but quite unsuited for symbolic matrices.
  *
+ *  @param det may be set to true to save a lot of space if one is only
+ *  interested in the diagonal elements (i.e. for calculating determinants).
+ *  The others are set to zero in this case.
  *  @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 matrix::gauss_elimination(const bool det)
 {
     ensure_if_modifiable();
+    const unsigned m = this->rows();
+    const unsigned n = this->cols();
+    GINAC_ASSERT(!det || n==m);
     int sign = 1;
-    ex piv;
-    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) {
-            piv = this->m[r2*col+r1] / this->m[r1*col+r1];
-            for (unsigned c=r1+1; c<col; ++c)
-                this->m[r2*col+c] -= piv * this->m[r1*col+c];
-            for (unsigned c=0; c<=r1; ++c)
-                this->m[r2*col+c] = _ex0();
+    
+    unsigned r0 = 0;
+    for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+        int indx = pivot(r0, r1, true);
+        if (indx == -1) {
+            sign = 0;
+            if (det)
+                return 0;  // leaves *this in a messy state
+        }
+        if (indx>=0) {
+            if (indx > 0)
+                sign = -sign;
+            for (unsigned r2=r0+1; r2<m; ++r2) {
+                ex piv = this->m[r2*n+r1] / this->m[r0*n+r1];
+                for (unsigned c=r1+1; c<n; ++c) {
+                    this->m[r2*n+c] -= piv * this->m[r0*n+c];
+                    if (!this->m[r2*n+c].info(info_flags::numeric))
+                        this->m[r2*n+c] = this->m[r2*n+c].normal();
+                }
+                // fill up left hand side with zeros
+                for (unsigned c=0; c<=r1; ++c)
+                    this->m[r2*n+c] = _ex0();
+            }
+            if (det) {
+                // save space by deleting no longer needed elements
+                for (unsigned c=r0+1; c<n; ++c)
+                    this->m[r0*n+c] = _ex0();
+            }
+            ++r0;
         }
     }
     
@@ -999,26 +967,46 @@ int matrix::gauss_elimination(void)
 }
 
 
-/** Perform the steps of division free elimination to bring the matrix
+/** Perform the steps of division free elimination to bring the m x n matrix
  *  into an upper echelon form.
  *
+ *  @param det may be set to true to save a lot of space if one is only
+ *  interested in the diagonal elements (i.e. for calculating determinants).
+ *  The others are set to zero in this case.
  *  @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::division_free_elimination(void)
+int matrix::division_free_elimination(const bool det)
 {
-    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[r1*col+r1]*this->m[r2*col+c] - this->m[r2*col+r1]*this->m[r1*col+c];
-            for (unsigned c=0; c<=r1; ++c)
-                this->m[r2*col+c] = _ex0();
+    const unsigned m = this->rows();
+    const unsigned n = this->cols();
+    GINAC_ASSERT(!det || n==m);
+    int sign = 1;
+    
+    unsigned r0 = 0;
+    for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+        int indx = pivot(r0, r1, true);
+        if (indx==-1) {
+            sign = 0;
+            if (det)
+                return 0;  // leaves *this in a messy state
+        }
+        if (indx>=0) {
+            if (indx>0)
+                sign = -sign;
+            for (unsigned r2=r0+1; r2<m; ++r2) {
+                for (unsigned c=r1+1; c<n; ++c)
+                    this->m[r2*n+c] = (this->m[r0*n+r1]*this->m[r2*n+c] - this->m[r2*n+r1]*this->m[r0*n+c]).expand();
+                // fill up left hand side with zeros
+                for (unsigned c=0; c<=r1; ++c)
+                    this->m[r2*n+c] = _ex0();
+            }
+            if (det) {
+                // save space by deleting no longer needed elements
+                for (unsigned c=r0+1; c<n; ++c)
+                    this->m[r0*n+c] = _ex0();
+            }
+            ++r0;
         }
     }
     
@@ -1032,11 +1020,11 @@ int matrix::division_free_elimination(void)
  *  is possible, since we know the divisor at each step.
  *  
  *  @param det may be set to true to save a lot of space if one is only
- *  interested in the last element (i.e. for calculating determinants), the
+ *  interested in the last element (i.e. for calculating determinants). The
  *  others are set to zero in this case.
  *  @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::fraction_free_elimination(bool det)
+int matrix::fraction_free_elimination(const bool det)
 {
     // Method:
     // (single-step fraction free elimination scheme, already known to Jordan)
@@ -1062,12 +1050,13 @@ int matrix::fraction_free_elimination(bool det)
     // and D{m[k+1](r,c)} by
     //     D{m[k-1](k-1,k-1)}.
     
-    GINAC_ASSERT(!det || row==col);
     ensure_if_modifiable();
-    if (rows()==1)
-        return 1;
-    
+    const unsigned m = this->rows();
+    const unsigned n = this->cols();
+    GINAC_ASSERT(!det || n==m);
     int sign = 1;
+    if (m==1)
+        return 1;
     ex divisor_n = 1;
     ex divisor_d = 1;
     ex dividend_n;
@@ -1081,62 +1070,70 @@ int matrix::fraction_free_elimination(bool det)
     // need GCDs) since the elements of *this might be unnormalized, which
     // makes things more complicated than they need to be.
     matrix tmp_n(*this);
-    matrix tmp_d(row,col);  // for denominators, if needed
+    matrix tmp_d(m,n);  // for denominators, if needed
     lst srl;  // symbol replacement list
-    exvector::iterator it = m.begin();
+    exvector::iterator it = this->m.begin();
     exvector::iterator tmp_n_it = tmp_n.m.begin();
     exvector::iterator tmp_d_it = tmp_d.m.begin();
-    for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it) {
+    for (; it!= this->m.end(); ++it, ++tmp_n_it, ++tmp_d_it) {
         (*tmp_n_it) = (*it).normal().to_rational(srl);
         (*tmp_d_it) = (*tmp_n_it).denom();
         (*tmp_n_it) = (*tmp_n_it).numer();
     }
     
-    for (unsigned r1=0; r1<row-1; ++r1) {
-        int indx = tmp_n.pivot(r1);
-        if (det && indx==-1)
-            return 0;  // FIXME: what to do if det is false, some day?
-        if (indx>0) {
-            sign = -sign;
-            // rows r1 and indx were swapped, so pivot matrix tmp_d:
-            for (unsigned c=0; c<col; ++c)
-                tmp_d.m[row*indx+c].swap(tmp_d.m[row*r1+c]);
+    unsigned r0 = 0;
+    for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+        int indx = tmp_n.pivot(r0, r1, true);
+        if (indx==-1) {
+            sign = 0;
+            if (det)
+                return 0;
         }
-        if (r1>0) {
-            divisor_n = tmp_n.m[(r1-1)*col+(r1-1)].expand();
-            divisor_d = tmp_d.m[(r1-1)*col+(r1-1)].expand();
-            // save space by deleting no longer needed elements:
-            if (det) {
-                for (unsigned c=0; c<col; ++c) {
-                    tmp_n.m[(r1-1)*col+c] = 0;
-                    tmp_d.m[(r1-1)*col+c] = 1;
+        if (indx>=0) {
+            if (indx>0) {
+                sign = -sign;
+                // tmp_n's rows r0 and indx were swapped, do the same in tmp_d:
+                for (unsigned c=r1; c<n; ++c)
+                    tmp_d.m[n*indx+c].swap(tmp_d.m[n*r0+c]);
+            }
+            for (unsigned r2=r0+1; r2<m; ++r2) {
+                for (unsigned c=r1+1; c<n; ++c) {
+                    dividend_n = (tmp_n.m[r0*n+r1]*tmp_n.m[r2*n+c]*
+                                  tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]
+                                 -tmp_n.m[r2*n+r1]*tmp_n.m[r0*n+c]*
+                                  tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+                    dividend_d = (tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]*
+                                  tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+                    bool check = divide(dividend_n, divisor_n,
+                                        tmp_n.m[r2*n+c], true);
+                    check &= divide(dividend_d, divisor_d,
+                                    tmp_d.m[r2*n+c], true);
+                    GINAC_ASSERT(check);
                 }
+                // fill up left hand side with zeros
+                for (unsigned c=0; c<=r1; ++c)
+                    tmp_n.m[r2*n+c] = _ex0();
             }
-        }
-        for (unsigned r2=r1+1; r2<row; ++r2) {
-            for (unsigned c=r1+1; c<col; ++c) {
-                dividend_n = (tmp_n.m[r1*col+r1]*tmp_n.m[r2*col+c]*
-                              tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]
-                             -tmp_n.m[r2*col+r1]*tmp_n.m[r1*col+c]*
-                              tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand();
-                dividend_d = (tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]*
-                              tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand();
-                bool check = divide(dividend_n, divisor_n,
-                                    tmp_n.m[r2*col+c],true);
-                check &= divide(dividend_d, divisor_d,
-                                tmp_d.m[r2*col+c],true);
-                GINAC_ASSERT(check);
+            if ((r1<n-1)&&(r0<m-1)) {
+                // compute next iteration's divisor
+                divisor_n = tmp_n.m[r0*n+r1].expand();
+                divisor_d = tmp_d.m[r0*n+r1].expand();
+                if (det) {
+                    // save space by deleting no longer needed elements
+                    for (unsigned c=0; c<n; ++c) {
+                        tmp_n.m[r0*n+c] = _ex0();
+                        tmp_d.m[r0*n+c] = _ex1();
+                    }
+                }
             }
-            // fill up left hand side.
-            for (unsigned c=0; c<=r1; ++c)
-                tmp_n.m[r2*col+c] = _ex0();
+            ++r0;
         }
     }
     // repopulate *this matrix:
-    it = m.begin();
+    it = this->m.begin();
     tmp_n_it = tmp_n.m.begin();
     tmp_d_it = tmp_d.m.begin();
-    for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
+    for (; it!= this->m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
         (*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl);
     
     return sign;
@@ -1149,68 +1146,72 @@ int matrix::fraction_free_elimination(bool det)
  *  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 ro is the row from where to begin
+ *  @param co is the column 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)
+int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
 {
     unsigned k = ro;
-    
-    if (symbolic) {  // search first non-zero
-        for (unsigned r=ro; r<row; ++r) {
-            if (!m[r*col+ro].expand().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 (symbolic) {
+        // search first non-zero element in column co beginning at row ro
+        while ((k<row) && (this->m[k*col+co].expand().is_zero()))
+            ++k;
+    } else {
+        // search largest element in column co beginning at row ro
+        GINAC_ASSERT(is_ex_of_type(this->m[k*col+co],numeric));
+        unsigned kmax = k+1;
+        numeric mmax = abs(ex_to_numeric(m[kmax*col+co]));
+        while (kmax<row) {
+            GINAC_ASSERT(is_ex_of_type(this->m[kmax*col+co],numeric));
+            numeric tmp = ex_to_numeric(this->m[kmax*col+co]);
+            if (abs(tmp) > mmax) {
+                mmax = tmp;
+                k = kmax;
             }
+            ++kmax;
         }
+        if (!mmax.is_zero())
+            k = kmax;
     }
-    if (m[k*col+ro].is_zero())
+    if (k==row)
+        // all elements in column co below row ro vanish
         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]);
-        }
-        return k;
-    }
-    return 0;
+    if (k==ro)
+        // matrix needs no pivoting
+        return 0;
+    // matrix needs pivoting, so swap rows k and ro
+    ensure_if_modifiable();
+    for (unsigned c=0; c<col; ++c)
+        m[k*col+c].swap(m[ro*col+c]);
+    
+    return k;
 }
 
 /** Convert list of lists to matrix. */
 ex lst_to_matrix(const ex &l)
 {
-	if (!is_ex_of_type(l, lst))
-		throw(std::invalid_argument("argument to lst_to_matrix() must be a lst"));
-
-	// Find number of rows and columns
-	unsigned rows = l.nops(), cols = 0, i, j;
-	for (i=0; i<rows; i++)
-		if (l.op(i).nops() > cols)
-			cols = l.op(i).nops();
-
-	// Allocate and fill matrix
-	matrix &m = *new matrix(rows, cols);
-	for (i=0; i<rows; i++)
-		for (j=0; j<cols; j++)
-			if (l.op(i).nops() > j)
-				m.set(i, j, l.op(i).op(j));
-			else
-				m.set(i, j, ex(0));
-	return m;
+    if (!is_ex_of_type(l, lst))
+        throw(std::invalid_argument("argument to lst_to_matrix() must be a lst"));
+    
+    // Find number of rows and columns
+    unsigned rows = l.nops(), cols = 0, i, j;
+    for (i=0; i<rows; i++)
+        if (l.op(i).nops() > cols)
+            cols = l.op(i).nops();
+    
+    // Allocate and fill matrix
+    matrix &m = *new matrix(rows, cols);
+    for (i=0; i<rows; i++)
+        for (j=0; j<cols; j++)
+            if (l.op(i).nops() > j)
+                m.set(i, j, l.op(i).op(j));
+            else
+                m.set(i, j, ex(0));
+    return m;
 }
 
 //////////
diff --git a/ginac/matrix.h b/ginac/matrix.h
index b737a38c..baaa8a5b 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -87,19 +87,18 @@ public:
     const ex & operator() (unsigned ro, unsigned co) const;
     matrix & set(unsigned ro, unsigned co, ex value);
     matrix transpose(void) const;
-    ex determinant(unsigned options = determinant_algo::automatic) const;
+    ex determinant(unsigned algo = determinant_algo::automatic) const;
     ex trace(void) const;
     ex charpoly(const symbol & lambda) const;
     matrix inverse(void) const;
-    matrix fraction_free_elim(const matrix & vars, const matrix & v) const;
-    matrix solve(const matrix & vars, const matrix & rhs) const;
-    matrix old_solve(const matrix & v) const;  // FIXME: may be removed
+    matrix solve(const matrix & vars, const matrix & rhs,
+                 unsigned algo = solve_algo::automatic) const;
 protected:
     ex determinant_minor(void) const;
-    int gauss_elimination(void);
-    int division_free_elimination(void);
-    int fraction_free_elimination(bool det = false);
-    int pivot(unsigned ro, bool symbolic=true);
+    int gauss_elimination(const bool det = false);
+    int division_free_elimination(const bool det = false);
+    int fraction_free_elimination(const bool det = false);
+    int pivot(unsigned ro, unsigned co, bool symbolic = true);
     
 // member variables
 protected:
diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 0df296eb..bcf140d6 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -153,7 +153,8 @@ numeric::numeric(int i) : basic(TINFO_numeric)
     // emphasizes efficiency:
     value = new ::cl_I((long) i);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -166,7 +167,8 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
     // emphasizes efficiency:
     value = new ::cl_I((unsigned long)i);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -176,7 +178,8 @@ numeric::numeric(long i) : basic(TINFO_numeric)
     debugmsg("numeric constructor from long",LOGLEVEL_CONSTRUCT);
     value = new ::cl_I(i);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -186,7 +189,8 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric)
     debugmsg("numeric constructor from ulong",LOGLEVEL_CONSTRUCT);
     value = new ::cl_I(i);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -201,7 +205,8 @@ numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
     value = new ::cl_I(numer);
     *value = *value / ::cl_I(denom);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -215,7 +220,8 @@ numeric::numeric(double d) : basic(TINFO_numeric)
     value = new cl_N;
     *value = cl_float(d, cl_default_float_format);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -286,7 +292,8 @@ numeric::numeric(const cl_N & z) : basic(TINFO_numeric)
     debugmsg("numeric constructor from cl_N", LOGLEVEL_CONSTRUCT);
     value = new ::cl_N(z);
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
@@ -329,7 +336,8 @@ numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l
         }
     }
     calchash();
-    setflag(status_flags::evaluated|
+    setflag(status_flags::evaluated |
+            status_flags::expanded |
             status_flags::hash_calculated);
 }
 
diff --git a/ginac/power.cpp b/ginac/power.cpp
index fc52a5f8..53cafc52 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -525,11 +525,10 @@ ex power::derivative(const symbol & s) const
 {
     if (exponent.info(info_flags::real)) {
         // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
-        //return mul(mul(exponent, power(basis, exponent - _ex1())), basis.diff(s));
         epvector newseq;
         newseq.reserve(2);
         newseq.push_back(expair(basis, exponent - _ex1()));
-        newseq.push_back(expair(basis.diff(s),_ex1()));
+        newseq.push_back(expair(basis.diff(s), _ex1()));
         return mul(newseq, exponent);
     } else {
         // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
diff --git a/ginac/version.h b/ginac/version.h
index 35801017..c49a8dcb 100644
--- a/ginac/version.h
+++ b/ginac/version.h
@@ -26,6 +26,6 @@
 /* Major, minor, and micro version number of the GiNaC library. */
 #define GINACLIB_MAJOR_VERSION 0
 #define GINACLIB_MINOR_VERSION 6
-#define GINACLIB_MICRO_VERSION 3
+#define GINACLIB_MICRO_VERSION 4
 
 #endif // ndef __GINAC_VERSION_H__