X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnumeric.cpp;h=40aa8e2ee90ec1b18e10a0e75a6d0084de81c81b;hp=c80ff3a7757e57685a63e7b583d78f959833871b;hb=5ea14a22768031a1cd4abce2926d5359c5d0c15f;hpb=7bc327c75aaa3118de14dfcee59bcf0fd40e3f4a

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index c80ff3a7..40aa8e2e 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -1962,24 +1962,34 @@ lanczos_coeffs::lanczos_coeffs()
 	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;
@@ -1988,17 +1998,25 @@ const numeric lgamma(const numeric &x)
 		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;
 	}
@@ -2006,6 +2024,12 @@ const numeric tgamma(const numeric &x)
 		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! */