- Made determinant_algo (in flags.h) really work.

- Fixed minor bugs in matrix.cpp.
- Small Documantation updates.

diff --git a/INSTALL b/INSTALL
index 3730e91c0489aaf8c329e6bbb316ecc0fbf2cc0e..10e0b58fce43fea90b9e5b091178dc965e7254b1 100644 (file)
--- a/INSTALL
+++ b/INSTALL
@@ -113,3 +113,4 @@ to compile, install and work properly:
  < 5.14.39    | `VERBOTEN' by license (please bite your favorite lawyer)
  < 5.14.39,40 | compiles but does not feel happy at all (inconsistent!)
    5.14.41    | tested on egcs 1.1.1, gcc 2.95.2: only minor weirdnesses
  < 5.14.39    | `VERBOTEN' by license (please bite your favorite lawyer)
  < 5.14.39,40 | compiles but does not feel happy at all (inconsistent!)
    5.14.41    | tested on egcs 1.1.1, gcc 2.95.2: only minor weirdnesses

+   5.14.44    | G__cpp_ginaccint.C needs manual fixes, doesn't work well

diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index c1fa6257599938bc57c015e363a6822bf876a9c2..3e702b8e963c174bdf37a896b91f89d7f5003d18 100644 (file)
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -31,15 +31,13 @@ static unsigned integdom_matrix_determinants(void)
     
     for (int size=3; size<20; ++size) {
         matrix A(size,size);
     
     for (int size=3; size<20; ++size) {
         matrix A(size,size);

-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
+        // populate one element in each row:
+        for (int r=0; r<size-1; ++r)

             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));

-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
+        // set the last row to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int c=0; c<size; ++c)

             A.set(size-1,c,A(0,c)-A(size-2,c));
             A.set(size-1,c,A(0,c)-A(size-2,c));

-        }

         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl


@@ -51,29 +49,30 @@ static unsigned integdom_matrix_determinants(void)
     return result;
 }
 
     return result;
 }
 

-/* determinants of some sparse symbolic matrices with multivariate rational
- * function coefficients. */
+/* determinants of some symbolic matrices with multivariate rational function
+ * coefficients. */

 static unsigned rational_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
 static unsigned rational_matrix_determinants(void)
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");

-
+    

     for (int size=3; size<8; ++size) {
         matrix A(size,size);
         for (int r=0; r<size-1; ++r) {
     for (int size=3; size<8; ++size) {
         matrix A(size,size);
         for (int r=0; r<size-1; ++r) {

-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 4, false, false, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, false, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
+            // populate one or two elements in each row:
+            for (int ec=0; ec<2; ++ec) {
+                ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
+                ex denom;
+                do {
+                    denom = sparse_tree(a, b, c, rand()%2, false, false, false);
+                } while (denom.is_zero());
+                A.set(r,unsigned(rand()%size),numer/denom);
+            }

         }
         }

+        // set the last row to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int co=0; co<size; ++co)
+            A.set(size-1,co,A(0,co)-A(size-2,co));

         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl


@@ -85,28 +84,29 @@ static unsigned rational_matrix_determinants(void)
     return result;
 }
 
     return result;
 }
 

-/* Some quite wild determinants with functions and stuff like that. */
-static unsigned wild_matrix_determinants(void)
+/* Some quite funny determinants with functions and stuff like that inside. */
+static unsigned funny_matrix_determinants(void)

 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
 {
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     

-    for (int size=3; size<6; ++size) {
+    for (int size=3; size<7; ++size) {

         matrix A(size,size);
         matrix A(size,size);

-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 3, true, true, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, true, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
+        for (int co=0; co<size-1; ++co) {
+            // populate one or two elements in each row:
+            for (int ec=0; ec<2; ++ec) {
+                ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
+                ex denom;
+                do {
+                    denom = sparse_tree(a, b, c, rand()%2, false, true, false);
+                } while (denom.is_zero());
+                A.set(unsigned(rand()%size),co,numer/denom);
+            }

         }
         }

+        // set the last column to a linear combination of two other lines
+        // to guarantee that the determinant is zero:
+        for (int ro=0; ro<size; ++ro)
+            A.set(ro,size-1,A(ro,0)-A(ro,size-2));

         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl


@@ -118,6 +118,41 @@ static unsigned wild_matrix_determinants(void)
     return result;
 }
 
     return result;
 }
 

+/* compare results from different determinant algorithms.*/
+static unsigned compare_matrix_determinants(void)
+{
+    unsigned result = 0;
+    symbol a("a");
+    
+    for (int size=2; size<6; ++size) {
+        matrix A(size,size);
+        for (int co=0; co<size; ++co) {
+            for (int ro=0; ro<size; ++ro) {
+                // populate some elements
+                ex elem = 0;
+                if (rand()%(size-1) == 0)
+                    elem = sparse_tree(a, a, a, rand()%3, false, true, false);
+                A.set(ro,co,elem);
+            }
+        }
+        ex det_gauss = A.determinant(determinant_algo::gauss);
+        ex det_laplace = A.determinant(determinant_algo::laplace);
+        ex det_bareiss = A.determinant(determinant_algo::bareiss);
+        if ((det_gauss-det_laplace).normal() != 0 ||
+            (det_bareiss-det_laplace).normal() != 0) {
+            clog << "Determinant of " << size << "x" << size << " matrix "
+                 << endl << A << endl
+                 << "is inconsistent between different algorithms:" << endl
+                 << "Gauss elimination:   " << det_gauss << endl
+                 << "Minor elimination:   " << det_laplace << endl
+                 << "Fraction-free elim.: " << det_bareiss << endl;
+            ++result;
+        }
+    }
+    
+    return result;
+}
+

 unsigned check_matrices(void)
 {
     unsigned result = 0;
 unsigned check_matrices(void)
 {
     unsigned result = 0;


@@ -127,7 +162,8 @@ unsigned check_matrices(void)
     
     result += integdom_matrix_determinants();  cout << '.' << flush;
     result += rational_matrix_determinants();  cout << '.' << flush;
     
     result += integdom_matrix_determinants();  cout << '.' << flush;
     result += rational_matrix_determinants();  cout << '.' << flush;

-    result += wild_matrix_determinants();  cout << '.' << flush;
+    result += funny_matrix_determinants();  cout << '.' << flush;
+    result += compare_matrix_determinants();  cout << '.' << flush;

     
     if (!result) {
         cout << " passed " << endl;
     
     if (!result) {
         cout << " passed " << endl;

diff --git a/check/exam_lsolve.cpp b/check/exam_lsolve.cpp
index bcbc25fee3959e1ff951c5fb78f313b86a832c87..3deddbc2c86ed8c9a4fe3bac3d8c586785c8fb15 100644 (file)
--- a/check/exam_lsolve.cpp
+++ b/check/exam_lsolve.cpp
@@ -32,7 +32,7 @@ static unsigned exam_lsolve1(void)
     eq = (3*x+5 == numeric(8));
     aux = lsolve(eq, x);
     if (aux != 1) {
     eq = (3*x+5 == numeric(8));
     aux = lsolve(eq, x);
     if (aux != 1) {

-        result++;
+        ++result;

         clog << "solution of 3*x+5==8 erroneously returned "
              << aux << endl;
     }
         clog << "solution of 3*x+5==8 erroneously returned "
              << aux << endl;
     }


@@ -60,7 +60,7 @@ static unsigned exam_lsolve2a(void)
     // It should have returned [x==(3+b^2)/(a+b),y==(3-a*b)/(a+b)]
     if (!normal(sol_x - (3+pow(b,2))/(a+b)).is_zero() ||
         !normal(sol_y - (3-a*b)/(a+b)).is_zero()) {
     // It should have returned [x==(3+b^2)/(a+b),y==(3-a*b)/(a+b)]
     if (!normal(sol_x - (3+pow(b,2))/(a+b)).is_zero() ||
         !normal(sol_y - (3-a*b)/(a+b)).is_zero()) {

-        result++;
+        ++result;

         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }


@@ -88,7 +88,7 @@ static unsigned exam_lsolve2b(void)
     // It should have returned [x==43/17,y==-10/17]
     if (!(sol_x - numeric(43,17)).is_zero() ||
         !(sol_y - numeric(-10,17)).is_zero()) {
     // It should have returned [x==43/17,y==-10/17]
     if (!(sol_x - numeric(43,17)).is_zero() ||
         !(sol_y - numeric(-10,17)).is_zero()) {

-        result++;
+        ++result;

         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }


@@ -116,7 +116,35 @@ static unsigned exam_lsolve2c(void)
     // It should have returned [x==-3/2*I,y==-1/2]
     if (!(sol_x - numeric(-3,2)*I).is_zero() ||
         !(sol_y - numeric(-1,2)).is_zero()) {
     // It should have returned [x==-3/2*I,y==-1/2]
     if (!(sol_x - numeric(-3,2)*I).is_zero() ||
         !(sol_y - numeric(-1,2)).is_zero()) {

-        result++;
+        ++result;
+        clog << "solution of the system " << eqns << " for " << vars
+             << " erroneously returned " << sol << endl;
+    }
+    
+    return result;
+}
+
+static unsigned exam_lsolve2S(void)
+{
+    // A degenerate example that went wrong in GiNaC 0.6.2.
+    unsigned result = 0;
+    symbol x("x"), y("y"), t("t");
+    lst eqns, vars;
+    ex sol;
+    
+    // Create the linear system [0*x+0*y==0,0*x+1*y==t]...
+    eqns.append(0*x+0*y==0).append(0*x+1*y==t);
+    // ...to be solved for [x,y]...
+    vars.append(x).append(y);
+    // ...and solve it:
+    sol = lsolve(eqns, vars);
+    ex sol_x = sol.op(0).rhs();  // rhs of solution for first variable (x)
+    ex sol_y = sol.op(1).rhs();  // rhs of solution for second variable (y)
+    
+    // It should have returned [x==x,y==t]
+    if (!(sol_x - x).is_zero() ||
+        !(sol_y - t).is_zero()) {
+        ++result;

         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }
         clog << "solution of the system " << eqns << " for " << vars
              << " erroneously returned " << sol << endl;
     }


@@ -135,6 +163,7 @@ unsigned exam_lsolve(void)
     result += exam_lsolve2a();  cout << '.' << flush;
     result += exam_lsolve2b();  cout << '.' << flush;
     result += exam_lsolve2c();  cout << '.' << flush;
     result += exam_lsolve2a();  cout << '.' << flush;
     result += exam_lsolve2b();  cout << '.' << flush;
     result += exam_lsolve2c();  cout << '.' << flush;

+    result += exam_lsolve2S();  cout << '.' << flush;

     
     if (!result) {
         cout << " passed " << endl;
     
     if (!result) {
         cout << " passed " << endl;

diff --git a/config.guess b/config.guess
index 499496436cd502d1fbd4861bd479191222eb825e..323f8c9c50d06043e9eab27672673dc149ea256f 100755 (executable)
--- a/config.guess
+++ b/config.guess
@@ -3,7 +3,7 @@
 #   Copyright (C) 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000
 #   Free Software Foundation, Inc.
 
 #   Copyright (C) 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000
 #   Free Software Foundation, Inc.
 

-version='2000-06-13'
+version='2000-07-24'

 
 # This file is free software; you can redistribute it and/or modify it
 # under the terms of the GNU General Public License as published by
 
 # This file is free software; you can redistribute it and/or modify it
 # under the terms of the GNU General Public License as published by


@@ -1088,6 +1088,17 @@ EOF
     DS/*:UNIX_System_V:*:*)
        echo ${UNAME_MACHINE}-${UNAME_SYSTEM}-${UNAME_RELEASE}
        exit 0 ;;
     DS/*:UNIX_System_V:*:*)
        echo ${UNAME_MACHINE}-${UNAME_SYSTEM}-${UNAME_RELEASE}
        exit 0 ;;

+    *:Plan9:*:*)
+       # "uname -m" is not consistent, so use $cputype instead. 386
+       # is converted to i386 for consistency with other x86
+       # operating systems.
+       if test "$cputype" = "386"; then
+           UNAME_MACHINE=i386
+       else
+           UNAME_MACHINE="$cputype"
+       fi
+       echo ${UNAME_MACHINE}-unknown-plan9
+       exit 0 ;;

 esac
 
 #echo '(No uname command or uname output not recognized.)' 1>&2
 esac
 
 #echo '(No uname command or uname output not recognized.)' 1>&2

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 536bf897a439267d05d4cc7f84001ffa382cd8f2..bc330c7d5aec5b7a4052dd01480f545c91178d79 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -2075,7 +2075,7 @@ negative real axis where the points on the axis itself belong to the
 upper part (i.e. continuous with quadrant II).  The inverse
 trigonometric and hyperbolic functions are not defined for complex
 arguments by the C++ standard, however.  Here, we follow the conventions
 upper part (i.e. continuous with quadrant II).  The inverse
 trigonometric and hyperbolic functions are not defined for complex
 arguments by the C++ standard, however.  Here, we follow the conventions

-used by CLN, which in turn follow the carefully structured definitions
+used by CLN, which in turn follow the carefully designed definitions

 in the Common Lisp standard.  Hopefully, future revisions of the C++
 standard incorporate these functions in the complex domain in a manner
 compatible with Common Lisp.
 in the Common Lisp standard.  Hopefully, future revisions of the C++
 standard incorporate these functions in the complex domain in a manner
 compatible with Common Lisp.

diff --git a/ginac/flags.h b/ginac/flags.h
index 77e3ec0117c7312f3f2bf28aeefc9dac97635652..41a890b11c34aa14c86cf13535d69bea5b824b7a 100644 (file)
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -39,9 +39,11 @@ public:
          };
 };
 
          };
 };
 

-class determinant_options {
+class determinant_algo {

 public:
 public:

-    enum { laplace,
+    enum { automatic,
+           gauss,
+           laplace,

            bareiss
          };
 };
            bareiss
          };
 };

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index e0d1e91edcaa0520ce5f94c017dff93e3e67a8a2..b560100f2437296cd6f752a372c35aeb7fcb8994 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -45,21 +45,21 @@ namespace GiNaC {
 // absolute value
 //////////
 
 // absolute value
 //////////
 

-static ex abs_evalf(const ex & x)
+static ex abs_evalf(const ex & arg)

 {
     BEGIN_TYPECHECK
 {
     BEGIN_TYPECHECK

-        TYPECHECK(x,numeric)
-    END_TYPECHECK(abs(x))
+        TYPECHECK(arg,numeric)
+    END_TYPECHECK(abs(arg))

     
     

-    return abs(ex_to_numeric(x));
+    return abs(ex_to_numeric(arg));

 }
 
 }
 

-static ex abs_eval(const ex & x)
+static ex abs_eval(const ex & arg)

 {
 {

-    if (is_ex_exactly_of_type(x, numeric))
-        return abs(ex_to_numeric(x));
+    if (is_ex_exactly_of_type(arg, numeric))
+        return abs(ex_to_numeric(arg));

     else
     else

-        return abs(x).hold();
+        return abs(arg).hold();

 }
 
 REGISTER_FUNCTION(abs, eval_func(abs_eval).
 }
 
 REGISTER_FUNCTION(abs, eval_func(abs_eval).


@@ -70,41 +70,41 @@ REGISTER_FUNCTION(abs, eval_func(abs_eval).
 // Complex sign
 //////////
 
 // Complex sign
 //////////
 

-static ex csgn_evalf(const ex & x)
+static ex csgn_evalf(const ex & arg)

 {
     BEGIN_TYPECHECK
 {
     BEGIN_TYPECHECK

-        TYPECHECK(x,numeric)
-    END_TYPECHECK(csgn(x))
+        TYPECHECK(arg,numeric)
+    END_TYPECHECK(csgn(arg))

     
     

-    return csgn(ex_to_numeric(x));
+    return csgn(ex_to_numeric(arg));

 }
 
 }
 

-static ex csgn_eval(const ex & x)
+static ex csgn_eval(const ex & arg)

 {
 {

-    if (is_ex_exactly_of_type(x, numeric))
-        return csgn(ex_to_numeric(x));
+    if (is_ex_exactly_of_type(arg, numeric))
+        return csgn(ex_to_numeric(arg));

     
     

-    else if (is_ex_exactly_of_type(x, mul)) {
-        numeric oc = ex_to_numeric(x.op(x.nops()-1));
+    else if (is_ex_exactly_of_type(arg, mul)) {
+        numeric oc = ex_to_numeric(arg.op(arg.nops()-1));

         if (oc.is_real()) {
             if (oc > 0)
                 // csgn(42*x) -> csgn(x)
         if (oc.is_real()) {
             if (oc > 0)
                 // csgn(42*x) -> csgn(x)

-                return csgn(x/oc).hold();
+                return csgn(arg/oc).hold();

             else
                 // csgn(-42*x) -> -csgn(x)
             else
                 // csgn(-42*x) -> -csgn(x)

-                return -csgn(x/oc).hold();
+                return -csgn(arg/oc).hold();

         }
         if (oc.real().is_zero()) {
             if (oc.imag() > 0)
                 // csgn(42*I*x) -> csgn(I*x)
         }
         if (oc.real().is_zero()) {
             if (oc.imag() > 0)
                 // csgn(42*I*x) -> csgn(I*x)

-                return csgn(I*x/oc).hold();
+                return csgn(I*arg/oc).hold();

             else
                 // csgn(-42*I*x) -> -csgn(I*x)
             else
                 // csgn(-42*I*x) -> -csgn(I*x)

-                return -csgn(I*x/oc).hold();
+                return -csgn(I*arg/oc).hold();

         }
        }
    
         }
        }
    

-    return csgn(x).hold();
+    return csgn(arg).hold();

 }
 
 static ex csgn_series(const ex & arg,
 }
 
 static ex csgn_series(const ex & arg,


@@ -113,13 +113,10 @@ static ex csgn_series(const ex & arg,
                       unsigned options)
 {
     const ex arg_pt = arg.subs(rel);
                       unsigned options)
 {
     const ex arg_pt = arg.subs(rel);

-    if (arg_pt.info(info_flags::numeric)) {
-        if (ex_to_numeric(arg_pt).real().is_zero())
-            throw (std::domain_error("csgn_series(): on imaginary axis"));
-        epvector seq;
-        seq.push_back(expair(csgn(arg_pt), _ex0()));
-        return pseries(rel,seq);
-    }
+    if (arg_pt.info(info_flags::numeric) &&
+        ex_to_numeric(arg_pt).real().is_zero())
+        throw (std::domain_error("csgn_series(): on imaginary axis"));
+    

     epvector seq;
     seq.push_back(expair(csgn(arg_pt), _ex0()));
     return pseries(rel,seq);
     epvector seq;
     seq.push_back(expair(csgn(arg_pt), _ex0()));
     return pseries(rel,seq);


@@ -129,6 +126,61 @@ REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
                         evalf_func(csgn_evalf).
                         series_func(csgn_series));
 
                         evalf_func(csgn_evalf).
                         series_func(csgn_series));
 

+
+//////////
+// Eta function: log(x*y) == log(x) + log(y) + eta(x,y).
+//////////
+
+static ex eta_evalf(const ex & x, const ex & y)
+{
+    BEGIN_TYPECHECK
+        TYPECHECK(x,numeric)
+        TYPECHECK(y,numeric)
+    END_TYPECHECK(eta(x,y))
+        
+    numeric xim = imag(ex_to_numeric(x));
+    numeric yim = imag(ex_to_numeric(y));
+    numeric xyim = imag(ex_to_numeric(x*y));
+    return evalf(I/4*Pi)*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1));
+}
+
+static ex eta_eval(const ex & x, const ex & y)
+{
+    if (is_ex_exactly_of_type(x, numeric) &&
+        is_ex_exactly_of_type(y, numeric)) {
+        // don't call eta_evalf here because it would call Pi.evalf()!
+        numeric xim = imag(ex_to_numeric(x));
+        numeric yim = imag(ex_to_numeric(y));
+        numeric xyim = imag(ex_to_numeric(x*y));
+        return (I/4)*Pi*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1));
+    }
+    
+    return eta(x,y).hold();
+}
+
+static ex eta_series(const ex & arg1,
+                     const ex & arg2,
+                     const relational & rel,
+                     int order,
+                     unsigned options)
+{
+    const ex arg1_pt = arg1.subs(rel);
+    const ex arg2_pt = arg2.subs(rel);
+    if (ex_to_numeric(arg1_pt).imag().is_zero() ||
+        ex_to_numeric(arg2_pt).imag().is_zero() ||
+        ex_to_numeric(arg1_pt*arg2_pt).imag().is_zero()) {
+        throw (std::domain_error("eta_series(): on discontinuity"));
+    }
+    epvector seq;
+    seq.push_back(expair(eta(arg1_pt,arg2_pt), _ex0()));
+    return pseries(rel,seq);
+}
+
+REGISTER_FUNCTION(eta, eval_func(eta_eval).
+                       evalf_func(eta_evalf).
+                       series_func(eta_series));
+
+

 //////////
 // dilogarithm
 //////////
 //////////
 // dilogarithm
 //////////

diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 9781d935c1a2e777a2dfbf9326fba1da8eae8d22..95ed15a1606404e458a41e3aab64801e21da1f2b 100644 (file)
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -36,6 +36,9 @@ DECLARE_FUNCTION_1P(abs)
 /** Complex sign. */
 DECLARE_FUNCTION_1P(csgn)
 
 /** Complex sign. */
 DECLARE_FUNCTION_1P(csgn)
 

+/** Eta function: log(a*b) == log(a) + log(b) + eta(a, b). */
+DECLARE_FUNCTION_2P(eta)
+

 /** Sine. */
 DECLARE_FUNCTION_1P(sin)
 
 /** Sine. */
 DECLARE_FUNCTION_1P(sin)
 


@@ -87,7 +90,7 @@ DECLARE_FUNCTION_1P(Li2)
 /** Trilogarithm. */
 DECLARE_FUNCTION_1P(Li3)
 
 /** Trilogarithm. */
 DECLARE_FUNCTION_1P(Li3)
 

-// overloading at work: we cannot use the macros
+// overloading at work: we cannot use the macros here

 /** Riemann's Zeta-function. */
 extern const unsigned function_index_zeta1;
 inline function zeta(const ex & p1) {
 /** Riemann's Zeta-function. */
 extern const unsigned function_index_zeta1;
 inline function zeta(const ex & p1) {


@@ -106,7 +109,7 @@ DECLARE_FUNCTION_1P(tgamma)
 /** Beta-function. */
 DECLARE_FUNCTION_2P(beta)
 
 /** Beta-function. */
 DECLARE_FUNCTION_2P(beta)
 

-// overloading at work: we cannot use the macros
+// overloading at work: we cannot use the macros here

 /** Psi-function (aka digamma-function). */
 extern const unsigned function_index_psi1;
 inline function psi(const ex & p1) {
 /** Psi-function (aka digamma-function). */
 extern const unsigned function_index_psi1;
 inline function psi(const ex & p1) {

diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index bf60c1e5884fe9616f51055e54a4aa827fecb364..aa1cab4c1f365f8edca64acbc80195d8804647a9 100644 (file)
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -63,11 +63,10 @@ static ex lgamma_eval(const ex & x)
         // trap integer arguments:
         if (x.info(info_flags::integer)) {
             // lgamma(n) -> log((n-1)!) for postitive n
         // trap integer arguments:
         if (x.info(info_flags::integer)) {
             // lgamma(n) -> log((n-1)!) for postitive n

-            if (x.info(info_flags::posint)) {
+            if (x.info(info_flags::posint))

                 return log(factorial(x.exadd(_ex_1())));
                 return log(factorial(x.exadd(_ex_1())));

-            } else {
+            else

                 throw (pole_error("lgamma_eval(): logarithmic pole",0));
                 throw (pole_error("lgamma_eval(): logarithmic pole",0));

-            }

         }
         //  lgamma_evalf should be called here once it becomes available
     }
         }
         //  lgamma_evalf should be called here once it becomes available
     }


@@ -98,15 +97,16 @@ static ex lgamma_series(const ex & arg,
     // from which follows
     //   series(lgamma(x),x==-m,order) ==
     //   series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order);
     // from which follows
     //   series(lgamma(x),x==-m,order) ==
     //   series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order);

-    // This, however, seems to fail utterly because you run into branch-cut
-    // problems.  Somebody ought to implement it some day using an asymptotic
-    // series for tgamma:

     const ex arg_pt = arg.subs(rel);
     if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
     // if we got here we have to care for a simple pole of tgamma(-m):
     const ex arg_pt = arg.subs(rel);
     if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
     // if we got here we have to care for a simple pole of tgamma(-m):

-    throw (std::overflow_error("lgamma_series: please implement my at the poles"));
-    return _ex0();  // not reached
+    numeric m = -ex_to_numeric(arg_pt);
+    ex recur;
+    for (numeric p; p<=m; ++p)
+        recur += log(arg+p);
+    cout << recur << endl;
+    return (lgamma(arg+m+_ex1())-recur).series(rel, order, options);

 }
 
 
 }
 
 


@@ -202,7 +202,7 @@ static ex tgamma_series(const ex & arg,
     ex ser_denom = _ex1();
     for (numeric p; p<=m; ++p)
         ser_denom *= arg+p;
     ex ser_denom = _ex1();
     for (numeric p; p<=m; ++p)
         ser_denom *= arg+p;

-    return (tgamma(arg+m+_ex1())/ser_denom).series(rel, order+1);
+    return (tgamma(arg+m+_ex1())/ser_denom).series(rel, order+1, options);

 }
 
 
 }
 
 


@@ -299,21 +299,21 @@ static ex beta_series(const ex & arg1,
         throw do_taylor();  // caught by function::series()
     // trap the case where arg1 is on a pole:
     if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
     // trap the case where arg1 is on a pole:
     if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive))

-        arg1_ser = tgamma(arg1+*s).series(rel,order);
+        arg1_ser = tgamma(arg1+*s).series(rel, order, options);

     else
         arg1_ser = tgamma(arg1).series(rel,order);
     // trap the case where arg2 is on a pole:
     if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive))
     else
         arg1_ser = tgamma(arg1).series(rel,order);
     // trap the case where arg2 is on a pole:
     if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive))

-        arg2_ser = tgamma(arg2+*s).series(rel,order);
+        arg2_ser = tgamma(arg2+*s).series(rel, order, options);

     else
         arg2_ser = tgamma(arg2).series(rel,order);
     // trap the case where arg1+arg2 is on a pole:
     if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive))
     else
         arg2_ser = tgamma(arg2).series(rel,order);
     // trap the case where arg1+arg2 is on a pole:
     if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive))

-        arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel,order);
+        arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel, order, options);

     else
         arg1arg2_ser = tgamma(arg2+arg1).series(rel,order);
     // compose the result (expanding all the terms):
     else
         arg1arg2_ser = tgamma(arg2+arg1).series(rel,order);
     // compose the result (expanding all the terms):

-    return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel,order).expand();
+    return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel, order, options).expand();

 }
 
 
 }
 
 


@@ -410,7 +410,7 @@ static ex psi1_series(const ex & arg,
     ex recur;
     for (numeric p; p<=m; ++p)
         recur += power(arg+p,_ex_1());
     ex recur;
     for (numeric p; p<=m; ++p)
         recur += power(arg+p,_ex_1());

-    return (psi(arg+m+_ex1())-recur).series(rel, order);
+    return (psi(arg+m+_ex1())-recur).series(rel, order, options);

 }
 
 const unsigned function_index_psi1 =
 }
 
 const unsigned function_index_psi1 =


@@ -536,7 +536,7 @@ static ex psi2_series(const ex & n,
     for (numeric p; p<=m; ++p)
         recur += power(arg+p,-n+_ex_1());
     recur *= factorial(n)*power(_ex_1(),n);
     for (numeric p; p<=m; ++p)
         recur += power(arg+p,-n+_ex_1());
     recur *= factorial(n)*power(_ex_1(),n);

-    return (psi(n, arg+m+_ex1())-recur).series(rel, order);
+    return (psi(n, arg+m+_ex1())-recur).series(rel, order, options);

 }
 
 const unsigned function_index_psi2 =
 }
 
 const unsigned function_index_psi2 =

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index feae32fee687341c7a92a5d8ea4e6693dc61dc88..52b4de36555dfbc10953e1fb8a3c882f7abb80b5 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -152,7 +152,7 @@ void matrix::archive(archive_node &n) const
     exvector::const_iterator i = m.begin(), iend = m.end();
     while (i != iend) {
         n.add_ex("m", *i);
     exvector::const_iterator i = m.begin(), iend = m.end();
     while (i != iend) {
         n.add_ex("m", *i);

-        i++;
+        ++i;

     }
 }
 
     }
 }
 


@@ -174,15 +174,13 @@ void matrix::print(std::ostream & os, unsigned upper_precedence) const
     os << "[[ ";
     for (unsigned r=0; r<row-1; ++r) {
         os << "[[";
     os << "[[ ";
     for (unsigned r=0; r<row-1; ++r) {
         os << "[[";

-        for (unsigned c=0; c<col-1; ++c) {
+        for (unsigned c=0; c<col-1; ++c)

             os << m[r*col+c] << ",";
             os << m[r*col+c] << ",";

-        }

         os << m[col*(r+1)-1] << "]], ";
     }
     os << "[[";
         os << m[col*(r+1)-1] << "]], ";
     }
     os << "[[";

-    for (unsigned c=0; c<col-1; ++c) {
+    for (unsigned c=0; c<col-1; ++c)

         os << m[(row-1)*col+c] << ",";
         os << m[(row-1)*col+c] << ",";

-    }

     os << m[row*col-1] << "]] ]]";
 }
 
     os << m[row*col-1] << "]] ]]";
 }
 


@@ -192,15 +190,13 @@ void matrix::printraw(std::ostream & os) const
     os << "matrix(" << row << "," << col <<",";
     for (unsigned r=0; r<row-1; ++r) {
         os << "(";
     os << "matrix(" << row << "," << col <<",";
     for (unsigned r=0; r<row-1; ++r) {
         os << "(";

-        for (unsigned c=0; c<col-1; ++c) {
+        for (unsigned c=0; c<col-1; ++c)

             os << m[r*col+c] << ",";
             os << m[r*col+c] << ",";

-        }

         os << m[col*(r-1)-1] << "),";
     }
     os << "(";
         os << m[col*(r-1)-1] << "),";
     }
     os << "(";

-    for (unsigned c=0; c<col-1; ++c) {
+    for (unsigned c=0; c<col-1; ++c)

         os << m[(row-1)*col+c] << ",";
         os << m[(row-1)*col+c] << ",";

-    }

     os << m[row*col-1] << "))";
 }
 
     os << m[row*col-1] << "))";
 }
 


@@ -219,6 +215,9 @@ ex matrix::op(int i) const
 /** returns matrix entry at position (i/col, i%col). */
 ex & matrix::let_op(int i)
 {
 /** returns matrix entry at position (i/col, i%col). */
 ex & matrix::let_op(int i)
 {

+    GINAC_ASSERT(i>=0);
+    GINAC_ASSERT(i<nops());
+    

     return m[i];
 }
 
     return m[i];
 }
 


@@ -226,9 +225,9 @@ ex & matrix::let_op(int i)
 ex matrix::expand(unsigned options) const
 {
     exvector tmp(row*col);
 ex matrix::expand(unsigned options) const
 {
     exvector tmp(row*col);

-    for (unsigned i=0; i<row*col; ++i) {
-        tmp[i]=m[i].expand(options);
-    }
+    for (unsigned i=0; i<row*col; ++i)
+        tmp[i] = m[i].expand(options);
+    

     return matrix(row, col, tmp);
 }
 
     return matrix(row, col, tmp);
 }
 


@@ -242,9 +241,9 @@ bool matrix::has(const ex & other) const
     if (is_equal(*other.bp)) return true;
     
     // search all the elements
     if (is_equal(*other.bp)) return true;
     
     // search all the elements

-    for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
+    for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r)

         if ((*r).has(other)) return true;
         if ((*r).has(other)) return true;

-    }
+    

     return false;
 }
 
     return false;
 }
 


@@ -263,12 +262,10 @@ ex matrix::eval(int level) const
     
     // eval() entry by entry
     exvector m2(row*col);
     
     // eval() entry by entry
     exvector m2(row*col);

-    --level;    
-    for (unsigned r=0; r<row; ++r) {
-        for (unsigned c=0; c<col; ++c) {
+    --level;
+    for (unsigned r=0; r<row; ++r)
+        for (unsigned c=0; c<col; ++c)

             m2[r*col+c] = m[r*col+c].eval(level);
             m2[r*col+c] = m[r*col+c].eval(level);

-        }
-    }

     
     return (new matrix(row, col, m2))->setflag(status_flags::dynallocated |
                                                status_flags::evaluated );
     
     return (new matrix(row, col, m2))->setflag(status_flags::dynallocated |
                                                status_flags::evaluated );


@@ -291,11 +288,10 @@ ex matrix::evalf(int level) const
     // evalf() entry by entry
     exvector m2(row*col);
     --level;
     // evalf() entry by entry
     exvector m2(row*col);
     --level;

-    for (unsigned r=0; r<row; ++r) {
-        for (unsigned c=0; c<col; ++c) {
+    for (unsigned r=0; r<row; ++r)
+        for (unsigned c=0; c<col; ++c)

             m2[r*col+c] = m[r*col+c].evalf(level);
             m2[r*col+c] = m[r*col+c].evalf(level);

-        }
-    }
+    

     return matrix(row, col, m2);
 }
 
     return matrix(row, col, m2);
 }
 


@@ -343,11 +339,9 @@ matrix matrix::add(const matrix & other) const
     exvector sum(this->m);
     exvector::iterator i;
     exvector::const_iterator ci;
     exvector sum(this->m);
     exvector::iterator i;
     exvector::const_iterator ci;

-    for (i=sum.begin(), ci=other.m.begin();
-         i!=sum.end();
-         ++i, ++ci) {
+    for (i=sum.begin(), ci=other.m.begin(); i!=sum.end(); ++i, ++ci)

         (*i) += (*ci);
         (*i) += (*ci);

-    }
+    

     return matrix(row,col,sum);
 }
 
     return matrix(row,col,sum);
 }
 


@@ -363,11 +357,9 @@ matrix matrix::sub(const matrix & other) const
     exvector dif(this->m);
     exvector::iterator i;
     exvector::const_iterator ci;
     exvector dif(this->m);
     exvector::iterator i;
     exvector::const_iterator ci;

-    for (i=dif.begin(), ci=other.m.begin();
-         i!=dif.end();
-         ++i, ++ci) {
+    for (i=dif.begin(), ci=other.m.begin(); i!=dif.end(); ++i, ++ci)

         (*i) -= (*ci);
         (*i) -= (*ci);

-    }
+    

     return matrix(row,col,dif);
 }
 
     return matrix(row,col,dif);
 }
 


@@ -438,7 +430,7 @@ matrix matrix::transpose(void) const
 
 /** Determinant of square matrix.  This routine doesn't actually calculate the
  *  determinant, it only implements some heuristics about which algorithm to
 
 /** Determinant of square matrix.  This routine doesn't actually calculate the
  *  determinant, it only implements some heuristics about which algorithm to

- *  call.  If all the elements of the matrix are elements of an integral domain
+ *  run.  If all the elements of the matrix are elements of an integral domain

  *  the determinant is also in that integral domain and the result is expanded
  *  only.  If one or more elements are from a quotient field the determinant is
  *  usually also in that quotient field and the result is normalized before it
  *  the determinant is also in that integral domain and the result is expanded
  *  only.  If one or more elements are from a quotient field the determinant is
  *  usually also in that quotient field and the result is normalized before it


@@ -446,85 +438,105 @@ matrix matrix::transpose(void) const
  *  [[a/(a-b),1],[b/(a-b),1]] is returned as unity.  (In this respect, it
  *  behaves like MapleV and unlike Mathematica.)
  *
  *  [[a/(a-b),1],[b/(a-b),1]] is returned as unity.  (In this respect, it
  *  behaves like MapleV and unlike Mathematica.)
  *

+ *  @param     algo allows to chose an algorithm

  *  @return    the determinant as a new expression
  *  @return    the determinant as a new expression

- *  @exception logic_error (matrix not square) */
-ex matrix::determinant(void) const
+ *  @exception logic_error (matrix not square)
+ *  @see       determinant_algo */
+ex matrix::determinant(unsigned algo) const

 {
     if (row!=col)
         throw (std::logic_error("matrix::determinant(): matrix not square"));
     GINAC_ASSERT(row*col==m.capacity());
     if (this->row==1)  // continuation would be pointless
         return m[0];
 {
     if (row!=col)
         throw (std::logic_error("matrix::determinant(): matrix not square"));
     GINAC_ASSERT(row*col==m.capacity());
     if (this->row==1)  // continuation would be pointless
         return m[0];

-    
-    // Gather some information about the matrix:
+        
+    // Gather some statistical information about this matrix:

     bool numeric_flag = true;
     bool normal_flag = false;
     bool numeric_flag = true;
     bool normal_flag = false;

-    unsigned sparse_count = 0;  // count non-zero elements
+    unsigned sparse_count = 0;  // counts non-zero elements

     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {

-        if (!(*r).is_zero())
+        lst srl;  // symbol replacement list
+        ex rtest = (*r).to_rational(srl);
+        if (!rtest.is_zero())

             ++sparse_count;
             ++sparse_count;

-        if (!(*r).info(info_flags::numeric))
+        if (!rtest.info(info_flags::numeric))

             numeric_flag = false;
             numeric_flag = false;

-        if ((*r).info(info_flags::rational_function) &&
-            !(*r).info(info_flags::crational_polynomial))
+        if (!rtest.info(info_flags::crational_polynomial) &&
+             rtest.info(info_flags::rational_function))

             normal_flag = true;
     }
     
             normal_flag = true;
     }
     

-    // Purely numeric matrix handled by Gauss elimination
-    if (numeric_flag) {
-        ex det = 1;
-        matrix tmp(*this);
-        int sign = tmp.gauss_elimination();
-        for (int d=0; d<row; ++d)
-            det *= tmp.m[d*col+d];
-        return (sign*det);
-    }
-    
-    // Does anybody know when a matrix is really sparse?
-    // Maybe <~row/2.2 nonzero elements average in a row?
-    if (5*sparse_count<=row*col) {
-        // copy *this:
-        matrix tmp(*this);
-        int sign;
-        sign = tmp.fraction_free_elimination(true);
-        if (normal_flag)
-            return (sign*tmp.m[row*col-1]).normal();
-        else
-            return (sign*tmp.m[row*col-1]).expand();
+    // Here is the heuristics in case this routine has to decide:
+    if (algo == determinant_algo::automatic) {
+        // Minor expansion is generally a good starting point:
+        algo = determinant_algo::laplace;
+        // Does anybody know when a matrix is really sparse?
+        // Maybe <~row/2.236 nonzero elements average in a row?
+        if (5*sparse_count<=row*col)
+            algo = determinant_algo::bareiss;
+        // Purely numeric matrix can be handled by Gauss elimination.
+        // This overrides any prior decisions.
+        if (numeric_flag)
+            algo = determinant_algo::gauss;

     }
     
     }
     

-    // Now come the minor expansion schemes.  We always develop such that the
-    // smallest minors (i.e, the trivial 1x1 ones) are on the rightmost column.
-    // For this to be efficient it turns out that the emptiest columns (i.e.
-    // the ones with most zeros) should be the ones on the right hand side.
-    // Therefore we presort the columns of the matrix:
-    typedef std::pair<unsigned,unsigned> uintpair;  // # of zeros, column
-    std::vector<uintpair> c_zeros;  // number of zeros in column
-    for (unsigned c=0; c<col; ++c) {
-        unsigned acc = 0;
-        for (unsigned r=0; r<row; ++r)
-            if (m[r*col+c].is_zero())
-                ++acc;
-        c_zeros.push_back(uintpair(acc,c));
-    }
-    sort(c_zeros.begin(),c_zeros.end());
-    // unfortunately std::vector<uintpair> can't be used for permutation_sign.
-    std::vector<unsigned> pre_sort;
-    for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
-        pre_sort.push_back(i->second);
-    int sign = permutation_sign(pre_sort);
-    exvector result(row*col);  // represents sorted matrix
-    unsigned c = 0;
-    for (std::vector<unsigned>::iterator i=pre_sort.begin();
-         i!=pre_sort.end();
-         ++i,++c) {
-        for (unsigned r=0; r<row; ++r)
-            result[r*col+c] = m[r*col+(*i)];
+    switch(algo) {
+        case determinant_algo::gauss: {
+            ex det = 1;
+            matrix tmp(*this);
+            int sign = tmp.gauss_elimination();
+            for (unsigned d=0; d<row; ++d)
+                det *= tmp.m[d*col+d];
+            if (normal_flag)
+                return (sign*det).normal();
+            else
+                return (sign*det).expand();
+        }
+        case determinant_algo::bareiss: {
+            matrix tmp(*this);
+            int sign;
+            sign = tmp.fraction_free_elimination(true);
+            if (normal_flag)
+                return (sign*tmp.m[row*col-1]).normal();
+            else
+                return (sign*tmp.m[row*col-1]).expand();
+        }
+        case determinant_algo::laplace:
+        default: {
+            // This is the minor expansion scheme.  We always develop such
+            // that the smallest minors (i.e, the trivial 1x1 ones) are on the
+            // rightmost column.  For this to be efficient it turns out that
+            // the emptiest columns (i.e. the ones with most zeros) should be
+            // the ones on the right hand side.  Therefore we presort the
+            // columns of the matrix:
+            typedef std::pair<unsigned,unsigned> uintpair;
+            std::vector<uintpair> c_zeros;  // number of zeros in column
+            for (unsigned c=0; c<col; ++c) {
+                unsigned acc = 0;
+                for (unsigned r=0; r<row; ++r)
+                    if (m[r*col+c].is_zero())
+                        ++acc;
+                c_zeros.push_back(uintpair(acc,c));
+            }
+            sort(c_zeros.begin(),c_zeros.end());
+            std::vector<unsigned> pre_sort;
+            for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
+                pre_sort.push_back(i->second);
+            int sign = permutation_sign(pre_sort);
+            exvector result(row*col);  // represents sorted matrix
+            unsigned c = 0;
+            for (std::vector<unsigned>::iterator i=pre_sort.begin();
+                 i!=pre_sort.end();
+                 ++i,++c) {
+                for (unsigned r=0; r<row; ++r)
+                    result[r*col+c] = m[r*col+(*i)];
+            }
+            
+            if (normal_flag)
+                return sign*matrix(row,col,result).determinant_minor().normal();
+            return sign*matrix(row,col,result).determinant_minor();
+        }

     }
     }

-    
-    if (normal_flag)
-        return sign*matrix(row,col,result).determinant_minor().normal();
-    return sign*matrix(row,col,result).determinant_minor();

 }
 
 
 }
 
 


@@ -647,16 +659,6 @@ matrix matrix::inverse(void) const
     return tmp;
 }
 
     return tmp;
 }
 

-// superfluous helper function, to be removed:
-void matrix::swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
-{
-    ensure_if_modifiable();
-    
-    ex tmp = (*this)(r1,c1);
-    set(r1,c1,(*this)(r2,c2));
-    set(r2,c2,tmp);
-}
-

 
 /** Solve a set of equations for an m x n matrix by fraction-free Gaussian
  *  elimination.  Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
 
 /** Solve a set of equations for an m x n matrix by fraction-free Gaussian
  *  elimination.  Based on algorithm 9.1 from 'Algorithms for Computer Algebra'


@@ -669,7 +671,7 @@ void matrix::swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
 matrix matrix::fraction_free_elim(const matrix & vars,
                                   const matrix & rhs) const
 {
 matrix matrix::fraction_free_elim(const matrix & vars,
                                   const matrix & rhs) const
 {

-    // FIXME: use implementation of matrix::fraction_free_elimination
+    // FIXME: use implementation of matrix::fraction_free_elimination instead!

     if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
         throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
     
     if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
         throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
     


@@ -687,29 +689,29 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned k=0; (k<n)&&(r<m); ++k) {
         // find a nonzero pivot
         unsigned p;
     for (unsigned k=0; (k<n)&&(r<m); ++k) {
         // find a nonzero pivot
         unsigned p;

-        for (p=r; (p<m)&&(a(p,k).is_zero()); ++p) {}
+        for (p=r; (p<m)&&(a.m[p*a.cols()+k].is_zero()); ++p) {}

         // pivot is in row p
         if (p<m) {
             if (p!=r) {
                 // swap rows p and r
                 for (unsigned j=k; j<n; ++j)
         // pivot is in row p
         if (p<m) {
             if (p!=r) {
                 // swap rows p and r
                 for (unsigned j=k; j<n; ++j)

-                    a.swap(p,j,r,j);
-                b.swap(p,0,r,0);
+                    a.m[p*a.cols()+j].swap(a.m[r*a.cols()+j]);
+                b.m[p*b.cols()].swap(b.m[r*b.cols()]);

                 // keep track of sign changes due to row exchange
                 sign *= -1;
             }
             for (unsigned i=r+1; i<m; ++i) {
                 for (unsigned j=k+1; j<n; ++j) {
                 // keep track of sign changes due to row exchange
                 sign *= -1;
             }
             for (unsigned i=r+1; i<m; ++i) {
                 for (unsigned j=k+1; j<n; ++j) {

-                    a.set(i,j,(a(r,k)*a(i,j)
-                              -a(r,j)*a(i,k))/divisor);
-                    a.set(i,j,a(i,j).normal());
+                    a.set(i,j,(a.m[r*a.cols()+k]*a.m[i*a.cols()+j]
+                              -a.m[r*a.cols()+j]*a.m[i*a.cols()+k])/divisor);
+                    a.set(i,j,a.m[i*a.cols()+j].normal());

                 }
                 }

-                b.set(i,0,(a(r,k)*b(i,0)
-                          -b(r,0)*a(i,k))/divisor);
-                b.set(i,0,b(i,0).normal());
-                a.set(i,k,0);
+                b.set(i,0,(a.m[r*a.cols()+k]*b.m[i*b.cols()]
+                          -b.m[r*b.cols()]*a.m[i*a.cols()+k])/divisor);
+                b.set(i,0,b.m[i*b.cols()].normal());
+                a.set(i,k,_ex0());

             }
             }

-            divisor = a(r,k);
+            divisor = a.m[r*a.cols()+k];

             ++r;
         }
     }
             ++r;
         }
     }


@@ -720,8 +722,8 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned r=0; r<m; ++r) {
         int zero_in_this_row=0;
         for (unsigned c=0; c<n; ++c) {
     for (unsigned r=0; r<m; ++r) {
         int zero_in_this_row=0;
         for (unsigned c=0; c<n; ++c) {

-            if (a(r,c).is_zero())
-               zero_in_this_row++;
+            if (a.m[r*a.cols()+c].is_zero())
+               ++zero_in_this_row;

             else
                 break;
         }
             else
                 break;
         }


@@ -739,26 +741,26 @@ matrix matrix::fraction_free_elim(const matrix & vars,
             first_non_zero++;
         if (first_non_zero>n) {
             // row consists only of zeroes, corresponding rhs must be 0 as well
             first_non_zero++;
         if (first_non_zero>n) {
             // row consists only of zeroes, corresponding rhs must be 0 as well

-            if (!b(r,0).is_zero()) {
+            if (!b.m[r*b.cols()].is_zero()) {

                 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
             }
         } else {
             // assign solutions for vars between first_non_zero+1 and
             // last_assigned_sol-1: free parameters
             for (unsigned c=first_non_zero; c<last_assigned_sol-1; ++c)
                 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
             }
         } else {
             // assign solutions for vars between first_non_zero+1 and
             // last_assigned_sol-1: free parameters
             for (unsigned c=first_non_zero; c<last_assigned_sol-1; ++c)

-                sol.set(c,0,vars(c,0));
-            ex e = b(r,0);
+                sol.set(c,0,vars.m[c*vars.cols()]);
+            ex e = b.m[r*b.cols()];

             for (unsigned c=first_non_zero; c<n; ++c)
             for (unsigned c=first_non_zero; c<n; ++c)

-                e -= a(r,c)*sol(c,0);
+                e -= a.m[r*a.cols()+c]*sol.m[c*sol.cols()];

             sol.set(first_non_zero-1,0,
             sol.set(first_non_zero-1,0,

-                    (e/a(r,first_non_zero-1)).normal());
+                    (e/(a.m[r*a.cols()+(first_non_zero-1)])).normal());

             last_assigned_sol = first_non_zero;
         }
     }
     // assign solutions for vars between 1 and
     // last_assigned_sol-1: free parameters
     for (unsigned c=0; c<last_assigned_sol-1; ++c)
             last_assigned_sol = first_non_zero;
         }
     }
     // assign solutions for vars between 1 and
     // last_assigned_sol-1: free parameters
     for (unsigned c=0; c<last_assigned_sol-1; ++c)

-        sol.set(c,0,vars(c,0));
+        sol.set(c,0,vars.m[c*vars.cols()]);

     
 #ifdef DO_GINAC_ASSERT
     // test solution with echelon matrix
     
 #ifdef DO_GINAC_ASSERT
     // test solution with echelon matrix


@@ -1159,7 +1161,7 @@ int matrix::pivot(unsigned ro, bool symbolic)
     
     if (symbolic) {  // search first non-zero
         for (unsigned r=ro; r<row; ++r) {
     
     if (symbolic) {  // search first non-zero
         for (unsigned r=ro; r<row; ++r) {

-            if (!m[r*col+ro].is_zero()) {
+            if (!m[r*col+ro].expand().is_zero()) {

                 k = r;
                 break;
             }
                 k = r;
                 break;
             }

diff --git a/ginac/matrix.h b/ginac/matrix.h
index cef3dc36667b3c08d5b4c2c4a7525f5a6efcec0e..b737a38c72fec2e74b118cf57aaf285a295bb6ce 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -27,11 +27,6 @@
 #include "basic.h"
 #include "ex.h"
 
 #include "basic.h"
 #include "ex.h"
 

-namespace std {
-    // forward declaration, so <stdexcept> need not be included:
-    class range_error;
-}
-

 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
 #endif // ndef NO_NAMESPACE_GINAC
 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
 #endif // ndef NO_NAMESPACE_GINAC


@@ -82,9 +77,9 @@ protected:
     
     // non-virtual functions in this class
 public:
     
     // non-virtual functions in this class
 public:

-    unsigned rows(void) const        //! get number of rows.
+    unsigned rows(void) const        //! Get number of rows.

         { return row; }
         { return row; }

-    unsigned cols(void) const        //! get number of columns.
+    unsigned cols(void) const        //! Get number of columns.

         { return col; }
     matrix add(const matrix & other) const;
     matrix sub(const matrix & other) const;
         { return col; }
     matrix add(const matrix & other) const;
     matrix sub(const matrix & other) const;


@@ -92,7 +87,7 @@ public:
     const ex & operator() (unsigned ro, unsigned co) const;
     matrix & set(unsigned ro, unsigned co, ex value);
     matrix transpose(void) const;
     const ex & operator() (unsigned ro, unsigned co) const;
     matrix & set(unsigned ro, unsigned co, ex value);
     matrix transpose(void) const;

-    ex determinant(void) const;
+    ex determinant(unsigned options = determinant_algo::automatic) const;

     ex trace(void) const;
     ex charpoly(const symbol & lambda) const;
     matrix inverse(void) const;
     ex trace(void) const;
     ex charpoly(const symbol & lambda) const;
     matrix inverse(void) const;


@@ -105,7 +100,6 @@ protected:
     int division_free_elimination(void);
     int fraction_free_elimination(bool det = false);
     int pivot(unsigned ro, bool symbolic=true);
     int division_free_elimination(void);
     int fraction_free_elimination(bool det = false);
     int pivot(unsigned ro, bool symbolic=true);

-    void swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2);

     
 // member variables
 protected:
     
 // member variables
 protected:


@@ -147,8 +141,8 @@ inline unsigned cols(const matrix & m)
 inline matrix transpose(const matrix & m)
 { return m.transpose(); }
 
 inline matrix transpose(const matrix & m)
 { return m.transpose(); }
 

-inline ex determinant(const matrix & m)
-{ return m.determinant(); }
+inline ex determinant(const matrix & m, unsigned options = determinant_algo::automatic)
+{ return m.determinant(options); }

 
 inline ex trace(const matrix & m)
 { return m.trace(); }
 
 inline ex trace(const matrix & m)
 { return m.trace(); }

diff --git a/ginac/numeric.h b/ginac/numeric.h
index d55e1d134e70b719dd5b6753f6af5506d6f80640..c5086f433cac0142f542799cf39c85c9750247ea 100644 (file)
--- a/ginac/numeric.h
+++ b/ginac/numeric.h
@@ -83,7 +83,6 @@ class numeric : public basic
     friend const numeric atanh(const numeric & x);
     friend const numeric Li2(const numeric & x);
     friend const numeric zeta(const numeric & x);
     friend const numeric atanh(const numeric & x);
     friend const numeric Li2(const numeric & x);
     friend const numeric zeta(const numeric & x);

-    // friend const numeric bernoulli(const numeric & n);

     friend const numeric fibonacci(const numeric & n);
     friend numeric abs(const numeric & x);
     friend numeric mod(const numeric & a, const numeric & b);
     friend const numeric fibonacci(const numeric & n);
     friend numeric abs(const numeric & x);
     friend numeric mod(const numeric & a, const numeric & b);