Any complex number can be (un)archived properly.

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.

diff --git a/check/Makefile.am b/check/Makefile.am
index 8d2e451e29768b04973ae3341fcb837ce2fdddcd..a96dcfe7e798fa7950d9d8e1cd856031d236329d 100644 (file)
--- 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 (file)
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 b592220b4263f7236e133f630de975cefc054061..e4a1e60ab59f5a51567fd2efc3aa9f2e21d9ac5c 100644 (file)
--- 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 69a787b5d985b5e174ea61a6f38c7b89f2bb3a2c..49951470475de6ccf9de61b97b9f76e03a259d88 100644 (file)
--- 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());
 }