From: Richard Kreckel Date: Thu, 2 Dec 1999 22:28:00 +0000 (+0000) Subject: - added the beta function to GiNaC X-Git-Tag: release_0-5-0~98 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=c76a85d9ad2a9beeac6e26cb5c24a1bdbca911ed - added the beta function to GiNaC - threw out ginsh's beta function - added a check for first polygamma function --- diff --git a/check/inifcns_consist.cpp b/check/inifcns_consist.cpp index 024dad79..25e34197 100644 --- a/check/inifcns_consist.cpp +++ b/check/inifcns_consist.cpp @@ -169,8 +169,8 @@ static unsigned inifcns_consist_trans(void) return result; } -/* Simple tests on the Gamma combinatorial function. We stuff in arguments - * where the result exists in closed form and check if it's ok. */ +/* Simple tests on the Gamma function. We stuff in arguments where the results + * exists in closed form and check if it's ok. */ static unsigned inifcns_consist_gamma(void) { unsigned result = 0; @@ -214,6 +214,27 @@ static unsigned inifcns_consist_gamma(void) return result; } +/* Simple tests on the Psi-function (aka polygamma-function). We stuff in + arguments where the result exists in closed form and check if it's ok. */ +static unsigned inifcns_consist_psi(void) +{ + unsigned result = 0; + symbol x; + ex e; + + // We check psi(1) and psi(1/2) implicitly by calculating the curious + // little identity gamma(1)'/gamma(1) - gamma(1/2)'/gamma(1/2) == 2*log(2). + e += (gamma(x).diff(x)/gamma(x)).subs(x==numeric(1)); + e -= (gamma(x).diff(x)/gamma(x)).subs(x==numeric(1,2)); + if (e!=2*log(2)) { + clog << "gamma(1)'/gamma(1) - gamma(1/2)'/gamma(1/2) erroneously returned " + << e << " instead of 2*log(2)" << endl; + ++result; + } + + return result; +} + /* Simple tests on the Riemann Zeta function. We stuff in arguments where the * result exists in closed form and check if it's ok. Of course, this checks * the Bernoulli numbers as a side effect. */ @@ -253,6 +274,7 @@ unsigned inifcns_consist(void) result += inifcns_consist_cos(); result += inifcns_consist_trans(); result += inifcns_consist_gamma(); + result += inifcns_consist_psi(); result += inifcns_consist_zeta(); if ( !result ) { diff --git a/ginac/inifcns.h b/ginac/inifcns.h index 84e0dbcb..a5a52700 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -88,7 +88,11 @@ DECLARE_FUNCTION_1P(zeta) /** Gamma-function. */ DECLARE_FUNCTION_1P(gamma) +/** Beta-function. */ +DECLARE_FUNCTION_2P(beta) + /** Psi-function (aka polygamma-function). */ +// overloading @ work: we cannot use the macros extern const unsigned function_index_psi1; inline function psi(ex const & p1) { return function(function_index_psi1, p1); diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index 2fd58cce..0d59eb30 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -1,7 +1,7 @@ /** @file inifcns_gamma.cpp * - * Implementation of Gamma-function, Polygamma-functions, and some related - * stuff. */ + * Implementation of Gamma-function, Beta-function, Polygamma-functions, and + * some related stuff. */ /* * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany @@ -53,7 +53,7 @@ static ex gamma_eval(ex const & x) if (x.info(info_flags::posint)) { return factorial(ex_to_numeric(x).sub(numONE())); } else { - return numZERO(); // Infinity. Throw? What? + throw (std::domain_error("gamma_eval(): simple pole")); } } // trap half integer arguments: @@ -106,6 +106,77 @@ 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); +////////// +// Beta-function +////////// + +static ex beta_eval(ex const & x, ex const & y) +{ + if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { + numeric nx(ex_to_numeric(x)); + numeric ny(ex_to_numeric(y)); + // treat all problematic x and y that may not be passed into gamma, + // because they would throw there although beta(x,y) is well-defined: + if (nx.is_real() && nx.is_integer() && + ny.is_real() && ny.is_integer()) { + if (nx.is_negative()) { + if (nx<=-ny) + return numMINUSONE().power(ny)*beta(1-x-y, y); + else + throw (std::domain_error("beta_eval(): simple pole")); + } + if (ny.is_negative()) { + if (ny<=-nx) + return numMINUSONE().power(nx)*beta(1-y-x, x); + else + throw (std::domain_error("beta_eval(): simple pole")); + } + return gamma(x)*gamma(y)/gamma(x+y); + } + // no problem in numerator, but denominator has pole: + if ((nx+ny).is_real() && + (nx+ny).is_integer() && + !(nx+ny).is_positive()) + return exZERO(); + return gamma(x)*gamma(y)/gamma(x+y); + } + return beta(x,y).hold(); +} + +static ex beta_evalf(ex const & x, ex const & y) +{ + BEGIN_TYPECHECK + TYPECHECK(x,numeric) + TYPECHECK(y,numeric) + END_TYPECHECK(beta(x,y)) + + return gamma(ex_to_numeric(x))*gamma(ex_to_numeric(y)) + / gamma(ex_to_numeric(x+y)); +} + +static ex beta_diff(ex const & x, ex const & y, unsigned diff_param) +{ + GINAC_ASSERT(diff_param<2); + ex retval; + + if (diff_param==0) // d/dx beta(x,y) + retval = (psi(x)-psi(x+y))*beta(x,y); + if (diff_param==1) // d/dy beta(x,y) + retval = (psi(y)-psi(x+y))*beta(x,y); + return retval; +} + +static ex beta_series(ex const & x, ex const & y, symbol const & s, ex const & point, int order) +{ + if (x.is_equal(s) && point.is_zero()) { + 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 beta function")); +} + +REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, beta_series); + ////////// // Psi-function (aka polygamma-function) ////////// diff --git a/ginsh/ginsh.1 b/ginsh/ginsh.1 index 32b0cec6..cf3b945c 100644 --- a/ginsh/ginsh.1 +++ b/ginsh/ginsh.1 @@ -216,9 +216,6 @@ detail here. Please refer to the GiNaC documentation. .PP .RS \" GINSH_FCN_HELP_START -.BI beta( expression ", " expression ) -\- beta function -.br .BI charpoly( matrix ", " symbol ) \- characteristic polynomial of a matrix .br diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index c49e7c0d..2e02b6cc 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -242,7 +242,6 @@ static void push(const ex &e) * Built-in functions */ -static ex f_beta(const exprseq &e) {return gamma(e[0])*gamma(e[1])/gamma(e[0]+e[1]);} static ex f_denom(const exprseq &e) {return e[0].denom();} static ex f_eval1(const exprseq &e) {return e[0].eval();} static ex f_evalf1(const exprseq &e) {return e[0].evalf();} @@ -483,7 +482,6 @@ struct fcn_init { }; static const fcn_init builtin_fcns[] = { - {"beta", fcn_desc(f_beta, 2)}, {"charpoly", fcn_desc(f_charpoly, 2)}, {"coeff", fcn_desc(f_coeff, 3)}, {"collect", fcn_desc(f_collect, 2)}, @@ -759,8 +757,10 @@ int main(int argc, char **argv) insert_fcn_help("atan", "inverse tangent function"); insert_fcn_help("atan2", "inverse tangent function with two arguments"); insert_fcn_help("atanh", "inverse hyperbolic tangent function"); + insert_fcn_help("beta", "beta function"); insert_fcn_help("cos", "cosine function"); insert_fcn_help("cosh", "hyperbolic cosine function"); + insert_fcn_help("psi", "polygamma function"); insert_fcn_help("sin", "sine function"); insert_fcn_help("sinh", "hyperbolic sine function"); insert_fcn_help("tan", "tangent function");