* 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-2008 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
}
}
+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();
* @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()));
}
-/** Arcustangent.
+/** Numeric arcustangent.
*
* @param x complex number
* @return atan(x)
- * @exception pole_error("atan(): logarithmic pole",0) */
+ * @exception pole_error("atan(): logarithmic pole",0) if x==I or x==-I. */
const numeric atan(const numeric &x)
{
if (!x.is_real() &&
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()));
}
-/** Arcustangent.
+/** Numeric arcustangent of two arguments, analytically continued in a suitable way.
*
- * @param x real number
- * @param y real number
- * @return atan(y/x) */
+ * @param y complex number
+ * @param x complex number
+ * @return -I*log((x+I*y)/sqrt(x^2+y^2)), which is equal to atan(y/x) if y and
+ * x are both real.
+ * @exception pole_error("atan(): logarithmic pole",0) if y/x==+I or y/x==-I. */
const numeric atan(const numeric &y, const numeric &x)
{
+ 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()));
- else
- throw std::invalid_argument("atan(): complex argument");
+ 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)))
+ // Do not "simplify" this to -I/2*log((x+I*y)/(x-I*y))) or likewise.
+ // The branch cuts are easily messed up.
+ const cln::cl_N aux_p = x.to_cl_N()+cln::complex(0,1)*y.to_cl_N();
+ if (cln::zerop(aux_p)) {
+ // x+I*y==0 => y/x==I, so this is a pole (we have x!=0).
+ throw pole_error("atan(): logarithmic pole",0);
+ }
+ const cln::cl_N aux_m = x.to_cl_N()-cln::complex(0,1)*y.to_cl_N();
+ if (cln::zerop(aux_m)) {
+ // x-I*y==0 => y/x==-I, so this is a pole (we have x!=0).
+ throw pole_error("atan(): logarithmic pole",0);
+ }
+ 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();
}
coeffs[3].swap(coeffs_120);
}
+static const 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;
+ 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.to_cl_N()-cln::cl_N(1)/2)*log(temp)
+ + (x-cln::cl_N(1)/2)*log(temp)
- temp
+ log(A);
return result;
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;
+ 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.to_cl_N()-cln::cl_N(1)/2)
+ = sqrt(cln::cl_I(2)*pi_val) * expt(temp, x - cln::cl_N(1)/2)
* exp(-temp) * A;
return result;
}
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! */
/** 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;
}
{
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 numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()) + b2,
+ cln::the<cln::cl_I>(b.to_cl_N())) - b2);
} 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;
+ r = numeric(rem_quo.remainder);
return rem_quo.quotient;
} else {
r = *_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;
}
git://www.ginac.de / ginac.git / blobdiff