* 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-2003 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
{
// 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. <cln/object.h>)
+ if (i < (1U<<cl_value_len-1))
value = cln::cl_I(i);
else
value = cln::cl_I((long) i);
{
// 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,
- // 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. <cln/object.h>)
+ if (i < (1U<<cl_value_len-1))
value = cln::cl_I(i);
else
value = cln::cl_I((unsigned long) i);
// 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
// 31.4E-1 --> 31.4e-1_<Digits>
// 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");
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);
}
!is_a<print_latex>(c)) {
cln::print_real(c.s, ourflags, x);
} else { // rational output in LaTeX context
+ if (x < 0)
+ c.s << "-";
c.s << "\\frac{";
- cln::print_real(c.s, ourflags, cln::numerator(cln::the<cln::cl_RA>(x)));
+ cln::print_real(c.s, ourflags, cln::abs(cln::numerator(cln::the<cln::cl_RA>(x))));
c.s << "}{";
cln::print_real(c.s, ourflags, cln::denominator(cln::the<cln::cl_RA>(x)));
c.s << '}';
std::ios::fmtflags oldflags = c.s.flags();
c.s.setf(std::ios::scientific);
- if (this->is_rational() && !this->is_integer()) {
+ int oldprec = c.s.precision();
+ if (is_a<print_csrc_double>(c))
+ c.s.precision(16);
+ else
+ c.s.precision(7);
+ if (is_a<print_csrc_cl_N>(c) && this->is_integer()) {
+ c.s << "cln::cl_I(\"";
+ const cln::cl_R r = cln::realpart(cln::the<cln::cl_N>(value));
+ print_real_number(c,r);
+ c.s << "\")";
+ } else if (this->is_rational() && !this->is_integer()) {
if (compare(_num0) > 0) {
c.s << "(";
if (is_a<print_csrc_cl_N>(c))
c.s << ")";
} else {
if (is_a<print_csrc_cl_N>(c))
- c.s << "cln::cl_F(\"" << evalf() << "\")";
+ c.s << "cln::cl_F(\"" << evalf() << "_" << Digits << "\")";
else
c.s << to_double();
}
c.s.flags(oldflags);
+ c.s.precision(oldprec);
} else {
const std::string par_open = is_a<print_latex>(c) ? "{(" : "(";
const std::string mul_sym = is_a<print_latex>(c) ? " " : "*";
const cln::cl_R r = cln::realpart(cln::the<cln::cl_N>(value));
const cln::cl_R i = cln::imagpart(cln::the<cln::cl_N>(value));
+ if (is_a<print_python_repr>(c))
+ c.s << class_name() << "('";
if (cln::zerop(i)) {
// case 1, real: x or -x
if ((precedence() <= level) && (!this->is_nonneg_integer())) {
} 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
c.s << par_close;
}
}
+ if (is_a<print_python_repr>(c))
+ c.s << "')";
}
}
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
{
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
// 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; i<n/2; ++i) {
- cln::cl_RA B = 0;
- long n = 8;
- long m = 5;
- long d1 = i;
- long d2 = 2*i-1;
- for (int j=i; j>0; --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;
}
- return results[n/2];
+ if (n<next_r)
+ return results[n/2-1];
+
+ 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/object.h>)
+ if (p < (1UL<<cl_value_len/2)) {
+ for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+ c = cln::exquo(c * ((p3-i) * p_2), (i-1)*k);
+ b = b + c*results[k-1];
+ }
+ } else {
+ for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+ c = cln::exquo((c * (p3-i)) * p_2, cln::cl_I(i-1)*k);
+ b = b + c*results[k-1];
+ }
+ }
+ results.push_back(-b/(p+1));
+ }
+ next_r = n+2;
+ return results[n/2-1];
}
* In general, mod(a,b) has the sign of b or is zero, and irem(a,b) has the
* sign of a or is zero.
*
- * @return remainder of a/b if both are integer, 0 otherwise. */
+ * @return remainder of a/b if both are integer, 0 otherwise.
+ * @exception overflow_error (division by zero) if b is zero. */
const numeric irem(const numeric &a, const numeric &b)
{
+ if (b.is_zero())
+ throw std::overflow_error("numeric::irem(): division by zero");
if (a.is_integer() && b.is_integer())
return cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
cln::the<cln::cl_I>(b.to_cl_N()));
* and irem(a,b) has the sign of a or is zero.
*
* @return remainder of a/b and quotient stored in q if both are integer,
- * 0 otherwise. */
+ * 0 otherwise.
+ * @exception overflow_error (division by zero) if b is zero. */
const numeric irem(const numeric &a, const numeric &b, numeric &q)
{
+ if (b.is_zero())
+ throw std::overflow_error("numeric::irem(): division by zero");
if (a.is_integer() && b.is_integer()) {
const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
cln::the<cln::cl_I>(b.to_cl_N()));
/** Numeric integer quotient.
* Equivalent to Maple's iquo as far as sign conventions are concerned.
*
- * @return truncated quotient of a/b if both are integer, 0 otherwise. */
+ * @return truncated quotient of a/b if both are integer, 0 otherwise.
+ * @exception overflow_error (division by zero) if b is zero. */
const numeric iquo(const numeric &a, const numeric &b)
{
+ if (b.is_zero())
+ throw std::overflow_error("numeric::iquo(): division by zero");
if (a.is_integer() && b.is_integer())
return cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
cln::the<cln::cl_I>(b.to_cl_N()));
* r == a - iquo(a,b,r)*b.
*
* @return truncated quotient of a/b and remainder stored in r if both are
- * integer, 0 otherwise. */
+ * integer, 0 otherwise.
+ * @exception overflow_error (division by zero) if b is zero. */
const numeric iquo(const numeric &a, const numeric &b, numeric &r)
{
+ if (b.is_zero())
+ throw std::overflow_error("numeric::iquo(): division by zero");
if (a.is_integer() && b.is_integer()) {
const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
cln::the<cln::cl_I>(b.to_cl_N()));
git://www.ginac.de / ginac.git / blobdiff