[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] lnpoch_pos (specfunc) does not handle large inputs properly

From: Reed A. Cartwright
Subject: [Bug-gsl] lnpoch_pos (specfunc) does not handle large inputs properly
Date: Fri, 20 Nov 2015 12:17:22 -0700

lnpoch_pos(a,x) calculates Log[Pochhammer[a,x]], where a >0 and a+x > 0

However, for some situations it does not calculate the proper value due to
improper conditionals.

lnpoch_pos(2**60,20) should return 831.77; however, it returns 0.

This is because the first if statement is true, which (eventually) results
in lnpoch_pos calculating the result via lngamma(a+x) - lngamma(a), which
will be wrong because a+x == a due to machine precision.

lnpoch_pos does appear to have a code (the stirling approximation) for this
situation.  This code path is used if `absx < 0.1*a && a > 15.0`.  However
it is not reached because earlier, `absx*log(GSL_MAX_DBL(a,2.0)) > 0.1` is
tested and found to be true.

I'm not sure why the earlier test exists, but it does seem like it needs to
be modified so that the stirling approximation will be used for appropriate

Reed A. Cartwright, PhD
Barrett Honors Faculty
Assistant Professor of Genomics, Evolution, and Bioinformatics
School of Life Sciences
Human and Comparative Genomics Laboratory
The Biodesign Institute
Arizona State University
Address: The Biodesign Institute, PO Box 875301, Tempe, AZ 85287-5301 USA
Packages: The Biodesign Institute, 1001 S. McAllister Ave, Tempe, AZ
85287-5301 USA
Office: Biodesign A-224A, 1-480-965-9949
Website: http://cartwrig.ht/

reply via email to

[Prev in Thread] Current Thread [Next in Thread]