* of special functions or implement the interface to the bignum package. */
/*
- * GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2019 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
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
+#ifdef HAVE_CONFIG_H
#include "config.h"
-
-#include <vector>
-#include <stdexcept>
-#include <string>
-#include <sstream>
-#include <limits>
+#endif
#include "numeric.h"
#include "ex.h"
#include "operators.h"
#include "archive.h"
-#include "tostring.h"
#include "utils.h"
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
// CLN should pollute the global namespace as little as possible. Hence, we
// include most of it here and include only the part needed for properly
// declaring cln::cl_number in numeric.h. This can only be safely done in
//////////
/** default ctor. Numerically it initializes to an integer zero. */
-numeric::numeric() : basic(&numeric::tinfo_static)
+numeric::numeric()
{
value = cln::cl_I(0);
setflag(status_flags::evaluated | status_flags::expanded);
// public
-numeric::numeric(int i) : basic(&numeric::tinfo_static)
+numeric::numeric(int i)
{
// 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
+ // first. This is due to the behavior of the cl_I-ctor, which
// emphasizes efficiency. However, if the integer is small enough
// we save space and dereferences by using an immediate type.
// (C.f. <cln/object.h>)
}
-numeric::numeric(unsigned int i) : basic(&numeric::tinfo_static)
+numeric::numeric(unsigned int 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
+ // first. This is due to the behavior of the cl_I-ctor, which
// emphasizes efficiency. However, if the integer is small enough
// we save space and dereferences by using an immediate type.
// (C.f. <cln/object.h>)
}
-numeric::numeric(long i) : basic(&numeric::tinfo_static)
+numeric::numeric(long i)
{
value = cln::cl_I(i);
setflag(status_flags::evaluated | status_flags::expanded);
}
-numeric::numeric(unsigned long i) : basic(&numeric::tinfo_static)
+numeric::numeric(unsigned long i)
{
value = cln::cl_I(i);
setflag(status_flags::evaluated | status_flags::expanded);
}
+numeric::numeric(long long i)
+{
+ value = cln::cl_I(i);
+ setflag(status_flags::evaluated | status_flags::expanded);
+}
+
+numeric::numeric(unsigned long long i)
+{
+ value = cln::cl_I(i);
+ setflag(status_flags::evaluated | status_flags::expanded);
+}
/** Constructor for rational numerics a/b.
*
* @exception overflow_error (division by zero) */
-numeric::numeric(long numer, long denom) : basic(&numeric::tinfo_static)
+numeric::numeric(long numer, long denom)
{
if (!denom)
throw std::overflow_error("division by zero");
}
-numeric::numeric(double d) : basic(&numeric::tinfo_static)
+numeric::numeric(double d)
{
// 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
/** ctor from C-style string. It also accepts complex numbers in GiNaC
* notation like "2+5*I". */
-numeric::numeric(const char *s) : basic(&numeric::tinfo_static)
+numeric::numeric(const char *s)
{
cln::cl_N ctorval = 0;
// parse complex numbers (functional but not completely safe, unfortunately
// E to lower case
term = term.replace(term.find("E"),1,"e");
// append _<Digits> to term
- term += "_" + ToString((unsigned)Digits);
+ term += "_" + std::to_string((unsigned)Digits);
// construct float using cln::cl_F(const char *) ctor.
if (imaginary)
ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_F(term.c_str()));
/** Ctor from CLN types. This is for the initiated user or internal use
* only. */
-numeric::numeric(const cln::cl_N &z) : basic(&numeric::tinfo_static)
+numeric::numeric(const cln::cl_N &z)
{
value = z;
setflag(status_flags::evaluated | status_flags::expanded);
// archiving
//////////
-numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+/**
+ * Construct a floating point number from sign, mantissa, and exponent
+ */
+static const cln::cl_F make_real_float(const cln::cl_idecoded_float& dec)
{
- cln::cl_N ctorval = 0;
+ cln::cl_F x = cln::cl_float(dec.mantissa, cln::default_float_format);
+ x = cln::scale_float(x, dec.exponent);
+ cln::cl_F sign = cln::cl_float(dec.sign, cln::default_float_format);
+ x = cln::float_sign(sign, x);
+ return x;
+}
+/**
+ * Read serialized floating point number
+ */
+static const cln::cl_F read_real_float(std::istream& s)
+{
+ cln::cl_idecoded_float dec;
+ s >> dec.sign >> dec.mantissa >> dec.exponent;
+ const cln::cl_F x = make_real_float(dec);
+ return x;
+}
+
+void numeric::read_archive(const archive_node &n, lst &sym_lst)
+{
+ inherited::read_archive(n, sym_lst);
+ value = 0;
+
// Read number as string
std::string str;
if (n.find_string("number", str)) {
std::istringstream s(str);
- cln::cl_idecoded_float re, im;
+ cln::cl_R re, im;
char c;
s.get(c);
switch (c) {
- case 'R': // Integer-decoded real number
- s >> re.sign >> re.mantissa >> re.exponent;
- ctorval = re.sign * re.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), re.exponent);
+ case 'R':
+ // real FP (floating point) number
+ re = read_real_float(s);
+ value = re;
break;
- case 'C': // Integer-decoded complex number
- s >> re.sign >> re.mantissa >> re.exponent;
- s >> im.sign >> im.mantissa >> im.exponent;
- ctorval = cln::complex(re.sign * re.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), re.exponent),
- im.sign * im.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), im.exponent));
+ case 'C':
+ // both real and imaginary part are FP numbers
+ re = read_real_float(s);
+ im = read_real_float(s);
+ value = cln::complex(re, im);
break;
- default: // Ordinary number
+ case 'H':
+ // real part is a rational number,
+ // imaginary part is a FP number
+ s >> re;
+ im = read_real_float(s);
+ value = cln::complex(re, im);
+ break;
+ case 'J':
+ // real part is a FP number,
+ // imaginary part is a rational number
+ re = read_real_float(s);
+ s >> im;
+ value = cln::complex(re, im);
+ break;
+ default:
+ // both real and imaginary parts are rational
s.putback(c);
- s >> ctorval;
+ s >> value;
break;
}
}
- value = ctorval;
setflag(status_flags::evaluated | status_flags::expanded);
}
+GINAC_BIND_UNARCHIVER(numeric);
+
+static void write_real_float(std::ostream& s, const cln::cl_R& n)
+{
+ const cln::cl_idecoded_float dec = cln::integer_decode_float(cln::the<cln::cl_F>(n));
+ s << dec.sign << ' ' << dec.mantissa << ' ' << dec.exponent;
+}
void numeric::archive(archive_node &n) const
{
inherited::archive(n);
// Write number as string
+
+ const cln::cl_R re = cln::realpart(value);
+ const cln::cl_R im = cln::imagpart(value);
+ const bool re_rationalp = cln::instanceof(re, cln::cl_RA_ring);
+ const bool im_rationalp = cln::instanceof(im, cln::cl_RA_ring);
+
+ // Non-rational numbers are written in an integer-decoded format
+ // to preserve the precision
std::ostringstream s;
- if (this->is_crational())
+ if (re_rationalp && im_rationalp)
s << value;
- else {
- // Non-rational numbers are written in an integer-decoded format
- // to preserve the precision
- if (this->is_real()) {
- cln::cl_idecoded_float re = cln::integer_decode_float(cln::the<cln::cl_F>(value));
- s << "R";
- s << re.sign << " " << re.mantissa << " " << re.exponent;
- } else {
- cln::cl_idecoded_float re = cln::integer_decode_float(cln::the<cln::cl_F>(cln::realpart(cln::the<cln::cl_N>(value))));
- cln::cl_idecoded_float im = cln::integer_decode_float(cln::the<cln::cl_F>(cln::imagpart(cln::the<cln::cl_N>(value))));
- s << "C";
- s << re.sign << " " << re.mantissa << " " << re.exponent << " ";
- s << im.sign << " " << im.mantissa << " " << im.exponent;
- }
+ else if (zerop(im)) {
+ // real FP (floating point) number
+ s << 'R';
+ write_real_float(s, re);
+ } else if (re_rationalp) {
+ s << 'H'; // just any unique character
+ // real part is a rational number,
+ // imaginary part is a FP number
+ s << re << ' ';
+ write_real_float(s, im);
+ } else if (im_rationalp) {
+ s << 'J';
+ // real part is a FP number,
+ // imaginary part is a rational number
+ write_real_float(s, re);
+ s << ' ' << im;
+ } else {
+ // both real and imaginary parts are floating point
+ s << 'C';
+ write_real_float(s, re);
+ s << ' ';
+ write_real_float(s, im);
}
n.add_string("number", s.str());
}
-DEFAULT_UNARCHIVE(numeric)
-
//////////
// functions overriding virtual functions from base classes
//////////
// Rational number
const cln::cl_I numer = cln::numerator(cln::the<cln::cl_RA>(x));
const cln::cl_I denom = cln::denominator(cln::the<cln::cl_RA>(x));
- if (cln::plusp(x) > 0) {
+ if (cln::plusp(x)) {
c.s << "(";
print_integer_csrc(c, numer);
} else {
}
}
+template<typename T1, typename T2>
+static inline bool coerce(T1& dst, const T2& arg);
+
+/**
+ * @brief Check if CLN integer can be converted into int
+ *
+ * @sa http://www.ginac.de/pipermail/cln-list/2006-October/000248.html
+ */
+template<>
+inline bool coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg)
+{
+ static const cln::cl_I cl_max_int =
+ (cln::cl_I)(long)(std::numeric_limits<int>::max());
+ static const cln::cl_I cl_min_int =
+ (cln::cl_I)(long)(std::numeric_limits<int>::min());
+ if ((arg >= cl_min_int) && (arg <= cl_max_int)) {
+ dst = cl_I_to_int(arg);
+ return true;
+ }
+ return false;
+}
+
+template<>
+inline bool coerce<unsigned int, cln::cl_I>(unsigned int& dst, const cln::cl_I& arg)
+{
+ static const cln::cl_I cl_max_uint =
+ (cln::cl_I)(unsigned long)(std::numeric_limits<unsigned int>::max());
+ if ((! minusp(arg)) && (arg <= cl_max_uint)) {
+ dst = cl_I_to_uint(arg);
+ return true;
+ }
+ return false;
+}
+
/** Helper function to print real number in C++ source format using cl_N types.
*
* @see numeric::print() */
{
if (cln::instanceof(x, cln::cl_I_ring)) {
- // Integer number
- c.s << "cln::cl_I(\"";
- print_real_number(c, x);
- c.s << "\")";
-
+ int dst;
+ // fixnum
+ if (coerce(dst, cln::the<cln::cl_I>(x))) {
+ // can be converted to native int
+ if (dst < 0)
+ c.s << '(' << dst << ')';
+ else
+ c.s << dst;
+ } else {
+ // bignum
+ c.s << "cln::cl_I(\"";
+ print_real_number(c, x);
+ c.s << "\")";
+ }
} else if (cln::instanceof(x, cln::cl_RA_ring)) {
// Rational number
case info_flags::numeric:
case info_flags::polynomial:
case info_flags::rational_function:
+ case info_flags::expanded:
return true;
case info_flags::real:
return is_real();
case info_flags::negative:
return is_negative();
case info_flags::nonnegative:
- return !is_negative();
+ return is_zero() || is_positive();
case info_flags::posint:
return is_pos_integer();
case info_flags::negint:
return is_odd();
case info_flags::prime:
return is_prime();
- case info_flags::algebraic:
- return !is_real();
}
return false;
}
/** Evaluation of numbers doesn't do anything at all. */
-ex numeric::eval(int level) const
+ex numeric::eval() const
{
- // Warning: if this is ever gonna do something, the ex ctors from all kinds
- // of numbers should be checking for status_flags::evaluated.
return this->hold();
}
* currently set. In case the object already was a floating point number the
* precision is trimmed to match the currently set default.
*
- * @param level ignored, only needed for overriding basic::evalf.
* @return an ex-handle to a numeric. */
-ex numeric::evalf(int level) const
+ex numeric::evalf() const
{
- // level can safely be discarded for numeric objects.
return numeric(cln::cl_float(1.0, cln::default_float_format) * value);
}
return other;
else if (&other==_num0_p)
return *this;
-
- return static_cast<const numeric &>((new numeric(value + other.value))->
- setflag(status_flags::dynallocated));
+
+ return dynallocate<numeric>(value + other.value);
}
// hack is supposed to keep the number of distinct numeric objects low.
if (&other==_num0_p || cln::zerop(other.value))
return *this;
-
- return static_cast<const numeric &>((new numeric(value - other.value))->
- setflag(status_flags::dynallocated));
+
+ return dynallocate<numeric>(value - other.value);
}
else if (&other==_num1_p)
return *this;
- return static_cast<const numeric &>((new numeric(value * other.value))->
- setflag(status_flags::dynallocated));
+ return dynallocate<numeric>(value * other.value);
}
return *this;
if (cln::zerop(cln::the<cln::cl_N>(other.value)))
throw std::overflow_error("division by zero");
- return static_cast<const numeric &>((new numeric(value / other.value))->
- setflag(status_flags::dynallocated));
+
+ return dynallocate<numeric>(value / other.value);
}
else
return *_num0_p;
}
- return static_cast<const numeric &>((new numeric(cln::expt(value, other.value)))->
- setflag(status_flags::dynallocated));
+
+ return dynallocate<numeric>(cln::expt(value, other.value));
}
if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_RA_ring)) {
const cln::cl_I s = cln::lcm(cln::denominator(r), cln::denominator(i));
return numeric(cln::complex(cln::numerator(r)*(cln::exquo(s,cln::denominator(r))),
- cln::numerator(i)*(cln::exquo(s,cln::denominator(i)))));
+ cln::numerator(i)*(cln::exquo(s,cln::denominator(i)))));
}
}
// at least one float encountered
* @return arbitrary precision numerical exp(x). */
const numeric exp(const numeric &x)
{
- return cln::exp(x.to_cl_N());
+ return numeric(cln::exp(x.to_cl_N()));
}
{
if (x.is_zero())
throw pole_error("log(): logarithmic pole",0);
- return cln::log(x.to_cl_N());
+ return numeric(cln::log(x.to_cl_N()));
}
* @return arbitrary precision numerical sin(x). */
const numeric sin(const numeric &x)
{
- return cln::sin(x.to_cl_N());
+ return numeric(cln::sin(x.to_cl_N()));
}
* @return arbitrary precision numerical cos(x). */
const numeric cos(const numeric &x)
{
- return cln::cos(x.to_cl_N());
+ return numeric(cln::cos(x.to_cl_N()));
}
* @return arbitrary precision numerical tan(x). */
const numeric tan(const numeric &x)
{
- return cln::tan(x.to_cl_N());
+ return numeric(cln::tan(x.to_cl_N()));
}
* @return arbitrary precision numerical asin(x). */
const numeric asin(const numeric &x)
{
- return cln::asin(x.to_cl_N());
+ return numeric(cln::asin(x.to_cl_N()));
}
* @return arbitrary precision numerical acos(x). */
const numeric acos(const numeric &x)
{
- return cln::acos(x.to_cl_N());
+ return numeric(cln::acos(x.to_cl_N()));
}
x.real().is_zero() &&
abs(x.imag()).is_equal(*_num1_p))
throw pole_error("atan(): logarithmic pole",0);
- return cln::atan(x.to_cl_N());
+ return numeric(cln::atan(x.to_cl_N()));
}
if (x.is_zero() && y.is_zero())
return *_num0_p;
if (x.is_real() && y.is_real())
- return cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
- cln::the<cln::cl_R>(y.to_cl_N()));
+ return numeric(cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
+ cln::the<cln::cl_R>(y.to_cl_N())));
// Compute -I*log((x+I*y)/sqrt(x^2+y^2))
// == -I*log((x+I*y)/sqrt((x+I*y)*(x-I*y)))
// x-I*y==0 => y/x==-I, so this is a pole (we have x!=0).
throw pole_error("atan(): logarithmic pole",0);
}
- return cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m));
+ return numeric(cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m)));
}
* @return arbitrary precision numerical sinh(x). */
const numeric sinh(const numeric &x)
{
- return cln::sinh(x.to_cl_N());
+ return numeric(cln::sinh(x.to_cl_N()));
}
* @return arbitrary precision numerical cosh(x). */
const numeric cosh(const numeric &x)
{
- return cln::cosh(x.to_cl_N());
+ return numeric(cln::cosh(x.to_cl_N()));
}
* @return arbitrary precision numerical tanh(x). */
const numeric tanh(const numeric &x)
{
- return cln::tanh(x.to_cl_N());
+ return numeric(cln::tanh(x.to_cl_N()));
}
* @return arbitrary precision numerical asinh(x). */
const numeric asinh(const numeric &x)
{
- return cln::asinh(x.to_cl_N());
+ return numeric(cln::asinh(x.to_cl_N()));
}
* @return arbitrary precision numerical acosh(x). */
const numeric acosh(const numeric &x)
{
- return cln::acosh(x.to_cl_N());
+ return numeric(cln::acosh(x.to_cl_N()));
}
* @return arbitrary precision numerical atanh(x). */
const numeric atanh(const numeric &x)
{
- return cln::atanh(x.to_cl_N());
+ return numeric(cln::atanh(x.to_cl_N()));
}
return Li2_series(x, prec);
}
+
/** Numeric evaluation of Dilogarithm. The domain is the entire complex plane,
* the branch cut lies along the positive real axis, starting at 1 and
* continuous with quadrant IV.
*
* @return arbitrary precision numerical Li2(x). */
-const numeric Li2(const numeric &x)
+const cln::cl_N Li2_(const cln::cl_N& value)
{
- if (x.is_zero())
- return *_num0_p;
+ if (zerop(value))
+ return 0;
// what is the desired float format?
// first guess: default format
cln::float_format_t prec = cln::default_float_format;
- const cln::cl_N value = x.to_cl_N();
// second guess: the argument's format
- if (!x.real().is_rational())
+ if (!instanceof(realpart(value), cln::cl_RA_ring))
prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
- else if (!x.imag().is_rational())
+ else if (!instanceof(imagpart(value), cln::cl_RA_ring))
prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
if (value==1) // may cause trouble with log(1-x)
- cln::zeta(2, prec)
- Li2_projection(cln::recip(value), prec));
else
- return Li2_projection(x.to_cl_N(), prec);
+ return Li2_projection(value, prec);
+}
+
+const numeric Li2(const numeric &x)
+{
+ const cln::cl_N x_ = x.to_cl_N();
+ if (zerop(x_))
+ return *_num0_p;
+ const cln::cl_N result = Li2_(x_);
+ return numeric(result);
}
if (x.is_real()) {
const int aux = (int)(cln::double_approx(cln::the<cln::cl_R>(x.to_cl_N())));
if (cln::zerop(x.to_cl_N()-aux))
- return cln::zeta(aux);
+ return numeric(cln::zeta(aux));
}
throw dunno();
}
std::vector<cln::cl_N> *current_vector;
};
-std::vector<cln::cl_N>* lanczos_coeffs::coeffs = 0;
+std::vector<cln::cl_N>* lanczos_coeffs::coeffs = nullptr;
bool lanczos_coeffs::sufficiently_accurate(int digits)
{ if (digits<=20) {
{
cln::cl_N A = (*current_vector)[0];
int size = current_vector->size();
- for (int i=1; i<size; ++i)
- A = A + (*current_vector)[i]/(x+cln::cl_I(-1+i));
+ for (int i=1; i<size; ++i)
+ A = A + (*current_vector)[i]/(x+cln::cl_I(-1+i));
return A;
}
coeffs[3].swap(coeffs_120);
}
+static cln::float_format_t guess_precision(const cln::cl_N& x)
+{
+ cln::float_format_t prec = cln::default_float_format;
+ if (!instanceof(realpart(x), cln::cl_RA_ring))
+ prec = cln::float_format(cln::the<cln::cl_F>(realpart(x)));
+ if (!instanceof(imagpart(x), cln::cl_RA_ring))
+ prec = cln::float_format(cln::the<cln::cl_F>(imagpart(x)));
+ return prec;
+}
/** The Gamma function.
* Use the Lanczos approximation. If the coefficients used here are not
* sufficiently many or sufficiently accurate, more can be calculated
* using the program doc/examples/lanczos.cpp. In that case, be sure to
* read the comments in that file. */
-const numeric lgamma(const numeric &x)
+const cln::cl_N lgamma(const cln::cl_N &x)
{
+ cln::float_format_t prec = guess_precision(x);
lanczos_coeffs lc;
- if (lc.sufficiently_accurate(Digits)) {
- cln::cl_N pi_val = cln::pi(cln::default_float_format);
- if (x.real() < 0.5)
- return log(pi_val) - log(sin(pi_val*x.to_cl_N()))
- - lgamma(numeric(1).sub(x));
- cln::cl_N A = lc.calc_lanczos_A(x.to_cl_N());
- cln::cl_N temp = x.to_cl_N() + lc.get_order() - cln::cl_N(1)/2;
- cln::cl_N result = log(cln::cl_I(2)*pi_val)/2
- + (x.to_cl_N()-cln::cl_N(1)/2)*log(temp)
- - temp
- + log(A);
- return result;
+ if (lc.sufficiently_accurate(prec)) {
+ cln::cl_N pi_val = cln::pi(prec);
+ if (realpart(x) < 0.5)
+ return cln::log(pi_val) - cln::log(sin(pi_val*x))
+ - lgamma(1 - x);
+ cln::cl_N A = lc.calc_lanczos_A(x);
+ cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2;
+ cln::cl_N result = log(cln::cl_I(2)*pi_val)/2
+ + (x-cln::cl_N(1)/2)*log(temp)
+ - temp
+ + log(A);
+ return result;
}
else
throw dunno();
}
-const numeric tgamma(const numeric &x)
+const numeric lgamma(const numeric &x)
+{
+ const cln::cl_N x_ = x.to_cl_N();
+ const cln::cl_N result = lgamma(x_);
+ return numeric(result);
+}
+
+const cln::cl_N tgamma(const cln::cl_N &x)
{
+ cln::float_format_t prec = guess_precision(x);
lanczos_coeffs lc;
- if (lc.sufficiently_accurate(Digits)) {
- cln::cl_N pi_val = cln::pi(cln::default_float_format);
- if (x.real() < 0.5)
- return pi_val/(sin(pi_val*x))/(tgamma(numeric(1).sub(x)).to_cl_N());
- cln::cl_N A = lc.calc_lanczos_A(x.to_cl_N());
- cln::cl_N temp = x.to_cl_N() + lc.get_order() - cln::cl_N(1)/2;
- cln::cl_N result
- = sqrt(cln::cl_I(2)*pi_val) * expt(temp, x.to_cl_N()-cln::cl_N(1)/2)
- * exp(-temp) * A;
- return result;
+ if (lc.sufficiently_accurate(prec)) {
+ cln::cl_N pi_val = cln::pi(prec);
+ if (realpart(x) < 0.5)
+ return pi_val/(cln::sin(pi_val*x))/tgamma(1 - x);
+ cln::cl_N A = lc.calc_lanczos_A(x);
+ cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2;
+ cln::cl_N result = sqrt(cln::cl_I(2)*pi_val)
+ * expt(temp, x - cln::cl_N(1)/2)
+ * exp(-temp) * A;
+ return result;
}
else
throw dunno();
}
+const numeric tgamma(const numeric &x)
+{
+ const cln::cl_N x_ = x.to_cl_N();
+ const cln::cl_N result = tgamma(x_);
+ return numeric(result);
+}
/** The psi function (aka polygamma function).
* This is only a stub! */
next_r = 4;
}
if (n<next_r)
- return results[n/2-1];
+ return numeric(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_I c = 1; // seed for binomial coefficients
cln::cl_RA b = cln::cl_RA(p-1)/-2;
// The CLN manual says: "The conversion from `unsigned int' works only
// if the argument is < 2^29" (This is for 32 Bit machines. More
c = cln::exquo((c * (p+3-2*k)) * (p/2-k+1), cln::cl_I(2*k-1)*k);
b = b + c*results[k-1];
}
- }
+ }
results.push_back(-b/(p+1));
}
next_r = n+2;
- return results[n/2-1];
+ return numeric(results[n/2-1]);
}
// F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
if (n.is_zero())
return *_num0_p;
- if (n.is_negative())
- if (n.is_even())
+ if (n.is_negative()) {
+ if (n.is_even()) {
return -fibonacci(-n);
- else
+ }
+ else {
return fibonacci(-n);
+ }
+ }
cln::cl_I u(0);
cln::cl_I v(1);
if (n.is_even())
// Here we don't use the squaring formula because one multiplication
// is cheaper than two squarings.
- return u * ((v << 1) - u);
+ return numeric(u * ((v << 1) - u));
else
- return cln::square(u) + cln::square(v);
+ return numeric(cln::square(u) + cln::square(v));
}
/** Absolute value. */
const numeric abs(const numeric& x)
{
- return cln::abs(x.to_cl_N());
+ return numeric(cln::abs(x.to_cl_N()));
}
const numeric mod(const numeric &a, const numeric &b)
{
if (a.is_integer() && b.is_integer())
- return cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
- cln::the<cln::cl_I>(b.to_cl_N()));
+ return numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
+ cln::the<cln::cl_I>(b.to_cl_N())));
else
return *_num0_p;
}
/** Modulus (in symmetric representation).
- * Equivalent to Maple's mods.
*
- * @return a mod b in the range [-iquo(abs(b)-1,2), iquo(abs(b),2)]. */
-const numeric smod(const numeric &a, const numeric &b)
-{
- if (a.is_integer() && b.is_integer()) {
- const cln::cl_I b2 = cln::ceiling1(cln::the<cln::cl_I>(b.to_cl_N()) >> 1) - 1;
- return cln::mod(cln::the<cln::cl_I>(a.to_cl_N()) + b2,
- cln::the<cln::cl_I>(b.to_cl_N())) - b2;
+ * @return a mod b in the range [-iquo(abs(b),2), iquo(abs(b),2)]. */
+const numeric smod(const numeric &a_, const numeric &b_)
+{
+ if (a_.is_integer() && b_.is_integer()) {
+ const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
+ const cln::cl_I b = cln::the<cln::cl_I>(b_.to_cl_N());
+ const cln::cl_I b2 = b >> 1;
+ const cln::cl_I m = cln::mod(a, b);
+ const cln::cl_I m_b = m - b;
+ const cln::cl_I ret = m > b2 ? m_b : m;
+ return numeric(ret);
} else
return *_num0_p;
}
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()));
+ return numeric(cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
+ cln::the<cln::cl_I>(b.to_cl_N())));
else
return *_num0_p;
}
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()));
- q = rem_quo.quotient;
- return rem_quo.remainder;
+ q = numeric(rem_quo.quotient);
+ return numeric(rem_quo.remainder);
} else {
q = *_num0_p;
return *_num0_p;
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()));
+ return numeric(cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
+ cln::the<cln::cl_I>(b.to_cl_N())));
else
return *_num0_p;
}
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()));
- r = rem_quo.remainder;
- return rem_quo.quotient;
+ r = numeric(rem_quo.remainder);
+ return numeric(rem_quo.quotient);
} else {
r = *_num0_p;
return *_num0_p;
const numeric gcd(const numeric &a, const numeric &b)
{
if (a.is_integer() && b.is_integer())
- return cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
- cln::the<cln::cl_I>(b.to_cl_N()));
+ return numeric(cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
+ cln::the<cln::cl_I>(b.to_cl_N())));
else
return *_num1_p;
}
const numeric lcm(const numeric &a, const numeric &b)
{
if (a.is_integer() && b.is_integer())
- return cln::lcm(cln::the<cln::cl_I>(a.to_cl_N()),
- cln::the<cln::cl_I>(b.to_cl_N()));
+ return numeric(cln::lcm(cln::the<cln::cl_I>(a.to_cl_N()),
+ cln::the<cln::cl_I>(b.to_cl_N())));
else
return a.mul(b);
}
* where imag(x)>0. */
const numeric sqrt(const numeric &x)
{
- return cln::sqrt(x.to_cl_N());
+ return numeric(cln::sqrt(x.to_cl_N()));
}
if (x.is_integer()) {
cln::cl_I root;
cln::isqrt(cln::the<cln::cl_I>(x.to_cl_N()), &root);
- return root;
+ return numeric(root);
} else
return *_num0_p;
}
cln::default_float_format = cln::float_format(prec);
// call registered callbacks
- std::vector<digits_changed_callback>::const_iterator it = callbacklist.begin(), end = callbacklist.end();
- for (; it != end; ++it) {
- (*it)(digitsdiff);
+ for (auto it : callbacklist) {
+ (it)(digitsdiff);
}
return *this;
git://www.ginac.de / ginac.git / blobdiff