X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=b1a4033cd73fe5c57f51dd0cb858bc843b8d2fcc;hp=995672e3c3372a845d3a250f1b95321a1a4eb02b;hb=27d6204effdef95a00af461fff98024e290dbaa7;hpb=c20f44225542fd9e4c0cda17660af6a276d3ed9f diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 995672e3..b1a4033c 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -40,6 +40,7 @@ #include "numeric.h" #include "ex.h" +#include "print.h" #include "archive.h" #include "debugmsg.h" #include "utils.h" @@ -72,8 +73,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(numeric, basic) // operator and helpers ////////// -// public - /** default ctor. Numerically it initializes to an integer zero. */ numeric::numeric() : basic(TINFO_numeric) { @@ -82,19 +81,13 @@ numeric::numeric() : basic(TINFO_numeric) setflag(status_flags::evaluated | status_flags::expanded); } -// protected - -/** For use by copy ctor and assignment operator. */ void numeric::copy(const numeric &other) { inherited::copy(other); value = other.value; } -void numeric::destroy(bool call_parent) -{ - if (call_parent) inherited::destroy(call_parent); -} +DEFAULT_DESTROY(numeric) ////////// // other ctors @@ -266,7 +259,6 @@ numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric) // archiving ////////// -/** Construct object from archive_node. */ numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { debugmsg("numeric ctor from archive_node", LOGLEVEL_CONSTRUCT); @@ -304,13 +296,6 @@ numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l setflag(status_flags::evaluated | status_flags::expanded); } -/** Unarchive the object. */ -ex numeric::unarchive(const archive_node &n, const lst &sym_lst) -{ - return (new numeric(n, sym_lst))->setflag(status_flags::dynallocated); -} - -/** Archive the object. */ void numeric::archive(archive_node &n) const { inherited::archive(n); @@ -348,8 +333,10 @@ void numeric::archive(archive_node &n) const #endif } +DEFAULT_UNARCHIVE(numeric) + ////////// -// functions overriding virtual functions from bases classes +// functions overriding virtual functions from base classes ////////// /** Helper function to print a real number in a nicer way than is CLN's @@ -359,146 +346,144 @@ void numeric::archive(archive_node &n) const * want to visibly distinguish from cl_LF. * * @see numeric::print() */ -static void print_real_number(std::ostream &os, const cln::cl_R &num) +static void print_real_number(const print_context & c, const cln::cl_R &x) { cln::cl_print_flags ourflags; - if (cln::instanceof(num, cln::cl_RA_ring)) { - // case 1: integer or rational, nothing special to do: - cln::print_real(os, ourflags, num); + if (cln::instanceof(x, cln::cl_RA_ring)) { + // case 1: integer or rational + if (cln::instanceof(x, cln::cl_I_ring) || + !is_a(c)) { + cln::print_real(c.s, ourflags, x); + } else { // rational output in LaTeX context + c.s << "\\frac{"; + cln::print_real(c.s, ourflags, cln::numerator(cln::the(x))); + c.s << "}{"; + cln::print_real(c.s, ourflags, cln::denominator(cln::the(x))); + c.s << '}'; + } } else { // case 2: float // make CLN believe this number has default_float_format, so it prints // 'E' as exponent marker instead of 'L': - ourflags.default_float_format = cln::float_format(cln::the(num)); - cln::print_real(os, ourflags, num); + ourflags.default_float_format = cln::float_format(cln::the(x)); + cln::print_real(c.s, ourflags, x); } - return; } /** This method adds to the output so it blends more consistently together * with the other routines and produces something compatible to ginsh input. * * @see print_real_number() */ -void numeric::print(std::ostream &os, unsigned upper_precedence) const +void numeric::print(const print_context & c, unsigned level) const { debugmsg("numeric print", LOGLEVEL_PRINT); - cln::cl_R r = cln::realpart(cln::the(value)); - cln::cl_R i = cln::imagpart(cln::the(value)); - if (cln::zerop(i)) { - // case 1, real: x or -x - if ((precedence<=upper_precedence) && (!this->is_nonneg_integer())) { - os << "("; - print_real_number(os, r); - os << ")"; + + if (is_a(c)) { + + c.s << std::string(level, ' ') << cln::the(value) + << " (" << class_name() << ")" + << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec + << std::endl; + + } else if (is_a(c)) { + + std::ios::fmtflags oldflags = c.s.flags(); + c.s.setf(std::ios::scientific); + if (this->is_rational() && !this->is_integer()) { + if (compare(_num0) > 0) { + c.s << "("; + if (is_a(c)) + c.s << "cln::cl_F(\"" << numer().evalf() << "\")"; + else + c.s << numer().to_double(); + } else { + c.s << "-("; + if (is_a(c)) + c.s << "cln::cl_F(\"" << -numer().evalf() << "\")"; + else + c.s << -numer().to_double(); + } + c.s << "/"; + if (is_a(c)) + c.s << "cln::cl_F(\"" << denom().evalf() << "\")"; + else + c.s << denom().to_double(); + c.s << ")"; } else { - print_real_number(os, r); + if (is_a(c)) + c.s << "cln::cl_F(\"" << evalf() << "\")"; + else + c.s << to_double(); } + c.s.flags(oldflags); + } else { - if (cln::zerop(r)) { - // case 2, imaginary: y*I or -y*I - if ((precedence<=upper_precedence) && (i < 0)) { - if (i == -1) { - os << "(-I)"; - } else { - os << "("; - print_real_number(os, i); - os << "*I)"; - } + const std::string par_open = is_a(c) ? "{(" : "("; + const std::string par_close = is_a(c) ? ")}" : ")"; + const std::string imag_sym = is_a(c) ? "i" : "I"; + const std::string mul_sym = is_a(c) ? " " : "*"; + const cln::cl_R r = cln::realpart(cln::the(value)); + const cln::cl_R i = cln::imagpart(cln::the(value)); + if (cln::zerop(i)) { + // case 1, real: x or -x + if ((precedence() <= level) && (!this->is_nonneg_integer())) { + c.s << par_open; + print_real_number(c, r); + c.s << par_close; } else { - if (i == 1) { - os << "I"; - } else { + print_real_number(c, r); + } + } else { + if (cln::zerop(r)) { + // case 2, imaginary: y*I or -y*I + if ((precedence() <= level) && (i < 0)) { if (i == -1) { - os << "-I"; + c.s << par_open+imag_sym+par_close; } else { - print_real_number(os, i); - os << "*I"; + c.s << par_open; + print_real_number(c, i); + c.s << mul_sym+imag_sym+par_close; } - } - } - } else { - // case 3, complex: x+y*I or x-y*I or -x+y*I or -x-y*I - if (precedence <= upper_precedence) - os << "("; - print_real_number(os, r); - if (i < 0) { - if (i == -1) { - os << "-I"; } else { - print_real_number(os, i); - os << "*I"; + if (i == 1) { + c.s << imag_sym; + } else { + if (i == -1) { + c.s << "-" << imag_sym; + } else { + print_real_number(c, i); + c.s << mul_sym+imag_sym; + } + } } } else { - if (i == 1) { - os << "+I"; + // case 3, complex: x+y*I or x-y*I or -x+y*I or -x-y*I + if (precedence() <= level) + c.s << par_open; + print_real_number(c, r); + if (i < 0) { + if (i == -1) { + c.s << "-"+imag_sym; + } else { + print_real_number(c, i); + c.s << mul_sym+imag_sym; + } } else { - os << "+"; - print_real_number(os, i); - os << "*I"; + if (i == 1) { + c.s << "+"+imag_sym; + } else { + c.s << "+"; + print_real_number(c, i); + c.s << mul_sym+imag_sym; + } } + if (precedence() <= level) + c.s << par_close; } - if (precedence <= upper_precedence) - os << ")"; } } } - -void numeric::printraw(std::ostream &os) const -{ - // The method printraw doesn't do much, it simply uses CLN's operator<<() - // for output, which is ugly but reliable. e.g: 2+2i - debugmsg("numeric printraw", LOGLEVEL_PRINT); - os << class_name() << "(" << cln::the(value) << ")"; -} - - -void numeric::printtree(std::ostream &os, unsigned indent) const -{ - debugmsg("numeric printtree", LOGLEVEL_PRINT); - os << std::string(indent,' ') << cln::the(value) - << " (numeric): " - << "hash=" << hashvalue - << " (0x" << std::hex << hashvalue << std::dec << ")" - << ", flags=" << flags << std::endl; -} - - -void numeric::printcsrc(std::ostream &os, unsigned type, unsigned upper_precedence) const -{ - debugmsg("numeric print csrc", LOGLEVEL_PRINT); - std::ios::fmtflags oldflags = os.flags(); - os.setf(std::ios::scientific); - if (this->is_rational() && !this->is_integer()) { - if (compare(_num0()) > 0) { - os << "("; - if (type == csrc_types::ctype_cl_N) - os << "cln::cl_F(\"" << numer().evalf() << "\")"; - else - os << numer().to_double(); - } else { - os << "-("; - if (type == csrc_types::ctype_cl_N) - os << "cln::cl_F(\"" << -numer().evalf() << "\")"; - else - os << -numer().to_double(); - } - os << "/"; - if (type == csrc_types::ctype_cl_N) - os << "cln::cl_F(\"" << denom().evalf() << "\")"; - else - os << denom().to_double(); - os << ")"; - } else { - if (type == csrc_types::ctype_cl_N) - os << "cln::cl_F(\"" << evalf() << "\")"; - else - os << to_double(); - } - os.flags(oldflags); -} - - bool numeric::info(unsigned inf) const { switch (inf) { @@ -552,9 +537,9 @@ bool numeric::info(unsigned inf) const * sign as a multiplicative factor. */ bool numeric::has(const ex &other) const { - if (!is_exactly_of_type(*other.bp, numeric)) + if (!is_ex_exactly_of_type(other, numeric)) return false; - const numeric &o = static_cast(const_cast(*other.bp)); + 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 @@ -598,8 +583,8 @@ ex numeric::evalf(int level) const int numeric::compare_same_type(const basic &other) const { - GINAC_ASSERT(is_exactly_of_type(other, numeric)); - const numeric &o = static_cast(const_cast(other)); + GINAC_ASSERT(is_exactly_a(other)); + const numeric &o = static_cast(other); return this->compare(o); } @@ -607,10 +592,10 @@ int numeric::compare_same_type(const basic &other) const bool numeric::is_equal_same_type(const basic &other) const { - GINAC_ASSERT(is_exactly_of_type(other,numeric)); - const numeric *o = static_cast(&other); + GINAC_ASSERT(is_exactly_a(other)); + const numeric &o = static_cast(other); - return this->is_equal(*o); + return this->is_equal(o); } @@ -641,10 +626,9 @@ unsigned numeric::calchash(void) const const numeric numeric::add(const numeric &other) const { // Efficiency shortcut: trap the neutral element by pointer. - static const numeric * _num0p = &_num0(); - if (this==_num0p) + if (this==_num0_p) return other; - else if (&other==_num0p) + else if (&other==_num0_p) return *this; return numeric(cln::the(value)+cln::the(other.value)); @@ -664,10 +648,9 @@ const numeric numeric::sub(const numeric &other) const const numeric numeric::mul(const numeric &other) const { // Efficiency shortcut: trap the neutral element by pointer. - static const numeric * _num1p = &_num1(); - if (this==_num1p) + if (this==_num1_p) return other; - else if (&other==_num1p) + else if (&other==_num1_p) return *this; return numeric(cln::the(value)*cln::the(other.value)); @@ -691,8 +674,7 @@ const numeric numeric::div(const numeric &other) const const numeric numeric::power(const numeric &other) const { // Efficiency shortcut: trap the neutral exponent by pointer. - static const numeric * _num1p = &_num1(); - if (&other==_num1p) + if (&other==_num1_p) return *this; if (cln::zerop(cln::the(value))) { @@ -703,7 +685,7 @@ const numeric numeric::power(const numeric &other) const else if (cln::minusp(cln::realpart(cln::the(other.value)))) throw std::overflow_error("numeric::eval(): division by zero"); else - return _num0(); + return _num0; } return numeric(cln::expt(cln::the(value),cln::the(other.value))); } @@ -712,10 +694,9 @@ const numeric numeric::power(const numeric &other) const const numeric &numeric::add_dyn(const numeric &other) const { // Efficiency shortcut: trap the neutral element by pointer. - static const numeric * _num0p = &_num0(); - if (this==_num0p) + if (this==_num0_p) return other; - else if (&other==_num0p) + else if (&other==_num0_p) return *this; return static_cast((new numeric(cln::the(value)+cln::the(other.value)))-> @@ -733,10 +714,9 @@ const numeric &numeric::sub_dyn(const numeric &other) const const numeric &numeric::mul_dyn(const numeric &other) const { // Efficiency shortcut: trap the neutral element by pointer. - static const numeric * _num1p = &_num1(); - if (this==_num1p) + if (this==_num1_p) return other; - else if (&other==_num1p) + else if (&other==_num1_p) return *this; return static_cast((new numeric(cln::the(value)*cln::the(other.value)))-> @@ -756,8 +736,7 @@ const numeric &numeric::div_dyn(const numeric &other) const const numeric &numeric::power_dyn(const numeric &other) const { // Efficiency shortcut: trap the neutral exponent by pointer. - static const numeric * _num1p=&_num1(); - if (&other==_num1p) + if (&other==_num1_p) return *this; if (cln::zerop(cln::the(value))) { @@ -768,7 +747,7 @@ const numeric &numeric::power_dyn(const numeric &other) const else if (cln::minusp(cln::realpart(cln::the(other.value)))) throw std::overflow_error("numeric::eval(): division by zero"); else - return _num0(); + return _num0; } return static_cast((new numeric(cln::expt(cln::the(value),cln::the(other.value))))-> setflag(status_flags::dynallocated)); @@ -961,13 +940,13 @@ bool numeric::is_real(void) const bool numeric::operator==(const numeric &other) const { - return equal(cln::the(value), cln::the(other.value)); + return cln::equal(cln::the(value), cln::the(other.value)); } bool numeric::operator!=(const numeric &other) const { - return !equal(cln::the(value), cln::the(other.value)); + return !cln::equal(cln::the(value), cln::the(other.value)); } @@ -1135,16 +1114,16 @@ const numeric numeric::numer(void) const const numeric numeric::denom(void) const { if (this->is_integer()) - return _num1(); + return _num1; - if (instanceof(value, cln::cl_RA_ring)) + 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))); if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_I_ring)) - return _num1(); + return _num1; 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)) @@ -1153,7 +1132,7 @@ const numeric numeric::denom(void) const return numeric(cln::lcm(cln::denominator(r), cln::denominator(i))); } // at least one float encountered - return _num1(); + return _num1; } @@ -1171,15 +1150,6 @@ int numeric::int_length(void) const return 0; } - -////////// -// static member variables -////////// - -// protected - -unsigned numeric::precedence = 30; - ////////// // global constants ////////// @@ -1266,7 +1236,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)) throw pole_error("atan(): logarithmic pole",0); return cln::atan(x.to_cl_N()); } @@ -1417,7 +1387,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; // what is the desired float format? // first guess: default format @@ -1456,10 +1426,7 @@ const numeric zeta(const numeric &x) if (cln::zerop(x.to_cl_N()-aux)) return cln::zeta(aux); } - std::clog << "zeta(" << x - << "): Does anybody know a good way to calculate this numerically?" - << std::endl; - return numeric(0); + throw dunno(); } @@ -1467,17 +1434,11 @@ const numeric zeta(const numeric &x) * This is only a stub! */ const numeric lgamma(const numeric &x) { - std::clog << "lgamma(" << x - << "): Does anybody know a good way to calculate this numerically?" - << std::endl; - return numeric(0); + throw dunno(); } const numeric tgamma(const numeric &x) { - std::clog << "tgamma(" << x - << "): Does anybody know a good way to calculate this numerically?" - << std::endl; - return numeric(0); + throw dunno(); } @@ -1485,10 +1446,7 @@ const numeric tgamma(const numeric &x) * This is only a stub! */ const numeric psi(const numeric &x) { - std::clog << "psi(" << x - << "): Does anybody know a good way to calculate this numerically?" - << std::endl; - return numeric(0); + throw dunno(); } @@ -1496,10 +1454,7 @@ const numeric psi(const numeric &x) * This is only a stub! */ const numeric psi(const numeric &n, const numeric &x) { - std::clog << "psi(" << n << "," << x - << "): Does anybody know a good way to calculate this numerically?" - << std::endl; - return numeric(0); + throw dunno(); } @@ -1523,8 +1478,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)) + return _num1; if (!n.is_nonneg_integer()) throw std::range_error("numeric::doublefactorial(): argument must be integer >= -1"); @@ -1541,12 +1496,12 @@ 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)!=-1) return numeric(cln::binomial(n.to_int(),k.to_int())); else - return _num0(); + return _num0; } else { - return _num_1().power(k)*binomial(k-n-_num1(),k); + return _num_1.power(k)*binomial(k-n-_num1,k); } } @@ -1597,16 +1552,16 @@ const numeric bernoulli(const numeric &nn) // we don't use it.) // the special cases not covered by the algorithm below - if (nn.is_equal(_num1())) - return _num_1_2(); + if (nn.is_equal(_num1)) + return _num_1_2; if (nn.is_odd()) - return _num0(); + return _num0; // store nonvanishing Bernoulli numbers here static std::vector< cln::cl_RA > results; static int highest_result = 0; // algorithm not applicable to B(0), so just store it - if (results.size()==0) + if (results.empty()) results.push_back(cln::cl_RA(1)); int n = nn.to_long(); @@ -1658,7 +1613,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; if (n.is_negative()) if (n.is_even()) return -fibonacci(-n); @@ -1710,7 +1665,7 @@ 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; } @@ -1725,7 +1680,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; } @@ -1741,7 +1696,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; } @@ -1760,8 +1715,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; + return _num0; } } @@ -1773,10 +1728,10 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q) const numeric iquo(const numeric &a, const numeric &b) { if (a.is_integer() && b.is_integer()) - return truncate1(cln::the(a.to_cl_N()), - cln::the(b.to_cl_N())); + return cln::truncate1(cln::the(a.to_cl_N()), + cln::the(b.to_cl_N())); else - return _num0(); + return _num0; } @@ -1794,8 +1749,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; + return _num0; } } @@ -1810,7 +1765,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; } @@ -1850,7 +1805,7 @@ const numeric isqrt(const numeric &x) cln::isqrt(cln::the(x.to_cl_N()), &root); return root; } else - return _num0(); + return _num0; }