#include <cln/cl_output.h>
#include <cln/cl_integer_io.h>
#include <cln/cl_integer_ring.h>
-#include <cln/cl_univpoly_integer.h>
#include <cln/cl_rational_io.h>
#include <cln/cl_rational_ring.h>
#include <cln/cl_lfloat_class.h>
#include <cl_output.h>
#include <cl_integer_io.h>
#include <cl_integer_ring.h>
-#include <cl_univpoly_integer.h>
#include <cl_rational_io.h>
#include <cl_rational_ring.h>
#include <cl_lfloat_class.h>
// 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<n as well, we store all results computed so far, which
- // gives us the same convenient runtime behaviour as method (1).
+ // previous ones, the computation is necessarily very expensive. There are
+ // several other ways of computing them, a particularly good one being
+ // cl_I s = 1;
+ // cl_I c = n+1;
+ // cl_RA Bern = 0;
+ // for (unsigned i=0; i<n; i++) {
+ // c = exquo(c*(i-n),(i+2));
+ // Bern = Bern + c*s/(i+2);
+ // s = s + expt_pos(cl_I(i+2),n);
+ // }
+ // return Bern;
+ //
+ // But if somebody works with the n'th Bernoulli number she is likely to
+ // also need all previous Bernoulli numbers we need a complete remember
+ // table and above divide and conquer algorithm is not suited to build one
+ // up. The code below is adapted from Pari's function bernvec().
+ //
+ // (There is an interesting relation with the tangent polynomials described
+ // in `Concrete Mathematics', which leads to a program twice as fast as our
+ // implementation below, but it requires storing one such polynomial in
+ // addition to the remember table. This doubles the memory footprint so
+ // we don't use it.)
// the special cases not covered by the algorithm below
- if (nn.is_zero())
- return _num1();
if (!nn.compare(_num1()))
return numeric(-1,2);
- if (!nn.compare(_num2()))
- return numeric(::cl_I(1)/::cl_I(6));
if (nn.is_odd())
return _num0();
- // let CLN set up the integer ring we are working in
- ::cl_univpoly_integer_ring myring = ::cl_find_univpoly_ring(::cl_I_ring);
-
- // store nonvanishing Bernoulli numbers and the last tangent poly here
+ // store nonvanishing Bernoulli numbers here
static vector< ::cl_RA > 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; i<n/2; ++i) {
+ ::cl_RA B = 0;
+ long n = 8;
+ long m = 5;
+ long d1 = i;
+ long d2 = 2*i-1;
+ for (int j=i; j>0; --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];
}
n.add_ex("point", point);
}
-
-/*
- * Functions overriding virtual functions from base classes
- */
+//////////
+// functions overriding virtual functions from bases classes
+//////////
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
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. */
{
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;
/** 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());
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);
}
* @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<symbol *>(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);
}
// 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();
nul.push_back(expair(Order(_ex1()), _ex0()));
return pseries(relational(var,point), nul);
}
-
+
// Series multiplication
epvector new_seq;
GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
const symbol *s = static_cast<symbol *>(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);
}
git://www.ginac.de / ginac.git / commitdiff