+// 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*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));
+ } 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);
+ return (doublefactorial(n.mul(_num2).sub(_num1)).div(pow(_num2,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));
+ return (pow(_num_2, n).div(doublefactorial(n.mul(_num2).sub(_num1))))*sqrt(Pi);
+ }
+ }
+ // tgamma_evalf should be called here once it becomes available
+ }
+
+ return tgamma(x).hold();
+}
+
+
+static ex tgamma_deriv(const ex & x, unsigned deriv_param)
+{
+ GINAC_ASSERT(deriv_param==0);
+
+ // d/dx tgamma(x) -> psi(x)*tgamma(x)
+ return psi(x)*tgamma(x);
+}
+
+
+static ex tgamma_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 use the recurrence relation
+ // tgamma(x) == tgamma(x+1) / x
+ // from which follows
+ // series(tgamma(x),x==-m,order) ==
+ // series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order+1);
+ 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 at -m:
+ const numeric m = -ex_to<numeric>(arg_pt);
+ ex ser_denom = _ex1;
+ for (numeric p; p<=m; ++p)
+ ser_denom *= arg+p;
+ return (tgamma(arg+m+_ex1)/ser_denom).series(rel, order+1, options);
+}
+
+
+REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval).
+ evalf_func(tgamma_evalf).
+ derivative_func(tgamma_deriv).
+ series_func(tgamma_series).
+ latex_name("\\Gamma"));
+
+
+//////////
+// beta-function