X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=278741f1e43aec63aa875dd56dd8092687ff58b6;hp=b6c95e01e0c9854ee3e2ab7f97156b6966afc573;hb=80ceed0c855ea06ff5acd9375ad2e0c5ca386373;hpb=8b9b3aa132366c534e5cc8c7a6116afabd885f7d diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index b6c95e01..278741f1 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-2004 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2006 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" @@ -73,7 +73,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(numeric, basic, ////////// /** default ctor. Numerically it initializes to an integer zero. */ -numeric::numeric() : basic(TINFO_numeric) +numeric::numeric() : basic(&numeric::tinfo_static) { value = cln::cl_I(0); setflag(status_flags::evaluated | status_flags::expanded); @@ -85,7 +85,7 @@ numeric::numeric() : basic(TINFO_numeric) // public -numeric::numeric(int i) : basic(TINFO_numeric) +numeric::numeric(int i) : basic(&numeric::tinfo_static) { // 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 @@ -100,14 +100,14 @@ numeric::numeric(int i) : basic(TINFO_numeric) } -numeric::numeric(unsigned int i) : basic(TINFO_numeric) +numeric::numeric(unsigned int i) : basic(&numeric::tinfo_static) { // 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 // 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)); @@ -115,14 +115,14 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric) } -numeric::numeric(long i) : basic(TINFO_numeric) +numeric::numeric(long i) : basic(&numeric::tinfo_static) { value = cln::cl_I(i); setflag(status_flags::evaluated | status_flags::expanded); } -numeric::numeric(unsigned long i) : basic(TINFO_numeric) +numeric::numeric(unsigned long i) : basic(&numeric::tinfo_static) { value = cln::cl_I(i); setflag(status_flags::evaluated | status_flags::expanded); @@ -132,7 +132,7 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric) /** Constructor for rational numerics a/b. * * @exception overflow_error (division by zero) */ -numeric::numeric(long numer, long denom) : basic(TINFO_numeric) +numeric::numeric(long numer, long denom) : basic(&numeric::tinfo_static) { if (!denom) throw std::overflow_error("division by zero"); @@ -141,7 +141,7 @@ numeric::numeric(long numer, long denom) : basic(TINFO_numeric) } -numeric::numeric(double d) : basic(TINFO_numeric) +numeric::numeric(double d) : basic(&numeric::tinfo_static) { // 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 @@ -153,7 +153,7 @@ numeric::numeric(double d) : basic(TINFO_numeric) /** ctor from C-style string. It also accepts complex numbers in GiNaC * notation like "2+5*I". */ -numeric::numeric(const char *s) : basic(TINFO_numeric) +numeric::numeric(const char *s) : basic(&numeric::tinfo_static) { cln::cl_N ctorval = 0; // parse complex numbers (functional but not completely safe, unfortunately @@ -232,12 +232,13 @@ numeric::numeric(const char *s) : basic(TINFO_numeric) /** Ctor from CLN types. This is for the initiated user or internal use * only. */ -numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric) +numeric::numeric(const cln::cl_N &z) : basic(&numeric::tinfo_static) { value = z; setflag(status_flags::evaluated | status_flags::expanded); } + ////////// // archiving ////////// @@ -602,6 +603,11 @@ bool numeric::info(unsigned inf) const return false; } +bool numeric::is_polynomial(const ex & var) const +{ + return true; +} + int numeric::degree(const ex & s) const { return 0; @@ -623,22 +629,29 @@ ex numeric::coeff(const ex & s, int n) const * results: (2+I).has(-2) -> true. But this is consistent, since we also * would like to have (-2+I).has(2) -> true and we want to think about the * sign as a multiplicative factor. */ -bool numeric::has(const ex &other) const +bool numeric::has(const ex &other, unsigned options) const { if (!is_exactly_a(other)) return false; 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; } @@ -760,7 +773,7 @@ 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(other.value,_num1.value)) + if (&other==_num1_p || cln::equal(other.value,_num1_p->value)) return *this; if (cln::zerop(value)) { @@ -771,7 +784,7 @@ const numeric numeric::power(const numeric &other) const 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(value, other.value)); } @@ -857,7 +870,7 @@ 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(other.value, _num1.value)) + if (&other==_num1_p || cln::equal(other.value, _num1_p->value)) return *this; if (cln::zerop(value)) { @@ -868,7 +881,7 @@ const numeric &numeric::power_dyn(const numeric &other) const 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(value, other.value)))-> setflag(status_flags::dynallocated)); @@ -919,6 +932,18 @@ const numeric numeric::inverse() const return numeric(cln::recip(value)); } +/** Return the step function of a numeric. The imaginary part of it is + * ignored because the step function is generally considered real but + * a numeric may develop a small imaginary part due to rounding errors. + */ +numeric numeric::step() const +{ cln::cl_R r = cln::realpart(value); + if(cln::zerop(r)) + return numeric(1,2); + if(cln::plusp(r)) + return 1; + return 0; +} /** Return the complex half-plane (left or right) in which the number lies. * csgn(x)==0 for x==0, csgn(x)==1 for Re(x)>0 or Re(x)=0 and Im(x)>0, @@ -1237,7 +1262,7 @@ 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))); @@ -1246,7 +1271,7 @@ const numeric numeric::denom() const 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)) @@ -1255,7 +1280,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; } @@ -1359,7 +1384,7 @@ 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()); } @@ -1510,7 +1535,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 @@ -1601,8 +1626,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"); @@ -1619,17 +1644,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(k+1)/gamma(n-k+1) or a suitable limit - throw std::range_error("numeric::binomial(): don´t know how to evaluate that."); + throw std::range_error("numeric::binomial(): don't know how to evaluate that."); } @@ -1681,9 +1706,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; @@ -1700,20 +1725,20 @@ 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; - const unsigned p3 = p+3; - const unsigned pm = p-2; - unsigned i, k, p_2; - // test if intermediate unsigned int can be represented by immediate - // objects by CLN (i.e. < 2^29 for 32 Bit machines, see ) + cln::cl_RA b = cln::cl_RA(p-1)/-2; + // The CLN manual says: "The conversion from `unsigned int' works only + // if the argument is < 2^29" (This is for 32 Bit machines. More + // generally, cl_value_len is the limiting exponent of 2. We must make + // sure that no intermediates are created which exceed this value. The + // largest intermediate is (p+3-2*k)*(p/2-k+1) <= (p^2+p)/2. if (p < (1UL<(a.to_cl_N()), cln::the(b.to_cl_N())); else - return _num0; + return *_num0_p; } @@ -1818,7 +1843,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; } @@ -1837,7 +1862,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; } @@ -1859,8 +1884,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; } } @@ -1878,7 +1903,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; } @@ -1899,8 +1924,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; } } @@ -1915,7 +1940,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; } @@ -1955,7 +1980,7 @@ const numeric isqrt(const numeric &x) cln::isqrt(cln::the(x.to_cl_N()), &root); return root; } else - return _num0; + return *_num0_p; } @@ -1991,14 +2016,25 @@ _numeric_digits::_numeric_digits() throw(std::runtime_error("I told you not to do instantiate me!")); too_late = true; cln::default_float_format = cln::float_format(17); + + // add callbacks for built-in functions + // like ... add_callback(Li_lookuptable); } /** Assign a native long to global Digits object. */ _numeric_digits& _numeric_digits::operator=(long prec) { + long digitsdiff = prec - digits; digits = prec; - cln::default_float_format = cln::float_format(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); + } + return *this; } @@ -2018,6 +2054,13 @@ void _numeric_digits::print(std::ostream &os) const } +/** Add a new callback function. */ +void _numeric_digits::add_callback(digits_changed_callback callback) +{ + callbacklist.push_back(callback); +} + + std::ostream& operator<<(std::ostream &os, const _numeric_digits &e) { e.print(os);