X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=f27fd3aca60d7d1b090a59aa6562e3bb2bc33620;hp=1b98040314e22f679649fd0386c28e8838805d38;hb=c12c8ec3c5cf0c75f061f6c52d04206277bbdcca;hpb=a297391b22c1a4ec3464d13191f8cb33332c9863 diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 1b980403..f27fd3ac 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-2016 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 @@ -24,21 +24,22 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ +#ifdef HAVE_CONFIG_H #include "config.h" - -#include -#include -#include -#include -#include +#endif #include "numeric.h" #include "ex.h" #include "operators.h" #include "archive.h" -#include "tostring.h" #include "utils.h" +#include +#include +#include +#include +#include + // 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 +74,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 +86,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 +107,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 +128,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 +145,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 +154,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 +166,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 @@ -223,7 +224,7 @@ numeric::numeric(const char *s) : basic(&numeric::tinfo_static) // E to lower case term = term.replace(term.find("E"),1,"e"); // append _ to term - term += "_" + ToString((unsigned)Digits); + term += "_" + std::to_string((unsigned)Digits); // construct float using cln::cl_F(const char *) ctor. if (imaginary) ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_F(term.c_str())); @@ -244,7 +245,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 +256,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 '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); } +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(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(value)); - s << "R"; - s << re.sign << " " << re.mantissa << " " << re.exponent; - } else { - cln::cl_idecoded_float re = cln::integer_decode_float(cln::the(cln::realpart(cln::the(value)))); - cln::cl_idecoded_float im = cln::integer_decode_float(cln::the(cln::imagpart(cln::the(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 ////////// @@ -440,20 +501,20 @@ static void print_real_cl_N(const print_context & c, const cln::cl_R & x) { if (cln::instanceof(x, cln::cl_I_ring)) { - int dst; - // fixnum - if (coerce(dst, cln::the(x))) { - // can be converted to native int - if (dst < 0) - c.s << "(-" << dst << ")"; - else - c.s << dst; - } else { - // bignum - c.s << "cln::cl_I(\""; - print_real_number(c, x); - c.s << "\")"; - } + int dst; + // fixnum + if (coerce(dst, cln::the(x))) { + // can be converted to native int + if (dst < 0) + c.s << "(-" << dst << ")"; + else + c.s << dst; + } else { + // bignum + c.s << "cln::cl_I(\""; + print_real_number(c, x); + c.s << "\")"; + } } else if (cln::instanceof(x, cln::cl_RA_ring)) { // Rational number @@ -640,7 +701,7 @@ bool numeric::info(unsigned inf) const case info_flags::negative: return is_negative(); case info_flags::nonnegative: - return !is_negative(); + return is_zero() || is_positive(); case info_flags::posint: return is_pos_integer(); case info_flags::negint: @@ -653,8 +714,6 @@ bool numeric::info(unsigned inf) const return is_odd(); case info_flags::prime: return is_prime(); - case info_flags::algebraic: - return !is_real(); } return false; } @@ -714,10 +773,8 @@ bool numeric::has(const ex &other, unsigned options) const /** Evaluation of numbers doesn't do anything at all. */ -ex numeric::eval(int level) const +ex numeric::eval() const { - // Warning: if this is ever gonna do something, the ex ctors from all kinds - // of numbers should be checking for status_flags::evaluated. return this->hold(); } @@ -727,11 +784,9 @@ ex numeric::eval(int level) const * currently set. In case the object already was a floating point number the * precision is trimmed to match the currently set default. * - * @param level ignored, only needed for overriding basic::evalf. * @return an ex-handle to a numeric. */ -ex numeric::evalf(int level) const +ex numeric::evalf() const { - // level can safely be discarded for numeric objects. return numeric(cln::cl_float(1.0, cln::default_float_format) * value); } @@ -868,9 +923,8 @@ const numeric &numeric::add_dyn(const numeric &other) const return other; else if (&other==_num0_p) return *this; - - return static_cast((new numeric(value + other.value))-> - setflag(status_flags::dynallocated)); + + return dynallocate(value + other.value); } @@ -884,9 +938,8 @@ const numeric &numeric::sub_dyn(const numeric &other) const // hack is supposed to keep the number of distinct numeric objects low. if (&other==_num0_p || cln::zerop(other.value)) return *this; - - return static_cast((new numeric(value - other.value))-> - setflag(status_flags::dynallocated)); + + return dynallocate(value - other.value); } @@ -903,8 +956,7 @@ const numeric &numeric::mul_dyn(const numeric &other) const else if (&other==_num1_p) return *this; - return static_cast((new numeric(value * other.value))-> - setflag(status_flags::dynallocated)); + return dynallocate(value * other.value); } @@ -922,8 +974,8 @@ const numeric &numeric::div_dyn(const numeric &other) const return *this; if (cln::zerop(cln::the(other.value))) throw std::overflow_error("division by zero"); - return static_cast((new numeric(value / other.value))-> - setflag(status_flags::dynallocated)); + + return dynallocate(value / other.value); } @@ -949,8 +1001,8 @@ const numeric &numeric::power_dyn(const numeric &other) const else return *_num0_p; } - return static_cast((new numeric(cln::expt(value, other.value)))-> - setflag(status_flags::dynallocated)); + + return dynallocate(cln::expt(value, other.value)); } @@ -1314,7 +1366,7 @@ const numeric numeric::numer() const if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_RA_ring)) { const cln::cl_I s = cln::lcm(cln::denominator(r), cln::denominator(i)); return numeric(cln::complex(cln::numerator(r)*(cln::exquo(s,cln::denominator(r))), - cln::numerator(i)*(cln::exquo(s,cln::denominator(i))))); + cln::numerator(i)*(cln::exquo(s,cln::denominator(i))))); } } // at least one float encountered @@ -1469,7 +1521,7 @@ const numeric atan(const numeric &y, const numeric &x) return *_num0_p; if (x.is_real() && y.is_real()) return numeric(cln::atan(cln::the(x.to_cl_N()), - cln::the(y.to_cl_N()))); + cln::the(y.to_cl_N()))); // Compute -I*log((x+I*y)/sqrt(x^2+y^2)) // == -I*log((x+I*y)/sqrt((x+I*y)*(x-I*y))) @@ -1713,8 +1765,8 @@ cln::cl_N lanczos_coeffs::calc_lanczos_A(const cln::cl_N &x) const { cln::cl_N A = (*current_vector)[0]; int size = current_vector->size(); - for (int i=1; i(b.to_cl_N()) >> 1) - 1; - return numeric(cln::mod(cln::the(a.to_cl_N()) + b2, - cln::the(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(a_.to_cl_N()); + const cln::cl_I b = cln::the(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; } @@ -2471,9 +2529,8 @@ _numeric_digits& _numeric_digits::operator=(long prec) cln::default_float_format = cln::float_format(prec); // call registered callbacks - std::vector::const_iterator it = callbacklist.begin(), end = callbacklist.end(); - for (; it != end; ++it) { - (*it)(digitsdiff); + for (auto it : callbacklist) { + (it)(digitsdiff); } return *this;