# underflow handling in besselK (PR#2179)

**maechler@stat.math.ethz.ch
**
maechler@stat.math.ethz.ch

*Fri, 18 Oct 2002 17:16:38 +0200 (MET DST)*

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