[GiNaC-devel] [PATCH] [BUGFIX]: Any complex number can be (un)archived properly.
Alexei Sheplyakov
varg at theor.jinr.ru
Wed Jul 30 18:32:46 CEST 2008
A complex number can have an exact real part and inexact (floating point)
imaginary part, and vice a versa. Handle these cases properly in the archiving
code instead of bailing out with a bizzare error message.
Thanks to Chris Bouchard for reporting the bug.
NOTE: this fix changes the format of GiNaC archives, so formally it breaks
the binary compatibility. However, "mixed" complex numbers (i.e. ones with
exact real part and inexact imaginary part) can not be archived with previous
versions of GiNaC, so there's nothing to break.
---
check/Makefile.am | 4 ++
check/numeric_archive.cpp | 48 ++++++++++++++++++
configure.ac | 4 +-
ginac/numeric.cpp | 118 ++++++++++++++++++++++++++++++++++-----------
4 files changed, 143 insertions(+), 31 deletions(-)
create mode 100644 check/numeric_archive.cpp
diff --git a/check/Makefile.am b/check/Makefile.am
index 8d2e451..a96dcfe 100644
--- a/check/Makefile.am
+++ b/check/Makefile.am
@@ -7,6 +7,7 @@ CHECKS = check_numeric \
EXAMS = exam_paranoia \
exam_heur_gcd \
+ exam_numeric_archive \
exam_numeric \
exam_powerlaws \
exam_inifcns \
@@ -71,6 +72,9 @@ exam_paranoia_LDADD = ../ginac/libginac.la
exam_heur_gcd_SOURCES = heur_gcd_bug.cpp
exam_heur_gcd_LDADD = ../ginac/libginac.la
+exam_numeric_archive_SOURCES = numeric_archive.cpp
+exam_numeric_archive_LDADD = ../ginac/libginac.la
+
exam_numeric_SOURCES = exam_numeric.cpp
exam_numeric_LDADD = ../ginac/libginac.la
diff --git a/check/numeric_archive.cpp b/check/numeric_archive.cpp
new file mode 100644
index 0000000..8f69a7b
--- /dev/null
+++ b/check/numeric_archive.cpp
@@ -0,0 +1,48 @@
+/**
+ * @file numeric_archive.cpp Check for a bug in numeric::archive
+ *
+ * numeric::archive used to fail if the real part of a complex number
+ * is a rational number and the imaginary part is a floating point one.
+ */
+#include <cln/cln.h>
+#include <iostream>
+#include <stdexcept>
+#include <vector>
+#include <iterator>
+#include <algorithm>
+#include <sstream>
+#include "ginac.h"
+using namespace cln;
+using namespace GiNaC;
+
+struct archive_unarchive_check
+{
+ cl_N operator()(const cl_N& n) const
+ {
+ ex e = numeric(n);
+ archive ar;
+ ar.archive_ex(e, "test");
+ lst l;
+ ex check = ar.unarchive_ex(l, "test");
+ if (!check.is_equal(e)) {
+ std::ostringstream s;
+ s << __FILE__ << ':' << __LINE__ << ": expected: " << e << ", got " << check;
+ throw std::logic_error(s.str());
+ }
+ return n;
+ }
+};
+
+int main(int argc, char** argv)
+{
+ const cl_I one(1);
+ std::cout << "checking if numeric::archive handles complex numbers properly" << std::endl;
+ const cl_R three_fp = cl_float(3.0, default_float_format);
+ std::vector<cl_N> numbers;
+ numbers.push_back(complex(one, three_fp));
+ numbers.push_back(complex(three_fp, one));
+ numbers.push_back(complex(three_fp, three_fp));
+ numbers.push_back(complex(one, one));
+ std::for_each(numbers.begin(), numbers.end(), archive_unarchive_check());
+ return 0;
+}
diff --git a/configure.ac b/configure.ac
index a0adef0..2518057 100644
--- a/configure.ac
+++ b/configure.ac
@@ -44,8 +44,8 @@ dnl The version number in newly created archives will be ARCHIVE_VERSION.
dnl Archives version (ARCHIVE_VERSION-ARCHIVE_AGE) thru ARCHIVE_VERSION can
dnl be read by this version of the GiNaC library.
-ARCHIVE_VERSION=2
-ARCHIVE_AGE=2
+ARCHIVE_VERSION=3
+ARCHIVE_AGE=3
AC_SUBST(ARCHIVE_VERSION)
AC_SUBST(ARCHIVE_AGE)
diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 502b0ea..32839e7 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -255,60 +255,120 @@ 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;
+}
+numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(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 '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;
- 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 '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: // Ordinary number
+ 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);
}
+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());
}
--
1.5.6
Best regards,
Alexei
--
All science is either physics or stamp collecting.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-devel/attachments/20080730/56432bf9/attachment.sig
More information about the GiNaC-devel
mailing list