X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=e621607a730cabb8878cd5b5cf7d2c0b0dd71b54;hp=72668220da202597a09fe11564fdc4517e5a05db;hb=c86cff51ac5a42f86387ef7bc767f1274137350b;hpb=0c8c12f5393a7ca5f78233af5ee81593184c35fb diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 72668220..e621607a 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -35,7 +35,6 @@ #include "ex.h" #include "print.h" #include "archive.h" -#include "debugmsg.h" #include "tostring.h" #include "utils.h" @@ -63,14 +62,12 @@ namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS(numeric, basic) ////////// -// default ctor, dtor, copy ctor assignment -// operator and helpers +// default ctor, dtor, copy ctor, assignment operator and helpers ////////// /** default ctor. Numerically it initializes to an integer zero. */ numeric::numeric() : basic(TINFO_numeric) { - debugmsg("numeric default ctor", LOGLEVEL_CONSTRUCT); value = cln::cl_I(0); setflag(status_flags::evaluated | status_flags::expanded); } @@ -91,7 +88,6 @@ DEFAULT_DESTROY(numeric) numeric::numeric(int i) : basic(TINFO_numeric) { - debugmsg("numeric ctor from int",LOGLEVEL_CONSTRUCT); // 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 // emphasizes efficiency. However, if the integer is small enough, @@ -107,7 +103,6 @@ numeric::numeric(int i) : basic(TINFO_numeric) numeric::numeric(unsigned int i) : basic(TINFO_numeric) { - debugmsg("numeric ctor from uint",LOGLEVEL_CONSTRUCT); // 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, @@ -123,7 +118,6 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric) numeric::numeric(long i) : basic(TINFO_numeric) { - debugmsg("numeric ctor from long",LOGLEVEL_CONSTRUCT); value = cln::cl_I(i); setflag(status_flags::evaluated | status_flags::expanded); } @@ -131,7 +125,6 @@ numeric::numeric(long i) : basic(TINFO_numeric) numeric::numeric(unsigned long i) : basic(TINFO_numeric) { - debugmsg("numeric ctor from ulong",LOGLEVEL_CONSTRUCT); value = cln::cl_I(i); setflag(status_flags::evaluated | status_flags::expanded); } @@ -141,7 +134,6 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric) * @exception overflow_error (division by zero) */ numeric::numeric(long numer, long denom) : basic(TINFO_numeric) { - debugmsg("numeric ctor from long/long",LOGLEVEL_CONSTRUCT); if (!denom) throw std::overflow_error("division by zero"); value = cln::cl_I(numer) / cln::cl_I(denom); @@ -151,7 +143,6 @@ numeric::numeric(long numer, long denom) : basic(TINFO_numeric) numeric::numeric(double d) : basic(TINFO_numeric) { - debugmsg("numeric ctor from double",LOGLEVEL_CONSTRUCT); // 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 // will not be promoted to cl_LF if overflow occurs: @@ -164,40 +155,47 @@ numeric::numeric(double d) : basic(TINFO_numeric) * notation like "2+5*I". */ numeric::numeric(const char *s) : basic(TINFO_numeric) { - debugmsg("numeric ctor from string",LOGLEVEL_CONSTRUCT); cln::cl_N ctorval = 0; // parse complex numbers (functional but not completely safe, unfortunately // std::string does not understand regexpese): // ss should represent a simple sum like 2+5*I - std::string ss(s); - // make it safe by adding explicit sign + std::string ss = s; + std::string::size_type delim; + + // make this implementation safe by adding explicit sign if (ss.at(0) != '+' && ss.at(0) != '-' && ss.at(0) != '#') ss = '+' + ss; - std::string::size_type delim; + + // We use 'E' as exponent marker in the output, but some people insist on + // writing 'e' at input, so let's substitute them right at the beginning: + while ((delim = ss.find("e"))!=std::string::npos) + ss.replace(delim,1,"E"); + + // main parser loop: do { // chop ss into terms from left to right std::string term; bool imaginary = false; delim = ss.find_first_of(std::string("+-"),1); // Do we have an exponent marker like "31.415E-1"? If so, hop on! - if ((delim != std::string::npos) && (ss.at(delim-1) == 'E')) + if (delim!=std::string::npos && ss.at(delim-1)=='E') delim = ss.find_first_of(std::string("+-"),delim+1); term = ss.substr(0,delim); - if (delim != std::string::npos) + if (delim!=std::string::npos) ss = ss.substr(delim); // is the term imaginary? - if (term.find("I") != std::string::npos) { + if (term.find("I")!=std::string::npos) { // erase 'I': - term = term.replace(term.find("I"),1,""); + term.erase(term.find("I"),1); // erase '*': - if (term.find("*") != std::string::npos) - term = term.replace(term.find("*"),1,""); + if (term.find("*")!=std::string::npos) + term.erase(term.find("*"),1); // correct for trivial +/-I without explicit factor on I: - if (term.size() == 1) - term += "1"; + if (term.size()==1) + term += '1'; imaginary = true; } - if (term.find(".") != std::string::npos) { + if (term.find('.')!=std::string::npos || term.find('E')!=std::string::npos) { // CLN's short type cl_SF is not very useful within the GiNaC // framework where we are mainly interested in the arbitrary // precision type cl_LF. Hence we go straight to the construction @@ -208,7 +206,7 @@ numeric::numeric(const char *s) : basic(TINFO_numeric) // 31.4E-1 --> 31.4e-1_ // and s on. // No exponent marker? Let's add a trivial one. - if (term.find("E") == std::string::npos) + if (term.find("E")==std::string::npos) term += "E0"; // E to lower case term = term.replace(term.find("E"),1,"e"); @@ -220,13 +218,13 @@ numeric::numeric(const char *s) : basic(TINFO_numeric) else ctorval = ctorval + cln::cl_F(term.c_str()); } else { - // not a floating point number... + // this is not a floating point number... if (imaginary) ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_R(term.c_str())); else ctorval = ctorval + cln::cl_R(term.c_str()); } - } while(delim != std::string::npos); + } while (delim != std::string::npos); value = ctorval; setflag(status_flags::evaluated | status_flags::expanded); } @@ -236,7 +234,6 @@ numeric::numeric(const char *s) : basic(TINFO_numeric) * only. */ numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric) { - debugmsg("numeric ctor from cl_N", LOGLEVEL_CONSTRUCT); value = z; setflag(status_flags::evaluated | status_flags::expanded); } @@ -247,7 +244,6 @@ numeric::numeric(const cln::cl_N &z) : basic(TINFO_numeric) numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { - debugmsg("numeric ctor from archive_node", LOGLEVEL_CONSTRUCT); cln::cl_N ctorval = 0; // Read number as string @@ -347,8 +343,6 @@ static void print_real_number(const print_context & c, const cln::cl_R &x) * @see print_real_number() */ void numeric::print(const print_context & c, unsigned level) const { - debugmsg("numeric print", LOGLEVEL_PRINT); - if (is_a(c)) { c.s << std::string(level, ' ') << cln::the(value) @@ -395,6 +389,8 @@ void numeric::print(const print_context & c, unsigned level) const 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 (is_a(c)) + c.s << class_name() << "('"; if (cln::zerop(i)) { // case 1, real: x or -x if ((precedence() <= level) && (!this->is_nonneg_integer())) { @@ -452,6 +448,8 @@ void numeric::print(const print_context & c, unsigned level) const c.s << par_close; } } + if (is_a(c)) + c.s << "')"; } } @@ -1490,7 +1488,7 @@ const numeric bernoulli(const numeric &nn) { if (!nn.is_integer() || nn.is_negative()) throw std::range_error("numeric::bernoulli(): argument must be integer >= 0"); - + // Method: // // The Bernoulli numbers are rational numbers that may be computed using @@ -1514,46 +1512,49 @@ const numeric bernoulli(const numeric &nn) // But if somebody works with the n'th Bernoulli number she is likely to // also need all previous Bernoulli numbers. So we need a complete remember // table and above divide and conquer algorithm is not suited to build one - // up. The code below is adapted from Pari's function bernvec(). + // up. The formula below accomplishes this. It is a modification of the + // defining formula above but the computation of the binomial coefficients + // is carried along in an inline fashion. It also honors the fact that + // B_n is zero when n is odd and greater than 1. // // (There is an interesting relation with the tangent polynomials described - // in `Concrete Mathematics', which leads to a program twice as fast as our - // implementation below, but it requires storing one such polynomial in + // in `Concrete Mathematics', which leads to a program a little faster as + // our implementation below, but it requires storing one such polynomial in // addition to the remember table. This doubles the memory footprint so // we don't use it.) - + + const unsigned n = nn.to_int(); + // the special cases not covered by the algorithm below - if (nn.is_equal(_num1)) - return _num_1_2; - if (nn.is_odd()) - return _num0; - + if (n & 1) + return (n==1) ? _num_1_2 : _num0; + if (!n) + return _num1; + // 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.empty()) - results.push_back(cln::cl_RA(1)); - - int n = nn.to_long(); - for (int i=highest_result; i0; --j) { - B = cln::cl_I(n*m) * (B+results[j]) / (d1*d2); - n += 4; - m += 2; - d1 -= 1; - d2 -= 2; + static unsigned next_r = 0; + + // algorithm not applicable to B(2), so just store it + if (!next_r) { + results.push_back(cln::recip(cln::cl_RA(6))); + next_r = 4; + } + for (unsigned p=next_r; p<=n; p+=2) { + cln::cl_I c = 1; + cln::cl_RA b = cln::cl_RA(1-p)/2; + const unsigned p3 = p+3; + const unsigned p2 = p+2; + const unsigned pm = p-2; + unsigned i, k; + for (i=2, k=0; i <= pm; i += 2, k++) { + c = cln::exquo(c * ((p3 - i)*(p2 - i)), (i - 1)*i); + b = b + c * results[k]; } - B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2)); - results.push_back(B); - ++highest_result; + results.push_back(-b / (p+1)); + next_r += 2; } - return results[n/2]; + return results[n/2 - 1]; } @@ -1835,7 +1836,6 @@ _numeric_digits::operator long() /** Append global Digits object to ostream. */ void _numeric_digits::print(std::ostream &os) const { - debugmsg("_numeric_digits print", LOGLEVEL_PRINT); os << digits; }