X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=50d873a9f6c5ac0da255553db40f64fb14a75057;hp=c80ff3a7757e57685a63e7b583d78f959833871b;hb=3f0b0165865bbb297901e9542fced88a0e32298e;hpb=7bc327c75aaa3118de14dfcee59bcf0fd40e3f4a diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index c80ff3a7..50d873a9 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -73,7 +73,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 +85,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 +106,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 +127,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 +144,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 +153,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 +165,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 @@ -244,7 +244,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 +255,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 '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 '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; - default: // Ordinary number + 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 ////////// @@ -1962,24 +2022,34 @@ lanczos_coeffs::lanczos_coeffs() coeffs[3].swap(coeffs_120); } +static const cln::float_format_t guess_precision(const cln::cl_N& x) +{ + cln::float_format_t prec = cln::default_float_format; + if (!instanceof(realpart(x), cln::cl_RA_ring)) + prec = cln::float_format(cln::the(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; + 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.to_cl_N()-cln::cl_N(1)/2)*log(temp) + + (x-cln::cl_N(1)/2)*log(temp) - temp + log(A); return result; @@ -1988,17 +2058,25 @@ const numeric lgamma(const numeric &x) 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; + 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.to_cl_N()-cln::cl_N(1)/2) + = sqrt(cln::cl_I(2)*pi_val) * expt(temp, x - cln::cl_N(1)/2) * exp(-temp) * A; return result; } @@ -2006,6 +2084,12 @@ const numeric tgamma(const numeric &x) 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! */ @@ -2137,7 +2221,7 @@ const numeric bernoulli(const numeric &nn) next_r = 4; } if (n(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; } @@ -2339,7 +2429,7 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r) const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); r = numeric(rem_quo.remainder); - return rem_quo.quotient; + return numeric(rem_quo.quotient); } else { r = *_num0_p; return *_num0_p;