m[r*n+c];
+ for (unsigned c=0; czero_in_last_row)(zero_in_this_row=n));
 zero_in_last_row = zero_in_this_row;
+ // Gather some statistical information about the augmented matrix:
+ bool numeric_flag = true;
+ for (exvector::const_iterator r=aug.m.begin(); r!=aug.m.end(); ++r) {
+ if (!(*r).info(info_flags::numeric))
+ numeric_flag = false;
}
#endif // def DO_GINAC_ASSERT
 // assemble solution
 matrix sol(n,1);
 unsigned last_assigned_sol = n+1;
 for (int r=m1; r>=0; r) {
 unsigned first_non_zero = 1;
 while ((first_non_zero<=n)&&(a(r,first_non_zero1).is_zero()))
 first_non_zero++;
 if (first_non_zero>n) {
 // row consists only of zeroes, corresponding rhs must be 0 as well
 if (!b.m[r*b.cols()].is_zero()) {
 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
 }
 } else {
 // assign solutions for vars between first_non_zero+1 and
 // last_assigned_sol1: free parameters
 for (unsigned c=first_non_zero; cm[r*col+c];
 for (unsigned c=0; c0; r) {
 for (unsigned i=r; irow==1)
 return m[0];
 if (this>row==2)
+ const unsigned n = this>cols();
+ if (n==1)
+ return m[0].expand();
+ if (n==2)
return (m[0]*m[3]m[2]*m[1]).expand();
 if (this>row==3)
+ if (n==3)
return (m[0]*m[4]*m[8]m[0]*m[5]*m[7]
m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
m[1]*m[5]*m[6]m[2]*m[4]*m[6]).expand();
@@ 875,8 +820,8 @@ ex matrix::determinant_minor(void) const
// This algorithm can best be understood by looking at a naive
// implementation of Laplaceexpansion, like this one:
// ex det;
 // matrix minorM(this>row1,this>col1);
 // for (unsigned r1=0; r1row; ++r1) {
+ // matrix minorM(this>rows()1,this>cols()1);
+ // for (unsigned r1=0; r1rows(); ++r1) {
// // shortcut if element(r1,0) vanishes
// if (m[r1*col].is_zero())
// continue;
@@ 906,10 +851,10 @@ ex matrix::determinant_minor(void) const
// Unique flipper counter for partitioning into minors
std::vector Pkey;
 Pkey.reserve(this>col);
+ Pkey.reserve(n);
// key for minor determinant (a subpartition of Pkey)
std::vector Mkey;
 Mkey.reserve(this>col1);
+ Mkey.reserve(n1);
// we store our subminors in maps, keys being the rows they arise from
typedef std::map,class ex> Rmap;
typedef std::map,class ex>::value_type Rmap_value;
@@ 917,34 +862,34 @@ ex matrix::determinant_minor(void) const
Rmap B;
ex det;
// initialize A with last column:
 for (unsigned r=0; rcol; ++r) {
+ for (unsigned r=0; rcol*r+this>col1]));
+ A.insert(Rmap_value(Pkey,m[n*(r+1)1]));
}
// proceed from right to left through matrix
 for (int c=this>col2; c>=0; c) {
+ for (int c=n2; c>=0; c) {
Pkey.erase(Pkey.begin(),Pkey.end()); // don't change capacity
Mkey.erase(Mkey.begin(),Mkey.end());
 for (unsigned i=0; icolc; ++i)
+ for (unsigned i=0; icolc; ++r) {
+ for (unsigned r=0; rcol+c].is_zero())
+ if (m[Pkey[r]*n+c].is_zero())
continue;
// create the sorted key for all possible minors
Mkey.erase(Mkey.begin(),Mkey.end());
 for (unsigned i=0; icolc; ++i)
+ for (unsigned i=0; icol+c]*A[Mkey];
+ det = m[Pkey[r]*n+c]*A[Mkey];
else
 det += m[Pkey[r]*this>col+c]*A[Mkey];
+ det += m[Pkey[r]*n+c]*A[Mkey];
}
// prevent buildup of deep nesting of expressions saves time:
det = det.expand();
@@ 952,13 +897,13 @@ ex matrix::determinant_minor(void) const
if (!det.is_zero())
B.insert(Rmap_value(Pkey,det));
// increment our strange flipper counter
 for (fc=this>colc; fc>0; fc) {
+ for (fc=nc; fc>0; fc) {
++Pkey[fc1];
if (Pkey[fc1]colc)
 for (unsigned j=fc; jcolc; ++j)
+ if (fcrows();
+ const unsigned n = this>cols();
+ GINAC_ASSERT(!det  n==m);
int sign = 1;
 ex piv;
 for (unsigned r1=0; r1 0)
 sign = sign;
 for (unsigned r2=r1+1; r2m[r2*col+r1] / this>m[r1*col+r1];
 for (unsigned c=r1+1; cm[r2*col+c] = piv * this>m[r1*col+c];
 for (unsigned c=0; c<=r1; ++c)
 this>m[r2*col+c] = _ex0();
+
+ unsigned r0 = 0;
+ for (unsigned r1=0; (r1=0) {
+ if (indx > 0)
+ sign = sign;
+ for (unsigned r2=r0+1; r2m[r2*n+r1] / this>m[r0*n+r1];
+ for (unsigned c=r1+1; cm[r2*n+c] = piv * this>m[r0*n+c];
+ if (!this>m[r2*n+c].info(info_flags::numeric))
+ this>m[r2*n+c] = this>m[r2*n+c].normal();
+ }
+ // fill up left hand side with zeros
+ for (unsigned c=0; c<=r1; ++c)
+ this>m[r2*n+c] = _ex0();
+ }
+ if (det) {
+ // save space by deleting no longer needed elements
+ for (unsigned c=r0+1; cm[r0*n+c] = _ex0();
+ }
+ ++r0;
}
}
@@ 999,26 +967,46 @@ int matrix::gauss_elimination(void)
}
/** Perform the steps of division free elimination to bring the matrix
+/** Perform the steps of division free elimination to bring the m x n matrix
* into an upper echelon form.
*
+ * @param det may be set to true to save a lot of space if one is only
+ * interested in the diagonal elements (i.e. for calculating determinants).
+ * The others are set to zero in this case.
* @return sign is 1 if an even number of rows was swapped, 1 if an odd
* number of rows was swapped and 0 if the matrix is singular. */
int matrix::division_free_elimination(void)
+int matrix::division_free_elimination(const bool det)
{
 int sign = 1;
ensure_if_modifiable();
 for (unsigned r1=0; r10)
 sign = sign;
 for (unsigned r2=r1+1; r2m[r2*col+c] = this>m[r1*col+r1]*this>m[r2*col+c]  this>m[r2*col+r1]*this>m[r1*col+c];
 for (unsigned c=0; c<=r1; ++c)
 this>m[r2*col+c] = _ex0();
+ const unsigned m = this>rows();
+ const unsigned n = this>cols();
+ GINAC_ASSERT(!det  n==m);
+ int sign = 1;
+
+ unsigned r0 = 0;
+ for (unsigned r1=0; (r1=0) {
+ if (indx>0)
+ sign = sign;
+ for (unsigned r2=r0+1; r2m[r2*n+c] = (this>m[r0*n+r1]*this>m[r2*n+c]  this>m[r2*n+r1]*this>m[r0*n+c]).expand();
+ // fill up left hand side with zeros
+ for (unsigned c=0; c<=r1; ++c)
+ this>m[r2*n+c] = _ex0();
+ }
+ if (det) {
+ // save space by deleting no longer needed elements
+ for (unsigned c=r0+1; cm[r0*n+c] = _ex0();
+ }
+ ++r0;
}
}
@@ 1032,11 +1020,11 @@ int matrix::division_free_elimination(void)
* is possible, since we know the divisor at each step.
*
* @param det may be set to true to save a lot of space if one is only
 * interested in the last element (i.e. for calculating determinants), the
+ * interested in the last element (i.e. for calculating determinants). The
* others are set to zero in this case.
* @return sign is 1 if an even number of rows was swapped, 1 if an odd
* number of rows was swapped and 0 if the matrix is singular. */
int matrix::fraction_free_elimination(bool det)
+int matrix::fraction_free_elimination(const bool det)
{
// Method:
// (singlestep fraction free elimination scheme, already known to Jordan)
@@ 1062,12 +1050,13 @@ int matrix::fraction_free_elimination(bool det)
// and D{m[k+1](r,c)} by
// D{m[k1](k1,k1)}.
 GINAC_ASSERT(!det  row==col);
ensure_if_modifiable();
 if (rows()==1)
 return 1;

+ const unsigned m = this>rows();
+ const unsigned n = this>cols();
+ GINAC_ASSERT(!det  n==m);
int sign = 1;
+ if (m==1)
+ return 1;
ex divisor_n = 1;
ex divisor_d = 1;
ex dividend_n;
@@ 1081,62 +1070,70 @@ int matrix::fraction_free_elimination(bool det)
// need GCDs) since the elements of *this might be unnormalized, which
// makes things more complicated than they need to be.
matrix tmp_n(*this);
 matrix tmp_d(row,col); // for denominators, if needed
+ matrix tmp_d(m,n); // for denominators, if needed
lst srl; // symbol replacement list
 exvector::iterator it = m.begin();
+ exvector::iterator it = this>m.begin();
exvector::iterator tmp_n_it = tmp_n.m.begin();
exvector::iterator tmp_d_it = tmp_d.m.begin();
 for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it) {
+ for (; it!= this>m.end(); ++it, ++tmp_n_it, ++tmp_d_it) {
(*tmp_n_it) = (*it).normal().to_rational(srl);
(*tmp_d_it) = (*tmp_n_it).denom();
(*tmp_n_it) = (*tmp_n_it).numer();
}
 for (unsigned r1=0; r10) {
 sign = sign;
 // rows r1 and indx were swapped, so pivot matrix tmp_d:
 for (unsigned c=0; c0) {
 divisor_n = tmp_n.m[(r11)*col+(r11)].expand();
 divisor_d = tmp_d.m[(r11)*col+(r11)].expand();
 // save space by deleting no longer needed elements:
 if (det) {
 for (unsigned c=0; c=0) {
+ if (indx>0) {
+ sign = sign;
+ // tmp_n's rows r0 and indx were swapped, do the same in tmp_d:
+ for (unsigned c=r1; cm.begin();
tmp_n_it = tmp_n.m.begin();
tmp_d_it = tmp_d.m.begin();
 for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
+ for (; it!= this>m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
(*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl);
return sign;
@@ 1149,68 +1146,72 @@ int matrix::fraction_free_elimination(bool det)
* where the element was found. With (symbolic==true) it does the same thing
* with the first nonzero element.
*
 * @param ro is the row to be inspected
+ * @param ro is the row from where to begin
+ * @param co is the column to be inspected
* @param symbolic signal if we want the first nonzero element to be pivoted
* (true) or the one with the largest absolute value (false).
* @return 0 if no interchange occured, 1 if all are zero (usually signaling
* a degeneracy) and positive integer k means that rows ro and k were swapped.
*/
int matrix::pivot(unsigned ro, bool symbolic)
+int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
{
unsigned k = ro;

 if (symbolic) { // search first nonzero
 for (unsigned r=ro; r maxn &&
 !tmp.is_zero()) {
 maxn = tmp;
 k = r;
+ if (symbolic) {
+ // search first nonzero element in column co beginning at row ro
+ while ((km[k*col+co].expand().is_zero()))
+ ++k;
+ } else {
+ // search largest element in column co beginning at row ro
+ GINAC_ASSERT(is_ex_of_type(this>m[k*col+co],numeric));
+ unsigned kmax = k+1;
+ numeric mmax = abs(ex_to_numeric(m[kmax*col+co]));
+ while (kmaxm[kmax*col+co],numeric));
+ numeric tmp = ex_to_numeric(this>m[kmax*col+co]);
+ if (abs(tmp) > mmax) {
+ mmax = tmp;
+ k = kmax;
}
+ ++kmax;
}
+ if (!mmax.is_zero())
+ k = kmax;
}
 if (m[k*col+ro].is_zero())
+ if (k==row)
+ // all elements in column co below row ro vanish
return 1;
 if (k!=ro) { // swap rows
 ensure_if_modifiable();
 for (unsigned c=0; c cols)
 cols = l.op(i).nops();

 // Allocate and fill matrix
 matrix &m = *new matrix(rows, cols);
 for (i=0; i j)
 m.set(i, j, l.op(i).op(j));
 else
 m.set(i, j, ex(0));
 return m;
+ if (!is_ex_of_type(l, lst))
+ throw(std::invalid_argument("argument to lst_to_matrix() must be a lst"));
+
+ // Find number of rows and columns
+ unsigned rows = l.nops(), cols = 0, i, j;
+ for (i=0; i cols)
+ cols = l.op(i).nops();
+
+ // Allocate and fill matrix
+ matrix &m = *new matrix(rows, cols);
+ for (i=0; i j)
+ m.set(i, j, l.op(i).op(j));
+ else
+ m.set(i, j, ex(0));
+ return m;
}
//////////
diff git a/ginac/matrix.h b/ginac/matrix.h
index b737a38c..baaa8a5b 100644
 a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ 87,19 +87,18 @@ public:
const ex & operator() (unsigned ro, unsigned co) const;
matrix & set(unsigned ro, unsigned co, ex value);
matrix transpose(void) const;
 ex determinant(unsigned options = determinant_algo::automatic) const;
+ ex determinant(unsigned algo = determinant_algo::automatic) const;
ex trace(void) const;
ex charpoly(const symbol & lambda) const;
matrix inverse(void) const;
 matrix fraction_free_elim(const matrix & vars, const matrix & v) const;
 matrix solve(const matrix & vars, const matrix & rhs) const;
 matrix old_solve(const matrix & v) const; // FIXME: may be removed
+ matrix solve(const matrix & vars, const matrix & rhs,
+ unsigned algo = solve_algo::automatic) const;
protected:
ex determinant_minor(void) const;
 int gauss_elimination(void);
 int division_free_elimination(void);
 int fraction_free_elimination(bool det = false);
 int pivot(unsigned ro, bool symbolic=true);
+ int gauss_elimination(const bool det = false);
+ int division_free_elimination(const bool det = false);
+ int fraction_free_elimination(const bool det = false);
+ int pivot(unsigned ro, unsigned co, bool symbolic = true);
// member variables
protected:
diff git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 0df296eb..bcf140d6 100644
 a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ 153,7 +153,8 @@ numeric::numeric(int i) : basic(TINFO_numeric)
// emphasizes efficiency:
value = new ::cl_I((long) i);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 166,7 +167,8 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
// emphasizes efficiency:
value = new ::cl_I((unsigned long)i);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 176,7 +178,8 @@ numeric::numeric(long i) : basic(TINFO_numeric)
debugmsg("numeric constructor from long",LOGLEVEL_CONSTRUCT);
value = new ::cl_I(i);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 186,7 +189,8 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric)
debugmsg("numeric constructor from ulong",LOGLEVEL_CONSTRUCT);
value = new ::cl_I(i);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 201,7 +205,8 @@ numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
value = new ::cl_I(numer);
*value = *value / ::cl_I(denom);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 215,7 +220,8 @@ numeric::numeric(double d) : basic(TINFO_numeric)
value = new cl_N;
*value = cl_float(d, cl_default_float_format);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 286,7 +292,8 @@ numeric::numeric(const cl_N & z) : basic(TINFO_numeric)
debugmsg("numeric constructor from cl_N", LOGLEVEL_CONSTRUCT);
value = new ::cl_N(z);
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
@@ 329,7 +336,8 @@ numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l
}
}
calchash();
 setflag(status_flags::evaluated
+ setflag(status_flags::evaluated 
+ status_flags::expanded 
status_flags::hash_calculated);
}
diff git a/ginac/power.cpp b/ginac/power.cpp
index fc52a5f8..53cafc52 100644
 a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ 525,11 +525,10 @@ ex power::derivative(const symbol & s) const
{
if (exponent.info(info_flags::real)) {
// D(b^r) = r * b^(r1) * D(b) (faster than the formula below)
 //return mul(mul(exponent, power(basis, exponent  _ex1())), basis.diff(s));
epvector newseq;
newseq.reserve(2);
newseq.push_back(expair(basis, exponent  _ex1()));
 newseq.push_back(expair(basis.diff(s),_ex1()));
+ newseq.push_back(expair(basis.diff(s), _ex1()));
return mul(newseq, exponent);
} else {
// D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
diff git a/ginac/version.h b/ginac/version.h
index 35801017..c49a8dcb 100644
 a/ginac/version.h
+++ b/ginac/version.h
@@ 26,6 +26,6 @@
/* Major, minor, and micro version number of the GiNaC library. */
#define GINACLIB_MAJOR_VERSION 0
#define GINACLIB_MINOR_VERSION 6
#define GINACLIB_MICRO_VERSION 3
+#define GINACLIB_MICRO_VERSION 4
#endif // ndef __GINAC_VERSION_H__

2.20.1