X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=26e012b9f75f60fc104897c7e8648500b69fa251;hb=eba3e6126d33fff5f1231906ffd8d8d65f14c953;hp=f9c4a2781be7cdaecb22cbb7f3b2aee20d41a998;hpb=13896261ee985f23e5b5648532e70f0cce704ede;p=ginac.git diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index f9c4a278..26e012b9 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-2003 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2005 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 @@ -21,7 +21,7 @@ * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ #include "config.h" @@ -107,7 +107,7 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric) // emphasizes efficiency. However, if the integer is small enough // we save space and dereferences by using an immediate type. // (C.f. ) - if (i < (1U << (cl_value_len-1))) + if (i < (1UL << (cl_value_len-1))) value = cln::cl_I(i); else value = cln::cl_I(static_cast(i)); @@ -238,6 +238,7 @@ numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric) setflag(status_flags::evaluated | status_flags::expanded); } + ////////// // archiving ////////// @@ -281,7 +282,7 @@ void numeric::archive(archive_node &n) const // Write number as string std::ostringstream s; if (this->is_crational()) - s << cln::the(value); + s << value; else { // Non-rational numbers are written in an integer-decoded format // to preserve the precision @@ -417,8 +418,8 @@ static void print_real_cl_N(const print_context & c, const cln::cl_R & x) void numeric::print_numeric(const print_context & c, const char *par_open, const char *par_close, const char *imag_sym, const char *mul_sym, unsigned level) const { - const cln::cl_R r = cln::realpart(cln::the(value)); - const cln::cl_R i = cln::imagpart(cln::the(value)); + const cln::cl_R r = cln::realpart(value); + const cln::cl_R i = cln::imagpart(value); if (cln::zerop(i)) { @@ -514,9 +515,9 @@ void numeric::do_print_csrc(const print_csrc & c, unsigned level) const else c.s << "float>("; - print_real_csrc(c, cln::realpart(cln::the(value))); + print_real_csrc(c, cln::realpart(value)); c.s << ","; - print_real_csrc(c, cln::imagpart(cln::the(value))); + print_real_csrc(c, cln::imagpart(value)); c.s << ")"; } @@ -535,16 +536,16 @@ void numeric::do_print_csrc_cl_N(const print_csrc_cl_N & c, unsigned level) cons // Complex number c.s << "cln::complex("; - print_real_cl_N(c, cln::realpart(cln::the(value))); + print_real_cl_N(c, cln::realpart(value)); c.s << ","; - print_real_cl_N(c, cln::imagpart(cln::the(value))); + print_real_cl_N(c, cln::imagpart(value)); c.s << ")"; } } void numeric::do_print_tree(const print_tree & c, unsigned level) const { - c.s << std::string(level, ' ') << cln::the(value) + c.s << std::string(level, ' ') << value << " (" << class_name() << ")" << " @" << this << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec << std::endl; @@ -630,15 +631,22 @@ bool numeric::has(const ex &other) const const numeric &o = ex_to(other); if (this->is_equal(o) || this->is_equal(-o)) return true; - if (o.imag().is_zero()) // e.g. scan for 3 in -3*I - return (this->real().is_equal(o) || this->imag().is_equal(o) || - this->real().is_equal(-o) || this->imag().is_equal(-o)); + if (o.imag().is_zero()) { // e.g. scan for 3 in -3*I + if (!this->real().is_equal(*_num0_p)) + if (this->real().is_equal(o) || this->real().is_equal(-o)) + return true; + if (!this->imag().is_equal(*_num0_p)) + if (this->imag().is_equal(o) || this->imag().is_equal(-o)) + return true; + return false; + } else { if (o.is_equal(I)) // e.g scan for I in 42*I return !this->is_real(); if (o.real().is_zero()) // e.g. scan for 2*I in 2*I+1 - return (this->real().has(o*I) || this->imag().has(o*I) || - this->real().has(-o*I) || this->imag().has(-o*I)); + if (!this->imag().is_equal(*_num0_p)) + if (this->imag().is_equal(o*I) || this->imag().is_equal(-o*I)) + return true; } return false; } @@ -663,8 +671,15 @@ ex numeric::eval(int level) const ex numeric::evalf(int level) const { // level can safely be discarded for numeric objects. - return numeric(cln::cl_float(1.0, cln::default_float_format) * - (cln::the(value))); + return numeric(cln::cl_float(1.0, cln::default_float_format) * value); +} + +ex numeric::conjugate() const +{ + if (is_real()) { + return *this; + } + return numeric(cln::conjugate(this->value)); } // protected @@ -694,7 +709,7 @@ unsigned numeric::calchash() const // equivalence relation on numbers). As a consequence, 3 and 3.0 share // the same hashvalue. That shouldn't really matter, though. setflag(status_flags::hash_calculated); - hashvalue = golden_ratio_hash(cln::equal_hashcode(cln::the(value))); + hashvalue = golden_ratio_hash(cln::equal_hashcode(value)); return hashvalue; } @@ -715,7 +730,7 @@ unsigned numeric::calchash() const * a numeric object. */ const numeric numeric::add(const numeric &other) const { - return numeric(cln::the(value)+cln::the(other.value)); + return numeric(value + other.value); } @@ -723,7 +738,7 @@ const numeric numeric::add(const numeric &other) const * result as a numeric object. */ const numeric numeric::sub(const numeric &other) const { - return numeric(cln::the(value)-cln::the(other.value)); + return numeric(value - other.value); } @@ -731,7 +746,7 @@ const numeric numeric::sub(const numeric &other) const * result as a numeric object. */ const numeric numeric::mul(const numeric &other) const { - return numeric(cln::the(value)*cln::the(other.value)); + return numeric(value * other.value); } @@ -741,9 +756,9 @@ const numeric numeric::mul(const numeric &other) const * @exception overflow_error (division by zero) */ const numeric numeric::div(const numeric &other) const { - if (cln::zerop(cln::the(other.value))) + if (cln::zerop(other.value)) throw std::overflow_error("numeric::div(): division by zero"); - return numeric(cln::the(value)/cln::the(other.value)); + return numeric(value / other.value); } @@ -753,20 +768,20 @@ const numeric numeric::power(const numeric &other) const { // Shortcut for efficiency and numeric stability (as in 1.0 exponent): // trap the neutral exponent. - if (&other==_num1_p || cln::equal(cln::the(other.value),cln::the(_num1.value))) + if (&other==_num1_p || cln::equal(other.value,_num1_p->value)) return *this; - if (cln::zerop(cln::the(value))) { - if (cln::zerop(cln::the(other.value))) + if (cln::zerop(value)) { + if (cln::zerop(other.value)) throw std::domain_error("numeric::eval(): pow(0,0) is undefined"); - else if (cln::zerop(cln::realpart(cln::the(other.value)))) + else if (cln::zerop(cln::realpart(other.value))) throw std::domain_error("numeric::eval(): pow(0,I) is undefined"); - else if (cln::minusp(cln::realpart(cln::the(other.value)))) + else if (cln::minusp(cln::realpart(other.value))) throw std::overflow_error("numeric::eval(): division by zero"); else - return _num0; + return *_num0_p; } - return numeric(cln::expt(cln::the(value),cln::the(other.value))); + return numeric(cln::expt(value, other.value)); } @@ -783,7 +798,7 @@ const numeric &numeric::add_dyn(const numeric &other) const else if (&other==_num0_p) return *this; - return static_cast((new numeric(cln::the(value)+cln::the(other.value)))-> + return static_cast((new numeric(value + other.value))-> setflag(status_flags::dynallocated)); } @@ -796,10 +811,10 @@ const numeric &numeric::sub_dyn(const numeric &other) const { // Efficiency shortcut: trap the neutral exponent (first by pointer). This // hack is supposed to keep the number of distinct numeric objects low. - if (&other==_num0_p || cln::zerop(cln::the(other.value))) + if (&other==_num0_p || cln::zerop(other.value)) return *this; - return static_cast((new numeric(cln::the(value)-cln::the(other.value)))-> + return static_cast((new numeric(value - other.value))-> setflag(status_flags::dynallocated)); } @@ -817,7 +832,7 @@ const numeric &numeric::mul_dyn(const numeric &other) const else if (&other==_num1_p) return *this; - return static_cast((new numeric(cln::the(value)*cln::the(other.value)))-> + return static_cast((new numeric(value * other.value))-> setflag(status_flags::dynallocated)); } @@ -836,7 +851,7 @@ 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(cln::the(value)/cln::the(other.value)))-> + return static_cast((new numeric(value / other.value))-> setflag(status_flags::dynallocated)); } @@ -850,20 +865,20 @@ const numeric &numeric::power_dyn(const numeric &other) const // Efficiency shortcut: trap the neutral exponent (first try by pointer, then // try harder, since calls to cln::expt() below may return amazing results for // floating point exponent 1.0). - if (&other==_num1_p || cln::equal(cln::the(other.value),cln::the(_num1.value))) + if (&other==_num1_p || cln::equal(other.value, _num1_p->value)) return *this; - if (cln::zerop(cln::the(value))) { - if (cln::zerop(cln::the(other.value))) + if (cln::zerop(value)) { + if (cln::zerop(other.value)) throw std::domain_error("numeric::eval(): pow(0,0) is undefined"); - else if (cln::zerop(cln::realpart(cln::the(other.value)))) + else if (cln::zerop(cln::realpart(other.value))) throw std::domain_error("numeric::eval(): pow(0,I) is undefined"); - else if (cln::minusp(cln::realpart(cln::the(other.value)))) + else if (cln::minusp(cln::realpart(other.value))) throw std::overflow_error("numeric::eval(): division by zero"); else - return _num0; + return *_num0_p; } - return static_cast((new numeric(cln::expt(cln::the(value),cln::the(other.value))))-> + return static_cast((new numeric(cln::expt(value, other.value)))-> setflag(status_flags::dynallocated)); } @@ -907,9 +922,9 @@ const numeric &numeric::operator=(const char * s) /** Inverse of a number. */ const numeric numeric::inverse() const { - if (cln::zerop(cln::the(value))) + if (cln::zerop(value)) throw std::overflow_error("numeric::inverse(): division by zero"); - return numeric(cln::recip(cln::the(value))); + return numeric(cln::recip(value)); } @@ -920,16 +935,16 @@ const numeric numeric::inverse() const * @see numeric::compare(const numeric &other) */ int numeric::csgn() const { - if (cln::zerop(cln::the(value))) + if (cln::zerop(value)) return 0; - cln::cl_R r = cln::realpart(cln::the(value)); + cln::cl_R r = cln::realpart(value); if (!cln::zerop(r)) { if (cln::plusp(r)) return 1; else return -1; } else { - if (cln::plusp(cln::imagpart(cln::the(value)))) + if (cln::plusp(cln::imagpart(value))) return 1; else return -1; @@ -953,25 +968,25 @@ int numeric::compare(const numeric &other) const return cln::compare(cln::the(value), cln::the(other.value)); else { // No, first cln::compare real parts... - cl_signean real_cmp = cln::compare(cln::realpart(cln::the(value)), cln::realpart(cln::the(other.value))); + cl_signean real_cmp = cln::compare(cln::realpart(value), cln::realpart(other.value)); if (real_cmp) return real_cmp; // ...and then the imaginary parts. - return cln::compare(cln::imagpart(cln::the(value)), cln::imagpart(cln::the(other.value))); + return cln::compare(cln::imagpart(value), cln::imagpart(other.value)); } } bool numeric::is_equal(const numeric &other) const { - return cln::equal(cln::the(value),cln::the(other.value)); + return cln::equal(value, other.value); } /** True if object is zero. */ bool numeric::is_zero() const { - return cln::zerop(cln::the(value)); + return cln::zerop(value); } @@ -1056,13 +1071,13 @@ bool numeric::is_real() const bool numeric::operator==(const numeric &other) const { - return cln::equal(cln::the(value), cln::the(other.value)); + return cln::equal(value, other.value); } bool numeric::operator!=(const numeric &other) const { - return !cln::equal(cln::the(value), cln::the(other.value)); + return !cln::equal(value, other.value); } @@ -1073,8 +1088,8 @@ bool numeric::is_cinteger() const if (cln::instanceof(value, cln::cl_I_ring)) return true; else if (!this->is_real()) { // complex case, handle n+m*I - if (cln::instanceof(cln::realpart(cln::the(value)), cln::cl_I_ring) && - cln::instanceof(cln::imagpart(cln::the(value)), cln::cl_I_ring)) + if (cln::instanceof(cln::realpart(value), cln::cl_I_ring) && + cln::instanceof(cln::imagpart(value), cln::cl_I_ring)) return true; } return false; @@ -1088,8 +1103,8 @@ bool numeric::is_crational() const if (cln::instanceof(value, cln::cl_RA_ring)) return true; else if (!this->is_real()) { // complex case, handle Q(i): - if (cln::instanceof(cln::realpart(cln::the(value)), cln::cl_RA_ring) && - cln::instanceof(cln::imagpart(cln::the(value)), cln::cl_RA_ring)) + if (cln::instanceof(cln::realpart(value), cln::cl_RA_ring) && + cln::instanceof(cln::imagpart(value), cln::cl_RA_ring)) return true; } return false; @@ -1165,7 +1180,7 @@ long numeric::to_long() const double numeric::to_double() const { GINAC_ASSERT(this->is_real()); - return cln::double_approx(cln::realpart(cln::the(value))); + return cln::double_approx(cln::realpart(value)); } @@ -1174,21 +1189,21 @@ double numeric::to_double() const */ cln::cl_N numeric::to_cl_N() const { - return cln::cl_N(cln::the(value)); + return value; } /** Real part of a number. */ const numeric numeric::real() const { - return numeric(cln::realpart(cln::the(value))); + return numeric(cln::realpart(value)); } /** Imaginary part of a number. */ const numeric numeric::imag() const { - return numeric(cln::imagpart(cln::the(value))); + return numeric(cln::imagpart(value)); } @@ -1205,8 +1220,8 @@ const numeric numeric::numer() const return numeric(cln::numerator(cln::the(value))); else if (!this->is_real()) { // complex case, handle Q(i): - const cln::cl_RA r = cln::the(cln::realpart(cln::the(value))); - const cln::cl_RA i = cln::the(cln::imagpart(cln::the(value))); + const cln::cl_RA r = cln::the(cln::realpart(value)); + const cln::cl_RA i = cln::the(cln::imagpart(value)); if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_I_ring)) return numeric(*this); if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_RA_ring)) @@ -1230,16 +1245,16 @@ const numeric numeric::numer() const const numeric numeric::denom() const { if (cln::instanceof(value, cln::cl_I_ring)) - return _num1; // integer case + return *_num1_p; // integer case if (cln::instanceof(value, cln::cl_RA_ring)) return numeric(cln::denominator(cln::the(value))); if (!this->is_real()) { // complex case, handle Q(i): - const cln::cl_RA r = cln::the(cln::realpart(cln::the(value))); - const cln::cl_RA i = cln::the(cln::imagpart(cln::the(value))); + const cln::cl_RA r = cln::the(cln::realpart(value)); + const cln::cl_RA i = cln::the(cln::imagpart(value)); if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_I_ring)) - return _num1; + return *_num1_p; if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_RA_ring)) return numeric(cln::denominator(i)); if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_I_ring)) @@ -1248,7 +1263,7 @@ const numeric numeric::denom() const return numeric(cln::lcm(cln::denominator(r), cln::denominator(i))); } // at least one float encountered - return _num1; + return *_num1_p; } @@ -1287,14 +1302,14 @@ const numeric exp(const numeric &x) /** Natural logarithm. * - * @param z complex number + * @param x complex number * @return arbitrary precision numerical log(x). * @exception pole_error("log(): logarithmic pole",0) */ -const numeric log(const numeric &z) +const numeric log(const numeric &x) { - if (z.is_zero()) + if (x.is_zero()) throw pole_error("log(): logarithmic pole",0); - return cln::log(z.to_cl_N()); + return cln::log(x.to_cl_N()); } @@ -1345,14 +1360,14 @@ const numeric acos(const numeric &x) /** Arcustangent. * - * @param z complex number - * @return atan(z) + * @param x complex number + * @return atan(x) * @exception pole_error("atan(): logarithmic pole",0) */ const numeric atan(const numeric &x) { if (!x.is_real() && x.real().is_zero() && - abs(x.imag()).is_equal(_num1)) + abs(x.imag()).is_equal(*_num1_p)) throw pole_error("atan(): logarithmic pole",0); return cln::atan(x.to_cl_N()); } @@ -1503,7 +1518,7 @@ static cln::cl_N Li2_projection(const cln::cl_N &x, const numeric Li2(const numeric &x) { if (x.is_zero()) - return _num0; + return *_num0_p; // what is the desired float format? // first guess: default format @@ -1515,7 +1530,7 @@ const numeric Li2(const numeric &x) else if (!x.imag().is_rational()) prec = cln::float_format(cln::the(cln::imagpart(value))); - if (cln::the(value)==1) // may cause trouble with log(1-x) + if (value==1) // may cause trouble with log(1-x) return cln::zeta(2, prec); if (cln::abs(value) > 1) @@ -1594,8 +1609,8 @@ const numeric factorial(const numeric &n) * @exception range_error (argument must be integer >= -1) */ const numeric doublefactorial(const numeric &n) { - if (n.is_equal(_num_1)) - return _num1; + if (n.is_equal(*_num_1_p)) + return *_num1_p; if (!n.is_nonneg_integer()) throw std::range_error("numeric::doublefactorial(): argument must be integer >= -1"); @@ -1612,17 +1627,17 @@ const numeric binomial(const numeric &n, const numeric &k) { if (n.is_integer() && k.is_integer()) { if (n.is_nonneg_integer()) { - if (k.compare(n)!=1 && k.compare(_num0)!=-1) + if (k.compare(n)!=1 && k.compare(*_num0_p)!=-1) return numeric(cln::binomial(n.to_int(),k.to_int())); else - return _num0; + return *_num0_p; } else { - return _num_1.power(k)*binomial(k-n-_num1,k); + return _num_1_p->power(k)*binomial(k-n-(*_num1_p),k); } } - // should really be gamma(n+1)/gamma(r+1)/gamma(n-r+1) or a suitable limit - throw std::range_error("numeric::binomial(): don´t know how to evaluate that."); + // should really be gamma(n+1)/gamma(k+1)/gamma(n-k+1) or a suitable limit + throw std::range_error("numeric::binomial(): don't know how to evaluate that."); } @@ -1674,9 +1689,9 @@ const numeric bernoulli(const numeric &nn) // the special cases not covered by the algorithm below if (n & 1) - return (n==1) ? _num_1_2 : _num0; + return (n==1) ? (*_num_1_2_p) : (*_num0_p); if (!n) - return _num1; + return *_num1_p; // store nonvanishing Bernoulli numbers here static std::vector< cln::cl_RA > results; @@ -1693,7 +1708,7 @@ const numeric bernoulli(const numeric &nn) results.reserve(n/2); for (unsigned p=next_r; p<=n; p+=2) { cln::cl_I c = 1; // seed for binonmial coefficients - cln::cl_RA b = cln::cl_RA(1-p)/2; + cln::cl_RA b = cln::cl_RA(p-1)/-2; const unsigned p3 = p+3; const unsigned pm = p-2; unsigned i, k, p_2; @@ -1744,7 +1759,7 @@ const numeric fibonacci(const numeric &n) // hence // F(2n+2) = F(n+1)*(2*F(n) + F(n+1)) if (n.is_zero()) - return _num0; + return *_num0_p; if (n.is_negative()) if (n.is_even()) return -fibonacci(-n); @@ -1796,14 +1811,14 @@ const numeric mod(const numeric &a, const numeric &b) return cln::mod(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); else - return _num0; + return *_num0_p; } /** Modulus (in symmetric representation). * Equivalent to Maple's mods. * - * @return a mod b in the range [-iquo(abs(m)-1,2), iquo(abs(m),2)]. */ + * @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()) { @@ -1811,7 +1826,7 @@ const numeric smod(const numeric &a, const numeric &b) return cln::mod(cln::the(a.to_cl_N()) + b2, cln::the(b.to_cl_N())) - b2; } else - return _num0; + return *_num0_p; } @@ -1830,7 +1845,7 @@ const numeric irem(const numeric &a, const numeric &b) return cln::rem(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); else - return _num0; + return *_num0_p; } @@ -1852,8 +1867,8 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q) q = rem_quo.quotient; return rem_quo.remainder; } else { - q = _num0; - return _num0; + q = *_num0_p; + return *_num0_p; } } @@ -1871,7 +1886,7 @@ const numeric iquo(const numeric &a, const numeric &b) return cln::truncate1(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); else - return _num0; + return *_num0_p; } @@ -1892,8 +1907,8 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r) r = rem_quo.remainder; return rem_quo.quotient; } else { - r = _num0; - return _num0; + r = *_num0_p; + return *_num0_p; } } @@ -1908,7 +1923,7 @@ const numeric gcd(const numeric &a, const numeric &b) return cln::gcd(cln::the(a.to_cl_N()), cln::the(b.to_cl_N())); else - return _num1; + return *_num1_p; } @@ -1927,16 +1942,16 @@ const numeric lcm(const numeric &a, const numeric &b) /** Numeric square root. - * If possible, sqrt(z) should respect squares of exact numbers, i.e. sqrt(4) + * If possible, sqrt(x) should respect squares of exact numbers, i.e. sqrt(4) * should return integer 2. * - * @param z numeric argument - * @return square root of z. Branch cut along negative real axis, the negative - * real axis itself where imag(z)==0 and real(z)<0 belongs to the upper part - * where imag(z)>0. */ -const numeric sqrt(const numeric &z) + * @param x numeric argument + * @return square root of x. Branch cut along negative real axis, the negative + * real axis itself where imag(x)==0 and real(x)<0 belongs to the upper part + * where imag(x)>0. */ +const numeric sqrt(const numeric &x) { - return cln::sqrt(z.to_cl_N()); + return cln::sqrt(x.to_cl_N()); } @@ -1948,7 +1963,7 @@ const numeric isqrt(const numeric &x) cln::isqrt(cln::the(x.to_cl_N()), &root); return root; } else - return _num0; + return *_num0_p; }