From af0c47009ca7a15af966430bdf1a72fe05c1c6f9 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Thu, 20 Apr 2000 04:16:23 +0000 Subject: [PATCH 1/1] - 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. --- ginac/basic.cpp | 3 + ginac/ex.h | 2 - ginac/expairseq.cpp | 13 +-- ginac/inifcns.cpp | 18 +++- ginac/inifcns_gamma.cpp | 18 ++-- ginac/inifcns_trans.cpp | 27 ++---- ginac/matrix.cpp | 202 ++++++++++++++++++++++++++-------------- ginac/matrix.h | 5 +- ginac/normal.cpp | 40 ++++---- ginac/pseries.cpp | 3 + 10 files changed, 194 insertions(+), 137 deletions(-) diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 5606f364..2a4a3e39 100644 --- 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 { + // 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; } diff --git a/ginac/ex.h b/ginac/ex.h index bf8d0d6f..b5e511e6 100644 --- 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)) -// typedef vector exvector; - #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 d4e13188..bdfc2ff7 100644 --- 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); } -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); @@ -563,7 +552,7 @@ ex expairseq::expand(unsigned options) const // 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); } diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index c03f4a15..f7d2864e 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -107,8 +107,24 @@ static ex csgn_eval(const ex & x) 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 diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index 85275f62..213d9d7f 100644 --- 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... * - * @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)) { @@ -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. - // 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 } diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 924d195e..81f92a93 100644 --- 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)) { + 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); @@ -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); - + // 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()) @@ -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); - } // the branch cut: + } // on the branch cut: const ex point = rel.rhs(); const symbol *s = static_cast(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). diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 0296c183..f76aa9d7 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -32,6 +32,7 @@ #include "debugmsg.h" #include "power.h" #include "symbol.h" +#include "normal.h" #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? - if (4*sparse_countrow==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; r1row; ++r1) { - // shortcut if element(r1,0) vanishes - if (m[r1*col].is_zero()) - continue; - // assemble the minor matrix - for (unsigned r=0; rrow==1) @@ -1057,29 +1015,110 @@ ex matrix::determinant_minor_dense(void) const 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 sigma(col); - for (unsigned i=0; i0) + 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; c0) 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; r2m[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; cm[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(); } diff --git a/ginac/matrix.h b/ginac/matrix.h index 65eebde9..a36276f4 100644 --- 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; - 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); diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 9e24b99b..cdaf8367 100644 --- 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 */ - 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) */ - 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 */ - 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 */ - 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 */ - 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] */ - 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] */ - 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] */ - 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 */ - 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; @@ -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 */ - 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 */ - 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 */ - 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 */ - 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() */ - 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 */ - 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