underflow handling in besselK (PR#2179)

[maechler@stat.math.ethz.ch]

>>>>> "Ben" == Ben Bolker <bolker@zoo.ufl.edu>
>>>>>     on Thu, 17 Oct 2002 19:58:33 +0200 (MET DST) writes:

    Ben> The besselK() function knows about overflows/underflows internally;
    Ben> there is a constant xmax_BESS_K in src/nmath/bessel.h (and referred to
    Ben> only in bessel_k.c), equal to 705.342, which is checked if expon.scaled is
    Ben> FALSE.  (The equivalent number for bessel_i.c is 709, defined as
    Ben> 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
    Ben> an overflow), but seems wrong for besselK (since the problem here is an
    Ben> *underflow*).  The answer may be as simple as returning 0 rather than +Inf
    Ben> in this case.
 
Yes, indeed.

    Ben> If this is appropriate, the following patch should do the job:
    Ben> 
    Ben> <...........>

"R-patched" now contains something like your patch, namely :

 	if(ex <= 0 || (*ize == 1 && ex > xmax_BESS_K)) {
+	    if(ex <= 0) {
 	    ML_ERROR(ME_RANGE);
-	    *ncalc = *nb;
 	    for(i=0; i < *nb; i++)
 		bk[i] = ML_POSINF;
+	    } else /* would only have underflow */
+		for(i=0; i < *nb; i++)
+		    bk[i] = 0.;
+	    *ncalc = *nb;
 	    return;
 	}

Thank you for the report!

Martin Maechler <maechler@stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._