[GiNaC-devel] [PATCH 5/8] tgamma, lgamma: convert to take cl_N as an argument; provide wrappers for GiNaC::numeric

[Alexei Sheplyakov]

Implicit conversion from cl_N to numeric considered harmful, part 5.

---
 ginac/numeric.cpp |   58 +++++++++++++++++++++++++++++++++++++---------------
 1 files changed, 41 insertions(+), 17 deletions(-)

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index a06be42..704854c 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! */
-- 
1.5.4.2


-- 
All science is either physics or stamp collecting.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
URL: <http://www.ginac.de/pipermail/ginac-devel/attachments/20080319/24b27309/attachment.sig>