- if (x.info(info_flags::numeric)) {
- // trap integer arguments:
- if (x.info(info_flags::integer)) {
- // gamma(n+1) -> n! for postitive n
- if (x.info(info_flags::posint)) {
- return factorial(ex_to_numeric(x).sub(_num1()));
- } else {
- throw (std::domain_error("gamma_eval(): simple pole"));
- }
- }
- // trap half integer arguments:
- if ((x*2).info(info_flags::integer)) {
- // trap positive x==(n+1/2)
- // gamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
- if ((x*_ex2()).info(info_flags::posint)) {
- numeric n = ex_to_numeric(x).sub(_num1_2());
- numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1()));
- coefficient = coefficient.div(pow(_num2(),n));
- return coefficient * pow(Pi,_ex1_2());
- } else {
- // trap negative x==(-n+1/2)
- // gamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
- numeric n = abs(ex_to_numeric(x).sub(_num1_2()));
- numeric coefficient = pow(_num_2(), n);
- coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));;
- return coefficient*power(Pi,_ex1_2());
- }
- }
- // gamma_evalf should be called here once it becomes available
- }
-
- return gamma(x).hold();
-}
+ if (x.info(info_flags::numeric)) {
+ // trap integer arguments:
+ if (x.info(info_flags::integer)) {
+ // lgamma(n) -> log((n-1)!) for postitive n
+ if (x.info(info_flags::posint))
+ return log(factorial(x + _ex_1));
+ else
+ throw (pole_error("lgamma_eval(): logarithmic pole",0));
+ }
+ // lgamma_evalf should be called here once it becomes available
+ }
+
+ return lgamma(x).hold();
+}
+
+
+static ex lgamma_deriv(const ex & x, unsigned deriv_param)
+{
+ GINAC_ASSERT(deriv_param==0);
+
+ // d/dx lgamma(x) -> psi(x)
+ return psi(x);
+}
+
+
+static ex lgamma_series(const ex & arg,
+ const relational & rel,
+ int order,
+ unsigned options)
+{
+ // method:
+ // Taylor series where there is no pole falls back to psi function
+ // evaluation.
+ // On a pole at -m we could use the recurrence relation
+ // lgamma(x) == lgamma(x+1)-log(x)
+ // from which follows
+ // series(lgamma(x),x==-m,order) ==
+ // series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order);
+ const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
+ if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
+ throw do_taylor(); // caught by function::series()
+ // if we got here we have to care for a simple pole of tgamma(-m):
+ numeric m = -ex_to<numeric>(arg_pt);
+ ex recur;
+ for (numeric p = 0; p<=m; ++p)
+ recur += log(arg+p);
+ return (lgamma(arg+m+_ex1)-recur).series(rel, order, options);
+}
+
+
+REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval).
+ evalf_func(lgamma_evalf).
+ derivative_func(lgamma_deriv).
+ series_func(lgamma_series).
+ latex_name("\\log \\Gamma"));
+
+
+//////////
+// true Gamma function
+//////////
+
+static ex tgamma_evalf(const ex & x)
+{
+ if (is_exactly_a<numeric>(x)) {
+ try {
+ return tgamma(ex_to<numeric>(x));
+ } catch (const dunno &e) { }
+ }
+
+ return tgamma(x).hold();
+}
+
+
+/** Evaluation of tgamma(x), the true Gamma function. Knows about integer
+ * arguments, half-integer arguments and that's it. Somebody ought to provide
+ * some good numerical evaluation some day...
+ *
+ * @exception pole_error("tgamma_eval(): simple pole",0) */
+static ex tgamma_eval(const ex & x)
+{
+ if (x.info(info_flags::numeric)) {
+ // trap integer arguments:
+ const numeric two_x = (*_num2_p)*ex_to<numeric>(x);
+ if (two_x.is_even()) {
+ // tgamma(n) -> (n-1)! for postitive n
+ if (two_x.is_positive()) {
+ return factorial(ex_to<numeric>(x).sub(*_num1_p));
+ } else {
+ throw (pole_error("tgamma_eval(): simple pole",1));
+ }
+ }
+ // trap half integer arguments:
+ if (two_x.is_integer()) {
+ // trap positive x==(n+1/2)
+ // tgamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
+ if (two_x.is_positive()) {
+ const numeric n = ex_to<numeric>(x).sub(*_num1_2_p);
+ return (doublefactorial(n.mul(*_num2_p).sub(*_num1_p)).div(pow(*_num2_p,n))) * sqrt(Pi);
+ } else {
+ // trap negative x==(-n+1/2)
+ // tgamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
+ const numeric n = abs(ex_to<numeric>(x).sub(*_num1_2_p));
+ return (pow(*_num_2_p, n).div(doublefactorial(n.mul(*_num2_p).sub(*_num1_p))))*sqrt(Pi);
+ }
+ }
+ // tgamma_evalf should be called here once it becomes available
+ }
+
+ return tgamma(x).hold();
+}