X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=5649353f308f68a0c26e715034b84977c5702979;hp=cec39c8c974da9b235ac1a4a89a6f4d950f06fbd;hb=94e5076a9ad0d142ee4604e9b078a10083efd24c;hpb=83a7ee99a947cbbf331018b803ad6be43a9ccd45 diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index cec39c8c..5649353f 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,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(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()); } @@ -1379,7 +1439,7 @@ const numeric I = numeric(cln::complex(cln::cl_I(0),cln::cl_I(1))); * @return arbitrary precision numerical exp(x). */ const numeric exp(const numeric &x) { - return cln::exp(x.to_cl_N()); + return numeric(cln::exp(x.to_cl_N())); } @@ -1392,7 +1452,7 @@ const numeric log(const numeric &x) { if (x.is_zero()) throw pole_error("log(): logarithmic pole",0); - return cln::log(x.to_cl_N()); + return numeric(cln::log(x.to_cl_N())); } @@ -1401,7 +1461,7 @@ const numeric log(const numeric &x) * @return arbitrary precision numerical sin(x). */ const numeric sin(const numeric &x) { - return cln::sin(x.to_cl_N()); + return numeric(cln::sin(x.to_cl_N())); } @@ -1410,7 +1470,7 @@ const numeric sin(const numeric &x) * @return arbitrary precision numerical cos(x). */ const numeric cos(const numeric &x) { - return cln::cos(x.to_cl_N()); + return numeric(cln::cos(x.to_cl_N())); } @@ -1419,7 +1479,7 @@ const numeric cos(const numeric &x) * @return arbitrary precision numerical tan(x). */ const numeric tan(const numeric &x) { - return cln::tan(x.to_cl_N()); + return numeric(cln::tan(x.to_cl_N())); } @@ -1428,7 +1488,7 @@ const numeric tan(const numeric &x) * @return arbitrary precision numerical asin(x). */ const numeric asin(const numeric &x) { - return cln::asin(x.to_cl_N()); + return numeric(cln::asin(x.to_cl_N())); } @@ -1437,7 +1497,7 @@ const numeric asin(const numeric &x) * @return arbitrary precision numerical acos(x). */ const numeric acos(const numeric &x) { - return cln::acos(x.to_cl_N()); + return numeric(cln::acos(x.to_cl_N())); } @@ -1452,7 +1512,7 @@ const numeric atan(const numeric &x) x.real().is_zero() && abs(x.imag()).is_equal(*_num1_p)) throw pole_error("atan(): logarithmic pole",0); - return cln::atan(x.to_cl_N()); + return numeric(cln::atan(x.to_cl_N())); } @@ -1468,8 +1528,8 @@ const numeric atan(const numeric &y, const numeric &x) if (x.is_zero() && y.is_zero()) return *_num0_p; if (x.is_real() && y.is_real()) - return cln::atan(cln::the(x.to_cl_N()), - cln::the(y.to_cl_N())); + return numeric(cln::atan(cln::the(x.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))) @@ -1485,7 +1545,7 @@ const numeric atan(const numeric &y, const numeric &x) // x-I*y==0 => y/x==-I, so this is a pole (we have x!=0). throw pole_error("atan(): logarithmic pole",0); } - return cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m)); + return numeric(cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m))); } @@ -1494,7 +1554,7 @@ const numeric atan(const numeric &y, const numeric &x) * @return arbitrary precision numerical sinh(x). */ const numeric sinh(const numeric &x) { - return cln::sinh(x.to_cl_N()); + return numeric(cln::sinh(x.to_cl_N())); } @@ -1503,7 +1563,7 @@ const numeric sinh(const numeric &x) * @return arbitrary precision numerical cosh(x). */ const numeric cosh(const numeric &x) { - return cln::cosh(x.to_cl_N()); + return numeric(cln::cosh(x.to_cl_N())); } @@ -1512,7 +1572,7 @@ const numeric cosh(const numeric &x) * @return arbitrary precision numerical tanh(x). */ const numeric tanh(const numeric &x) { - return cln::tanh(x.to_cl_N()); + return numeric(cln::tanh(x.to_cl_N())); } @@ -1521,7 +1581,7 @@ const numeric tanh(const numeric &x) * @return arbitrary precision numerical asinh(x). */ const numeric asinh(const numeric &x) { - return cln::asinh(x.to_cl_N()); + return numeric(cln::asinh(x.to_cl_N())); } @@ -1530,7 +1590,7 @@ const numeric asinh(const numeric &x) * @return arbitrary precision numerical acosh(x). */ const numeric acosh(const numeric &x) { - return cln::acosh(x.to_cl_N()); + return numeric(cln::acosh(x.to_cl_N())); } @@ -1539,7 +1599,7 @@ const numeric acosh(const numeric &x) * @return arbitrary precision numerical atanh(x). */ const numeric atanh(const numeric &x) { - return cln::atanh(x.to_cl_N()); + return numeric(cln::atanh(x.to_cl_N())); } @@ -1611,24 +1671,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 +1700,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 +1725,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(); } @@ -1953,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; @@ -1979,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; } @@ -1997,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! */ @@ -2128,7 +2221,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; } @@ -2248,8 +2341,8 @@ 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 numeric(cln::mod(cln::the(a.to_cl_N()) + b2, + cln::the(b.to_cl_N())) - b2); } else return *_num0_p; } @@ -2267,8 +2360,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 +2382,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 +2401,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 +2422,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 +2438,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 +2452,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 +2469,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 +2479,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; }