X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=a3872b96650b9843505126cce338757ade54169a;hp=a368491ad55c69363fa932daf9135087d712e01f;hb=725021581cc862520c1f04b253ecb86f28032f69;hpb=ec397a253020039a0b0ea4ec5e91dde24100366a

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index a368491a..a3872b96 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -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
@@ -90,10 +90,10 @@ numeric::numeric(int i) : basic(TINFO_numeric)
 {
 	// 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);
@@ -105,10 +105,10 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 {
 	// 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);
@@ -168,8 +168,8 @@ numeric::numeric(const char *s) : basic(TINFO_numeric)
 
 	// 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 = ss.replace(delim,1,'E');
+	while ((delim = ss.find("e"))!=std::string::npos)
+		ss.replace(delim,1,"E");
 
 	// main parser loop:
 	do {
@@ -184,12 +184,12 @@ numeric::numeric(const char *s) : basic(TINFO_numeric)
 		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';
@@ -206,10 +206,10 @@ 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');
+			term = term.replace(term.find("E"),1,"e");
 			// append _<Digits> to term
 			term += "_" + ToString((unsigned)Digits);
 			// construct float using cln::cl_F(const char *) ctor.
@@ -322,8 +322,10 @@ static void print_real_number(const print_context & c, const cln::cl_R &x)
 		    !is_a<print_latex>(c)) {
 			cln::print_real(c.s, ourflags, x);
 		} else {  // rational output in LaTeX context
+			if (x < 0)
+				c.s << "-";
 			c.s << "\\frac{";
-			cln::print_real(c.s, ourflags, cln::numerator(cln::the<cln::cl_RA>(x)));
+			cln::print_real(c.s, ourflags, cln::abs(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 << '}';
@@ -354,7 +356,17 @@ void numeric::print(const print_context & c, unsigned level) const
 
 		std::ios::fmtflags oldflags = c.s.flags();
 		c.s.setf(std::ios::scientific);
-		if (this->is_rational() && !this->is_integer()) {
+		int oldprec = c.s.precision();
+		if (is_a<print_csrc_double>(c))
+			c.s.precision(16);
+		else
+			c.s.precision(7);
+		if (is_a<print_csrc_cl_N>(c) && this->is_integer()) {
+			c.s << "cln::cl_I(\"";
+			const cln::cl_R r = cln::realpart(cln::the<cln::cl_N>(value));
+			print_real_number(c,r);
+			c.s << "\")";
+		} else if (this->is_rational() && !this->is_integer()) {
 			if (compare(_num0) > 0) {
 				c.s << "(";
 				if (is_a<print_csrc_cl_N>(c))
@@ -376,11 +388,12 @@ void numeric::print(const print_context & c, unsigned level) const
 			c.s << ")";
 		} else {
 			if (is_a<print_csrc_cl_N>(c))
-				c.s << "cln::cl_F(\"" << evalf() << "\")";
+				c.s << "cln::cl_F(\"" << evalf() << "_" << Digits << "\")";
 			else
 				c.s << to_double();
 		}
 		c.s.flags(oldflags);
+		c.s.precision(oldprec);
 
 	} else {
 		const std::string par_open  = is_a<print_latex>(c) ? "{(" : "(";
@@ -389,6 +402,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())) {
@@ -401,25 +416,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
@@ -446,6 +455,8 @@ void numeric::print(const print_context & c, unsigned level) const
 					c.s << par_close;
 			}
 		}
+		if (is_a<print_python_repr>(c))
+			c.s << "')";
 	}
 }
 
@@ -494,6 +505,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
@@ -1484,7 +1510,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
@@ -1508,46 +1534,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;
+	}
+	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));
 	}
-	return results[n/2];
+	next_r = n+2;
+	return results[n/2-1];
 }
 
 
@@ -1654,9 +1695,12 @@ 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()));
@@ -1671,9 +1715,12 @@ 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()));
@@ -1689,9 +1736,12 @@ 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()));
@@ -1705,9 +1755,12 @@ 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()));