Mercurial > octave-dspies

changeset 17708:f10b7a578e2c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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	Sun Oct 20 23:43:26 2013 +0200
+++ b/libinterp/corefcn/mappers.cc	Tue Jul 02 14:17:33 2013 +0100
@@ -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	Sun Oct 20 23:43:26 2013 +0200
+++ b/liboctave/numeric/lo-specfun.cc	Tue Jul 02 14:17:33 2013 +0100
@@ -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;
 }