Yup. Does this diff fix it for you?
On 2/6/14, Daniel Dickman <didick...@gmail.com> wrote:
> I think I recently ran into a similar issue but I suspect the root cause
> might be the same. I think the floorl function is wrong for numbers slightly
> larger than -1 to numbers slightly below 0. In this range floorl returns -0
> instead of -1.
>
>> On Feb 5, 2014, at 3:57 AM, David Coppa <dco...@openbsd.org> wrote:
>>
>>
>> Hi!
>>
>> I hit this problem while working on updating math/R from version
>> 2.15.3 to the latest version (3.0.2).
>>
>> It started happening since upstream switched from double functions
>> to C99 long double functions (expl, fabsl, ...), during the R-3
>> development cycle.
>>
>> Take the following reduced test-case, adapted from what R's code
>> does:
>>
>> ---8<---
>>
>> #include <stdio.h>
>> #include <stdlib.h>
>> #include <math.h>
>>
>> int main(void) {
>> double theta = 1;
>> long double lambda, pr, pr2;
>>
>> lambda = (0.5*theta);
>> pr = exp(-lambda);
>> pr2 = expl(-lambda);
>>
>> printf("theta == %g, pr == %Lg, pr2 == %Lg\n", theta, pr, pr2);
>> exit(0);
>> }
>>
>> ---8<---
>>
>> This produces the following output on Linux (x86_64):
>>
>> theta == 1, pr == 0.606531, pr2 == 0.606531
>>
>> While on OpenBSD -current amd64:
>>
>> theta == 1, pr == 0.606531, pr2 == nan
>>
>> And indeed R-3's testsuite fails with the error message
>> "NaNs produced":
>>
>> Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
>>> stopifnot(all.equal(p00, exp(-lam/2)),
>> + all.equal(p.0, exp(-lam/2)))
>> Error: all.equal(p.0, exp(-lam/2)) is not TRUE
>> Execution halted
>>
>> Is this a bug in our expl() ?
>>
>> Ciao,
>> David
>>
>
Index: s_floorl.c
===================================================================
RCS file: /cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -r1.2 s_floorl.c
--- s_floorl.c 25 Jul 2011 16:20:09 -0000 1.2
+++ s_floorl.c 7 Feb 2014 07:01:59 -0000
@@ -35,10 +35,12 @@
jj0 = (se&0x7fff)-0x3fff;
if(jj0<31) {
if(jj0<0) { /* raise inexact if x != 0 */
- if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
- if(sx==0) {se=0;i0=i1=0;}
- else if(((se&0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ if(huge+x>0.0) {
+ if(sx==0) {
+ return 0.0L;
+ } else if(((se&0x7fff)|i0|i1)!=0) {
+ return -1.0L;
+ }
}
} else {
i = (0x7fffffff)>>jj0;
