> Date: Wed, 5 Feb 2014 01:57:33 -0700
> From: David Coppa <dco...@openbsd.org>
>
>
> 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() ?
Oh, btw, the quad-precision code used on sparc64 gets this right. So
the bug is probably somewhere in src/lib/limb/src/ld80.
