define lgamma and tgamma taking cl_N as an argument.

[ginac.git] / ginac / numeric.cpp

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index c80ff3a7757e57685a63e7b583d78f959833871b..69a787b5d985b5e174ea61a6f38c7b89f2bb3a2c 100644 (file)
--- 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! */
@@ -2137,7 +2161,7 @@ const numeric bernoulli(const numeric &nn)
                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) {
@@ -2162,7 +2186,7 @@ const numeric bernoulli(const numeric &nn)
                results.push_back(-b/(p+1));
        }
        next_r = n+2;
-       return results[n/2-1];
+       return numeric(results[n/2-1]);
 }
 
 
@@ -2219,9 +2243,9 @@ const numeric fibonacci(const numeric &n)
        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)); 
 }
 
 
@@ -2339,7 +2363,7 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r)
                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 = numeric(rem_quo.remainder);
-               return rem_quo.quotient;
+               return numeric(rem_quo.quotient);
        } else {
                r = *_num0_p;
                return *_num0_p;