X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=e621607a730cabb8878cd5b5cf7d2c0b0dd71b54;hp=ef4e3e929ce5732f854bfae80d3f72f9f55e2b34;hb=c86cff51ac5a42f86387ef7bc767f1274137350b;hpb=bfb395b4e24e659fc46abc2d2ad4305168d1c6f2

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index ef4e3e92..e621607a 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -29,20 +29,13 @@
 #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,7 +88,6 @@ 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, 
@@ -113,7 +103,6 @@ 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, 
@@ -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,19 +297,13 @@ 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)
 
 //////////
-// functions overriding virtual functions from bases classes
+// functions overriding virtual functions from base classes
 //////////
 
 /** Helper function to print a real number in a nicer way than is CLN's
@@ -346,20 +313,28 @@ DEFAULT_UNARCHIVE(numeric)
  *  want to visibly distinguish from cl_LF.
  *
  *  @see numeric::print() */
-static void print_real_number(std::ostream &os, const cln::cl_R &num)
+static void print_real_number(const print_context & c, const cln::cl_R &x)
 {
 	cln::cl_print_flags ourflags;
-	if (cln::instanceof(num, cln::cl_RA_ring)) {
-		// case 1: integer or rational, nothing special to do:
-		cln::print_real(os, ourflags, num);
+	if (cln::instanceof(x, cln::cl_RA_ring)) {
+		// case 1: integer or rational
+		if (cln::instanceof(x, cln::cl_I_ring) ||
+		    !is_a<print_latex>(c)) {
+			cln::print_real(c.s, ourflags, x);
+		} else {  // rational output in LaTeX context
+			c.s << "\\frac{";
+			cln::print_real(c.s, ourflags, cln::numerator(cln::the<cln::cl_RA>(x)));
+			c.s << "}{";
+			cln::print_real(c.s, ourflags, cln::denominator(cln::the<cln::cl_RA>(x)));
+			c.s << '}';
+		}
 	} else {
 		// case 2: float
 		// make CLN believe this number has default_float_format, so it prints
 		// 'E' as exponent marker instead of 'L':
-		ourflags.default_float_format = cln::float_format(cln::the<cln::cl_F>(num));
-		cln::print_real(os, ourflags, num);
+		ourflags.default_float_format = cln::float_format(cln::the<cln::cl_F>(x));
+		cln::print_real(c.s, ourflags, x);
 	}
-	return;
 }
 
 /** This method adds to the output so it blends more consistently together
@@ -368,41 +343,39 @@ static void print_real_number(std::ostream &os, const cln::cl_R &num)
  *  @see print_real_number() */
 void numeric::print(const print_context & c, unsigned level) const
 {
-	debugmsg("numeric print", LOGLEVEL_PRINT);
-
-	if (is_of_type(c, print_tree)) {
+	if (is_a<print_tree>(c)) {
 
 		c.s << std::string(level, ' ') << cln::the<cln::cl_N>(value)
 		    << " (" << class_name() << ")"
 		    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 		    << std::endl;
 
-	} else if (is_of_type(c, print_csrc)) {
+	} else if (is_a<print_csrc>(c)) {
 
 		std::ios::fmtflags oldflags = c.s.flags();
 		c.s.setf(std::ios::scientific);
 		if (this->is_rational() && !this->is_integer()) {
-			if (compare(_num0()) > 0) {
+			if (compare(_num0) > 0) {
 				c.s << "(";
-				if (is_of_type(c, print_csrc_cl_N))
+				if (is_a<print_csrc_cl_N>(c))
 					c.s << "cln::cl_F(\"" << numer().evalf() << "\")";
 				else
 					c.s << numer().to_double();
 			} else {
 				c.s << "-(";
-				if (is_of_type(c, print_csrc_cl_N))
+				if (is_a<print_csrc_cl_N>(c))
 					c.s << "cln::cl_F(\"" << -numer().evalf() << "\")";
 				else
 					c.s << -numer().to_double();
 			}
 			c.s << "/";
-			if (is_of_type(c, print_csrc_cl_N))
+			if (is_a<print_csrc_cl_N>(c))
 				c.s << "cln::cl_F(\"" << denom().evalf() << "\")";
 			else
 				c.s << denom().to_double();
 			c.s << ")";
 		} else {
-			if (is_of_type(c, print_csrc_cl_N))
+			if (is_a<print_csrc_cl_N>(c))
 				c.s << "cln::cl_F(\"" << evalf() << "\")";
 			else
 				c.s << to_double();
@@ -410,30 +383,32 @@ void numeric::print(const print_context & c, unsigned level) const
 		c.s.flags(oldflags);
 
 	} else {
-		const std::string par_open  = is_of_type(c, print_latex) ? "{(" : "(";
-		const std::string par_close = is_of_type(c, print_latex) ? ")}" : ")";
-		const std::string imag_sym  = is_of_type(c, print_latex) ? "i" : "I";
-		const std::string mul_sym   = is_of_type(c, print_latex) ? " " : "*";
+		const std::string par_open  = is_a<print_latex>(c) ? "{(" : "(";
+		const std::string par_close = is_a<print_latex>(c) ? ")}" : ")";
+		const std::string imag_sym  = is_a<print_latex>(c) ? "i" : "I";
+		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())) {
+			if ((precedence() <= level) && (!this->is_nonneg_integer())) {
 				c.s << par_open;
-				print_real_number(c.s, r);
+				print_real_number(c, r);
 				c.s << par_close;
 			} else {
-				print_real_number(c.s, r);
+				print_real_number(c, r);
 			}
 		} else {
 			if (cln::zerop(r)) {
 				// case 2, imaginary:  y*I  or  -y*I
-				if ((precedence <= level) && (i < 0)) {
+				if ((precedence() <= level) && (i < 0)) {
 					if (i == -1) {
 						c.s << par_open+imag_sym+par_close;
 					} else {
 						c.s << par_open;
-						print_real_number(c.s, i);
+						print_real_number(c, i);
 						c.s << mul_sym+imag_sym+par_close;
 					}
 				} else {
@@ -443,21 +418,21 @@ void numeric::print(const print_context & c, unsigned level) const
 						if (i == -1) {
 							c.s << "-" << imag_sym;
 						} else {
-							print_real_number(c.s, i);
+							print_real_number(c, i);
 							c.s << mul_sym+imag_sym;
 						}
 					}
 				}
 			} else {
 				// case 3, complex:  x+y*I  or  x-y*I  or  -x+y*I  or  -x-y*I
-				if (precedence <= level)
+				if (precedence() <= level)
 					c.s << par_open;
-				print_real_number(c.s, r);
+				print_real_number(c, r);
 				if (i < 0) {
 					if (i == -1) {
 						c.s << "-"+imag_sym;
 					} else {
-						print_real_number(c.s, i);
+						print_real_number(c, i);
 						c.s << mul_sym+imag_sym;
 					}
 				} else {
@@ -465,14 +440,16 @@ void numeric::print(const print_context & c, unsigned level) const
 						c.s << "+"+imag_sym;
 					} else {
 						c.s << "+";
-						print_real_number(c.s, i);
+						print_real_number(c, i);
 						c.s << mul_sym+imag_sym;
 					}
 				}
-				if (precedence <= level)
+				if (precedence() <= level)
 					c.s << par_close;
 			}
 		}
+		if (is_a<print_python_repr>(c))
+			c.s << "')";
 	}
 }
 
@@ -529,9 +506,9 @@ bool numeric::info(unsigned inf) const
  *  sign as a multiplicative factor. */
 bool numeric::has(const ex &other) const
 {
-	if (!is_exactly_of_type(*other.bp, numeric))
+	if (!is_ex_exactly_of_type(other, numeric))
 		return false;
-	const numeric &o = static_cast<numeric &>(const_cast<basic &>(*other.bp));
+	const numeric &o = ex_to<numeric>(other);
 	if (this->is_equal(o) || this->is_equal(-o))
 		return true;
 	if (o.imag().is_zero())  // e.g. scan for 3 in -3*I
@@ -575,8 +552,8 @@ ex numeric::evalf(int level) const
 
 int numeric::compare_same_type(const basic &other) const
 {
-	GINAC_ASSERT(is_exactly_of_type(other, numeric));
-	const numeric &o = static_cast<numeric &>(const_cast<basic &>(other));
+	GINAC_ASSERT(is_exactly_a<numeric>(other));
+	const numeric &o = static_cast<const numeric &>(other);
 	
 	return this->compare(o);
 }
@@ -584,10 +561,10 @@ 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));
-	const numeric *o = static_cast<const numeric *>(&other);
+	GINAC_ASSERT(is_exactly_a<numeric>(other));
+	const numeric &o = static_cast<const numeric &>(other);
 	
-	return this->is_equal(*o);
+	return this->is_equal(o);
 }
 
 
@@ -618,10 +595,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));
@@ -641,10 +617,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));
@@ -668,8 +643,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))) {
@@ -680,7 +654,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)));
 }
@@ -689,10 +663,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)))->
@@ -710,10 +683,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)))->
@@ -733,8 +705,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))) {
@@ -745,7 +716,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));
@@ -1112,7 +1083,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)));
@@ -1121,7 +1092,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))
@@ -1130,7 +1101,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;
 }
 
 
@@ -1148,15 +1119,6 @@ int numeric::int_length(void) const
 		return 0;
 }
 
-
-//////////
-// static member variables
-//////////
-
-// protected
-
-unsigned numeric::precedence = 30;
-
 //////////
 // global constants
 //////////
@@ -1243,7 +1205,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());
 }
@@ -1394,7 +1356,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
@@ -1433,10 +1395,7 @@ const numeric zeta(const numeric &x)
 		if (cln::zerop(x.to_cl_N()-aux))
 			return cln::zeta(aux);
 	}
-	std::clog << "zeta(" << x
-			  << "): Does anybody know a good way to calculate this numerically?"
-			  << std::endl;
-	return numeric(0);
+	throw dunno();
 }
 
 
@@ -1444,17 +1403,11 @@ const numeric zeta(const numeric &x)
  *  This is only a stub! */
 const numeric lgamma(const numeric &x)
 {
-	std::clog << "lgamma(" << x
-	          << "): Does anybody know a good way to calculate this numerically?"
-	          << std::endl;
-	return numeric(0);
+	throw dunno();
 }
 const numeric tgamma(const numeric &x)
 {
-	std::clog << "tgamma(" << x
-	          << "): Does anybody know a good way to calculate this numerically?"
-	          << std::endl;
-	return numeric(0);
+	throw dunno();
 }
 
 
@@ -1462,10 +1415,7 @@ const numeric tgamma(const numeric &x)
  *  This is only a stub! */
 const numeric psi(const numeric &x)
 {
-	std::clog << "psi(" << x
-	          << "): Does anybody know a good way to calculate this numerically?"
-	          << std::endl;
-	return numeric(0);
+	throw dunno();
 }
 
 
@@ -1473,10 +1423,7 @@ const numeric psi(const numeric &x)
  *  This is only a stub! */
 const numeric psi(const numeric &n, const numeric &x)
 {
-	std::clog << "psi(" << n << "," << x
-	          << "): Does anybody know a good way to calculate this numerically?"
-	          << std::endl;
-	return numeric(0);
+	throw dunno();
 }
 
 
@@ -1500,8 +1447,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");
@@ -1518,12 +1465,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);
 		}
 	}
 	
@@ -1541,7 +1488,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
@@ -1565,46 +1512,49 @@ 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.size()==0)
-		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;
+	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;
+	}
+	for (unsigned p=next_r; p<=n;  p+=2) {
+		cln::cl_I  c = 1;
+		cln::cl_RA b = cln::cl_RA(1-p)/2;
+		const unsigned p3 = p+3;
+		const unsigned p2 = p+2;
+		const unsigned pm = p-2;
+		unsigned i, k;
+		for (i=2, k=0; i <= pm; i += 2, k++) {
+			c = cln::exquo(c * ((p3 - i)*(p2 - i)), (i - 1)*i);
+			b = b + c * results[k];
 		}
-		B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2));
-		results.push_back(B);
-		++highest_result;
+		results.push_back(-b / (p+1));
+		next_r += 2;
 	}
-	return results[n/2];
+	return results[n/2 - 1];
 }
 
 
@@ -1635,7 +1585,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);
@@ -1687,7 +1637,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;
 }
 
 
@@ -1702,7 +1652,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;
 }
 
 
@@ -1718,7 +1668,7 @@ const numeric irem(const numeric &a, const numeric &b)
 		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;
 }
 
 
@@ -1737,8 +1687,8 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q)
 		q = rem_quo.quotient;
 		return rem_quo.remainder;
 	} else {
-		q = _num0();
-		return _num0();
+		q = _num0;
+		return _num0;
 	}
 }
 
@@ -1753,7 +1703,7 @@ const numeric iquo(const numeric &a, const numeric &b)
 		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;
 }
 
 
@@ -1771,8 +1721,8 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r)
 		r = rem_quo.remainder;
 		return rem_quo.quotient;
 	} else {
-		r = _num0();
-		return _num0();
+		r = _num0;
+		return _num0;
 	}
 }
 
@@ -1787,7 +1737,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;
 }
 
 
@@ -1827,7 +1777,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;
 }
 
 
@@ -1886,7 +1836,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;
 }