Mercurial > hg > octave-nkf
changeset 17708:f10b7a578e2c
Correct return values of gamma() (see Numerical, item 3 on Projects page).
* liboctave/numeric/lo-specfun.cc(xgamma): Check for special cases of
NaN, -Inf, -integer and return NaN. Check for special case of 0 and
return -Inf.
* libinterp/corefcn/mappers.cc(Fgamma): Add %!tests for exceptional values.
author | Craig Hudson <c_hudson_phd@hotmail.com> |
---|---|
date | Tue, 02 Jul 2013 14:17:33 +0100 |
parents | e26b47e9285e |
children | 5415a9cd61d4 |
files | libinterp/corefcn/mappers.cc liboctave/numeric/lo-specfun.cc |
diffstat | 2 files changed, 26 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/mappers.cc +++ b/libinterp/corefcn/mappers.cc @@ -1056,8 +1056,9 @@ %! assert (gamma (x), v, sqrt (eps ("single"))); %!test -%! x = [-1, 0, 1, Inf]; -%! v = [Inf, Inf, 1, Inf]; +%! ## Test exceptional values +%! x = [-Inf, -1, -0, 0, 1, Inf, NaN]; +%! v = [NaN, NaN, -Inf, Inf, 1, Inf, NaN]; %! assert (gamma (x), v); %! assert (gamma (single (x)), single (v));
--- a/liboctave/numeric/lo-specfun.cc +++ b/liboctave/numeric/lo-specfun.cc @@ -366,16 +366,22 @@ { double result; - if (xisnan (x)) - result = x; - else if ((x <= 0 && D_NINT (x) == x) || xisinf (x)) + if (xisnan (x) || (x < 0 && (xisinf (x) || D_NINT (x) == x))) + result = octave_NaN; + else if (x == 0 && xnegative_sign (x)) + result = -octave_Inf; + else if (x == 0 || xisinf (x)) result = octave_Inf; else + { #if defined (HAVE_TGAMMA) - result = tgamma (x); + result = tgamma (x); #else - F77_XFCN (xdgamma, XDGAMMA, (x, result)); + F77_XFCN (xdgamma, XDGAMMA, (x, result)); #endif + if (xisinf (result) && ((int)floor (x) % 2)) + result = -octave_Inf; + } return result; } @@ -431,16 +437,22 @@ { float result; - if (xisnan (x)) - result = x; - else if ((x <= 0 && D_NINT (x) == x) || xisinf (x)) + if (xisnan (x) || (x < 0 && (xisinf (x) || D_NINT (x) == x))) + result = octave_Float_NaN; + else if (x == 0 && xnegative_sign (x)) + result = -octave_Float_Inf; + else if (x == 0 || xisinf (x)) result = octave_Float_Inf; else -#if defined (HAVE_TGAMMAF) - result = tgammaf (x); + { +#if defined (HAVE_TGAMMA) + result = tgammaf (x); #else - F77_XFCN (xgamma, XGAMMA, (x, result)); + F77_XFCN (xgamma, XGAMMA, (x, result)); #endif + if (xisinf (result) && ((int)floor (x) % 2)) + result = -octave_Float_Inf; + } return result; }