- expairseq.cpp: moved expairseq::to_rational to...

- ...normal.cpp.
- included new method matrix::determinant_bareiss, which doesn't rely on GCDs.
- log_series now honors the branch cut.

diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index 5606f3640e7d4d0fb9952ebefc830ec79478db3c..2a4a3e39b9e3629b79b9e2e3cfd7a5307428eaa6 100644 (file)
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -205,6 +205,9 @@ bool basic::info(unsigned inf) const
 /** Number of operands/members. */
 unsigned basic::nops() const
 {
 /** Number of operands/members. */
 unsigned basic::nops() const
 {

+    // iterating from 0 to nops() on atomic objects should be an empty loop,
+    // and accessing their elements is a range error.  Container objects should
+    // override this.

     return 0;
 }
 
     return 0;
 }
 

diff --git a/ginac/ex.h b/ginac/ex.h
index bf8d0d6f45edb04246c8f4594b7d4d35620253b0..b5e511e6e51f9330904e63908456148d77ec48b8 100644 (file)
--- a/ginac/ex.h
+++ b/ginac/ex.h
@@ -42,8 +42,6 @@ class lst;
 // are defined in utils.h and not visible from outside.
 extern const ex & _ex0(void);     //  single ex(numeric(0))
 
 // are defined in utils.h and not visible from outside.
 extern const ex & _ex0(void);     //  single ex(numeric(0))
 

-// typedef vector<ex> exvector;
-

 #define INLINE_EX_CONSTRUCTORS
 
 /** Lightweight wrapper for GiNaC's symbolic objects.  Basically all it does is
 #define INLINE_EX_CONSTRUCTORS
 
 /** Lightweight wrapper for GiNaC's symbolic objects.  Basically all it does is

diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index d4e13188c21c5cff02bd8b1ae2d66e370d49e165..bdfc2ff73f97f814bb4b1e87422a60f9645c67e5 100644 (file)
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -361,17 +361,6 @@ ex expairseq::normal(lst &sym_lst, lst &repl_lst, int level) const
     return n.bp->basic::normal(sym_lst,repl_lst,level);
 }
 
     return n.bp->basic::normal(sym_lst,repl_lst,level);
 }
 

-ex expairseq::to_rational(lst &repl_lst) const
-{
-    epvector s;
-    s.reserve(seq.size());
-    for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
-        s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst),
-                                                  (*it).coeff));
-    }
-    return thisexpairseq(s, overall_coeff);
-}
-

 ex expairseq::subs(const lst & ls, const lst & lr) const
 {
     epvector * vp=subschildren(ls,lr);
 ex expairseq::subs(const lst & ls, const lst & lr) const
 {
     epvector * vp=subschildren(ls,lr);


@@ -563,7 +552,7 @@ ex expairseq::expand(unsigned options) const
 
 // protected
 
 
 // protected
 

-ex expairseq::thisexpairseq(const epvector & v,const ex & oc) const
+ex expairseq::thisexpairseq(const epvector & v, const ex & oc) const

 {
     return expairseq(v,oc);
 }
 {
     return expairseq(v,oc);
 }

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index c03f4a15a6e29b39d4dea7d5269711fc5ce5f3a4..f7d2864ef9ae64035eae269d140b6c36b6d41169 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -107,8 +107,24 @@ static ex csgn_eval(const ex & x)
     return csgn(x).hold();
 }
 
     return csgn(x).hold();
 }
 

+static ex csgn_series(const ex & x, const relational & rel, int order)
+{
+    const ex x_pt = x.subs(rel);
+    if (x_pt.info(info_flags::numeric)) {
+        if (ex_to_numeric(x_pt).real().is_zero())
+            throw (std::domain_error("csgn_series(): on imaginary axis"));
+        epvector seq;
+        seq.push_back(expair(csgn(x_pt), _ex0()));
+        return pseries(rel,seq);
+    }
+    epvector seq;
+    seq.push_back(expair(csgn(x_pt), _ex0()));
+    return pseries(rel,seq);
+}
+

 REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
 REGISTER_FUNCTION(csgn, eval_func(csgn_eval).

-                        evalf_func(csgn_evalf));
+                        evalf_func(csgn_evalf).
+                        series_func(csgn_series));

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

diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index 85275f62aabec258d301d794bcc764a36f92ce21..213d9d7ff04ea70fd902b39e9bd4b4a282256209 100644 (file)
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -56,7 +56,7 @@ static ex lgamma_evalf(const ex & x)
  *  Knows about integer arguments and that's it.  Somebody ought to provide
  *  some good numerical evaluation some day...
  *
  *  Knows about integer arguments and that's it.  Somebody ought to provide
  *  some good numerical evaluation some day...
  *

- *  @exception std::domain_error("lgamma_eval(): simple pole") */
+ *  @exception std::domain_error("lgamma_eval(): logarithmic pole") */

 static ex lgamma_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
 static ex lgamma_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {


@@ -90,20 +90,20 @@ static ex lgamma_series(const ex & x, const relational & rel, int order)
     // method:
     // Taylor series where there is no pole falls back to psi function
     // evaluation.
     // method:
     // Taylor series where there is no pole falls back to psi function
     // evaluation.

-    // On a pole at -m use the recurrence relation
+    // On a pole at -m we could use the recurrence relation

     //   lgamma(x) == lgamma(x+1)-log(x)
     // from which follows
     //   series(lgamma(x),x==-m,order) ==
     //   lgamma(x) == lgamma(x+1)-log(x)
     // from which follows
     //   series(lgamma(x),x==-m,order) ==

-    //   series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order+1);
+    //   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 x_pt = x.subs(rel);
     if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
     const ex x_pt = x.subs(rel);
     if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()

-    // if we got here we have to care for a simple pole at -m:
-    numeric m = -ex_to_numeric(x_pt);
-    ex ser_sub = _ex0();
-    for (numeric p; p<=m; ++p)
-        ser_sub += log(x+p);
-    return (lgamma(x+m+_ex1())-ser_sub).series(rel, order);
+    // if we got here we have to care for a simple pole of tgamma(-m):
+    throw (std::domain_error("lgamma_series: please implemnt my at the poles"));
+    return _ex0();  // not reached

 }
 
 
 }
 
 

diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index 924d195e7f44187eb8c6e427fa4e65166bfc46b3..81f92a93aca908aa6f698e79dd1a4380d5998d61 100644 (file)
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -109,16 +109,16 @@ static ex log_evalf(const ex & x)
 static ex log_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
 static ex log_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {

+        if (x.is_equal(_ex0()))  // log(0) -> infinity
+            throw(std::domain_error("log_eval(): log(0)"));
+        if (x.info(info_flags::real) && x.info(info_flags::negative))
+            return (log(-x)+I*Pi);

         if (x.is_equal(_ex1()))  // log(1) -> 0
             return _ex0();
         if (x.is_equal(_ex1()))  // log(1) -> 0
             return _ex0();

-        if (x.is_equal(_ex_1())) // log(-1) -> I*Pi
-            return (I*Pi);        

         if (x.is_equal(I))       // log(I) -> Pi*I/2
             return (Pi*I*_num1_2());
         if (x.is_equal(-I))      // log(-I) -> -Pi*I/2
             return (Pi*I*_num_1_2());
         if (x.is_equal(I))       // log(I) -> Pi*I/2
             return (Pi*I*_num1_2());
         if (x.is_equal(-I))      // log(-I) -> -Pi*I/2
             return (Pi*I*_num_1_2());

-        if (x.is_equal(_ex0()))  // log(0) -> infinity
-            throw(std::domain_error("log_eval(): log(0)"));

         // log(float)
         if (!x.info(info_flags::crational))
             return log_evalf(x);
         // log(float)
         if (!x.info(info_flags::crational))
             return log_evalf(x);


@@ -139,12 +139,12 @@ static ex log_eval(const ex & x)
 static ex log_deriv(const ex & x, unsigned deriv_param)
 {
     GINAC_ASSERT(deriv_param==0);
 static ex log_deriv(const ex & x, unsigned deriv_param)
 {
     GINAC_ASSERT(deriv_param==0);

-
+    

     // d/dx log(x) -> 1/x
     return power(x, _ex_1());
 }
 
     // d/dx log(x) -> 1/x
     return power(x, _ex_1());
 }
 

-/*static ex log_series(const ex &x, const relational &rel, int order)
+static ex log_series(const ex &x, const relational &rel, int order)

 {
     const ex x_pt = x.subs(rel);
     if (!x_pt.info(info_flags::negative) && !x_pt.is_zero())
 {
     const ex x_pt = x.subs(rel);
     if (!x_pt.info(info_flags::negative) && !x_pt.is_zero())


@@ -154,27 +154,16 @@ static ex log_deriv(const ex & x, unsigned deriv_param)
         epvector seq;
         seq.push_back(expair(log(x), _ex0()));
         return pseries(rel, seq);
         epvector seq;
         seq.push_back(expair(log(x), _ex0()));
         return pseries(rel, seq);

-    } // the branch cut:
+    } // on the branch cut:

     const ex point = rel.rhs();
     const symbol *s = static_cast<symbol *>(rel.lhs().bp);
     const symbol foo;
     // compute the formal series:
     ex replx = series(log(x),*s==foo,order).subs(foo==point);
     epvector seq;
     const ex point = rel.rhs();
     const symbol *s = static_cast<symbol *>(rel.lhs().bp);
     const symbol foo;
     // compute the formal series:
     ex replx = series(log(x),*s==foo,order).subs(foo==point);
     epvector seq;

-    // FIXME: this is probably off by 2 or so:

     seq.push_back(expair(-I*csgn(x*I)*Pi,_ex0()));
     seq.push_back(expair(Order(_ex1()),order));
     seq.push_back(expair(-I*csgn(x*I)*Pi,_ex0()));
     seq.push_back(expair(Order(_ex1()),order));

-    return series(replx + pseries(rel, seq),rel,order);
-}*/
-
-static ex log_series(const ex &x, const relational &r, int order)
-{
-       if (x.subs(r).is_zero()) {
-               epvector seq;
-               seq.push_back(expair(log(x), _ex0()));
-               return pseries(r, seq);
-       } else
-               throw do_taylor();
+    return series(replx - I*Pi + pseries(rel, seq),rel,order);

 }
 
 REGISTER_FUNCTION(log, eval_func(log_eval).
 }
 
 REGISTER_FUNCTION(log, eval_func(log_eval).

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 0296c18340df9f16b746298df6c61577ba5c050b..f76aa9d768f0de53419c0285ebf9fc9485226596 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -32,6 +32,7 @@
 #include "debugmsg.h"
 #include "power.h"
 #include "symbol.h"
 #include "debugmsg.h"
 #include "power.h"
 #include "symbol.h"

+#include "normal.h"

 
 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
 
 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {


@@ -475,15 +476,8 @@ ex matrix::determinant(void) const
         return determinant_numeric();
     
     // Does anybody really know when a matrix is sparse?
         return determinant_numeric();
     
     // Does anybody really know when a matrix is sparse?

-    if (4*sparse_count<row*col) {  // < row/2 non-zero elements average in row
-        matrix M(*this);
-        int sign = M.fraction_free_elimination();
-        if (!sign)
-            return _ex0();
-        if (normal_flag)
-            return sign * M(row-1,col-1).normal();
-        else
-            return sign * M(row-1,col-1).expand();
+    if (4*sparse_count<row*col) {  // < row/2 nonzero elements average in a row
+        return determinant_bareiss(normal_flag);

     }
     
     // Now come the minor expansion schemes.  We always develop such that the
     }
     
     // Now come the minor expansion schemes.  We always develop such that the


@@ -516,8 +510,8 @@ ex matrix::determinant(void) const
     }
     
     if (normal_flag)
     }
     
     if (normal_flag)

-        return sign*matrix(row,col,result).determinant_minor_dense().normal();
-    return sign*matrix(row,col,result).determinant_minor_dense();
+        return sign*matrix(row,col,result).determinant_minor().normal();
+    return sign*matrix(row,col,result).determinant_minor();

 }
 
 
 }
 
 


@@ -570,8 +564,7 @@ ex matrix::charpoly(const symbol & lambda) const
     
     // The pure numeric case is traditionally rather common.  Hence, it is
     // trapped and we use Leverrier's algorithm which goes as row^3 for
     
     // The pure numeric case is traditionally rather common.  Hence, it is
     // trapped and we use Leverrier's algorithm which goes as row^3 for

-    // every coefficient.  The expensive section is the matrix multiplication,
-    // maybe this can be sped up even more?
+    // every coefficient.  The expensive section is the matrix multiplication.

     if (numeric_flag) {
         matrix B(*this);
         ex c = B.trace();
     if (numeric_flag) {
         matrix B(*this);
         ex c = B.trace();


@@ -676,7 +669,7 @@ ex matrix::ffe_get(unsigned r, unsigned c) const
 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: implement a Sasaki-Murao scheme which avoids division at all!
+    // FIXME: use implementation of matrix::fraction_free_elim

     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"));
     


@@ -883,6 +876,7 @@ ex matrix::determinant_numeric(void) const
     ex det = _ex1();
     ex piv;
     
     ex det = _ex1();
     ex piv;
     

+    // standard Gauss method:

     for (unsigned r1=0; r1<row; ++r1) {
         int indx = tmp.pivot(r1);
         if (indx == -1)
     for (unsigned r1=0; r1<row; ++r1) {
         int indx = tmp.pivot(r1);
         if (indx == -1)


@@ -902,42 +896,6 @@ ex matrix::determinant_numeric(void) const
 }
 
 
 }
 
 

-ex matrix::determinant_minor_sparse(void) const
-{
-    // for small matrices the algorithm does not make any sense:
-    if (this->row==1)
-        return m[0];
-    if (this->row==2)
-        return (m[0]*m[3]-m[2]*m[1]).expand();
-    if (this->row==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();
-    
-    ex det;
-    matrix minorM(this->row-1,this->col-1);
-    for (unsigned r1=0; r1<this->row; ++r1) {
-        // shortcut if element(r1,0) vanishes
-        if (m[r1*col].is_zero())
-            continue;
-        // assemble the minor matrix
-        for (unsigned r=0; r<minorM.rows(); ++r) {
-            for (unsigned c=0; c<minorM.cols(); ++c) {
-                if (r<r1)
-                    minorM.set(r,c,m[r*col+c+1]);
-                else
-                    minorM.set(r,c,m[(r+1)*col+c+1]);
-            }
-        }
-        // recurse down and care for sign:
-        if (r1%2)
-            det -= m[r1*col] * minorM.determinant_minor_sparse();
-        else
-            det += m[r1*col] * minorM.determinant_minor_sparse();
-    }
-    return det.expand();
-}
-

 /** Recursive determinant for small matrices having at least one symbolic
  *  entry.  The basic algorithm, known as Laplace-expansion, is enhanced by
  *  some bookkeeping to avoid calculation of the same submatrices ("minors")
 /** Recursive determinant for small matrices having at least one symbolic
  *  entry.  The basic algorithm, known as Laplace-expansion, is enhanced by
  *  some bookkeeping to avoid calculation of the same submatrices ("minors")


@@ -948,7 +906,7 @@ ex matrix::determinant_minor_sparse(void) const
  *
  *  @return the determinant as a new expression (in expanded form)
  *  @see matrix::determinant() */
  *
  *  @return the determinant as a new expression (in expanded form)
  *  @see matrix::determinant() */

-ex matrix::determinant_minor_dense(void) const
+ex matrix::determinant_minor(void) const

 {
     // for small matrices the algorithm does not make any sense:
     if (this->row==1)
 {
     // for small matrices the algorithm does not make any sense:
     if (this->row==1)


@@ -1057,29 +1015,110 @@ ex matrix::determinant_minor_dense(void) const
     return det;
 }
 
     return det;
 }
 

+/** Helper function to divide rational functions, as needed in any Bareiss
+ *  elimination scheme over a quotient field.
+ *  
+ *  @see divide() */
+bool rat_divide(const ex & a, const ex & b, ex & q, bool check_args = true)
+{
+    q = _ex0();
+    if (b.is_zero())
+        throw(std::overflow_error("rat_divide(): division by zero"));
+    if (a.is_zero())
+        return true;
+    if (is_ex_exactly_of_type(b, numeric)) {
+        q = a / b;
+        return true;
+    } else if (is_ex_exactly_of_type(a, numeric))
+        return false;
+    ex a_n = a.numer();
+    ex a_d = a.denom();
+    ex b_n = b.numer();
+    ex b_d = b.denom();
+    ex n;  // new numerator
+    ex d;  // new denominator
+    bool check = true;
+    check &= divide(a_n, b_n, n, check_args);
+    check &= divide(a_d, b_d, d, check_args);
+    q = n/d;
+    return check;
+}
+

 
 

-/** Determinant built by application of the full permutation group.  This
+/** Determinant computed by using fraction free elimination.  This

  *  routine is only called internally by matrix::determinant().
  *  routine is only called internally by matrix::determinant().

- *  NOTE: it is probably inefficient in all cases and may be eliminated. */
-ex matrix::determinant_perm(void) const
+ *
+ *  @param normalize may be set to false only in integral domains. */
+ex matrix::determinant_bareiss(bool normalize) const

 {
 {

-    if (rows()==1)  // speed things up
+    if (rows()==1)

         return m[0];
     
         return m[0];
     

-    ex det;
-    ex term;
-    vector<unsigned> sigma(col);
-    for (unsigned i=0; i<col; ++i)
-        sigma[i]=i;
+    int sign = 1;
+    ex divisor = 1;
+    ex dividend;

     
     

-    do {
-        term = (*this)(sigma[0],0);
-        for (unsigned i=1; i<col; ++i)
-            term *= (*this)(sigma[i],i);
-        det += permutation_sign(sigma)*term;
-    } while (next_permutation(sigma.begin(), sigma.end()));
+    // we populate a tmp matrix to subsequently operate on, it should
+    // be normalized even though this algorithm doesn't need GCDs since
+    // the elements of *this might be unnormalized, which complicates
+    // things:
+    matrix tmp(*this);
+    exvector::const_iterator i = m.begin();
+    exvector::iterator ti = tmp.m.begin();
+    for (; i!= m.end(); ++i, ++ti) {
+        if (normalize)
+            (*ti) = (*i).normal();
+        else
+            (*ti) = (*i);
+    }

     
     

-    return det;
+    for (unsigned r1=0; r1<row-1; ++r1) {
+        int indx = tmp.pivot(r1);
+        if (indx==-1)
+            return _ex0();
+        if (indx>0)
+            sign = -sign;
+        if (r1>0) {
+            divisor = tmp.m[(r1-1)*col+(r1-1)].expand();
+            // delete the elements we don't need anymore:
+            for (unsigned c=0; c<col; ++c)
+                tmp.m[(r1-1)*col+c] = _ex0();
+        }
+        for (unsigned r2=r1+1; r2<row; ++r2) {
+            for (unsigned c=r1+1; c<col; ++c) {
+                lst srl;  // symbol replacement list for .to_rational()
+                dividend = (tmp.m[r1*tmp.col+r1]*tmp.m[r2*tmp.col+c]
+                           -tmp.m[r2*tmp.col+r1]*tmp.m[r1*tmp.col+c]).expand();
+                if (normalize) {
+#ifdef DO_GINAC_ASSERT
+                    GINAC_ASSERT(rat_divide(dividend.to_rational(srl),
+                                            divisor.to_rational(srl),
+                                            tmp.m[r2*tmp.col+c],true));
+#else
+                    rat_divide(dividend.to_rational(srl),
+                               divisor.to_rational(srl),
+                               tmp.m[r2*tmp.col+c],false);
+#endif                    
+                }
+                else {
+#ifdef DO_GINAC_ASSERT
+                    GINAC_ASSERT(divide(dividend.to_rational(srl),
+                                        divisor.to_rational(srl),
+                                        tmp.m[r2*tmp.col+c],true));
+#else
+                    divide(dividend.to_rational(srl),
+                           divisor.to_rational(srl),
+                           tmp.m[r2*tmp.col+c],false);
+#endif                    
+                }
+                tmp.m[r2*tmp.col+c] = tmp.m[r2*tmp.col+c].subs(srl);
+            }
+            for (unsigned c=0; c<=r1; ++c)
+                tmp.m[r2*tmp.col+c] = _ex0();
+        }
+    }
+    
+    return sign*tmp.m[tmp.row*tmp.col-1];

 }
 
 
 }
 
 


@@ -1144,10 +1183,19 @@ int matrix::division_free_elimination(void)
  *  number of rows was swapped and 0 if the matrix is singular. */
 int matrix::fraction_free_elimination(void)
 {
  *  number of rows was swapped and 0 if the matrix is singular. */
 int matrix::fraction_free_elimination(void)
 {

+    ensure_if_modifiable();
+    
+    // first normal all elements:
+    for (exvector::iterator i=m.begin(); i!=m.end(); ++i)
+        (*i) = (*i).normal();
+    
+    // FIXME: this is unfinished, once matrix::determinant_bareiss is
+    // bulletproof, some code ought to be copy from there to here.

     int sign = 1;
     ex divisor = 1;
     int sign = 1;
     ex divisor = 1;

+    ex dividend;
+    lst srl;      // symbol replacement list for .to_rational()

     
     

-    ensure_if_modifiable();

     for (unsigned r1=0; r1<row-1; ++r1) {
         int indx = pivot(r1);
         if (indx==-1)
     for (unsigned r1=0; r1<row-1; ++r1) {
         int indx = pivot(r1);
         if (indx==-1)


@@ -1155,10 +1203,22 @@ int matrix::fraction_free_elimination(void)
         if (indx>0)
             sign = -sign;
         if (r1>0)
         if (indx>0)
             sign = -sign;
         if (r1>0)

-            divisor = this->m[(r1-1)*col + (r1-1)];
+            divisor = this->m[(r1-1)*col+(r1-1)].expand();

         for (unsigned r2=r1+1; r2<row; ++r2) {
         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])/divisor).normal();
+            for (unsigned c=r1+1; c<col; ++c) {
+                dividend = (this->m[r1*col+r1]*this->m[r2*col+c]
+                           -this->m[r2*col+r1]*this->m[r1*col+c]).expand();
+#ifdef DO_GINAC_ASSERT
+                GINAC_ASSERT(divide(dividend.to_rational(srl),
+                                    divisor.to_rational(srl),
+                                    this->m[r2*col+c]));
+#else
+                divide(dividend.to_rational(srl),
+                       divisor.to_rational(srl),
+                       this->m[r2*col+c]);
+#endif // DO_GINAC_ASSERT
+                this->m[r2*col+c] = this->m[r2*col+c].subs(srl);
+            }

             for (unsigned c=0; c<=r1; ++c)
                 this->m[r2*col+c] = _ex0();
         }
             for (unsigned c=0; c<=r1; ++c)
                 this->m[r2*col+c] = _ex0();
         }

diff --git a/ginac/matrix.h b/ginac/matrix.h
index 65eebde91697a504f7c03f72e205cbe32cd0272c..a36276f4d332b06f57afdc6747e7a5f54697cf72 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -96,9 +96,8 @@ public:
     matrix old_solve(const matrix & v) const;  // FIXME: may be removed
 protected:
     ex determinant_numeric(void) const;
     matrix old_solve(const matrix & v) const;  // FIXME: may be removed
 protected:
     ex determinant_numeric(void) const;

-    ex determinant_minor_sparse(void) const;
-    ex determinant_minor_dense(void) const;
-    ex determinant_perm(void) const;
+    ex determinant_minor(void) const;
+    ex determinant_bareiss(bool normalize=true) const;

     int gauss_elimination(void);
     int fraction_free_elimination(void);
     int division_free_elimination(void);
     int gauss_elimination(void);
     int fraction_free_elimination(void);
     int division_free_elimination(void);

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 9e24b99bc12e68bb2a5cdcbe263d4637bd9c64ea..cdaf8367f7cdfe68a98ff44e594f96ee03fb57bd 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -93,7 +93,6 @@ static struct _stat_print {
  *  @param e  expression to search
  *  @param x  pointer to first symbol found (returned)
  *  @return "false" if no symbol was found, "true" otherwise */
  *  @param e  expression to search
  *  @param x  pointer to first symbol found (returned)
  *  @return "false" if no symbol was found, "true" otherwise */

-

 static bool get_first_symbol(const ex &e, const symbol *&x)
 {
     if (is_ex_exactly_of_type(e, symbol)) {
 static bool get_first_symbol(const ex &e, const symbol *&x)
 {
     if (is_ex_exactly_of_type(e, symbol)) {


@@ -186,7 +185,6 @@ static void collect_symbols(const ex &e, sym_desc_vec &v)
  *  @param a  first multivariate polynomial
  *  @param b  second multivariate polynomial
  *  @param v  vector of sym_desc structs (filled in) */
  *  @param a  first multivariate polynomial
  *  @param b  second multivariate polynomial
  *  @param v  vector of sym_desc structs (filled in) */

-

 static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 {
     collect_symbols(a.eval(), v);   // eval() to expand assigned symbols
 static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 {
     collect_symbols(a.eval(), v);   // eval() to expand assigned symbols


@@ -247,7 +245,6 @@ static numeric lcmcoeff(const ex &e, const numeric &l)
  *
  *  @param e  multivariate polynomial (need not be expanded)
  *  @return LCM of denominators of coefficients */
  *
  *  @param e  multivariate polynomial (need not be expanded)
  *  @return LCM of denominators of coefficients */

-

 static numeric lcm_of_coefficients_denominators(const ex &e)
 {
     return lcmcoeff(e, _num1());
 static numeric lcm_of_coefficients_denominators(const ex &e)
 {
     return lcmcoeff(e, _num1());


@@ -258,7 +255,6 @@ static numeric lcm_of_coefficients_denominators(const ex &e)
  *
  *  @param e  multivariate polynomial (need not be expanded)
  *  @param lcm  LCM to multiply in */
  *
  *  @param e  multivariate polynomial (need not be expanded)
  *  @param lcm  LCM to multiply in */

-

 static ex multiply_lcm(const ex &e, const numeric &lcm)
 {
        if (is_ex_exactly_of_type(e, mul)) {
 static ex multiply_lcm(const ex &e, const numeric &lcm)
 {
        if (is_ex_exactly_of_type(e, mul)) {


@@ -288,7 +284,6 @@ static ex multiply_lcm(const ex &e, const numeric &lcm)
  *
  *  @param e  expanded polynomial
  *  @return integer content */
  *
  *  @param e  expanded polynomial
  *  @return integer content */

-

 numeric ex::integer_content(void) const
 {
     GINAC_ASSERT(bp!=0);
 numeric ex::integer_content(void) const
 {
     GINAC_ASSERT(bp!=0);


@@ -349,7 +344,6 @@ numeric mul::integer_content(void) const
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return quotient of a and b in Q[x] */
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return quotient of a and b in Q[x] */

-

 ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
     if (b.is_zero())
 ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
     if (b.is_zero())


@@ -400,7 +394,6 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return remainder of a(x) and b(x) in Q[x] */
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return remainder of a(x) and b(x) in Q[x] */

-

 ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
     if (b.is_zero())
 ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
     if (b.is_zero())


@@ -452,7 +445,6 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return pseudo-remainder of a(x) and b(x) in Z[x] */
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return pseudo-remainder of a(x) and b(x) in Z[x] */

-

 ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
     if (b.is_zero())
 ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
     if (b.is_zero())


@@ -506,14 +498,13 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
  *         coefficients (defaults to "true")
  *  @return "true" when exact division succeeds (quotient returned in q),
  *          "false" otherwise */
  *         coefficients (defaults to "true")
  *  @return "true" when exact division succeeds (quotient returned in q),
  *          "false" otherwise */

-

 bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 {
     q = _ex0();
     if (b.is_zero())
         throw(std::overflow_error("divide: division by zero"));
 bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 {
     q = _ex0();
     if (b.is_zero())
         throw(std::overflow_error("divide: division by zero"));

-       if (a.is_zero())
-               return true;
+    if (a.is_zero())
+        return true;

     if (is_ex_exactly_of_type(b, numeric)) {
         q = a / b;
         return true;
     if (is_ex_exactly_of_type(b, numeric)) {
         q = a / b;
         return true;


@@ -816,7 +807,6 @@ ex ex::primpart(const symbol &x) const
  *  @param x  variable in which to compute the primitive part
  *  @param c  previously computed content part
  *  @return primitive part */
  *  @param x  variable in which to compute the primitive part
  *  @param c  previously computed content part
  *  @return primitive part */

-

 ex ex::primpart(const symbol &x, const ex &c) const
 {
     if (is_zero())
 ex ex::primpart(const symbol &x, const ex &c) const
 {
     if (is_zero())


@@ -1039,7 +1029,6 @@ static ex red_gcd(const ex &a, const ex &b, const symbol *x)
  *  @param x  pointer to symbol (main variable) in which to compute the GCD in
  *  @return the GCD as a new expression
  *  @see gcd */
  *  @param x  pointer to symbol (main variable) in which to compute the GCD in
  *  @return the GCD as a new expression
  *  @see gcd */

-

 static ex sr_gcd(const ex &a, const ex &b, const symbol *x)
 {
 //clog << "sr_gcd(" << a << "," << b << ")\n";
 static ex sr_gcd(const ex &a, const ex &b, const symbol *x)
 {
 //clog << "sr_gcd(" << a << "," << b << ")\n";


@@ -1115,7 +1104,6 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x)
  *  @param e  expanded multivariate polynomial
  *  @return maximum coefficient
  *  @see heur_gcd */
  *  @param e  expanded multivariate polynomial
  *  @return maximum coefficient
  *  @see heur_gcd */

-

 numeric ex::max_coefficient(void) const
 {
     GINAC_ASSERT(bp!=0);
 numeric ex::max_coefficient(void) const
 {
     GINAC_ASSERT(bp!=0);


@@ -1171,7 +1159,6 @@ numeric mul::max_coefficient(void) const
  *  @param xi  modulus
  *  @return mapped polynomial
  *  @see heur_gcd */
  *  @param xi  modulus
  *  @return mapped polynomial
  *  @see heur_gcd */

-

 ex ex::smod(const numeric &xi) const
 {
     GINAC_ASSERT(bp!=0);
 ex ex::smod(const numeric &xi) const
 {
     GINAC_ASSERT(bp!=0);


@@ -1259,7 +1246,6 @@ class gcdheu_failed {};
  *  @return the GCD as a new expression
  *  @see gcd
  *  @exception gcdheu_failed() */
  *  @return the GCD as a new expression
  *  @see gcd
  *  @exception gcdheu_failed() */

-

 static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
 {
 //clog << "heur_gcd(" << a << "," << b << ")\n";
 static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
 {
 //clog << "heur_gcd(" << a << "," << b << ")\n";


@@ -1347,7 +1333,6 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return the GCD as a new expression */
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return the GCD as a new expression */

-

 ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
 {
 //clog << "gcd(" << a << "," << b << ")\n";
 ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
 {
 //clog << "gcd(" << a << "," << b << ")\n";


@@ -1682,7 +1667,7 @@ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
     for (unsigned i=0; i<repl_lst.nops(); i++)
         if (repl_lst.op(i).is_equal(e))
             return sym_lst.op(i);
     for (unsigned i=0; i<repl_lst.nops(); i++)
         if (repl_lst.op(i).is_equal(e))
             return sym_lst.op(i);

-
+    

     // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't contain symbols from the sym_lst
        // because subs() is not recursive
     // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't contain symbols from the sym_lst
        // because subs() is not recursive


@@ -1704,7 +1689,7 @@ static ex replace_with_symbol(const ex &e, lst &repl_lst)
     for (unsigned i=0; i<repl_lst.nops(); i++)
         if (repl_lst.op(i).op(1).is_equal(e))
             return repl_lst.op(i).op(0);
     for (unsigned i=0; i<repl_lst.nops(); i++)
         if (repl_lst.op(i).op(1).is_equal(e))
             return repl_lst.op(i).op(0);

-
+    

     // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't contain symbols from the sym_lst
        // because subs() is not recursive
     // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't contain symbols from the sym_lst
        // because subs() is not recursive


@@ -2079,7 +2064,8 @@ ex numeric::to_rational(lst &repl_lst) const
         if (!is_integer())
             return replace_with_symbol(*this, repl_lst);
     } else { // complex
         if (!is_integer())
             return replace_with_symbol(*this, repl_lst);
     } else { // complex

-        numeric re = real(), im = imag();
+        numeric re = real();
+        numeric im = imag();

         ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
         ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
         return re_ex + im_ex * replace_with_symbol(I, repl_lst);
         ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
         ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
         return re_ex + im_ex * replace_with_symbol(I, repl_lst);


@@ -2100,6 +2086,20 @@ ex power::to_rational(lst &repl_lst) const
 }
 
 
 }
 
 

+/** Implementation of ex::to_rational() for expairseqs.
+ *  @see ex::to_rational */
+ex expairseq::to_rational(lst &repl_lst) const
+{
+    epvector s;
+    s.reserve(seq.size());
+    for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
+        s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst),
+                                                  (*it).coeff));
+    }
+    return thisexpairseq(s, overall_coeff);
+}
+
+

 /** Rationalization of non-rational functions.
  *  This function converts a general expression to a rational polynomial
  *  by replacing all non-rational subexpressions (like non-rational numbers,
 /** Rationalization of non-rational functions.
  *  This function converts a general expression to a rational polynomial
  *  by replacing all non-rational subexpressions (like non-rational numbers,

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 277cef18a10253504f6eaf0e708718fba53cb50e..268041325aebd1ce11cdff7e40bcfe5871e8df0a 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -165,6 +165,9 @@ basic *pseries::duplicate() const
 void pseries::print(ostream &os, unsigned upper_precedence) const
 {
     debugmsg("pseries print", LOGLEVEL_PRINT);
 void pseries::print(ostream &os, unsigned upper_precedence) const
 {
     debugmsg("pseries print", LOGLEVEL_PRINT);

+    // This could be made better, since series expansion at x==1 might print
+    // -1+2*x+Order((-1+x)^2) instead of 1+2*(-1+x)+Order((-1+x)^2), which is
+    // correct but can be rather confusing.

     convert_to_poly().print(os, upper_precedence);
 }
 
     convert_to_poly().print(os, upper_precedence);
 }