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.