X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=4cfd7e4a6949da20f6470d23e6753c356e4f57e9;hp=795c76099172dba5d5284fd6fc08c4183428fc3d;hb=HEAD;hpb=f282e2b3a3873619006b6ffd3fa7fe65010c1299

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index 795c7609..900cac17 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-2011 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2024 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
@@ -24,15 +24,10 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
 #include "numeric.h"
 #include "ex.h"
 #include "operators.h"
 #include "archive.h"
-#include "tostring.h"
 #include "utils.h"
 
 #include <limits>
@@ -90,7 +85,7 @@ numeric::numeric()
 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>)
@@ -111,7 +106,7 @@ numeric::numeric(int i)
 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>)
@@ -142,6 +137,17 @@ numeric::numeric(unsigned long 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.
  *
@@ -225,7 +231,7 @@ numeric::numeric(const char *s)
 			// 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()));
@@ -443,7 +449,7 @@ static void print_real_csrc(const print_context & c, const cln::cl_R & x)
 		// 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 {
@@ -467,7 +473,7 @@ 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
+ * @sa https://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)
@@ -502,20 +508,20 @@ static void print_real_cl_N(const print_context & c, const cln::cl_R & x)
 {
 	if (cln::instanceof(x, cln::cl_I_ring)) {
 
-    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 << "\")";
-    }
+		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
@@ -715,8 +721,6 @@ bool numeric::info(unsigned inf) const
 			return is_odd();
 		case info_flags::prime:
 			return is_prime();
-		case info_flags::algebraic:
-			return !is_real();
 	}
 	return false;
 }
@@ -776,10 +780,8 @@ bool numeric::has(const ex &other, unsigned options) const
 
 
 /** 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();
 }
 
@@ -789,11 +791,9 @@ ex numeric::eval(int level) const
  *  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);
 }
 
@@ -930,9 +930,8 @@ const numeric &numeric::add_dyn(const numeric &other) const
 		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);
 }
 
 
@@ -946,9 +945,8 @@ const numeric &numeric::sub_dyn(const numeric &other) const
 	// 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);
 }
 
 
@@ -965,8 +963,7 @@ const numeric &numeric::mul_dyn(const numeric &other) const
 	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);
 }
 
 
@@ -984,8 +981,8 @@ const numeric &numeric::div_dyn(const numeric &other) const
 		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);
 }
 
 
@@ -1011,8 +1008,8 @@ const numeric &numeric::power_dyn(const numeric &other) const
 		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));
 }
 
 
@@ -1376,7 +1373,7 @@ const numeric numeric::numer() const
 		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
@@ -1531,7 +1528,7 @@ const numeric atan(const numeric &y, const numeric &x)
 		return *_num0_p;
 	if (x.is_real() && y.is_real())
 		return numeric(cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
-		                 cln::the<cln::cl_R>(y.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)))
@@ -1749,7 +1746,7 @@ class lanczos_coeffs
 		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) {
@@ -1775,8 +1772,8 @@ cln::cl_N lanczos_coeffs::calc_lanczos_A(const cln::cl_N &x) const
 {
 	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;
 }
 
@@ -2024,7 +2021,7 @@ lanczos_coeffs::lanczos_coeffs()
 	coeffs[3].swap(coeffs_120);
 }
 
-static const cln::float_format_t guess_precision(const cln::cl_N& x)
+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))
@@ -2050,11 +2047,11 @@ const cln::cl_N lgamma(const cln::cl_N &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;
+		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();
@@ -2077,10 +2074,10 @@ const cln::cl_N tgamma(const cln::cl_N &x)
 			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;
+		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();
@@ -2141,8 +2138,10 @@ const numeric doublefactorial(const numeric &n)
 
 /** 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. */
+ *  objects from n distinct objects.  If n is a negative integer, the formula
+ *  binomial(n,k) == (-1)^k*binomial(k-n-1,k) (if k>=0)
+ *  binomial(n,k) == (-1)^(n-k)*binomial(-k-1,n-k) (otherwise)
+ *  is used to compute the result. */
 const numeric binomial(const numeric &n, const numeric &k)
 {
 	if (n.is_integer() && k.is_integer()) {
@@ -2152,7 +2151,10 @@ const numeric binomial(const numeric &n, const numeric &k)
 			else
 				return *_num0_p;
 		} else {
-			return _num_1_p->power(k)*binomial(k-n-(*_num1_p),k);
+			if (k.is_nonneg_integer())
+				return _num_1_p->power(k)*binomial(k-n-(*_num1_p), k);
+			else
+				return _num_1_p->power(n-k)*binomial(-k-(*_num1_p), n-k);
 		}
 	}
 	
@@ -2227,7 +2229,7 @@ const numeric bernoulli(const numeric &nn)
 
 	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
@@ -2244,7 +2246,7 @@ const numeric bernoulli(const numeric &nn)
 				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;
@@ -2539,9 +2541,8 @@ _numeric_digits& _numeric_digits::operator=(long prec)
 	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;