X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=1d4de9f5b416d4e997f6d590f0c96ef7edb57e4c;hb=e98841136efa88c951edafc0cd43ba1343f20b5b;hp=5c8f65061819933dd37751ccf89cf5281628a67b;hpb=3b73d3c2d2da6c37895aab34f4eeffe4a156141e;p=ginac.git

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 5c8f6506..1d4de9f5 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -30,6 +30,7 @@
 #include <stdexcept>
 #include <string>
 #include <sstream>
+#include <limits>
 
 #include "numeric.h"
 #include "ex.h"
@@ -94,10 +95,10 @@ numeric::numeric(int i) : basic(TINFO_numeric)
 	// 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))
+	if (i < (1L << (cl_value_len-1)) && i >= -(1L << (cl_value_len-1)))
 		value = cln::cl_I(i);
 	else
-		value = cln::cl_I((long) i);
+		value = cln::cl_I(static_cast<long>(i));
 	setflag(status_flags::evaluated | status_flags::expanded);
 }
 
@@ -109,10 +110,10 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 	// 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))
+	if (i < (1U << (cl_value_len-1)))
 		value = cln::cl_I(i);
 	else
-		value = cln::cl_I((unsigned long) i);
+		value = cln::cl_I(static_cast<unsigned long>(i));
 	setflag(status_flags::evaluated | status_flags::expanded);
 }
 
@@ -130,7 +131,8 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric)
 	setflag(status_flags::evaluated | status_flags::expanded);
 }
 
-/** Ctor for rational numerics a/b.
+
+/** Constructor for rational numerics a/b.
  *
  *  @exception overflow_error (division by zero) */
 numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
@@ -243,7 +245,7 @@ numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric)
 // archiving
 //////////
 
-numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 {
 	cln::cl_N ctorval = 0;
 
@@ -314,7 +316,7 @@ DEFAULT_UNARCHIVE(numeric)
  *  want to visibly distinguish from cl_LF.
  *
  *  @see numeric::print() */
-static void print_real_number(const print_context & c, const cln::cl_R &x)
+static void print_real_number(const print_context & c, const cln::cl_R & x)
 {
 	cln::cl_print_flags ourflags;
 	if (cln::instanceof(x, cln::cl_RA_ring)) {
@@ -340,6 +342,82 @@ static void print_real_number(const print_context & c, const cln::cl_R &x)
 	}
 }
 
+/** Helper function to print integer number in C++ source format.
+ *
+ *  @see numeric::print() */
+static void print_integer_csrc(const print_context & c, const cln::cl_I & x)
+{
+	// Print small numbers in compact float format, but larger numbers in
+	// scientific format
+	const int max_cln_int = 536870911; // 2^29-1
+	if (x >= cln::cl_I(-max_cln_int) && x <= cln::cl_I(max_cln_int))
+		c.s << cln::cl_I_to_int(x) << ".0";
+	else
+		c.s << cln::double_approx(x);
+}
+
+/** Helper function to print real number in C++ source format.
+ *
+ *  @see numeric::print() */
+static void print_real_csrc(const print_context & c, const cln::cl_R & x)
+{
+	if (cln::instanceof(x, cln::cl_I_ring)) {
+
+		// Integer number
+		print_integer_csrc(c, cln::the<cln::cl_I>(x));
+
+	} else if (cln::instanceof(x, cln::cl_RA_ring)) {
+
+		// Rational number
+		const cln::cl_I numer = cln::numerator(cln::the<cln::cl_RA>(x));
+		const cln::cl_I denom = cln::denominator(cln::the<cln::cl_RA>(x));
+		if (cln::plusp(x) > 0) {
+			c.s << "(";
+			print_integer_csrc(c, numer);
+		} else {
+			c.s << "-(";
+			print_integer_csrc(c, -numer);
+		}
+		c.s << "/";
+		print_integer_csrc(c, denom);
+		c.s << ")";
+
+	} else {
+
+		// Anything else
+		c.s << cln::double_approx(x);
+	}
+}
+
+/** Helper function to print real number in C++ source format using cl_N types.
+ *
+ *  @see numeric::print() */
+static void print_real_cl_N(const print_context & c, const cln::cl_R & x)
+{
+	if (cln::instanceof(x, cln::cl_I_ring)) {
+
+		// Integer number
+		c.s << "cln::cl_I(\"";
+		print_real_number(c, x);
+		c.s << "\")";
+
+	} else if (cln::instanceof(x, cln::cl_RA_ring)) {
+
+		// Rational number
+		cln::cl_print_flags ourflags;
+		c.s << "cln::cl_RA(\"";
+		cln::print_rational(c.s, ourflags, cln::the<cln::cl_RA>(x));
+		c.s << "\")";
+
+	} else {
+
+		// Anything else
+		c.s << "cln::cl_F(\"";
+		print_real_number(c, cln::cl_float(1.0, cln::default_float_format) * x);
+		c.s << "_" << Digits << "\")";
+	}
+}
+
 /** This method adds to the output so it blends more consistently together
  *  with the other routines and produces something compatible to ginsh input.
  *  
@@ -353,56 +431,69 @@ void numeric::print(const print_context & c, unsigned level) const
 		    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 		    << std::endl;
 
+	} else if (is_a<print_csrc_cl_N>(c)) {
+
+		// CLN output
+		if (this->is_real()) {
+
+			// Real number
+			print_real_cl_N(c, cln::the<cln::cl_R>(value));
+
+		} else {
+
+			// Complex number
+			c.s << "cln::complex(";
+			print_real_cl_N(c, cln::realpart(cln::the<cln::cl_N>(value)));
+			c.s << ",";
+			print_real_cl_N(c, cln::imagpart(cln::the<cln::cl_N>(value)));
+			c.s << ")";
+		}
+
 	} else if (is_a<print_csrc>(c)) {
 
+		// C++ source output
 		std::ios::fmtflags oldflags = c.s.flags();
 		c.s.setf(std::ios::scientific);
 		int oldprec = c.s.precision();
+
+		// Set precision
 		if (is_a<print_csrc_double>(c))
-			c.s.precision(16);
+			c.s.precision(std::numeric_limits<double>::digits10 + 1);
 		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))
-					c.s << "cln::cl_F(\"" << numer().evalf() << "\")";
-				else
-					c.s << numer().to_double();
-			} else {
-				c.s << "-(";
-				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_a<print_csrc_cl_N>(c))
-				c.s << "cln::cl_F(\"" << denom().evalf() << "\")";
-			else
-				c.s << denom().to_double();
-			c.s << ")";
+			c.s.precision(std::numeric_limits<float>::digits10 + 1);
+
+		if (this->is_real()) {
+
+			// Real number
+			print_real_csrc(c, cln::the<cln::cl_R>(value));
+
 		} else {
-			if (is_a<print_csrc_cl_N>(c))
-				c.s << "cln::cl_F(\"" << evalf() << "_" << Digits << "\")";
+
+			// Complex number
+			c.s << "std::complex<";
+			if (is_a<print_csrc_double>(c))
+				c.s << "double>(";
 			else
-				c.s << to_double();
+				c.s << "float>(";
+
+			print_real_csrc(c, cln::realpart(cln::the<cln::cl_N>(value)));
+			c.s << ",";
+			print_real_csrc(c, cln::imagpart(cln::the<cln::cl_N>(value)));
+			c.s << ")";
 		}
+
 		c.s.flags(oldflags);
 		c.s.precision(oldprec);
 
 	} else {
+
 		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)) {
@@ -667,8 +758,9 @@ const numeric numeric::div(const numeric &other) const
  *  returns result as a numeric object. */
 const numeric numeric::power(const numeric &other) const
 {
-	// Efficiency shortcut: trap the neutral exponent by pointer.
-	if (&other==_num1_p)
+	// Shortcut for efficiency and numeric stability (as in 1.0 exponent):
+	// trap the neutral exponent.
+	if (&other==_num1_p || cln::equal(cln::the<cln::cl_N>(other.value),cln::the<cln::cl_N>(_num1.value)))
 		return *this;
 	
 	if (cln::zerop(cln::the<cln::cl_N>(value))) {
@@ -691,7 +783,8 @@ const numeric numeric::power(const numeric &other) const
  *  an ex object, where the result would end up on the heap anyways. */
 const numeric &numeric::add_dyn(const numeric &other) const
 {
-	// Efficiency shortcut: trap the neutral element by pointer.
+	// Efficiency shortcut: trap the neutral element by pointer.  This hack
+	// is supposed to keep the number of distinct numeric objects low.
 	if (this==_num0_p)
 		return other;
 	else if (&other==_num0_p)
@@ -751,8 +844,10 @@ const numeric &numeric::div_dyn(const numeric &other) const
  *  heap anyways. */
 const numeric &numeric::power_dyn(const numeric &other) const
 {
-	// Efficiency shortcut: trap the neutral exponent by pointer.
-	if (&other==_num1_p)
+	// Efficiency shortcut: trap the neutral exponent (first try by pointer, then
+	// try harder, since calls to cln::expt() below may return amazing results for
+	// floating point exponent 1.0).
+	if (&other==_num1_p || cln::equal(cln::the<cln::cl_N>(other.value),cln::the<cln::cl_N>(_num1.value)))
 		return *this;
 	
 	if (cln::zerop(cln::the<cln::cl_N>(value))) {