X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=000e15f8eb53acbb71277692a103dc28ca95b949;hb=bd05a471b5e57cdb65458b8d0af54f1743831103;hp=72668220da202597a09fe11564fdc4517e5a05db;hpb=0c8c12f5393a7ca5f78233af5ee81593184c35fb;p=ginac.git diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 72668220..000e15f8 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-2001 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2002 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 @@ -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,13 +88,12 @@ 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, - // i.e. satisfies cl_immediate_p(), we save space and dereferences by - // using an immediate type: - if (cln::cl_immediate_p(i)) + // emphasizes efficiency. However, if the integer is small enough + // we save space and dereferences by using an immediate type. + // (C.f. ) + if (i < (1U<) + if (i < (1U< 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())) { @@ -407,25 +403,19 @@ void numeric::print(const print_context & c, unsigned level) const } else { if (cln::zerop(r)) { // case 2, imaginary: y*I or -y*I - if ((precedence() <= level) && (i < 0)) { - if (i == -1) { - c.s << par_open+imag_sym+par_close; - } else { + if (i==1) + c.s << imag_sym; + else { + if (precedence()<=level) c.s << par_open; + if (i == -1) + c.s << "-" << imag_sym; + else { print_real_number(c, i); - c.s << mul_sym+imag_sym+par_close; - } - } else { - 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; - } + c.s << mul_sym+imag_sym; } + if (precedence()<=level) + c.s << par_close; } } else { // case 3, complex: x+y*I or x-y*I or -x+y*I or -x-y*I @@ -452,6 +442,8 @@ void numeric::print(const print_context & c, unsigned level) const c.s << par_close; } } + if (is_a(c)) + c.s << "')"; } } @@ -500,6 +492,21 @@ bool numeric::info(unsigned inf) const return false; } +int numeric::degree(const ex & s) const +{ + return 0; +} + +int numeric::ldegree(const ex & s) const +{ + return 0; +} + +ex numeric::coeff(const ex & s, int n) const +{ + return n==0 ? *this : _ex0; +} + /** Disassemble real part and imaginary part to scan for the occurrence of a * single number. Also handles the imaginary unit. It ignores the sign on * both this and the argument, which may lead to what might appear as funny @@ -1490,7 +1497,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 +1521,61 @@ 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; - } - B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2)); - results.push_back(B); - ++highest_result; + 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; + } + if (n) + if (p < (1UL<