From 8d16780dbd9f4da9397e638aca213745589818c0 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Fri, 19 Nov 1999 01:07:01 +0000 Subject: [PATCH] - fixed differentiation of gamma(x) - preliminary support for polygamma-functions psi(n,x) --- ginac/inifcns.h | 3 +++ ginac/inifcns_gamma.cpp | 48 ++++++++++++++++++++++++++++++++++++----- ginac/inifcns_trans.cpp | 6 +++--- ginac/numeric.cpp | 12 +++++++++-- 4 files changed, 59 insertions(+), 10 deletions(-) diff --git a/ginac/inifcns.h b/ginac/inifcns.h index 42f38389..0ccde8de 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -82,6 +82,9 @@ DECLARE_FUNCTION_1P(Li3) /** Gamma function. */ DECLARE_FUNCTION_1P(gamma) +/** Psi function (aka polygamma-function) */ +DECLARE_FUNCTION_2P(psi) + /** Factorial function. */ DECLARE_FUNCTION_1P(factorial) diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index fd3c32ef..7bd6669b 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -43,7 +43,7 @@ namespace GiNaC { * @exception fail_numeric("complex_infinity") or something similar... */ static ex gamma_eval(ex const & x) { - if ( x.info(info_flags::numeric) ) { + if (x.info(info_flags::numeric)) { // trap integer arguments: if ( x.info(info_flags::integer) ) { @@ -62,14 +62,14 @@ static ex gamma_eval(ex const & x) numeric n = ex_to_numeric(x).sub(numHALF()); numeric coefficient = doublefactorial(n.mul(numTWO()).sub(numONE())); coefficient = coefficient.div(numTWO().power(n)); - return coefficient * power(Pi,numHALF()); + return coefficient * pow(Pi,numHALF()); } 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(numHALF())); numeric coefficient = numeric(-2).power(n); coefficient = coefficient.div(doublefactorial(n.mul(numTWO()).sub(numONE())));; - return coefficient * power(Pi,numHALF()); + return coefficient*sqrt(Pi); } } } @@ -89,14 +89,14 @@ static ex gamma_diff(ex const & x, unsigned diff_param) { ASSERT(diff_param==0); - return power(x, -1); // FIXME + return psi(exZERO(),x)*gamma(x); } static ex gamma_series(ex const & x, symbol const & s, ex const & point, int order) { // FIXME: Only handle one special case for now... if (x.is_equal(s) && point.is_zero()) { - ex e = 1 / s - EulerGamma + s * (power(Pi, 2) / 12 + power(EulerGamma, 2) / 2) + Order(power(s, 2)); + ex e = 1 / s - EulerGamma + s * (pow(Pi, 2) / 12 + pow(EulerGamma, 2) / 2) + Order(pow(s, 2)); return e.series(s, point, order); } else throw(std::logic_error("don't know the series expansion of this particular gamma function")); @@ -104,4 +104,42 @@ static ex gamma_series(ex const & x, symbol const & s, ex const & point, int ord REGISTER_FUNCTION(gamma, gamma_eval, gamma_evalf, gamma_diff, gamma_series); +////////// +// psi function (aka polygamma function) +////////// + +/** Evaluation of polygamma-function psi(n,x). + * Somebody ought to provide some good numerical evaluation some day... */ +static ex psi_eval(ex const & n, ex const & x) +{ + if (n.info(info_flags::numeric) && x.info(info_flags::numeric)) { + // do some stuff... + } + return psi(n, x).hold(); +} + +static ex psi_evalf(ex const & n, ex const & x) +{ + BEGIN_TYPECHECK + TYPECHECK(n,numeric) + TYPECHECK(x,numeric) + END_TYPECHECK(psi(n,x)) + + return psi(ex_to_numeric(n), ex_to_numeric(x)); +} + +static ex psi_diff(ex const & n, ex const & x, unsigned diff_param) +{ + ASSERT(diff_param==0); + + return psi(n+1, x); +} + +static ex psi_series(ex const & n, ex const & x, symbol const & s, ex const & point, int order) +{ + throw(std::logic_error("Nobody told me how to series expand the psi function. :-(")); +} + +REGISTER_FUNCTION(psi, psi_eval, psi_evalf, psi_diff, psi_series); + } // namespace GiNaC diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 40968672..ce979686 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -465,13 +465,13 @@ static ex atan2_eval(ex const & y, ex const & x) static ex atan2_diff(ex const & y, ex const & x, unsigned diff_param) { ASSERT(diff_param<2); - + if (diff_param==0) { // d/dy atan(y,x) - return power(x*(1+y*y/(x*x)),-1); + return pow(x*(1+y*y/(x*x)),-1); } // d/dx atan(y,x) - return -y*power(x*x+y*y,-1); + return -y*pow(x*x+y*y,-1); } REGISTER_FUNCTION(atan2, atan2_eval, atan2_evalf, atan2_diff, NULL); diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 165ce594..ef278417 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -1063,10 +1063,18 @@ numeric atanh(numeric const & x) } /** The gamma function. - * stub stub stub stub stub stub! */ + * This is only a stub! */ numeric gamma(numeric const & x) { - clog << "gamma(): Nobody expects the Spanish inquisition" << endl; + clog << "gamma(): Does anybody know good way to calculate this numerically?" << endl; + return numeric(0); +} + +/** The psi function (aka polygamma function). + * This is only a stub! */ +numeric psi(numeric const & n, numeric const & x) +{ + clog << "psi(): Does anybody know good way to calculate this numerically?" << endl; return numeric(0); } -- 2.44.0