X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=e45fe9c98320dc3636154e93ed8d690e1b4e6a54;hp=99251af5342de885f18e43e421032933d0e1efc1;hb=cd22e73d44e3320898f62a0accdbbe005b33d3e5;hpb=58aaaf83c55c5784e1d94c2a794af96d7769165b diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 99251af5..e45fe9c9 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': // 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 '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; - default: // Ordinary number + 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 '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(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))) @@ -1611,24 +1663,24 @@ static cln::cl_N Li2_projection(const cln::cl_N &x, return Li2_series(x, prec); } + /** Numeric evaluation of Dilogarithm. The domain is the entire complex plane, * the branch cut lies along the positive real axis, starting at 1 and * continuous with quadrant IV. * * @return arbitrary precision numerical Li2(x). */ -const numeric Li2(const numeric &x) +const cln::cl_N Li2_(const cln::cl_N& value) { - if (x.is_zero()) - return *_num0_p; + if (zerop(value)) + return 0; // what is the desired float format? // first guess: default format cln::float_format_t prec = cln::default_float_format; - const cln::cl_N value = x.to_cl_N(); // second guess: the argument's format - if (!x.real().is_rational()) + if (!instanceof(realpart(value), cln::cl_RA_ring)) prec = cln::float_format(cln::the(cln::realpart(value))); - else if (!x.imag().is_rational()) + else if (!instanceof(imagpart(value), cln::cl_RA_ring)) prec = cln::float_format(cln::the(cln::imagpart(value))); if (value==1) // may cause trouble with log(1-x) @@ -1640,7 +1692,16 @@ const numeric Li2(const numeric &x) - cln::zeta(2, prec) - Li2_projection(cln::recip(value), prec)); else - return Li2_projection(x.to_cl_N(), prec); + return Li2_projection(value, prec); +} + +const numeric Li2(const numeric &x) +{ + const cln::cl_N x_ = x.to_cl_N(); + if (zerop(x_)) + return *_num0_p; + const cln::cl_N result = Li2_(x_); + return numeric(result); } @@ -1656,7 +1717,7 @@ const numeric zeta(const numeric &x) if (x.is_real()) { const int aux = (int)(cln::double_approx(cln::the(x.to_cl_N()))); if (cln::zerop(x.to_cl_N()-aux)) - return cln::zeta(aux); + return numeric(cln::zeta(aux)); } throw dunno(); } @@ -1678,7 +1739,7 @@ class lanczos_coeffs std::vector *current_vector; }; -std::vector* lanczos_coeffs::coeffs = 0; +std::vector* lanczos_coeffs::coeffs = nullptr; bool lanczos_coeffs::sufficiently_accurate(int digits) { if (digits<=20) { @@ -1704,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(realpart(x))); + if (!instanceof(imagpart(x), cln::cl_RA_ring)) + prec = cln::float_format(cln::the(imagpart(x))); + return prec; +} /** The Gamma function. * Use the Lanczos approximation. If the coefficients used here are not * sufficiently many or sufficiently accurate, more can be calculated * using the program doc/examples/lanczos.cpp. In that case, be sure to * read the comments in that file. */ -const numeric lgamma(const numeric &x) +const cln::cl_N lgamma(const cln::cl_N &x) { + cln::float_format_t prec = guess_precision(x); lanczos_coeffs lc; - if (lc.sufficiently_accurate(Digits)) { - cln::cl_N pi_val = cln::pi(cln::default_float_format); - if (x.real() < 0.5) - return log(pi_val) - log(sin(pi_val*x.to_cl_N())) - - lgamma(numeric(1).sub(x)); - cln::cl_N A = lc.calc_lanczos_A(x.to_cl_N()); - cln::cl_N temp = x.to_cl_N() + lc.get_order() - cln::cl_N(1)/2; - cln::cl_N result = log(cln::cl_I(2)*pi_val)/2 - + (x.to_cl_N()-cln::cl_N(1)/2)*log(temp) - - temp - + log(A); - return result; + if (lc.sufficiently_accurate(prec)) { + cln::cl_N pi_val = cln::pi(prec); + if (realpart(x) < 0.5) + return cln::log(pi_val) - cln::log(sin(pi_val*x)) + - lgamma(1 - x); + cln::cl_N A = lc.calc_lanczos_A(x); + cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2; + cln::cl_N result = log(cln::cl_I(2)*pi_val)/2 + + (x-cln::cl_N(1)/2)*log(temp) + - temp + + log(A); + return result; } else throw dunno(); } -const numeric tgamma(const numeric &x) +const numeric lgamma(const numeric &x) { + const cln::cl_N x_ = x.to_cl_N(); + const cln::cl_N result = lgamma(x_); + return numeric(result); +} + +const cln::cl_N tgamma(const cln::cl_N &x) +{ + cln::float_format_t prec = guess_precision(x); lanczos_coeffs lc; - if (lc.sufficiently_accurate(Digits)) { - cln::cl_N pi_val = cln::pi(cln::default_float_format); - if (x.real() < 0.5) - return pi_val/(sin(pi_val*x))/(tgamma(numeric(1).sub(x)).to_cl_N()); - cln::cl_N A = lc.calc_lanczos_A(x.to_cl_N()); - cln::cl_N temp = x.to_cl_N() + lc.get_order() - cln::cl_N(1)/2; - cln::cl_N result - = sqrt(cln::cl_I(2)*pi_val) * expt(temp, x.to_cl_N()-cln::cl_N(1)/2) - * exp(-temp) * A; - return result; + if (lc.sufficiently_accurate(prec)) { + cln::cl_N pi_val = cln::pi(prec); + if (realpart(x) < 0.5) + return pi_val/(cln::sin(pi_val*x))/tgamma(1 - x); + cln::cl_N A = lc.calc_lanczos_A(x); + cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2; + cln::cl_N result = sqrt(cln::cl_I(2)*pi_val) + * expt(temp, x - cln::cl_N(1)/2) + * exp(-temp) * A; + return result; } else throw dunno(); } +const numeric tgamma(const numeric &x) +{ + const cln::cl_N x_ = x.to_cl_N(); + const cln::cl_N result = tgamma(x_); + return numeric(result); +} /** The psi function (aka polygamma function). * This is only a stub! */ @@ -2128,7 +2213,7 @@ const numeric bernoulli(const numeric &nn) next_r = 4; } if (n(a.to_cl_N()), - cln::the(b.to_cl_N())); + return numeric(cln::mod(cln::the(a.to_cl_N()), + cln::the(b.to_cl_N()))); else return *_num0_p; } /** 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(b.to_cl_N()) >> 1) - 1; - return 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; } @@ -2267,8 +2358,8 @@ const numeric irem(const numeric &a, const numeric &b) if (b.is_zero()) throw std::overflow_error("numeric::irem(): division by zero"); if (a.is_integer() && b.is_integer()) - return cln::rem(cln::the(a.to_cl_N()), - cln::the(b.to_cl_N())); + return numeric(cln::rem(cln::the(a.to_cl_N()), + cln::the(b.to_cl_N()))); else return *_num0_p; } @@ -2289,8 +2380,8 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q) if (a.is_integer() && b.is_integer()) { const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); - q = rem_quo.quotient; - return rem_quo.remainder; + q = numeric(rem_quo.quotient); + return numeric(rem_quo.remainder); } else { q = *_num0_p; return *_num0_p; @@ -2308,8 +2399,8 @@ const numeric iquo(const numeric &a, const numeric &b) if (b.is_zero()) throw std::overflow_error("numeric::iquo(): division by zero"); if (a.is_integer() && b.is_integer()) - return cln::truncate1(cln::the(a.to_cl_N()), - cln::the(b.to_cl_N())); + return numeric(cln::truncate1(cln::the(a.to_cl_N()), + cln::the(b.to_cl_N()))); else return *_num0_p; } @@ -2329,8 +2420,8 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r) if (a.is_integer() && b.is_integer()) { const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); - r = rem_quo.remainder; - return rem_quo.quotient; + r = numeric(rem_quo.remainder); + return numeric(rem_quo.quotient); } else { r = *_num0_p; return *_num0_p; @@ -2345,8 +2436,8 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r) const numeric gcd(const numeric &a, const numeric &b) { if (a.is_integer() && b.is_integer()) - return cln::gcd(cln::the(a.to_cl_N()), - cln::the(b.to_cl_N())); + return numeric(cln::gcd(cln::the(a.to_cl_N()), + cln::the(b.to_cl_N()))); else return *_num1_p; } @@ -2359,8 +2450,8 @@ const numeric gcd(const numeric &a, const numeric &b) const numeric lcm(const numeric &a, const numeric &b) { if (a.is_integer() && b.is_integer()) - return cln::lcm(cln::the(a.to_cl_N()), - cln::the(b.to_cl_N())); + return numeric(cln::lcm(cln::the(a.to_cl_N()), + cln::the(b.to_cl_N()))); else return a.mul(b); } @@ -2376,7 +2467,7 @@ const numeric lcm(const numeric &a, const numeric &b) * where imag(x)>0. */ const numeric sqrt(const numeric &x) { - return cln::sqrt(x.to_cl_N()); + return numeric(cln::sqrt(x.to_cl_N())); } @@ -2386,7 +2477,7 @@ const numeric isqrt(const numeric &x) if (x.is_integer()) { cln::cl_I root; cln::isqrt(cln::the(x.to_cl_N()), &root); - return root; + return numeric(root); } else return *_num0_p; } @@ -2438,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;