From 3c47d2bc7b421bd504404a8a664349252dddda3a Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Fri, 9 Jun 2000 20:36:11 +0000 Subject: [PATCH] - pseries::print(): Make it print more nicely. - bernoulli(): A slightly slower but only half as memory consuming algorithm ripped from Pari. --- ginac/numeric.cpp | 104 ++++++++++------------ ginac/pseries.cpp | 215 ++++++++++++++++++++++++++-------------------- ginac/pseries.h | 4 +- 3 files changed, 165 insertions(+), 158 deletions(-) diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 3ce95876..f1a073aa 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -51,7 +51,6 @@ #include #include #include -#include #include #include #include @@ -65,7 +64,6 @@ #include #include #include -#include #include #include #include @@ -1490,80 +1488,64 @@ const numeric bernoulli(const numeric & nn) // The Bernoulli numbers are rational numbers that may be computed using // the relation // - // B_n = - 1/(n+1) * sum_{k=0}^{n-1}(binomial(n+1,k)*B_k) (1) + // B_n = - 1/(n+1) * sum_{k=0}^{n-1}(binomial(n+1,k)*B_k) // // with B(0) = 1. Since the n'th Bernoulli number depends on all the - // previous ones, the computation is necessarily very expensive. - // We can do better by computing the tangent numbers, sometimes called - // zag numbers. They are integers which speeds up computation - // considerably. Their relation to the Bernoulli numbers is given by - // - // B_n == T_{n-1} * (-)^(n/2-1) * n / (4^n-2^n) (2) - // - // for even n>=2. We can calculate the tangent numbers as the tangent - // polynomials T_n(x) evaluated at x == 0. The T_n(x) satisfy the - // recurrence relation - // - // T_{n+1}(x) == (1+x^2)*T_n(x)' - // - // with starting value T_0(x) = x, by which they will be computed. The - // n'th tangent polynomial has degree n+1. - // - // Method (2) is about 2-10 times faster than method (1) when it is - // implemented in CLN's efficient univariate polynomials and not much more - // memory consuming. Since any application that needs B_n is likely to - // need B_k, for k results; static int highest_result = 0; - static ::cl_UP_I T_n = myring->create(1); + // algorithm not applicable to B(0), so just store it + if (results.size()==0) + results.push_back(::cl_RA(1)); - // index of results vector - int m = (nn.to_int()-4)/2; - // look if the result is precomputed - if (m < highest_result) - return numeric(results[m]); - // reserve the whole chunk we'll need - if (results.capacity() < (unsigned)(m+1)) - results.reserve(m+1); - - // create the generating polynomial T_gen = 1+x^2 - ::cl_UP_I T_gen = myring->create(2); - ::set_coeff(T_gen, 0, cl_I(1)); - ::set_coeff(T_gen, 2, cl_I(1)); - ::finalize(T_gen); - // T_n will be iterated, start with T_1(x) == 1+x^2 == T_gen - if (highest_result==0) - T_n = T_gen; - // iterate to the n'th tangent polynomial - for (int n=highest_result; n<=m; ++n) { - // recur tangent polynomial twice - for (int i=0; i<2; ++i) - T_n = ::deriv(T_n)*T_gen; - // fetch T_{n-1} - ::cl_I T = ::coeff(T_n,0); - // compute B_n from the T_{n-1}: - ::cl_RA B = T * (n%2 ? ::cl_I(1) : ::cl_I(-1)); - B = B * 2*(n+2) / ((::cl_I(1)<<4*(n+2))-(::cl_I(1)<<2*(n+2))); + int n = nn.to_long(); + for (int i=highest_result; i0; --j) { + B = cl_I(n*m) * (B+results[j]) / (d1*d2); + n += 4; + m += 2; + d1 -= 1; + d2 -= 2; + } + B = (1 - ((B+1)/(2*i+3))) / (cl_I(1)<<(2*i+2)); results.push_back(B); ++highest_result; } - return numeric(results[m]); + return results[n/2]; } diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 58f8d629..9ed6aab3 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -151,10 +151,9 @@ void pseries::archive(archive_node &n) const n.add_ex("point", point); } - -/* - * Functions overriding virtual functions from base classes - */ +////////// +// functions overriding virtual functions from bases classes +////////// basic *pseries::duplicate() const { @@ -165,20 +164,46 @@ basic *pseries::duplicate() const 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); + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) { + // print a sign, if needed + if (i!=seq.begin()) + os << '+'; + if (!is_order_function(i->rest)) { + // print 'rest', i.e. the expansion coefficient + if (i->rest.info(info_flags::numeric) && + i->rest.info(info_flags::positive)) { + os << i->rest; + } else + os << "(" << i->rest << ')'; + // print 'coeff', something like (x-1)^42 + if (!i->coeff.is_zero()) { + os << '*'; + if (!point.is_zero()) + os << '(' << var-point << ')'; + else + os << var; + if (i->coeff.compare(_ex1())) { + os << '^'; + if (i->coeff.info(info_flags::negative)) + os << '(' << i->coeff << ')'; + else + os << i->coeff; + } + } + } else { + os << Order(power(var-point,i->coeff)); + } + } } void pseries::printraw(ostream &os) const { - debugmsg("pseries printraw", LOGLEVEL_PRINT); - os << "pseries(" << var << ";" << point << ";"; - for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) { - os << "(" << (*i).rest << "," << (*i).coeff << "),"; - } - os << ")"; + debugmsg("pseries printraw", LOGLEVEL_PRINT); + os << "pseries(" << var << ";" << point << ";"; + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) { + os << "(" << (*i).rest << "," << (*i).coeff << "),"; + } + os << ")"; } void pseries::printtree(ostream & os, unsigned indent) const @@ -264,37 +289,37 @@ int pseries::ldegree(const symbol &s) const ex pseries::coeff(const symbol &s, int n) const { if (var.is_equal(s)) { - if (seq.size() == 0) - return _ex0(); - - // Binary search in sequence for given power - numeric looking_for = numeric(n); - int lo = 0, hi = seq.size() - 1; - while (lo <= hi) { - int mid = (lo + hi) / 2; - GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric)); - int cmp = ex_to_numeric(seq[mid].coeff).compare(looking_for); - switch (cmp) { - case -1: - lo = mid + 1; - break; - case 0: - return seq[mid].rest; - case 1: - hi = mid - 1; - break; - default: - throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1")); - } - } - return _ex0(); + if (seq.size() == 0) + return _ex0(); + + // Binary search in sequence for given power + numeric looking_for = numeric(n); + int lo = 0, hi = seq.size() - 1; + while (lo <= hi) { + int mid = (lo + hi) / 2; + GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric)); + int cmp = ex_to_numeric(seq[mid].coeff).compare(looking_for); + switch (cmp) { + case -1: + lo = mid + 1; + break; + case 0: + return seq[mid].rest; + case 1: + hi = mid - 1; + break; + default: + throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1")); + } + } + return _ex0(); } else return convert_to_poly().coeff(s, n); } ex pseries::collect(const symbol &s) const { - return *this; + return *this; } /** Evaluate coefficients. */ @@ -302,9 +327,9 @@ ex pseries::eval(int level) const { if (level == 1) return this->hold(); - - if (level == -max_recursion_level) - throw (std::runtime_error("pseries::eval(): recursion limit exceeded")); + + if (level == -max_recursion_level) + throw (std::runtime_error("pseries::eval(): recursion limit exceeded")); // Construct a new series with evaluated coefficients epvector new_seq; @@ -320,12 +345,12 @@ ex pseries::eval(int level) const /** Evaluate coefficients numerically. */ ex pseries::evalf(int level) const { - if (level == 1) - return *this; - - if (level == -max_recursion_level) - throw (std::runtime_error("pseries::evalf(): recursion limit exceeded")); - + if (level == 1) + return *this; + + if (level == -max_recursion_level) + throw (std::runtime_error("pseries::evalf(): recursion limit exceeded")); + // Construct a new series with evaluated coefficients epvector new_seq; new_seq.reserve(seq.size()); @@ -339,21 +364,21 @@ ex pseries::evalf(int level) const ex pseries::subs(const lst & ls, const lst & lr) const { - // If expansion variable is being substituted, convert the series to a - // polynomial and do the substitution there because the result might - // no longer be a power series - if (ls.has(var)) - return convert_to_poly(true).subs(ls, lr); - - // Otherwise construct a new series with substituted coefficients and - // expansion point - epvector new_seq; - new_seq.reserve(seq.size()); - epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { - new_seq.push_back(expair(it->rest.subs(ls, lr), it->coeff)); - it++; - } + // If expansion variable is being substituted, convert the series to a + // polynomial and do the substitution there because the result might + // no longer be a power series + if (ls.has(var)) + return convert_to_poly(true).subs(ls, lr); + + // Otherwise construct a new series with substituted coefficients and + // expansion point + epvector new_seq; + new_seq.reserve(seq.size()); + epvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + new_seq.push_back(expair(it->rest.subs(ls, lr), it->coeff)); + it++; + } return (new pseries(relational(var,point.subs(ls, lr)), new_seq))->setflag(status_flags::dynallocated); } @@ -450,21 +475,21 @@ ex basic::series(const relational & r, int order) const * @see ex::series */ ex symbol::series(const relational & r, int order) const { - epvector seq; + epvector seq; const ex point = r.rhs(); GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol)); const symbol *s = static_cast(r.lhs().bp); - if (this->is_equal(*s)) { - if (order > 0 && !point.is_zero()) - seq.push_back(expair(point, _ex0())); - if (order > 1) - seq.push_back(expair(_ex1(), _ex1())); - else - seq.push_back(expair(Order(_ex1()), numeric(order))); - } else - seq.push_back(expair(*this, _ex0())); - return pseries(r, seq); + if (this->is_equal(*s)) { + if (order > 0 && !point.is_zero()) + seq.push_back(expair(point, _ex0())); + if (order > 1) + seq.push_back(expair(_ex1(), _ex1())); + else + seq.push_back(expair(Order(_ex1()), numeric(order))); + } else + seq.push_back(expair(*this, _ex0())); + return pseries(r, seq); } @@ -551,7 +576,7 @@ ex add::series(const relational & r, int order) const // Get first term from overall_coeff acc = overall_coeff.series(r, order); - + // Add remaining terms epvector::const_iterator it = seq.begin(); epvector::const_iterator itend = seq.end(); @@ -607,7 +632,7 @@ ex pseries::mul_series(const pseries &other) const nul.push_back(expair(Order(_ex1()), _ex0())); return pseries(relational(var,point), nul); } - + // Series multiplication epvector new_seq; @@ -761,25 +786,25 @@ ex pseries::series(const relational & r, int order) const GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol)); const symbol *s = static_cast(r.lhs().bp); - if (var.is_equal(*s) && point.is_equal(p)) { - if (order > degree(*s)) - return *this; - else { - epvector new_seq; - epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { - int o = ex_to_numeric(it->coeff).to_int(); - if (o >= order) { - new_seq.push_back(expair(Order(_ex1()), o)); - break; - } - new_seq.push_back(*it); - it++; - } - return pseries(r, new_seq); - } - } else - return convert_to_poly().series(r, order); + if (var.is_equal(*s) && point.is_equal(p)) { + if (order > degree(*s)) + return *this; + else { + epvector new_seq; + epvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + int o = ex_to_numeric(it->coeff).to_int(); + if (o >= order) { + new_seq.push_back(expair(Order(_ex1()), o)); + break; + } + new_seq.push_back(*it); + it++; + } + return pseries(r, new_seq); + } + } else + return convert_to_poly().series(r, order); } diff --git a/ginac/pseries.h b/ginac/pseries.h index 9346c255..a0a494ee 100644 --- a/ginac/pseries.h +++ b/ginac/pseries.h @@ -107,7 +107,7 @@ extern const type_info & typeid_pseries; * @see is_ex_of_type */ inline const pseries &ex_to_pseries(const ex &e) { - return static_cast(*e.bp); + return static_cast(*e.bp); } /** Convert the pseries object embedded in an expression to an ordinary @@ -120,7 +120,7 @@ inline const pseries &ex_to_pseries(const ex &e) * @see pseries::convert_to_poly */ inline ex series_to_poly(const ex &e) { - return (static_cast(*e.bp).convert_to_poly(true)); + return (static_cast(*e.bp).convert_to_poly(true)); } #ifndef NO_NAMESPACE_GINAC -- 2.44.0