author Richard Kreckel Fri, 9 Jun 2000 20:36:11 +0000 (20:36 +0000) committer Richard Kreckel Fri, 9 Jun 2000 20:36:11 +0000 (20:36 +0000)
- bernoulli(): A slightly slower but only half as memory consuming algorithm
ripped from Pari.

 ginac/numeric.cpp patch | blob | history ginac/pseries.cpp patch | blob | history ginac/pseries.h patch | blob | history

index 3ce95876d7d8547e8b9fbdbec9b9185881eb77fd..f1a073aa61ebf37af51a4170ac3d3cbcfec5cdd0 100644 (file)
@@ -51,7 +51,6 @@
#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>
@@ -65,7 +64,6 @@
#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>
@@ -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<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];
}

@@ -151,10 +151,9 @@ void pseries::archive(archive_node &n) const
}

-
-/*
- *  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<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);
}

@@ -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);
-
+
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<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);
}

@@ -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<const pseries &>(*e.bp);
+    return static_cast<const pseries &>(*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<const pseries &>(*e.bp).convert_to_poly(true));
+    return (static_cast<const pseries &>(*e.bp).convert_to_poly(true));
}

#ifndef NO_NAMESPACE_GINAC