underflow handling in besselK (PR#2179)

[bolker@zoo.ufl.edu]

  The besselK() function knows about overflows/underflows internally;
there is a constant xmax_BESS_K in src/nmath/bessel.h (and referred to
only in bessel_k.c), equal to 705.342, which is checked if expon.scaled is
FALSE.  (The equivalent number for bessel_i.c is 709, defined as
exparg_BESS in bessel.h.) However, besselK(x) silently returns +Inf if
x>705.342.  This behavior is reasonable for besselI (where the problem is
an overflow), but seems wrong for besselK (since the problem here is an
*underflow*).  The answer may be as simple as returning 0 rather than +Inf
in this case.
 
  If this is appropriate, the following patch should do the job:

*** bessel_k.c	Thu Oct 17 13:55:55 2002
--- bessel_k.c.dist	Thu Oct 17 13:48:53 2002
***************
*** 216,224 ****
  	    ML_ERROR(ME_RANGE);
  	    *ncalc = *nb;
  	    for(i=0; i < *nb; i++)
! 		if (ex<=0)
! 	   	   bk[i] = ML_POSINF;
! 	        else bk[i]=0.0;
  	    return;
  	}
  	k = 0;
--- 216,222 ----
  	    ML_ERROR(ME_RANGE);
  	    *ncalc = *nb;
  	    for(i=0; i < *nb; i++)
! 		bk[i] = ML_POSINF;
  	    return;
  	}
  	k = 0;

  (By the way, the test for x<=0 seems inappropriate; it looks like x<0 is 
caught elsewhere and that this is really a test for x==0.)

  [Apologies if this is not considered a bug, but my judgment is that it
is -- couldn't be bothered with a round of sending it to r-devel first.]

  Ben Bolker

--please do not edit the information below--

Version:
 platform = i686-pc-linux-gnu
 arch = i686
 os = linux-gnu
 system = i686, linux-gnu
 status = 
 major = 1
 minor = 6.0
 year = 2002
 month = 10
 day = 01
 language = R

Search Path:
 .GlobalEnv, package:ctest, Autoloads, package:base

-- 
318 Carr Hall                                bolker@zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._