]> www.ginac.de Git - ginac.git/blobdiff - ginac/numeric.cpp
* irem(), iquo(): throw an exception, when second argument vanishes.
[ginac.git] / ginac / numeric.cpp
index bdd8903c72068f4bb254cb277273af5f3be75342..858cbbd578ee0ba2fd332876249ddef43e8c026a 100644 (file)
@@ -7,7 +7,7 @@
  *  of special functions or implement the interface to the bignum package. */
 
 /*
- *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
 #include <vector>
 #include <stdexcept>
 #include <string>
-
-#if defined(HAVE_SSTREAM)
 #include <sstream>
-#elif defined(HAVE_STRSTREAM)
-#include <strstream>
-#else
-#error Need either sstream or strstream
-#endif
 
 #include "numeric.h"
 #include "ex.h"
 #include "print.h"
 #include "archive.h"
-#include "debugmsg.h"
+#include "tostring.h"
 #include "utils.h"
 
 // CLN should pollute the global namespace as little as possible.  Hence, we
@@ -69,14 +62,12 @@ namespace GiNaC {
 GINAC_IMPLEMENT_REGISTERED_CLASS(numeric, basic)
 
 //////////
-// default ctor, dtor, copy ctor assignment
-// operator and helpers
+// default ctor, dtor, copy ctor, assignment operator and helpers
 //////////
 
 /** default ctor. Numerically it initializes to an integer zero. */
 numeric::numeric() : basic(TINFO_numeric)
 {
-       debugmsg("numeric default ctor", LOGLEVEL_CONSTRUCT);
        value = cln::cl_I(0);
        setflag(status_flags::evaluated | status_flags::expanded);
 }
@@ -97,13 +88,12 @@ DEFAULT_DESTROY(numeric)
 
 numeric::numeric(int i) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from int",LOGLEVEL_CONSTRUCT);
        // Not the whole int-range is available if we don't cast to long
        // first.  This is due to the behaviour of the cl_I-ctor, which
-       // emphasizes efficiency.  However, if the integer is small enough
-       // i.e. satisfies cl_immediate_p(), we save space and dereferences by
-       // using an immediate type:
-       if (cln::cl_immediate_p(i))
+       // emphasizes efficiency.  However, if the integer is small enough
+       // we save space and dereferences by using an immediate type.
+       // (C.f. <cln/object.h>)
+       if (i < (1U<<cl_value_len-1))
                value = cln::cl_I(i);
        else
                value = cln::cl_I((long) i);
@@ -113,13 +103,12 @@ numeric::numeric(int i) : basic(TINFO_numeric)
 
 numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from uint",LOGLEVEL_CONSTRUCT);
        // Not the whole uint-range is available if we don't cast to ulong
        // first.  This is due to the behaviour of the cl_I-ctor, which
-       // emphasizes efficiency.  However, if the integer is small enough
-       // i.e. satisfies cl_immediate_p(), we save space and dereferences by
-       // using an immediate type:
-       if (cln::cl_immediate_p(i))
+       // emphasizes efficiency.  However, if the integer is small enough
+       // we save space and dereferences by using an immediate type.
+       // (C.f. <cln/object.h>)
+       if (i < (1U<<cl_value_len-1))
                value = cln::cl_I(i);
        else
                value = cln::cl_I((unsigned long) i);
@@ -129,7 +118,6 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 
 numeric::numeric(long i) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from long",LOGLEVEL_CONSTRUCT);
        value = cln::cl_I(i);
        setflag(status_flags::evaluated | status_flags::expanded);
 }
@@ -137,7 +125,6 @@ numeric::numeric(long i) : basic(TINFO_numeric)
 
 numeric::numeric(unsigned long i) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from ulong",LOGLEVEL_CONSTRUCT);
        value = cln::cl_I(i);
        setflag(status_flags::evaluated | status_flags::expanded);
 }
@@ -147,7 +134,6 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric)
  *  @exception overflow_error (division by zero) */
 numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from long/long",LOGLEVEL_CONSTRUCT);
        if (!denom)
                throw std::overflow_error("division by zero");
        value = cln::cl_I(numer) / cln::cl_I(denom);
@@ -157,7 +143,6 @@ numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
 
 numeric::numeric(double d) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from double",LOGLEVEL_CONSTRUCT);
        // We really want to explicitly use the type cl_LF instead of the
        // more general cl_F, since that would give us a cl_DF only which
        // will not be promoted to cl_LF if overflow occurs:
@@ -170,40 +155,47 @@ numeric::numeric(double d) : basic(TINFO_numeric)
  *  notation like "2+5*I". */
 numeric::numeric(const char *s) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from string",LOGLEVEL_CONSTRUCT);
        cln::cl_N ctorval = 0;
        // parse complex numbers (functional but not completely safe, unfortunately
        // std::string does not understand regexpese):
        // ss should represent a simple sum like 2+5*I
-       std::string ss(s);
-       // make it safe by adding explicit sign
+       std::string ss = s;
+       std::string::size_type delim;
+
+       // make this implementation safe by adding explicit sign
        if (ss.at(0) != '+' && ss.at(0) != '-' && ss.at(0) != '#')
                ss = '+' + ss;
-       std::string::size_type delim;
+
+       // We use 'E' as exponent marker in the output, but some people insist on
+       // writing 'e' at input, so let's substitute them right at the beginning:
+       while ((delim = ss.find("e"))!=std::string::npos)
+               ss.replace(delim,1,"E");
+
+       // main parser loop:
        do {
                // chop ss into terms from left to right
                std::string term;
                bool imaginary = false;
                delim = ss.find_first_of(std::string("+-"),1);
                // Do we have an exponent marker like "31.415E-1"?  If so, hop on!
-               if ((delim != std::string::npos) && (ss.at(delim-1) == 'E'))
+               if (delim!=std::string::npos && ss.at(delim-1)=='E')
                        delim = ss.find_first_of(std::string("+-"),delim+1);
                term = ss.substr(0,delim);
-               if (delim != std::string::npos)
+               if (delim!=std::string::npos)
                        ss = ss.substr(delim);
                // is the term imaginary?
-               if (term.find("I") != std::string::npos) {
+               if (term.find("I")!=std::string::npos) {
                        // erase 'I':
-                       term = term.replace(term.find("I"),1,"");
+                       term.erase(term.find("I"),1);
                        // erase '*':
-                       if (term.find("*") != std::string::npos)
-                               term = term.replace(term.find("*"),1,"");
+                       if (term.find("*")!=std::string::npos)
+                               term.erase(term.find("*"),1);
                        // correct for trivial +/-I without explicit factor on I:
-                       if (term.size() == 1)
-                               term += "1";
+                       if (term.size()==1)
+                               term += '1';
                        imaginary = true;
                }
-               if (term.find(".") != std::string::npos) {
+               if (term.find('.')!=std::string::npos || term.find('E')!=std::string::npos) {
                        // CLN's short type cl_SF is not very useful within the GiNaC
                        // framework where we are mainly interested in the arbitrary
                        // precision type cl_LF.  Hence we go straight to the construction
@@ -214,33 +206,25 @@ numeric::numeric(const char *s) : basic(TINFO_numeric)
                        // 31.4E-1   -->   31.4e-1_<Digits>
                        // and s on.
                        // No exponent marker?  Let's add a trivial one.
-                       if (term.find("E") == std::string::npos)
+                       if (term.find("E")==std::string::npos)
                                term += "E0";
                        // E to lower case
                        term = term.replace(term.find("E"),1,"e");
                        // append _<Digits> to term
-#if defined(HAVE_SSTREAM)
-                       std::ostringstream buf;
-                       buf << unsigned(Digits) << std::ends;
-                       term += "_" + buf.str();
-#else
-                       char buf[14];
-                       std::ostrstream(buf,sizeof(buf)) << unsigned(Digits) << std::ends;
-                       term += "_" + std::string(buf);
-#endif
+                       term += "_" + ToString((unsigned)Digits);
                        // construct float using cln::cl_F(const char *) ctor.
                        if (imaginary)
                                ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_F(term.c_str()));
                        else
                                ctorval = ctorval + cln::cl_F(term.c_str());
                } else {
-                       // not a floating point number...
+                       // this is not a floating point number...
                        if (imaginary)
                                ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_R(term.c_str()));
                        else
                                ctorval = ctorval + cln::cl_R(term.c_str());
                }
-       } while(delim != std::string::npos);
+       } while (delim != std::string::npos);
        value = ctorval;
        setflag(status_flags::evaluated | status_flags::expanded);
 }
@@ -250,7 +234,6 @@ numeric::numeric(const char *s) : basic(TINFO_numeric)
  *  only. */
 numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric)
 {
-       debugmsg("numeric ctor from cl_N", LOGLEVEL_CONSTRUCT);
        value = z;
        setflag(status_flags::evaluated | status_flags::expanded);
 }
@@ -261,17 +244,12 @@ numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric)
 
 numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
-       debugmsg("numeric ctor from archive_node", LOGLEVEL_CONSTRUCT);
        cln::cl_N ctorval = 0;
 
        // Read number as string
        std::string str;
        if (n.find_string("number", str)) {
-#ifdef HAVE_SSTREAM
                std::istringstream s(str);
-#else
-               std::istrstream s(str.c_str(), str.size() + 1);
-#endif
                cln::cl_idecoded_float re, im;
                char c;
                s.get(c);
@@ -301,12 +279,7 @@ void numeric::archive(archive_node &n) const
        inherited::archive(n);
 
        // Write number as string
-#ifdef HAVE_SSTREAM
        std::ostringstream s;
-#else
-       char buf[1024];
-       std::ostrstream s(buf, 1024);
-#endif
        if (this->is_crational())
                s << cln::the<cln::cl_N>(value);
        else {
@@ -324,13 +297,7 @@ void numeric::archive(archive_node &n) const
                        s << im.sign << " " << im.mantissa << " " << im.exponent;
                }
        }
-#ifdef HAVE_SSTREAM
        n.add_string("number", s.str());
-#else
-       s << ends;
-       std::string str(buf);
-       n.add_string("number", str);
-#endif
 }
 
 DEFAULT_UNARCHIVE(numeric)
@@ -376,8 +343,6 @@ static void print_real_number(const print_context & c, const cln::cl_R &x)
  *  @see print_real_number() */
 void numeric::print(const print_context & c, unsigned level) const
 {
-       debugmsg("numeric print", LOGLEVEL_PRINT);
-
        if (is_a<print_tree>(c)) {
 
                c.s << std::string(level, ' ') << cln::the<cln::cl_N>(value)
@@ -389,8 +354,13 @@ void numeric::print(const print_context & c, unsigned level) const
 
                std::ios::fmtflags oldflags = c.s.flags();
                c.s.setf(std::ios::scientific);
+               int oldprec = c.s.precision();
+               if (is_a<print_csrc_double>(c))
+                       c.s.precision(16);
+               else
+                       c.s.precision(7);
                if (this->is_rational() && !this->is_integer()) {
-                       if (compare(_num0()) > 0) {
+                       if (compare(_num0) > 0) {
                                c.s << "(";
                                if (is_a<print_csrc_cl_N>(c))
                                        c.s << "cln::cl_F(\"" << numer().evalf() << "\")";
@@ -416,6 +386,7 @@ void numeric::print(const print_context & c, unsigned level) const
                                c.s << to_double();
                }
                c.s.flags(oldflags);
+               c.s.precision(oldprec);
 
        } else {
                const std::string par_open  = is_a<print_latex>(c) ? "{(" : "(";
@@ -424,6 +395,8 @@ void numeric::print(const print_context & c, unsigned level) const
                const std::string mul_sym   = is_a<print_latex>(c) ? " " : "*";
                const cln::cl_R r = cln::realpart(cln::the<cln::cl_N>(value));
                const cln::cl_R i = cln::imagpart(cln::the<cln::cl_N>(value));
+               if (is_a<print_python_repr>(c))
+                       c.s << class_name() << "('";
                if (cln::zerop(i)) {
                        // case 1, real:  x  or  -x
                        if ((precedence() <= level) && (!this->is_nonneg_integer())) {
@@ -436,25 +409,19 @@ void numeric::print(const print_context & c, unsigned level) const
                } else {
                        if (cln::zerop(r)) {
                                // case 2, imaginary:  y*I  or  -y*I
-                               if ((precedence() <= level) && (i < 0)) {
-                                       if (i == -1) {
-                                               c.s << par_open+imag_sym+par_close;
-                                       } else {
+                               if (i==1)
+                                       c.s << imag_sym;
+                               else {
+                                       if (precedence()<=level)
                                                c.s << par_open;
+                                       if (i == -1)
+                                               c.s << "-" << imag_sym;
+                                       else {
                                                print_real_number(c, i);
-                                               c.s << mul_sym+imag_sym+par_close;
-                                       }
-                               } else {
-                                       if (i == 1) {
-                                               c.s << imag_sym;
-                                       } else {
-                                               if (i == -1) {
-                                                       c.s << "-" << imag_sym;
-                                               } else {
-                                                       print_real_number(c, i);
-                                                       c.s << mul_sym+imag_sym;
-                                               }
+                                               c.s << mul_sym+imag_sym;
                                        }
+                                       if (precedence()<=level)
+                                               c.s << par_close;
                                }
                        } else {
                                // case 3, complex:  x+y*I  or  x-y*I  or  -x+y*I  or  -x-y*I
@@ -481,6 +448,8 @@ void numeric::print(const print_context & c, unsigned level) const
                                        c.s << par_close;
                        }
                }
+               if (is_a<print_python_repr>(c))
+                       c.s << "')";
        }
 }
 
@@ -529,6 +498,21 @@ bool numeric::info(unsigned inf) const
        return false;
 }
 
+int numeric::degree(const ex & s) const
+{
+       return 0;
+}
+
+int numeric::ldegree(const ex & s) const
+{
+       return 0;
+}
+
+ex numeric::coeff(const ex & s, int n) const
+{
+       return n==0 ? *this : _ex0;
+}
+
 /** Disassemble real part and imaginary part to scan for the occurrence of a
  *  single number.  Also handles the imaginary unit.  It ignores the sign on
  *  both this and the argument, which may lead to what might appear as funny
@@ -583,7 +567,7 @@ ex numeric::evalf(int level) const
 
 int numeric::compare_same_type(const basic &other) const
 {
-       GINAC_ASSERT(is_exactly_of_type(other, numeric));
+       GINAC_ASSERT(is_exactly_a<numeric>(other));
        const numeric &o = static_cast<const numeric &>(other);
        
        return this->compare(o);
@@ -592,7 +576,7 @@ int numeric::compare_same_type(const basic &other) const
 
 bool numeric::is_equal_same_type(const basic &other) const
 {
-       GINAC_ASSERT(is_exactly_of_type(other,numeric));
+       GINAC_ASSERT(is_exactly_a<numeric>(other));
        const numeric &o = static_cast<const numeric &>(other);
        
        return this->is_equal(o);
@@ -626,10 +610,9 @@ unsigned numeric::calchash(void) const
 const numeric numeric::add(const numeric &other) const
 {
        // Efficiency shortcut: trap the neutral element by pointer.
-       static const numeric * _num0p = &_num0();
-       if (this==_num0p)
+       if (this==_num0_p)
                return other;
-       else if (&other==_num0p)
+       else if (&other==_num0_p)
                return *this;
        
        return numeric(cln::the<cln::cl_N>(value)+cln::the<cln::cl_N>(other.value));
@@ -649,10 +632,9 @@ const numeric numeric::sub(const numeric &other) const
 const numeric numeric::mul(const numeric &other) const
 {
        // Efficiency shortcut: trap the neutral element by pointer.
-       static const numeric * _num1p = &_num1();
-       if (this==_num1p)
+       if (this==_num1_p)
                return other;
-       else if (&other==_num1p)
+       else if (&other==_num1_p)
                return *this;
        
        return numeric(cln::the<cln::cl_N>(value)*cln::the<cln::cl_N>(other.value));
@@ -676,8 +658,7 @@ const numeric numeric::div(const numeric &other) const
 const numeric numeric::power(const numeric &other) const
 {
        // Efficiency shortcut: trap the neutral exponent by pointer.
-       static const numeric * _num1p = &_num1();
-       if (&other==_num1p)
+       if (&other==_num1_p)
                return *this;
        
        if (cln::zerop(cln::the<cln::cl_N>(value))) {
@@ -688,7 +669,7 @@ const numeric numeric::power(const numeric &other) const
                else if (cln::minusp(cln::realpart(cln::the<cln::cl_N>(other.value))))
                        throw std::overflow_error("numeric::eval(): division by zero");
                else
-                       return _num0();
+                       return _num0;
        }
        return numeric(cln::expt(cln::the<cln::cl_N>(value),cln::the<cln::cl_N>(other.value)));
 }
@@ -697,10 +678,9 @@ const numeric numeric::power(const numeric &other) const
 const numeric &numeric::add_dyn(const numeric &other) const
 {
        // Efficiency shortcut: trap the neutral element by pointer.
-       static const numeric * _num0p = &_num0();
-       if (this==_num0p)
+       if (this==_num0_p)
                return other;
-       else if (&other==_num0p)
+       else if (&other==_num0_p)
                return *this;
        
        return static_cast<const numeric &>((new numeric(cln::the<cln::cl_N>(value)+cln::the<cln::cl_N>(other.value)))->
@@ -718,10 +698,9 @@ const numeric &numeric::sub_dyn(const numeric &other) const
 const numeric &numeric::mul_dyn(const numeric &other) const
 {
        // Efficiency shortcut: trap the neutral element by pointer.
-       static const numeric * _num1p = &_num1();
-       if (this==_num1p)
+       if (this==_num1_p)
                return other;
-       else if (&other==_num1p)
+       else if (&other==_num1_p)
                return *this;
        
        return static_cast<const numeric &>((new numeric(cln::the<cln::cl_N>(value)*cln::the<cln::cl_N>(other.value)))->
@@ -741,8 +720,7 @@ const numeric &numeric::div_dyn(const numeric &other) const
 const numeric &numeric::power_dyn(const numeric &other) const
 {
        // Efficiency shortcut: trap the neutral exponent by pointer.
-       static const numeric * _num1p=&_num1();
-       if (&other==_num1p)
+       if (&other==_num1_p)
                return *this;
        
        if (cln::zerop(cln::the<cln::cl_N>(value))) {
@@ -753,7 +731,7 @@ const numeric &numeric::power_dyn(const numeric &other) const
                else if (cln::minusp(cln::realpart(cln::the<cln::cl_N>(other.value))))
                        throw std::overflow_error("numeric::eval(): division by zero");
                else
-                       return _num0();
+                       return _num0;
        }
        return static_cast<const numeric &>((new numeric(cln::expt(cln::the<cln::cl_N>(value),cln::the<cln::cl_N>(other.value))))->
                                             setflag(status_flags::dynallocated));
@@ -1120,7 +1098,7 @@ const numeric numeric::numer(void) const
 const numeric numeric::denom(void) const
 {
        if (this->is_integer())
-               return _num1();
+               return _num1;
        
        if (cln::instanceof(value, cln::cl_RA_ring))
                return numeric(cln::denominator(cln::the<cln::cl_RA>(value)));
@@ -1129,7 +1107,7 @@ const numeric numeric::denom(void) const
                const cln::cl_RA r = cln::the<cln::cl_RA>(cln::realpart(cln::the<cln::cl_N>(value)));
                const cln::cl_RA i = cln::the<cln::cl_RA>(cln::imagpart(cln::the<cln::cl_N>(value)));
                if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_I_ring))
-                       return _num1();
+                       return _num1;
                if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_RA_ring))
                        return numeric(cln::denominator(i));
                if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_I_ring))
@@ -1138,7 +1116,7 @@ const numeric numeric::denom(void) const
                        return numeric(cln::lcm(cln::denominator(r), cln::denominator(i)));
        }
        // at least one float encountered
-       return _num1();
+       return _num1;
 }
 
 
@@ -1242,7 +1220,7 @@ const numeric atan(const numeric &x)
 {
        if (!x.is_real() &&
            x.real().is_zero() &&
-           abs(x.imag()).is_equal(_num1()))
+           abs(x.imag()).is_equal(_num1))
                throw pole_error("atan(): logarithmic pole",0);
        return cln::atan(x.to_cl_N());
 }
@@ -1393,7 +1371,7 @@ static cln::cl_N Li2_projection(const cln::cl_N &x,
 const numeric Li2(const numeric &x)
 {
        if (x.is_zero())
-               return _num0();
+               return _num0;
        
        // what is the desired float format?
        // first guess: default format
@@ -1484,8 +1462,8 @@ const numeric factorial(const numeric &n)
  *  @exception range_error (argument must be integer >= -1) */
 const numeric doublefactorial(const numeric &n)
 {
-       if (n.is_equal(_num_1()))
-               return _num1();
+       if (n.is_equal(_num_1))
+               return _num1;
        
        if (!n.is_nonneg_integer())
                throw std::range_error("numeric::doublefactorial(): argument must be integer >= -1");
@@ -1502,12 +1480,12 @@ const numeric binomial(const numeric &n, const numeric &k)
 {
        if (n.is_integer() && k.is_integer()) {
                if (n.is_nonneg_integer()) {
-                       if (k.compare(n)!=1 && k.compare(_num0())!=-1)
+                       if (k.compare(n)!=1 && k.compare(_num0)!=-1)
                                return numeric(cln::binomial(n.to_int(),k.to_int()));
                        else
-                               return _num0();
+                               return _num0;
                } else {
-                       return _num_1().power(k)*binomial(k-n-_num1(),k);
+                       return _num_1.power(k)*binomial(k-n-_num1,k);
                }
        }
        
@@ -1525,7 +1503,7 @@ const numeric bernoulli(const numeric &nn)
 {
        if (!nn.is_integer() || nn.is_negative())
                throw std::range_error("numeric::bernoulli(): argument must be integer >= 0");
-       
+
        // Method:
        //
        // The Bernoulli numbers are rational numbers that may be computed using
@@ -1549,46 +1527,61 @@ const numeric bernoulli(const numeric &nn)
        // But if somebody works with the n'th Bernoulli number she is likely to
        // also need all previous Bernoulli numbers. So 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().
+       // up.  The formula below accomplishes this.  It is a modification of the
+       // defining formula above but the computation of the binomial coefficients
+       // is carried along in an inline fashion.  It also honors the fact that
+       // B_n is zero when n is odd and greater than 1.
        // 
        // (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
+       // in `Concrete Mathematics', which leads to a program a little faster 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.)
-       
+
+       const unsigned n = nn.to_int();
+
        // the special cases not covered by the algorithm below
-       if (nn.is_equal(_num1()))
-               return _num_1_2();
-       if (nn.is_odd())
-               return _num0();
-       
+       if (n & 1)
+               return (n==1) ? _num_1_2 : _num0;
+       if (!n)
+                return _num1;
+
        // store nonvanishing Bernoulli numbers here
        static std::vector< cln::cl_RA > results;
-       static int highest_result = 0;
-       // algorithm not applicable to B(0), so just store it
-       if (results.empty())
-               results.push_back(cln::cl_RA(1));
-       
-       int n = nn.to_long();
-       for (int i=highest_result; i<n/2; ++i) {
-               cln::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 = cln::cl_I(n*m) * (B+results[j]) / (d1*d2);
-                       n += 4;
-                       m += 2;
-                       d1 -= 1;
-                       d2 -= 2;
-               }
-               B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2));
-               results.push_back(B);
-               ++highest_result;
+       static unsigned next_r = 0;
+
+       // algorithm not applicable to B(2), so just store it
+       if (!next_r) {
+               results.push_back(cln::recip(cln::cl_RA(6)));
+               next_r = 4;
        }
-       return results[n/2];
+       if (n<next_r)
+               return results[n/2-1];
+
+       results.reserve(n/2);
+       for (unsigned p=next_r; p<=n;  p+=2) {
+               cln::cl_I  c = 1;  // seed for binonmial coefficients
+               cln::cl_RA b = cln::cl_RA(1-p)/2;
+               const unsigned p3 = p+3;
+               const unsigned pm = p-2;
+               unsigned i, k, p_2;
+               // test if intermediate unsigned int can be represented by immediate
+               // objects by CLN (i.e. < 2^29 for 32 Bit machines, see <cln/object.h>)
+               if (p < (1UL<<cl_value_len/2)) {
+                       for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+                               c = cln::exquo(c * ((p3-i) * p_2), (i-1)*k);
+                               b = b + c*results[k-1];
+                       }
+               } else {
+                       for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+                               c = cln::exquo((c * (p3-i)) * p_2, cln::cl_I(i-1)*k);
+                               b = b + c*results[k-1];
+                       }
+               }
+               results.push_back(-b/(p+1));
+       }
+       next_r = n+2;
+       return results[n/2-1];
 }
 
 
@@ -1619,7 +1612,7 @@ const numeric fibonacci(const numeric &n)
        // hence
        //      F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
        if (n.is_zero())
-               return _num0();
+               return _num0;
        if (n.is_negative())
                if (n.is_even())
                        return -fibonacci(-n);
@@ -1671,7 +1664,7 @@ const numeric mod(const numeric &a, const numeric &b)
                return cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
                                cln::the<cln::cl_I>(b.to_cl_N()));
        else
-               return _num0();
+               return _num0;
 }
 
 
@@ -1686,7 +1679,7 @@ const numeric smod(const numeric &a, const numeric &b)
                return cln::mod(cln::the<cln::cl_I>(a.to_cl_N()) + b2,
                                cln::the<cln::cl_I>(b.to_cl_N())) - b2;
        } else
-               return _num0();
+               return _num0;
 }
 
 
@@ -1695,14 +1688,17 @@ const numeric smod(const numeric &a, const numeric &b)
  *  In general, mod(a,b) has the sign of b or is zero, and irem(a,b) has the
  *  sign of a or is zero.
  *
- *  @return remainder of a/b if both are integer, 0 otherwise. */
+ *  @return remainder of a/b if both are integer, 0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric irem(const numeric &a, const numeric &b)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::irem(): division by zero");
        if (a.is_integer() && b.is_integer())
                return cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
                                cln::the<cln::cl_I>(b.to_cl_N()));
        else
-               return _num0();
+               return _num0;
 }
 
 
@@ -1712,17 +1708,20 @@ const numeric irem(const numeric &a, const numeric &b)
  *  and irem(a,b) has the sign of a or is zero.  
  *
  *  @return remainder of a/b and quotient stored in q if both are integer,
- *  0 otherwise. */
+ *  0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric irem(const numeric &a, const numeric &b, numeric &q)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::irem(): division by zero");
        if (a.is_integer() && b.is_integer()) {
                const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
                                                               cln::the<cln::cl_I>(b.to_cl_N()));
                q = rem_quo.quotient;
                return rem_quo.remainder;
        } else {
-               q = _num0();
-               return _num0();
+               q = _num0;
+               return _num0;
        }
 }
 
@@ -1730,14 +1729,17 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q)
 /** Numeric integer quotient.
  *  Equivalent to Maple's iquo as far as sign conventions are concerned.
  *  
- *  @return truncated quotient of a/b if both are integer, 0 otherwise. */
+ *  @return truncated quotient of a/b if both are integer, 0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric iquo(const numeric &a, const numeric &b)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::iquo(): division by zero");
        if (a.is_integer() && b.is_integer())
                return cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
                                  cln::the<cln::cl_I>(b.to_cl_N()));
        else
-               return _num0();
+               return _num0;
 }
 
 
@@ -1746,17 +1748,20 @@ const numeric iquo(const numeric &a, const numeric &b)
  *  r == a - iquo(a,b,r)*b.
  *
  *  @return truncated quotient of a/b and remainder stored in r if both are
- *  integer, 0 otherwise. */
+ *  integer, 0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric iquo(const numeric &a, const numeric &b, numeric &r)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::iquo(): division by zero");
        if (a.is_integer() && b.is_integer()) {
                const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
                                                               cln::the<cln::cl_I>(b.to_cl_N()));
                r = rem_quo.remainder;
                return rem_quo.quotient;
        } else {
-               r = _num0();
-               return _num0();
+               r = _num0;
+               return _num0;
        }
 }
 
@@ -1771,7 +1776,7 @@ const numeric gcd(const numeric &a, const numeric &b)
                return cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
                                cln::the<cln::cl_I>(b.to_cl_N()));
        else
-               return _num1();
+               return _num1;
 }
 
 
@@ -1811,7 +1816,7 @@ const numeric isqrt(const numeric &x)
                cln::isqrt(cln::the<cln::cl_I>(x.to_cl_N()), &root);
                return root;
        } else
-               return _num0();
+               return _num0;
 }
 
 
@@ -1870,7 +1875,6 @@ _numeric_digits::operator long()
 /** Append global Digits object to ostream. */
 void _numeric_digits::print(std::ostream &os) const
 {
-       debugmsg("_numeric_digits print", LOGLEVEL_PRINT);
        os << digits;
 }