X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_gamma.cpp;h=ef0169b5c9ea1539b5c92c10499baac02c521e9b;hp=d180065dbce9e743706097e086278fa9a51e630d;hb=a465898a5d582beb8bdf47e9684ddd931acebc75;hpb=c28015b35e3d6ac132a040032b28c79143a36d1f diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index d180065d..ef0169b5 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -31,11 +31,10 @@ #include "power.h" #include "relational.h" #include "symbol.h" +#include "symmetry.h" #include "utils.h" -#ifndef NO_NAMESPACE_GINAC namespace GiNaC { -#endif // ndef NO_NAMESPACE_GINAC ////////// // Logarithm of Gamma function @@ -43,11 +42,13 @@ namespace GiNaC { static ex lgamma_evalf(const ex & x) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - END_TYPECHECK(lgamma(x)) + if (is_exactly_a(x)) { + try { + return lgamma(ex_to(x)); + } catch (const dunno &e) { } + } - return lgamma(ex_to_numeric(x)); + return lgamma(x).hold(); } @@ -63,7 +64,7 @@ static ex lgamma_eval(const ex & x) if (x.info(info_flags::integer)) { // lgamma(n) -> log((n-1)!) for postitive n if (x.info(info_flags::posint)) - return log(factorial(x.exadd(_ex_1()))); + return log(factorial(x + _ex_1())); else throw (pole_error("lgamma_eval(): logarithmic pole",0)); } @@ -100,9 +101,9 @@ static ex lgamma_series(const ex & arg, 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); + numeric m = -ex_to(arg_pt); ex recur; - for (numeric p; p<=m; ++p) + for (numeric p = 0; p<=m; ++p) recur += log(arg+p); return (lgamma(arg+m+_ex1())-recur).series(rel, order, options); } @@ -111,7 +112,8 @@ static ex lgamma_series(const ex & arg, REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval). evalf_func(lgamma_evalf). derivative_func(lgamma_deriv). - series_func(lgamma_series)); + series_func(lgamma_series). + latex_name("\\log \\Gamma")); ////////// @@ -120,11 +122,13 @@ REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval). static ex tgamma_evalf(const ex & x) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - END_TYPECHECK(tgamma(x)) + if (is_exactly_a(x)) { + try { + return tgamma(ex_to(x)); + } catch (const dunno &e) { } + } - return tgamma(ex_to_numeric(x)); + return tgamma(x).hold(); } @@ -137,30 +141,27 @@ static ex tgamma_eval(const ex & x) { if (x.info(info_flags::numeric)) { // trap integer arguments: - if (x.info(info_flags::integer)) { + const numeric two_x = _num2()*ex_to(x); + if (two_x.is_even()) { // tgamma(n) -> (n-1)! for postitive n - if (x.info(info_flags::posint)) { - return factorial(ex_to_numeric(x).sub(_num1())); + if (two_x.is_positive()) { + return factorial(ex_to(x).sub(_num1())); } else { throw (pole_error("tgamma_eval(): simple pole",1)); } } // trap half integer arguments: - if ((x*2).info(info_flags::integer)) { + 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 ((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()); + if (two_x.is_positive()) { + const numeric n = ex_to(x).sub(_num1_2()); + return (doublefactorial(n.mul(_num2()).sub(_num1())).div(pow(_num2(),n))) * pow(Pi,_ex1_2()); } else { // trap negative x==(-n+1/2) // tgamma(-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()); + const numeric n = abs(ex_to(x).sub(_num1_2())); + return (pow(_num_2(), n).div(doublefactorial(n.mul(_num2()).sub(_num1()))))*power(Pi,_ex1_2()); } } // tgamma_evalf should be called here once it becomes available @@ -196,7 +197,7 @@ static ex tgamma_series(const ex & arg, 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: - numeric m = -ex_to_numeric(arg_pt); + const numeric m = -ex_to(arg_pt); ex ser_denom = _ex1(); for (numeric p; p<=m; ++p) ser_denom *= arg+p; @@ -207,7 +208,8 @@ static ex tgamma_series(const ex & arg, REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval). evalf_func(tgamma_evalf). derivative_func(tgamma_deriv). - series_func(tgamma_series)); + series_func(tgamma_series). + latex_name("\\Gamma")); ////////// @@ -216,12 +218,13 @@ REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval). static ex beta_evalf(const ex & x, const ex & y) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - TYPECHECK(y,numeric) - END_TYPECHECK(beta(x,y)) + if (is_exactly_a(x) && is_exactly_a(y)) { + try { + return tgamma(ex_to(x))*tgamma(ex_to(y))/tgamma(ex_to(x+y)); + } catch (const dunno &e) { } + } - return tgamma(ex_to_numeric(x))*tgamma(ex_to_numeric(y))/tgamma(ex_to_numeric(x+y)); + return beta(x,y).hold(); } @@ -231,8 +234,8 @@ static ex beta_eval(const ex & x, const ex & y) // treat all problematic x and y that may not be passed into tgamma, // because they would throw there although beta(x,y) is well-defined // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y) - numeric nx(ex_to_numeric(x)); - numeric ny(ex_to_numeric(y)); + const numeric nx = ex_to(x); + const numeric ny = ex_to(y); if (nx.is_real() && nx.is_integer() && ny.is_real() && ny.is_integer()) { if (nx.is_negative()) { @@ -251,11 +254,10 @@ static ex beta_eval(const ex & x, const ex & y) } // no problem in numerator, but denominator has pole: if ((nx+ny).is_real() && - (nx+ny).is_integer() && - !(nx+ny).is_positive()) + (nx+ny).is_integer() && + !(nx+ny).is_positive()) return _ex0(); - // everything is ok: - return tgamma(x)*tgamma(y)/tgamma(x+y); + // beta_evalf should be called here once it becomes available } return beta(x,y).hold(); @@ -290,24 +292,24 @@ static ex beta_series(const ex & arg1, const ex arg1_pt = arg1.subs(rel); const ex arg2_pt = arg2.subs(rel); GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol)); - const symbol *s = static_cast(rel.lhs().bp); + const symbol &s = ex_to(rel.lhs()); ex arg1_ser, arg2_ser, arg1arg2_ser; if ((!arg1_pt.info(info_flags::integer) || arg1_pt.info(info_flags::positive)) && (!arg2_pt.info(info_flags::integer) || arg2_pt.info(info_flags::positive))) throw do_taylor(); // caught by function::series() // trap the case where arg1 is on a pole: if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive)) - arg1_ser = tgamma(arg1+*s).series(rel, order, options); + arg1_ser = tgamma(arg1+s).series(rel, order, options); else arg1_ser = tgamma(arg1).series(rel,order); // trap the case where arg2 is on a pole: if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive)) - arg2_ser = tgamma(arg2+*s).series(rel, order, options); + arg2_ser = tgamma(arg2+s).series(rel, order, options); else arg2_ser = tgamma(arg2).series(rel,order); // trap the case where arg1+arg2 is on a pole: if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive)) - arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel, order, options); + arg1arg2_ser = tgamma(arg2+arg1+s).series(rel, order, options); else arg1arg2_ser = tgamma(arg2+arg1).series(rel,order); // compose the result (expanding all the terms): @@ -318,7 +320,9 @@ static ex beta_series(const ex & arg1, REGISTER_FUNCTION(beta, eval_func(beta_eval). evalf_func(beta_evalf). derivative_func(beta_deriv). - series_func(beta_series)); + series_func(beta_series). + latex_name("\\mbox{B}"). + set_symmetry(sy_symm(0, 1))); ////////// @@ -327,11 +331,13 @@ REGISTER_FUNCTION(beta, eval_func(beta_eval). static ex psi1_evalf(const ex & x) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - END_TYPECHECK(psi(x)) + if (is_exactly_a(x)) { + try { + return psi(ex_to(x)); + } catch (const dunno &e) { } + } - return psi(ex_to_numeric(x)); + return psi(x).hold(); } /** Evaluation of digamma-function psi(x). @@ -339,13 +345,13 @@ static ex psi1_evalf(const ex & x) static ex psi1_eval(const ex & x) { if (x.info(info_flags::numeric)) { - numeric nx = ex_to_numeric(x); + const numeric nx = ex_to(x); if (nx.is_integer()) { // integer case if (nx.is_positive()) { // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - Euler - numeric rat(0); - for (numeric i(nx+_num_1()); i.is_positive(); --i) + numeric rat = 0; + for (numeric i(nx+_num_1()); i>0; --i) rat += i.inverse(); return rat-Euler; } else { @@ -357,18 +363,18 @@ static ex psi1_eval(const ex & x) // half integer case if (nx.is_positive()) { // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - Euler - 2log(2) - numeric rat(0); - for (numeric i((nx+_num_1())*_num2()); i.is_positive(); i-=_num2()) - rat += _num2()*i.inverse(); - return rat-Euler-_ex2()*log(_ex2()); + numeric rat = 0; + for (numeric i = (nx+_num_1())*_num2(); i>0; i-=_num2()) + rat += _num2()*i.inverse(); + return rat-Euler-_ex2()*log(_ex2()); } else { // use the recurrence relation // psi(-m-1/2) == psi(-m-1/2+1) - 1 / (-m-1/2) // to relate psi(-m-1/2) to psi(1/2): // psi(-m-1/2) == psi(1/2) + r // where r == ((-1/2)^(-1) + ... + (-m-1/2)^(-1)) - numeric recur(0); - for (numeric p(nx); p<0; ++p) + numeric recur = 0; + for (numeric p = nx; p<0; ++p) recur -= pow(p, _num_1()); return recur+psi(_ex1_2()); } @@ -404,7 +410,7 @@ static ex psi1_series(const ex & arg, 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: - numeric m = -ex_to_numeric(arg_pt); + const numeric m = -ex_to(arg_pt); ex recur; for (numeric p; p<=m; ++p) recur += power(arg+p,_ex_1()); @@ -417,6 +423,7 @@ const unsigned function_index_psi1 = evalf_func(psi1_evalf). derivative_func(psi1_deriv). series_func(psi1_series). + latex_name("\\psi"). overloaded(2)); ////////// @@ -425,12 +432,13 @@ const unsigned function_index_psi1 = static ex psi2_evalf(const ex & n, const ex & x) { - BEGIN_TYPECHECK - TYPECHECK(n,numeric) - TYPECHECK(x,numeric) - END_TYPECHECK(psi(n,x)) + if (is_exactly_a(n) && is_exactly_a(x)) { + try { + return psi(ex_to(n),ex_to(x)); + } catch (const dunno &e) { } + } - return psi(ex_to_numeric(n), ex_to_numeric(x)); + return psi(n,x).hold(); } /** Evaluation of polygamma-function psi(n,x). @@ -445,8 +453,8 @@ static ex psi2_eval(const ex & n, const ex & x) return log(tgamma(x)); if (n.info(info_flags::numeric) && n.info(info_flags::posint) && x.info(info_flags::numeric)) { - numeric nn = ex_to_numeric(n); - numeric nx = ex_to_numeric(x); + const numeric nn = ex_to(n); + const numeric nx = ex_to(x); if (nx.is_integer()) { // integer case if (nx.is_equal(_num1())) @@ -458,8 +466,8 @@ static ex psi2_eval(const ex & n, const ex & x) // to relate psi(n,m) to psi(n,1): // psi(n,m) == psi(n,1) + r // where r == (-)^n * n! * (1^(-n-1) + ... + (m-1)^(-n-1)) - numeric recur(0); - for (numeric p(1); p(arg_pt); ex recur; for (numeric p; p<=m; ++p) recur += power(arg+p,-n+_ex_1()); @@ -543,9 +551,8 @@ const unsigned function_index_psi2 = evalf_func(psi2_evalf). derivative_func(psi2_deriv). series_func(psi2_series). + latex_name("\\psi"). overloaded(2)); -#ifndef NO_NAMESPACE_GINAC } // namespace GiNaC -#endif // ndef NO_NAMESPACE_GINAC