Mercurial > hg > octave-nkf
changeset 19552:c6437824681c stable
improve Matlab compatibility for gamma function (bug #43551)
* lo-specfun.cc (xgamma): Return Inf instead of NaN for negative
integer arguments.
* mappers.cc: Fix tests for gamma.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 20 Nov 2014 15:48:35 -0500 |
parents | 043440fa7006 |
children | 998628b7963a c490eac28bbb |
files | libinterp/corefcn/mappers.cc liboctave/numeric/lo-specfun.cc |
diffstat | 2 files changed, 17 insertions(+), 15 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/mappers.cc +++ b/libinterp/corefcn/mappers.cc @@ -1058,7 +1058,7 @@ %!test %! ## Test exceptional values %! x = [-Inf, -1, -0, 0, 1, Inf, NaN]; -%! v = [NaN, NaN, -Inf, Inf, 1, Inf, NaN]; +%! v = [Inf, Inf, -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,12 +366,15 @@ { double result; - if (xisnan (x) || (x < 0 && (xisinf (x) || D_NINT (x) == x))) + // Special cases for (near) compatibility with Matlab instead of + // tgamma. Matlab does not have -0. + + if (x == 0) + result = xnegative_sign (x) ? -octave_Inf : octave_Inf; + else if ((x < 0 && D_NINT (x) == x) || xisinf (x)) + result = octave_Inf; + else if (xisnan (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) @@ -379,8 +382,6 @@ #else F77_XFCN (xdgamma, XDGAMMA, (x, result)); #endif - if (xisinf (result) && (static_cast<int> (gnulib::floor (x)) % 2)) - result = -octave_Inf; } return result; @@ -437,12 +438,15 @@ { float result; - if (xisnan (x) || (x < 0 && (xisinf (x) || D_NINT (x) == x))) + // Special cases for (near) compatibility with Matlab instead of + // tgamma. Matlab does not have -0. + + if (x == 0) + result = xnegative_sign (x) ? -octave_Float_Inf : octave_Float_Inf; + else if ((x < 0 && D_NINT (x) == x) || xisinf (x)) + result = octave_Float_Inf; + else if (xisnan (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_TGAMMA) @@ -450,8 +454,6 @@ #else F77_XFCN (xgamma, XGAMMA, (x, result)); #endif - if (xisinf (result) && (static_cast<int> (gnulib::floor (x)) % 2)) - result = -octave_Float_Inf; } return result;