gamma(x) and lgamma(x) no longer warn when they correctly return Inf …

…or zero git-svn-id: https://svn.r-project.org/R/trunk@75839 00db46b3-68df-0310-9c12-caf00c1e9a41
wch · Dec 12, 2018 · 891224c · 891224c
1 parent 1e1a15f
commit 891224c
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 8 deletions.
diff --git a/doc/NEWS.Rd b/doc/NEWS.Rd
@@ -237,6 +237,10 @@
       \code{UNPROTECT_PTR}, which is not safe to mix with \code{UNPROTECT}
       (and with \code{PROTECT_WITH_INDEX}). Intended for use in parsers
       only.
+
+      \item \code{gamma(x)} and \code{lgamma(x)} no longer warn when
+      correctly returning \code{Inf} or underflowing to zero.  This
+      helps maximum likelihood and similar computations.
     }
   }
 
@@ -279,7 +283,7 @@
 
       \item \command{configure --with-blas} (without specifying a value)
       includes OpenBLAS in its search (before ATLAS and a generic BLAS).
-      
+
       \item The \command{configure} macro \samp{MAKEINFO} has been
       updated to \samp{TEXI2ANY}.
     }

diff --git a/src/nmath/gamma.c b/src/nmath/gamma.c
@@ -1,8 +1,8 @@
 /*
  *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 2000-2018 The R Core Team
+ *  Copyright (C) 2002-2018 The R Foundation
  *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000-2013 The R Core Team
- *  Copyright (C) 2002-2004 The R Foundation
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -182,12 +182,12 @@ double gammafn(double x)
 	/* gamma(x) for	 y = |x| > 10. */
 
 	if (x > xmax) {			/* Overflow */
-	    ML_ERROR(ME_RANGE, "gammafn");
+	    // No warning: +Inf is the best answer
 	    return ML_POSINF;
 	}
 
 	if (x < xmin) {			/* Underflow */
-	    ML_ERROR(ME_UNDERFLOW, "gammafn");
+	    // No warning: 0 is the best answer
 	    return 0.;
 	}
 

diff --git a/src/nmath/lgamma.c b/src/nmath/lgamma.c
@@ -1,7 +1,7 @@
 /*
  *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 2000-2018 The R Core Team
  *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000-2012 The R Core Team
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -72,7 +72,7 @@ double lgammafn_sign(double x, int *sgn)
 	*sgn = -1;
 
     if (x <= 0 && x == trunc(x)) { /* Negative integer argument */
-	ML_ERROR(ME_RANGE, "lgamma");
+	// No warning: this is the best answer; was  ML_ERROR(ME_RANGE, "lgamma");
 	return ML_POSINF;/* +Inf, since lgamma(x) = log|gamma(x)| */
     }
 
@@ -84,7 +84,7 @@ double lgammafn_sign(double x, int *sgn)
       ELSE  y = |x| > 10 ---------------------- */
 
     if (y > xmax) {
-	ML_ERROR(ME_RANGE, "lgamma");
+	// No warning: +Inf is the best answer
 	return ML_POSINF;
     }
 

diff --git a/tests/reg-tests-1d.R b/tests/reg-tests-1d.R
@@ -2415,6 +2415,17 @@ for (i in 1:2) { if (i == 1) { x <- i; rm(i) }}
 stopifnot(x == 1)
 
 
+## gamma & lgamma should not warn for correct limit cases:
+stopifnot(exprs = {
+    lgamma(0:-10) == Inf
+    gamma(-180.5) == 0
+    gamma(c(200,Inf)) == Inf
+    lgamma(c(10^(306:310), Inf)) == Inf
+})
+## had  "Warning message:  value out of range in 'lgamma' "  for ever
+
+
+
 ## keep at end
 rbind(last =  proc.time() - .pt,
       total = proc.time())