* of special functions or implement the interface to the bignum package. */
/*
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2000 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
#include <vector>
#include <stdexcept>
+#include <string>
+#include <strstream> //!!
#include "numeric.h"
#include "ex.h"
#include "config.h"
+#include "archive.h"
#include "debugmsg.h"
#include "utils.h"
// CLN should not pollute the global namespace, hence we include it here
-// instead of in some header file where it would propagate to other parts:
+// instead of in some header file where it would propagate to other parts.
+// Also, we only need a subset of CLN, so we don't include the complete cln.h:
#ifdef HAVE_CLN_CLN_H
-#include <CLN/cln.h>
+#include <cln/cl_integer_io.h>
+#include <cln/cl_integer_ring.h>
+#include <cln/cl_rational_io.h>
+#include <cln/cl_rational_ring.h>
+#include <cln/cl_lfloat_class.h>
+#include <cln/cl_lfloat_io.h>
+#include <cln/cl_real_io.h>
+#include <cln/cl_real_ring.h>
+#include <cln/cl_complex_io.h>
+#include <cln/cl_complex_ring.h>
+#include <cln/cl_numtheory.h>
#else
-#include <cln.h>
+#include <cl_integer_io.h>
+#include <cl_integer_ring.h>
+#include <cl_rational_io.h>
+#include <cl_rational_ring.h>
+#include <cl_lfloat_class.h>
+#include <cl_lfloat_io.h>
+#include <cl_real_io.h>
+#include <cl_real_ring.h>
+#include <cl_complex_io.h>
+#include <cl_complex_ring.h>
+#include <cl_numtheory.h>
#endif
#ifndef NO_GINAC_NAMESPACE
// linker has no problems finding text symbols for numerator or denominator
//#define SANE_LINKER
+GINAC_IMPLEMENT_REGISTERED_CLASS(numeric, basic)
+
//////////
// default constructor, destructor, copy constructor assignment
// operator and helpers
destroy(0);
}
-numeric::numeric(numeric const & other)
+numeric::numeric(const numeric & other)
{
debugmsg("numeric copy constructor", LOGLEVEL_CONSTRUCT);
copy(other);
}
-numeric const & numeric::operator=(numeric const & other)
+const numeric & numeric::operator=(const numeric & other)
{
debugmsg("numeric operator=", LOGLEVEL_ASSIGNMENT);
if (this != &other) {
// protected
-void numeric::copy(numeric const & other)
+void numeric::copy(const numeric & other)
{
basic::copy(other);
value = new cl_N(*other.value);
status_flags::hash_calculated);
}
+//////////
+// archiving
+//////////
+
+/** Construct object from archive_node. */
+numeric::numeric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+{
+ debugmsg("numeric constructor from archive_node", LOGLEVEL_CONSTRUCT);
+ value = new cl_N;
+#if 0 //!!
+ // This is how it should be implemented but we have no istringstream here...
+ string str;
+ if (n.find_string("number", str)) {
+ istringstream s(str);
+ s >> *value;
+ }
+#else
+ // Workaround for the above: read from strstream
+ string str;
+ if (n.find_string("number", str)) {
+ istrstream f(str.c_str(), str.size() + 1);
+ f >> *value;
+ }
+#endif
+ calchash();
+ setflag(status_flags::evaluated|
+ status_flags::hash_calculated);
+}
+
+/** 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);
+#if 0 //!!
+ // This is how it should be implemented but we have no ostringstream here...
+ ostringstream s;
+ s << *value;
+ n.add_string("number", s.str());
+#else
+ // Workaround for the above: write to strstream
+ char buf[1024];
+ ostrstream f(buf, 1024);
+ f << *value << ends;
+ string str(buf);
+ n.add_string("number", str);
+#endif
+}
+
//////////
// functions overriding virtual functions from bases classes
//////////
return new numeric(*this);
}
-// The method printraw doesn't do much, it simply uses CLN's operator<<() for
-// output, which is ugly but reliable. Examples:
-// 2+2i
-void numeric::printraw(ostream & os) const
-{
- debugmsg("numeric printraw", LOGLEVEL_PRINT);
- os << "numeric(" << *value << ")";
-}
-
-// The method print adds to the output so it blends more consistently together
-// with the other routines and produces something compatible to Maple input.
void numeric::print(ostream & os, unsigned upper_precedence) const
{
+ // The method print adds to the output so it blends more consistently
+ // together with the other routines and produces something compatible to
+ // ginsh input.
debugmsg("numeric print", LOGLEVEL_PRINT);
if (is_real()) {
// case 1, real: x or -x
}
}
+
+void numeric::printraw(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 << "numeric(" << *value << ")";
+}
void numeric::printtree(ostream & os, unsigned indent) const
{
debugmsg("numeric printtree", LOGLEVEL_PRINT);
int numeric::compare_same_type(basic const & other) const
{
GINAC_ASSERT(is_exactly_of_type(other, numeric));
- numeric const & o = static_cast<numeric &>(const_cast<basic &>(other));
+ const numeric & o = static_cast<numeric &>(const_cast<basic &>(other));
if (*value == *o.value) {
return 0;
bool numeric::is_equal_same_type(basic const & other) const
{
GINAC_ASSERT(is_exactly_of_type(other,numeric));
- numeric const *o = static_cast<numeric const *>(&other);
+ const numeric *o = static_cast<const numeric *>(&other);
return is_equal(*o);
}
/** Numerical addition method. Adds argument to *this and returns result as
* a new numeric object. */
-numeric numeric::add(numeric const & other) const
+numeric numeric::add(const numeric & other) const
{
return numeric((*value)+(*other.value));
}
/** Numerical subtraction method. Subtracts argument from *this and returns
* result as a new numeric object. */
-numeric numeric::sub(numeric const & other) const
+numeric numeric::sub(const numeric & other) const
{
return numeric((*value)-(*other.value));
}
/** Numerical multiplication method. Multiplies *this and argument and returns
* result as a new numeric object. */
-numeric numeric::mul(numeric const & other) const
+numeric numeric::mul(const numeric & other) const
{
static const numeric * _num1p=&_num1();
if (this==_num1p) {
* a new numeric object.
*
* @exception overflow_error (division by zero) */
-numeric numeric::div(numeric const & other) const
+numeric numeric::div(const numeric & other) const
{
if (::zerop(*other.value))
throw (std::overflow_error("division by zero"));
return numeric((*value)/(*other.value));
}
-numeric numeric::power(numeric const & other) const
+numeric numeric::power(const numeric & other) const
{
static const numeric * _num1p=&_num1();
if (&other==_num1p)
return numeric(::recip(*value)); // -> CLN
}
-numeric const & numeric::add_dyn(numeric const & other) const
+const numeric & numeric::add_dyn(const numeric & other) const
{
- return static_cast<numeric const &>((new numeric((*value)+(*other.value)))->
+ return static_cast<const numeric &>((new numeric((*value)+(*other.value)))->
setflag(status_flags::dynallocated));
}
-numeric const & numeric::sub_dyn(numeric const & other) const
+const numeric & numeric::sub_dyn(const numeric & other) const
{
- return static_cast<numeric const &>((new numeric((*value)-(*other.value)))->
+ return static_cast<const numeric &>((new numeric((*value)-(*other.value)))->
setflag(status_flags::dynallocated));
}
-numeric const & numeric::mul_dyn(numeric const & other) const
+const numeric & numeric::mul_dyn(const numeric & other) const
{
static const numeric * _num1p=&_num1();
if (this==_num1p) {
} else if (&other==_num1p) {
return *this;
}
- return static_cast<numeric const &>((new numeric((*value)*(*other.value)))->
+ return static_cast<const numeric &>((new numeric((*value)*(*other.value)))->
setflag(status_flags::dynallocated));
}
-numeric const & numeric::div_dyn(numeric const & other) const
+const numeric & numeric::div_dyn(const numeric & other) const
{
if (::zerop(*other.value))
throw (std::overflow_error("division by zero"));
- return static_cast<numeric const &>((new numeric((*value)/(*other.value)))->
+ return static_cast<const numeric &>((new numeric((*value)/(*other.value)))->
setflag(status_flags::dynallocated));
}
-numeric const & numeric::power_dyn(numeric const & other) const
+const numeric & numeric::power_dyn(const numeric & other) const
{
static const numeric * _num1p=&_num1();
if (&other==_num1p)
return *this;
if (::zerop(*value) && other.is_real() && ::minusp(realpart(*other.value)))
throw (std::overflow_error("division by zero"));
- return static_cast<numeric const &>((new numeric(::expt(*value,*other.value)))->
+ return static_cast<const numeric &>((new numeric(::expt(*value,*other.value)))->
setflag(status_flags::dynallocated));
}
-numeric const & numeric::operator=(int i)
+const numeric & numeric::operator=(int i)
{
return operator=(numeric(i));
}
-numeric const & numeric::operator=(unsigned int i)
+const numeric & numeric::operator=(unsigned int i)
{
return operator=(numeric(i));
}
-numeric const & numeric::operator=(long i)
+const numeric & numeric::operator=(long i)
{
return operator=(numeric(i));
}
-numeric const & numeric::operator=(unsigned long i)
+const numeric & numeric::operator=(unsigned long i)
{
return operator=(numeric(i));
}
-numeric const & numeric::operator=(double d)
+const numeric & numeric::operator=(double d)
{
return operator=(numeric(d));
}
-numeric const & numeric::operator=(char const * s)
+const numeric & numeric::operator=(char const * s)
{
return operator=(numeric(s));
}
* csgn(x)==0 for x==0, csgn(x)==1 for Re(x)>0 or Re(x)=0 and Im(x)>0,
* csgn(x)==-1 for Re(x)<0 or Re(x)=0 and Im(x)<0.
*
- * @see numeric::compare(numeric const & other) */
+ * @see numeric::compare(const numeric & other) */
int numeric::csgn(void) const
{
if (is_zero())
*
* @return csgn(*this-other)
* @see numeric::csgn(void) */
-int numeric::compare(numeric const & other) const
+int numeric::compare(const numeric & other) const
{
// Comparing two real numbers?
if (is_real() && other.is_real())
}
}
-bool numeric::is_equal(numeric const & other) const
+bool numeric::is_equal(const numeric & other) const
{
return (*value == *other.value);
}
return ::instanceof(*value, cl_R_ring); // -> CLN
}
-bool numeric::operator==(numeric const & other) const
+bool numeric::operator==(const numeric & other) const
{
return (*value == *other.value); // -> CLN
}
-bool numeric::operator!=(numeric const & other) const
+bool numeric::operator!=(const numeric & other) const
{
return (*value != *other.value); // -> CLN
}
/** Numerical comparison: less.
*
* @exception invalid_argument (complex inequality) */
-bool numeric::operator<(numeric const & other) const
+bool numeric::operator<(const numeric & other) const
{
if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) < The(cl_R)(*other.value)); // -> CLN
/** Numerical comparison: less or equal.
*
* @exception invalid_argument (complex inequality) */
-bool numeric::operator<=(numeric const & other) const
+bool numeric::operator<=(const numeric & other) const
{
if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) <= The(cl_R)(*other.value)); // -> CLN
/** Numerical comparison: greater.
*
* @exception invalid_argument (complex inequality) */
-bool numeric::operator>(numeric const & other) const
+bool numeric::operator>(const numeric & other) const
{
if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) > The(cl_R)(*other.value)); // -> CLN
/** Numerical comparison: greater or equal.
*
* @exception invalid_argument (complex inequality) */
-bool numeric::operator>=(numeric const & other) const
+bool numeric::operator>=(const numeric & other) const
{
if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) >= The(cl_R)(*other.value)); // -> CLN
/** Exponential function.
*
* @return arbitrary precision numerical exp(x). */
-numeric exp(numeric const & x)
+numeric exp(const numeric & x)
{
return ::exp(*x.value); // -> CLN
}
* @param z complex number
* @return arbitrary precision numerical log(x).
* @exception overflow_error (logarithmic singularity) */
-numeric log(numeric const & z)
+numeric log(const numeric & z)
{
if (z.is_zero())
throw (std::overflow_error("log(): logarithmic singularity"));
/** Numeric sine (trigonometric function).
*
* @return arbitrary precision numerical sin(x). */
-numeric sin(numeric const & x)
+numeric sin(const numeric & x)
{
return ::sin(*x.value); // -> CLN
}
/** Numeric cosine (trigonometric function).
*
* @return arbitrary precision numerical cos(x). */
-numeric cos(numeric const & x)
+numeric cos(const numeric & x)
{
return ::cos(*x.value); // -> CLN
}
/** Numeric tangent (trigonometric function).
*
* @return arbitrary precision numerical tan(x). */
-numeric tan(numeric const & x)
+numeric tan(const numeric & x)
{
return ::tan(*x.value); // -> CLN
}
/** Numeric inverse sine (trigonometric function).
*
* @return arbitrary precision numerical asin(x). */
-numeric asin(numeric const & x)
+numeric asin(const numeric & x)
{
return ::asin(*x.value); // -> CLN
}
/** Numeric inverse cosine (trigonometric function).
*
* @return arbitrary precision numerical acos(x). */
-numeric acos(numeric const & x)
+numeric acos(const numeric & x)
{
return ::acos(*x.value); // -> CLN
}
* @param z complex number
* @return atan(z)
* @exception overflow_error (logarithmic singularity) */
-numeric atan(numeric const & x)
+numeric atan(const numeric & x)
{
if (!x.is_real() &&
x.real().is_zero() &&
* @param x real number
* @param y real number
* @return atan(y/x) */
-numeric atan(numeric const & y, numeric const & x)
+numeric atan(const numeric & y, const numeric & x)
{
if (x.is_real() && y.is_real())
return ::atan(realpart(*x.value), realpart(*y.value)); // -> CLN
/** Numeric hyperbolic sine (trigonometric function).
*
* @return arbitrary precision numerical sinh(x). */
-numeric sinh(numeric const & x)
+numeric sinh(const numeric & x)
{
return ::sinh(*x.value); // -> CLN
}
/** Numeric hyperbolic cosine (trigonometric function).
*
* @return arbitrary precision numerical cosh(x). */
-numeric cosh(numeric const & x)
+numeric cosh(const numeric & x)
{
return ::cosh(*x.value); // -> CLN
}
/** Numeric hyperbolic tangent (trigonometric function).
*
* @return arbitrary precision numerical tanh(x). */
-numeric tanh(numeric const & x)
+numeric tanh(const numeric & x)
{
return ::tanh(*x.value); // -> CLN
}
/** Numeric inverse hyperbolic sine (trigonometric function).
*
* @return arbitrary precision numerical asinh(x). */
-numeric asinh(numeric const & x)
+numeric asinh(const numeric & x)
{
return ::asinh(*x.value); // -> CLN
}
/** Numeric inverse hyperbolic cosine (trigonometric function).
*
* @return arbitrary precision numerical acosh(x). */
-numeric acosh(numeric const & x)
+numeric acosh(const numeric & x)
{
return ::acosh(*x.value); // -> CLN
}
/** Numeric inverse hyperbolic tangent (trigonometric function).
*
* @return arbitrary precision numerical atanh(x). */
-numeric atanh(numeric const & x)
+numeric atanh(const numeric & x)
{
return ::atanh(*x.value); // -> CLN
}
/** Numeric evaluation of Riemann's Zeta function. Currently works only for
* integer arguments. */
-numeric zeta(numeric const & x)
+numeric zeta(const numeric & x)
{
// A dirty hack to allow for things like zeta(3.0), since CLN currently
// only knows about integer arguments and zeta(3).evalf() automatically
/** The gamma function.
* This is only a stub! */
-numeric gamma(numeric const & x)
+numeric gamma(const numeric & x)
{
clog << "gamma(" << x
<< "): Does anybody know good way to calculate this numerically?"
/** The psi function (aka polygamma function).
* This is only a stub! */
-numeric psi(numeric const & x)
+numeric psi(const numeric & x)
{
clog << "psi(" << x
<< "): Does anybody know good way to calculate this numerically?"
/** The psi functions (aka polygamma functions).
* This is only a stub! */
-numeric psi(numeric const & n, numeric const & x)
+numeric psi(const numeric & n, const numeric & x)
{
clog << "psi(" << n << "," << x
<< "): Does anybody know good way to calculate this numerically?"
/** Factorial combinatorial function.
*
* @exception range_error (argument must be integer >= 0) */
-numeric factorial(numeric const & nn)
+numeric factorial(const numeric & nn)
{
if (!nn.is_nonneg_integer())
throw (std::range_error("numeric::factorial(): argument must be integer >= 0"));
* useful in cases, like for exact results of Gamma(n+1/2) for instance.)
*
* @param n integer argument >= -1
- * @return n!! == n * (n-2) * (n-4) * ... * ({1|2}) with 0!! == 1 == (-1)!!
+ * @return n!! == n * (n-2) * (n-4) * ... * ({1|2}) with 0!! == (-1)!! == 1
* @exception range_error (argument must be integer >= -1) */
-numeric doublefactorial(numeric const & nn)
+numeric doublefactorial(const numeric & nn)
{
- // META-NOTE: The whole shit here will become obsolete and may be moved
- // out once CLN learns about double factorial, which should be as soon as
- // 1.0.3 rolls out!
-
- // We store the results separately for even and odd arguments. This has
- // the advantage that we don't have to compute any even result at all if
- // the function is always called with odd arguments and vice versa. There
- // is no tradeoff involved in this, it is guaranteed to save time as well
- // as memory. (If this is not enough justification consider the Gamma
- // function of half integer arguments: it only needs odd doublefactorials.)
- static vector<numeric> evenresults;
- static int highest_evenresult = -1;
- static vector<numeric> oddresults;
- static int highest_oddresult = -1;
-
if (nn == numeric(-1)) {
return _num1();
}
if (!nn.is_nonneg_integer()) {
throw (std::range_error("numeric::doublefactorial(): argument must be integer >= -1"));
}
- if (nn.is_even()) {
- int n = nn.div(_num2()).to_int();
- if (n <= highest_evenresult) {
- return evenresults[n];
- }
- if (evenresults.capacity() < (unsigned)(n+1)) {
- evenresults.reserve(n+1);
- }
- if (highest_evenresult < 0) {
- evenresults.push_back(_num1());
- highest_evenresult=0;
- }
- for (int i=highest_evenresult+1; i<=n; i++) {
- evenresults.push_back(numeric(evenresults[i-1].mul(numeric(i*2))));
- }
- highest_evenresult=n;
- return evenresults[n];
- } else {
- int n = nn.sub(_num1()).div(_num2()).to_int();
- if (n <= highest_oddresult) {
- return oddresults[n];
- }
- if (oddresults.capacity() < (unsigned)n) {
- oddresults.reserve(n+1);
- }
- if (highest_oddresult < 0) {
- oddresults.push_back(_num1());
- highest_oddresult=0;
- }
- for (int i=highest_oddresult+1; i<=n; i++) {
- oddresults.push_back(numeric(oddresults[i-1].mul(numeric(i*2+1))));
- }
- highest_oddresult=n;
- return oddresults[n];
- }
+ return numeric(::doublefactorial(nn.to_int())); // -> CLN
}
/** The Binomial coefficients. It computes the binomial coefficients. For
* integer n and k and positive n this is the number of ways of choosing k
* objects from n distinct objects. If n is negative, the formula
* binomial(n,k) == (-1)^k*binomial(k-n-1,k) is used to compute the result. */
-numeric binomial(numeric const & n, numeric const & k)
+numeric binomial(const numeric & n, const numeric & k)
{
if (n.is_integer() && k.is_integer()) {
if (n.is_nonneg_integer()) {
*
* @return the nth Bernoulli number (a rational number).
* @exception range_error (argument must be integer >= 0) */
-numeric bernoulli(numeric const & nn)
+numeric bernoulli(const numeric & nn)
{
if (!nn.is_integer() || nn.is_negative())
throw (std::range_error("numeric::bernoulli(): argument must be integer >= 0"));
}
/** Absolute value. */
-numeric abs(numeric const & x)
+numeric abs(const numeric & x)
{
return ::abs(*x.value); // -> CLN
}
*
* @return a mod b in the range [0,abs(b)-1] with sign of b if both are
* integer, 0 otherwise. */
-numeric mod(numeric const & a, numeric const & b)
+numeric mod(const numeric & a, const numeric & b)
{
if (a.is_integer() && b.is_integer())
return ::mod(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
* Equivalent to Maple's mods.
*
* @return a mod b in the range [-iquo(abs(m)-1,2), iquo(abs(m),2)]. */
-numeric smod(numeric const & a, numeric const & b)
+numeric smod(const numeric & a, const numeric & b)
{
// FIXME: Should this become a member function?
if (a.is_integer() && b.is_integer()) {
* sign of a or is zero.
*
* @return remainder of a/b if both are integer, 0 otherwise. */
-numeric irem(numeric const & a, numeric const & b)
+numeric irem(const numeric & a, const numeric & b)
{
if (a.is_integer() && b.is_integer())
return ::rem(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
*
* @return remainder of a/b and quotient stored in q if both are integer,
* 0 otherwise. */
-numeric irem(numeric const & a, numeric const & b, numeric & q)
+numeric irem(const numeric & a, const numeric & b, numeric & q)
{
if (a.is_integer() && b.is_integer()) { // -> CLN
cl_I_div_t rem_quo = truncate2(The(cl_I)(*a.value), The(cl_I)(*b.value));
* Equivalent to Maple's iquo as far as sign conventions are concerned.
*
* @return truncated quotient of a/b if both are integer, 0 otherwise. */
-numeric iquo(numeric const & a, numeric const & b)
+numeric iquo(const numeric & a, const numeric & b)
{
if (a.is_integer() && b.is_integer())
return truncate1(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
*
* @return truncated quotient of a/b and remainder stored in r if both are
* integer, 0 otherwise. */
-numeric iquo(numeric const & a, numeric const & b, numeric & r)
+numeric iquo(const numeric & a, const numeric & b, numeric & r)
{
if (a.is_integer() && b.is_integer()) { // -> CLN
cl_I_div_t rem_quo = truncate2(The(cl_I)(*a.value), The(cl_I)(*b.value));
* @return square root of z. Branch cut along negative real axis, the negative
* real axis itself where imag(z)==0 and real(z)<0 belongs to the upper part
* where imag(z)>0. */
-numeric sqrt(numeric const & z)
+numeric sqrt(const numeric & z)
{
return ::sqrt(*z.value); // -> CLN
}
/** Integer numeric square root. */
-numeric isqrt(numeric const & x)
+numeric isqrt(const numeric & x)
{
if (x.is_integer()) {
cl_I root;
*
* @return The GCD of two numbers if both are integer, a numerical 1
* if they are not. */
-numeric gcd(numeric const & a, numeric const & b)
+numeric gcd(const numeric & a, const numeric & b)
{
if (a.is_integer() && b.is_integer())
return ::gcd(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
*
* @return The LCM of two numbers if both are integer, the product of those
* two numbers if they are not. */
-numeric lcm(numeric const & a, numeric const & b)
+numeric lcm(const numeric & a, const numeric & b)
{
if (a.is_integer() && b.is_integer())
return ::lcm(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
{
assert(!too_late);
too_late = true;
- cl_default_float_format = cl_float_format(17);
+ cl_default_float_format = cl_float_format(17);
}
_numeric_digits& _numeric_digits::operator=(long prec)
git://www.ginac.de / ginac.git / blobdiff