X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_gamma.cpp;h=46c8ea9e00e717652b188a0979f18224bc639fa0;hp=baed308d366fba1a571e24989deca9f4c6c818a5;hb=94ead6345f31939c9c9eff70ce88d74cf36a8004;hpb=22abfbe8c78e339188096a5bf749a7c2d4f0a368

diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index baed308d..46c8ea9e 100644
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -4,7 +4,7 @@
  *  some related stuff. */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -21,9 +21,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <vector>
-#include <stdexcept>
-
 #include "inifcns.h"
 #include "constant.h"
 #include "pseries.h"
@@ -35,6 +32,9 @@
 #include "symmetry.h"
 #include "utils.h"
 
+#include <stdexcept>
+#include <vector>
+
 namespace GiNaC {
 
 //////////
@@ -54,8 +54,7 @@ static ex lgamma_evalf(const ex & x)
 
 
 /** Evaluation of lgamma(x), the natural logarithm of the Gamma function.
- *  Knows about integer arguments and that's it.  Somebody ought to provide
- *  some good numerical evaluation some day...
+ *  Handles integer arguments as a special case.
  *
  *  @exception GiNaC::pole_error("lgamma_eval(): logarithmic pole",0) */
 static ex lgamma_eval(const ex & x)
@@ -69,7 +68,8 @@ static ex lgamma_eval(const ex & x)
 			else
 				throw (pole_error("lgamma_eval(): logarithmic pole",0));
 		}
-		//  lgamma_evalf should be called here once it becomes available
+		if (!ex_to<numeric>(x).is_rational())
+			return lgamma(ex_to<numeric>(x));
 	}
 	
 	return lgamma(x).hold();
@@ -110,10 +110,26 @@ static ex lgamma_series(const ex & arg,
 }
 
 
+static ex lgamma_conjugate(const ex & x)
+{
+	// conjugate(lgamma(x))==lgamma(conjugate(x)) unless on the branch cut
+	// which runs along the negative real axis.
+	if (x.info(info_flags::positive)) {
+		return lgamma(x);
+	}
+	if (is_exactly_a<numeric>(x) &&
+	    !x.imag_part().is_zero()) {
+		return lgamma(x.conjugate());
+	}
+	return conjugate_function(lgamma(x)).hold();
+}
+
+
 REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval).
                           evalf_func(lgamma_evalf).
                           derivative_func(lgamma_deriv).
                           series_func(lgamma_series).
+                          conjugate_func(lgamma_conjugate).
                           latex_name("\\log \\Gamma"));
 
 
@@ -165,7 +181,8 @@ static ex tgamma_eval(const ex & x)
 				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
+		if (!ex_to<numeric>(x).is_rational())
+			return tgamma(ex_to<numeric>(x));
 	}
 	
 	return tgamma(x).hold();
@@ -206,10 +223,18 @@ static ex tgamma_series(const ex & arg,
 }
 
 
+static ex tgamma_conjugate(const ex & x)
+{
+	// conjugate(tgamma(x))==tgamma(conjugate(x))
+	return tgamma(x.conjugate());
+}
+
+
 REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval).
                           evalf_func(tgamma_evalf).
                           derivative_func(tgamma_deriv).
                           series_func(tgamma_series).
+                          conjugate_func(tgamma_conjugate).
                           latex_name("\\Gamma"));
 
 
@@ -221,7 +246,7 @@ static ex beta_evalf(const ex & x, const ex & y)
 {
 	if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y)) {
 		try {
-			return tgamma(ex_to<numeric>(x))*tgamma(ex_to<numeric>(y))/tgamma(ex_to<numeric>(x+y));
+			return exp(lgamma(ex_to<numeric>(x))+lgamma(ex_to<numeric>(y))-lgamma(ex_to<numeric>(x+y)));
 		} catch (const dunno &e) { }
 	}
 	
@@ -262,7 +287,8 @@ static ex beta_eval(const ex & x, const ex & y)
 		    (nx+ny).is_integer() &&
 		   !(nx+ny).is_positive())
 			 return _ex0;
-		// beta_evalf should be called here once it becomes available
+		if (!ex_to<numeric>(x).is_rational() || !ex_to<numeric>(x).is_rational())
+			return evalf(beta(x, y).hold());
 	}
 	
 	return beta(x,y).hold();
@@ -326,8 +352,8 @@ REGISTER_FUNCTION(beta, eval_func(beta_eval).
                         evalf_func(beta_evalf).
                         derivative_func(beta_deriv).
                         series_func(beta_series).
-                        latex_name("\\mbox{B}").
-						set_symmetry(sy_symm(0, 1)));
+                        latex_name("\\mathrm{B}").
+                        set_symmetry(sy_symm(0, 1)));
 
 
 //////////