X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=4c11ac4bf63dd1d4010794a87de82c0078cc20c9;hp=69a787b5d985b5e174ea61a6f38c7b89f2bb3a2c;hb=866b3eb23253a272788d8791b1ed023e63674d50;hpb=a7693a0f710b49494f95ce5a4a0953752e69c7f9

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 69a787b5..4c11ac4b 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-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2010 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
@@ -26,12 +26,6 @@
 
 #include "config.h"
 
-#include <vector>
-#include <stdexcept>
-#include <string>
-#include <sstream>
-#include <limits>
-
 #include "numeric.h"
 #include "ex.h"
 #include "operators.h"
@@ -39,6 +33,12 @@
 #include "tostring.h"
 #include "utils.h"
 
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
 // CLN should pollute the global namespace as little as possible.  Hence, we
 // include most of it here and include only the part needed for properly
 // declaring cln::cl_number in numeric.h.  This can only be safely done in
@@ -73,7 +73,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(numeric, basic,
 //////////
 
 /** default ctor. Numerically it initializes to an integer zero. */
-numeric::numeric() : basic(&numeric::tinfo_static)
+numeric::numeric()
 {
 	value = cln::cl_I(0);
 	setflag(status_flags::evaluated | status_flags::expanded);
@@ -85,7 +85,7 @@ numeric::numeric() : basic(&numeric::tinfo_static)
 
 // public
 
-numeric::numeric(int i) : basic(&numeric::tinfo_static)
+numeric::numeric(int i)
 {
 	// 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
@@ -106,7 +106,7 @@ numeric::numeric(int i) : basic(&numeric::tinfo_static)
 }
 
 
-numeric::numeric(unsigned int i) : basic(&numeric::tinfo_static)
+numeric::numeric(unsigned int i)
 {
 	// 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
@@ -127,14 +127,14 @@ numeric::numeric(unsigned int i) : basic(&numeric::tinfo_static)
 }
 
 
-numeric::numeric(long i) : basic(&numeric::tinfo_static)
+numeric::numeric(long i)
 {
 	value = cln::cl_I(i);
 	setflag(status_flags::evaluated | status_flags::expanded);
 }
 
 
-numeric::numeric(unsigned long i) : basic(&numeric::tinfo_static)
+numeric::numeric(unsigned long i)
 {
 	value = cln::cl_I(i);
 	setflag(status_flags::evaluated | status_flags::expanded);
@@ -144,7 +144,7 @@ numeric::numeric(unsigned long i) : basic(&numeric::tinfo_static)
 /** Constructor for rational numerics a/b.
  *
  *  @exception overflow_error (division by zero) */
-numeric::numeric(long numer, long denom) : basic(&numeric::tinfo_static)
+numeric::numeric(long numer, long denom)
 {
 	if (!denom)
 		throw std::overflow_error("division by zero");
@@ -153,7 +153,7 @@ numeric::numeric(long numer, long denom) : basic(&numeric::tinfo_static)
 }
 
 
-numeric::numeric(double d) : basic(&numeric::tinfo_static)
+numeric::numeric(double d)
 {
 	// 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
@@ -165,7 +165,7 @@ numeric::numeric(double d) : basic(&numeric::tinfo_static)
 
 /** ctor from C-style string.  It also accepts complex numbers in GiNaC
  *  notation like "2+5*I". */
-numeric::numeric(const char *s) : basic(&numeric::tinfo_static)
+numeric::numeric(const char *s)
 {
 	cln::cl_N ctorval = 0;
 	// parse complex numbers (functional but not completely safe, unfortunately
@@ -244,7 +244,7 @@ numeric::numeric(const char *s) : basic(&numeric::tinfo_static)
 
 /** Ctor from CLN types.  This is for the initiated user or internal use
  *  only. */
-numeric::numeric(const cln::cl_N &z) : basic(&numeric::tinfo_static)
+numeric::numeric(const cln::cl_N &z)
 {
 	value = z;
 	setflag(status_flags::evaluated | status_flags::expanded);
@@ -255,66 +255,126 @@ numeric::numeric(const cln::cl_N &z) : basic(&numeric::tinfo_static)
 // archiving
 //////////
 
-numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+/** 
+ * Construct a floating point number from sign, mantissa, and exponent 
+ */
+static const cln::cl_F make_real_float(const cln::cl_idecoded_float& dec)
 {
-	cln::cl_N ctorval = 0;
+	cln::cl_F x = cln::cl_float(dec.mantissa, cln::default_float_format);
+	x = cln::scale_float(x, dec.exponent);
+	cln::cl_F sign = cln::cl_float(dec.sign, cln::default_float_format);
+	x = cln::float_sign(sign, x);
+	return x;
+}
+
+/** 
+ * Read serialized floating point number 
+ */
+static const cln::cl_F read_real_float(std::istream& s)
+{
+	cln::cl_idecoded_float dec;
+	s >> dec.sign >> dec.mantissa >> dec.exponent;
+	const cln::cl_F x = make_real_float(dec);
+	return x;
+}
 
+void numeric::read_archive(const archive_node &n, lst &sym_lst)
+{
+	inherited::read_archive(n, sym_lst);
+	value = 0;
+	
 	// Read number as string
 	std::string str;
 	if (n.find_string("number", str)) {
 		std::istringstream s(str);
-		cln::cl_idecoded_float re, im;
+		cln::cl_R re, im;
 		char c;
 		s.get(c);
 		switch (c) {
-			case 'R':    // Integer-decoded real number
-				s >> re.sign >> re.mantissa >> re.exponent;
-				ctorval = re.sign * re.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), re.exponent);
+			case 'R':
+				// real FP (floating point) number
+				re = read_real_float(s);
+				value = re;
+				break;
+			case 'C':
+				// both real and imaginary part are FP numbers
+				re = read_real_float(s);
+				im = read_real_float(s); 
+				value = cln::complex(re, im);
 				break;
-			case 'C':    // Integer-decoded complex number
-				s >> re.sign >> re.mantissa >> re.exponent;
-				s >> im.sign >> im.mantissa >> im.exponent;
-				ctorval = cln::complex(re.sign * re.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), re.exponent),
-				                       im.sign * im.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), im.exponent));
+			case 'H':
+				// real part is a rational number,
+				// imaginary part is a FP number
+				s >> re;
+				im = read_real_float(s);
+				value = cln::complex(re, im);
 				break;
-			default:    // Ordinary number
+			case 'J':
+				// real part is a FP number,
+				// imaginary part is a rational number
+				re = read_real_float(s);
+				s >> im;
+				value = cln::complex(re, im);
+				break;
+			default:
+				// both real and imaginary parts are rational
 				s.putback(c);
-				s >> ctorval;
+				s >> value;
 				break;
 		}
 	}
-	value = ctorval;
 	setflag(status_flags::evaluated | status_flags::expanded);
 }
+GINAC_BIND_UNARCHIVER(numeric);
+
+static void write_real_float(std::ostream& s, const cln::cl_R& n)
+{
+	const cln::cl_idecoded_float dec = cln::integer_decode_float(cln::the<cln::cl_F>(n));
+	s << dec.sign << ' ' << dec.mantissa << ' ' << dec.exponent;
+}
 
 void numeric::archive(archive_node &n) const
 {
 	inherited::archive(n);
 
 	// Write number as string
+	
+	const cln::cl_R re = cln::realpart(value);
+	const cln::cl_R im = cln::imagpart(value);
+	const bool re_rationalp = cln::instanceof(re, cln::cl_RA_ring);
+	const bool im_rationalp = cln::instanceof(im, cln::cl_RA_ring);
+
+	// Non-rational numbers are written in an integer-decoded format
+	// to preserve the precision
 	std::ostringstream s;
-	if (this->is_crational())
+	if (re_rationalp && im_rationalp)
 		s << value;
-	else {
-		// Non-rational numbers are written in an integer-decoded format
-		// to preserve the precision
-		if (this->is_real()) {
-			cln::cl_idecoded_float re = cln::integer_decode_float(cln::the<cln::cl_F>(value));
-			s << "R";
-			s << re.sign << " " << re.mantissa << " " << re.exponent;
-		} else {
-			cln::cl_idecoded_float re = cln::integer_decode_float(cln::the<cln::cl_F>(cln::realpart(cln::the<cln::cl_N>(value))));
-			cln::cl_idecoded_float im = cln::integer_decode_float(cln::the<cln::cl_F>(cln::imagpart(cln::the<cln::cl_N>(value))));
-			s << "C";
-			s << re.sign << " " << re.mantissa << " " << re.exponent << " ";
-			s << im.sign << " " << im.mantissa << " " << im.exponent;
-		}
+	else if (zerop(im)) {
+		// real FP (floating point) number
+		s << 'R';
+		write_real_float(s, re);
+	} else if (re_rationalp) {
+		s << 'H'; // just any unique character
+		// real part is a rational number,
+		// imaginary part is a FP number
+		s << re << ' ';
+		write_real_float(s, im);
+	} else if (im_rationalp) {
+		s << 'J';
+		// real part is a FP number,
+		// imaginary part is a rational number
+		write_real_float(s, re);
+		s << ' ' << im;
+	} else  {
+		// both real and imaginary parts are floating point
+		s << 'C';
+		write_real_float(s, re);
+		s << ' ';
+		write_real_float(s, im);
 	}
 	n.add_string("number", s.str());
 }
 
-DEFAULT_UNARCHIVE(numeric)
-
 //////////
 // functions overriding virtual functions from base classes
 //////////
@@ -2218,11 +2278,14 @@ const numeric fibonacci(const numeric &n)
 	//      F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
 	if (n.is_zero())
 		return *_num0_p;
-	if (n.is_negative())
-		if (n.is_even())
+	if (n.is_negative()) {
+		if (n.is_even()) {
 			return -fibonacci(-n);
-		else
+		}
+		else {
 			return fibonacci(-n);
+		}
+	}
 	
 	cln::cl_I u(0);
 	cln::cl_I v(1);
@@ -2274,15 +2337,18 @@ const numeric mod(const numeric &a, const numeric &b)
 
 
 /** Modulus (in symmetric representation).
- *  Equivalent to Maple's mods.
  *
- *  @return a mod b in the range [-iquo(abs(b)-1,2), iquo(abs(b),2)]. */
-const numeric smod(const numeric &a, const numeric &b)
-{
-	if (a.is_integer() && b.is_integer()) {
-		const cln::cl_I b2 = cln::ceiling1(cln::the<cln::cl_I>(b.to_cl_N()) >> 1) - 1;
-		return numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()) + b2,
-		                cln::the<cln::cl_I>(b.to_cl_N())) - b2);
+ *  @return a mod b in the range [-iquo(abs(b),2), iquo(abs(b),2)]. */
+const numeric smod(const numeric &a_, const numeric &b_)
+{
+	if (a_.is_integer() && b_.is_integer()) {
+		const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
+		const cln::cl_I b = cln::the<cln::cl_I>(b_.to_cl_N());
+		const cln::cl_I b2 = b >> 1;
+		const cln::cl_I m = cln::mod(a, b);
+		const cln::cl_I m_b = m - b;
+		const cln::cl_I ret = m > b2 ? m_b : m;
+		return numeric(ret);
 	} else
 		return *_num0_p;
 }