- ...normal.cpp.
- included new method matrix::determinant_bareiss, which doesn't rely on GCDs.
- log_series now honors the branch cut.
/** 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;
}
// 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
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);
// 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 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).
- evalf_func(csgn_evalf));
+ evalf_func(csgn_evalf).
+ series_func(csgn_series));
//////////
// dilogarithm
* 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)) {
// 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) ==
- // 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()
- // 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
}
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(_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(_ex0())) // log(0) -> infinity
- throw(std::domain_error("log_eval(): log(0)"));
// log(float)
if (!x.info(info_flags::crational))
return log_evalf(x);
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());
}
-/*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())
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;
- // 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));
- 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).
#include "debugmsg.h"
#include "power.h"
#include "symbol.h"
+#include "normal.h"
#ifndef NO_NAMESPACE_GINAC
namespace GiNaC {
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
}
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();
}
// 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();
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"));
ex det = _ex1();
ex piv;
+ // standard Gauss method:
for (unsigned r1=0; r1<row; ++r1) {
int indx = tmp.pivot(r1);
if (indx == -1)
}
-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")
*
* @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)
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().
- * 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];
- 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];
}
* 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;
+ 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)
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 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();
}
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);
* @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)) {
* @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
*
* @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());
*
* @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)) {
*
* @param e expanded polynomial
* @return integer content */
-
numeric ex::integer_content(void) const
{
GINAC_ASSERT(bp!=0);
* @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())
* @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())
* @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())
* 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"));
- 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;
* @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())
* @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";
* @param e expanded multivariate polynomial
* @return maximum coefficient
* @see heur_gcd */
-
numeric ex::max_coefficient(void) const
{
GINAC_ASSERT(bp!=0);
* @param xi modulus
* @return mapped polynomial
* @see heur_gcd */
-
ex ex::smod(const numeric &xi) const
{
GINAC_ASSERT(bp!=0);
* @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";
* @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";
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
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
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);
}
+/** 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,
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);
}