exp() / expl() on Linux and OpenBSD (expl() bug?)

[David Coppa-4]

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

[Mark Kettenis]

> Date: Wed, 5 Feb 2014 01:57:33 -0700
> From: David Coppa <[hidden email]>
>
> 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() ?

Yes.

[Mark Kettenis]

> Date: Wed, 5 Feb 2014 01:57:33 -0700
> From: David Coppa <[hidden email]>
>
>
> 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.

[Donovan Watteau]

David Coppa wrote:

> 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

FWIW, it looks even stranger on loongson:

$ cc -o expl expl.c -O2 -pipe -lm
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ ./expl
theta == 1, pr == 0.606531, pr2 == 0.606531
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944

$ sysctl kern.version
kern.version=OpenBSD 5.5-beta (GENERIC) #106: Mon Feb  3 01:47:15 MST 2014
    [hidden email]:/usr/src/sys/arch/loongson/compile/GENERIC

[Daniel Dickman]

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 <[hidden email]> 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
>

[Martynas Venckus-4]

Yup.    Does this diff fix it for you?

On 2/6/14, Daniel Dickman <[hidden email]> 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 <[hidden email]> 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
>>
>

[David Coppa]

On Fri, Feb 7, 2014 at 8:07 AM, Martynas Venckus <[hidden email]> wrote:
> Yup.    Does this diff fix it for you?

Yeah! It works.

And R-3's testsuite is also happy now.

Thanks a lot! And thanks to Daniel too...

ok dcoppa@ to commit it, obviously

Ciao,
David

> On 2/6/14, Daniel Dickman <[hidden email]> 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 <[hidden email]> 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
>>>
>>

[Mark Kettenis]

> Date: Thu, 6 Feb 2014 23:07:58 -0800
> From: Martynas Venckus <[hidden email]>
>
> Yup.    Does this diff fix it for you?

Here's a diff that sticks a bit closer to the original code.  It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.

Index: s_floorl.c
===================================================================
RCS file: /home/cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_floorl.c
--- s_floorl.c 25 Jul 2011 16:20:09 -0000 1.2
+++ s_floorl.c 7 Feb 2014 14:43:19 -0000
@@ -38,7 +38,7 @@ floorl(long double x)
  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;}
+ { se=0xbfff;i0=0x80000000;i1=0;}
  }
     } else {
  i = (0x7fffffff)>>jj0;

[David Coppa]

On Fri, Feb 7, 2014 at 3:49 PM, Mark Kettenis <[hidden email]> wrote:

>> Date: Thu, 6 Feb 2014 23:07:58 -0800
>> From: Martynas Venckus <[hidden email]>
>>
>> Yup.    Does this diff fix it for you?
>
> Here's a diff that sticks a bit closer to the original code.  It's
> equivalent to your diff, and admittedly purely a matter of taste which
> version to prefer.
>
> Index: s_floorl.c
> ===================================================================
> RCS file: /home/cvs/src/lib/libm/src/ld80/s_floorl.c,v
> retrieving revision 1.2
> diff -u -p -r1.2 s_floorl.c
> --- s_floorl.c  25 Jul 2011 16:20:09 -0000      1.2
> +++ s_floorl.c  7 Feb 2014 14:43:19 -0000
> @@ -38,7 +38,7 @@ floorl(long double x)
>                 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;}
> +                       { se=0xbfff;i0=0x80000000;i1=0;}
>                 }
>             } else {
>                 i = (0x7fffffff)>>jj0;
>

Just tested, and this one works too, of course...

I'd be happy if this or the other one could be committed.

ciao,
David

[David Coppa]

On Fri, Feb 7, 2014 at 4:03 PM, David Coppa <[hidden email]> wrote:

> On Fri, Feb 7, 2014 at 3:49 PM, Mark Kettenis <[hidden email]> wrote:
>>> Date: Thu, 6 Feb 2014 23:07:58 -0800
>>> From: Martynas Venckus <[hidden email]>
>>>
>>> Yup.    Does this diff fix it for you?
>>
>> Here's a diff that sticks a bit closer to the original code.  It's
>> equivalent to your diff, and admittedly purely a matter of taste which
>> version to prefer.
>>
>> Index: s_floorl.c
>> ===================================================================
>> RCS file: /home/cvs/src/lib/libm/src/ld80/s_floorl.c,v
>> retrieving revision 1.2
>> diff -u -p -r1.2 s_floorl.c
>> --- s_floorl.c  25 Jul 2011 16:20:09 -0000      1.2
>> +++ s_floorl.c  7 Feb 2014 14:43:19 -0000
>> @@ -38,7 +38,7 @@ floorl(long double x)
>>                 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;}
>> +                       { se=0xbfff;i0=0x80000000;i1=0;}
>>                 }
>>             } else {
>>                 i = (0x7fffffff)>>jj0;
>>
>
> Just tested, and this one works too, of course...
>
> I'd be happy if this or the other one could be committed.

Ping

[Martynas-2]

> Here's a diff that sticks a bit closer to the original code.  It's
> equivalent to your diff, and admittedly purely a matter of taste which
> version to prefer.

I prefer my version better.  It's not '93 anymore and compilers are
able to convert 0.0L and -1.0L precisely, otherwise we have a huge
problem.  There's no need to obfuscate here by manually converting
to floating point representation.

Here's a diff which also includes same fix for ld128.  OK?

Index: ld128/s_floorl.c
===================================================================
RCS file: /cvs/src/lib/libm/src/ld128/s_floorl.c,v
retrieving revision 1.1
diff -u -r1.1 s_floorl.c
--- ld128/s_floorl.c 6 Jul 2011 00:02:42 -0000 1.1
+++ ld128/s_floorl.c 11 Feb 2014 05:24:15 -0000
@@ -34,10 +34,11 @@
  jj0 = ((i0>>48)&0x7fff)-0x3fff;
  if(jj0<48) {
     if(jj0<0) { /* raise inexact if x != 0 */
- if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
-    if(i0>=0) {i0=i1=0;}
+ if(huge+x>0.0) {
+    if(i0>=0)
+ return 0.0L;
     else if(((i0&0x7fffffffffffffffLL)|i1)!=0)
- { i0=0xbfff000000000000ULL;i1=0;}
+ return -1.0L;
  }
     } else {
  i = (0x0000ffffffffffffULL)>>jj0;
Index: ld80/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
--- ld80/s_floorl.c 25 Jul 2011 16:20:09 -0000 1.2
+++ ld80/s_floorl.c 11 Feb 2014 05:24:01 -0000
@@ -35,10 +35,11 @@
  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;}
+ if(huge+x>0.0) {
+    if(sx==0)
+ return 0.0L;
     else if(((se&0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ return -1.0L;
  }
     } else {
  i = (0x7fffffff)>>jj0;

[David Coppa]

On Tue, Feb 11, 2014 at 6:35 AM, Martynas Venckus
<[hidden email]> wrote:

>> Here's a diff that sticks a bit closer to the original code.  It's
>> equivalent to your diff, and admittedly purely a matter of taste which
>> version to prefer.
>
> I prefer my version better.  It's not '93 anymore and compilers are
> able to convert 0.0L and -1.0L precisely, otherwise we have a huge
> problem.  There's no need to obfuscate here by manually converting
> to floating point representation.
>
> Here's a diff which also includes same fix for ld128.  OK?

OK with me.

Ciao,
David

> Index: ld128/s_floorl.c
> ===================================================================
> RCS file: /cvs/src/lib/libm/src/ld128/s_floorl.c,v
> retrieving revision 1.1
> diff -u -r1.1 s_floorl.c
> --- ld128/s_floorl.c    6 Jul 2011 00:02:42 -0000       1.1
> +++ ld128/s_floorl.c    11 Feb 2014 05:24:15 -0000
> @@ -34,10 +34,11 @@
>         jj0 = ((i0>>48)&0x7fff)-0x3fff;
>         if(jj0<48) {
>             if(jj0<0) {         /* raise inexact if x != 0 */
> -               if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
> -                   if(i0>=0) {i0=i1=0;}
> +               if(huge+x>0.0) {
> +                   if(i0>=0)
> +                       return 0.0L;
>                     else if(((i0&0x7fffffffffffffffLL)|i1)!=0)
> -                       { i0=0xbfff000000000000ULL;i1=0;}
> +                       return -1.0L;
>                 }
>             } else {
>                 i = (0x0000ffffffffffffULL)>>jj0;
> Index: ld80/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
> --- ld80/s_floorl.c     25 Jul 2011 16:20:09 -0000      1.2
> +++ ld80/s_floorl.c     11 Feb 2014 05:24:01 -0000
> @@ -35,10 +35,11 @@
>         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;}
> +               if(huge+x>0.0) {
> +                   if(sx==0)
> +                       return 0.0L;
>                     else if(((se&0x7fff)|i0|i1)!=0)
> -                       { se=0xbfff;i0=i1=0;}
> +                       return -1.0L;
>                 }
>             } else {
>                 i = (0x7fffffff)>>jj0;
>

[Mark Kettenis]

> Date: Mon, 10 Feb 2014 22:35:02 -0700 (MST)
> From: Martynas Venckus <[hidden email]>
>
> > Here's a diff that sticks a bit closer to the original code.  It's
> > equivalent to your diff, and admittedly purely a matter of taste which
> > version to prefer.
>
> I prefer my version better.  It's not '93 anymore and compilers are
> able to convert 0.0L and -1.0L precisely, otherwise we have a huge
> problem.  There's no need to obfuscate here by manually converting
> to floating point representation.
>
> Here's a diff which also includes same fix for ld128.  OK?

Fair enough.  Please commit.

> Index: ld128/s_floorl.c
> ===================================================================
> RCS file: /cvs/src/lib/libm/src/ld128/s_floorl.c,v
> retrieving revision 1.1
> diff -u -r1.1 s_floorl.c
> --- ld128/s_floorl.c 6 Jul 2011 00:02:42 -0000 1.1
> +++ ld128/s_floorl.c 11 Feb 2014 05:24:15 -0000
> @@ -34,10 +34,11 @@
>   jj0 = ((i0>>48)&0x7fff)-0x3fff;
>   if(jj0<48) {
>      if(jj0<0) { /* raise inexact if x != 0 */
> - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
> -    if(i0>=0) {i0=i1=0;}
> + if(huge+x>0.0) {
> +    if(i0>=0)
> + return 0.0L;
>      else if(((i0&0x7fffffffffffffffLL)|i1)!=0)
> - { i0=0xbfff000000000000ULL;i1=0;}
> + return -1.0L;
>   }
>      } else {
>   i = (0x0000ffffffffffffULL)>>jj0;
> Index: ld80/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
> --- ld80/s_floorl.c 25 Jul 2011 16:20:09 -0000 1.2
> +++ ld80/s_floorl.c 11 Feb 2014 05:24:01 -0000
> @@ -35,10 +35,11 @@
>   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;}
> + if(huge+x>0.0) {
> +    if(sx==0)
> + return 0.0L;
>      else if(((se&0x7fff)|i0|i1)!=0)
> - { se=0xbfff;i0=i1=0;}
> + return -1.0L;
>   }
>      } else {
>   i = (0x7fffffff)>>jj0;
>
>

[Donovan Watteau]

On Thu, 6 Feb 2014, Donovan Watteau wrote:

> David Coppa wrote:
> > 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
>
> FWIW, it looks even stranger on loongson:
>
> $ cc -o expl expl.c -O2 -pipe -lm
> $ ./expl
> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
> $ ./expl
> theta == 1, pr == 0.606531, pr2 == 0.606531
> $ ./expl
> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
>
> $ sysctl kern.version
> kern.version=OpenBSD 5.5-beta (GENERIC) #106: Mon Feb  3 01:47:15 MST 2014
>     [hidden email]:/usr/src/sys/arch/loongson/compile/GENERIC

A fix has been committed, but there's still a problem on loongson with
libm updated:

$ ls -l /usr/lib/libm.so.*
-r--r--r--  1 root  bin  926033 Feb 12 12:17 /usr/lib/libm.so.9.0
$ cc -o expl expl.c -O2 -pipe -lm
$ for in in 1 2 3 4 5 6 ; do ./expl ; done
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531

[Martynas Venckus-4]

On 2/12/14, Donovan Watteau <[hidden email]> wrote:

> On Thu, 6 Feb 2014, Donovan Watteau wrote:
>> David Coppa wrote:
>> > 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
>>
>> FWIW, it looks even stranger on loongson:
>>
>> $ cc -o expl expl.c -O2 -pipe -lm
>> $ ./expl
>> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
>> $ ./expl
>> theta == 1, pr == 0.606531, pr2 == 0.606531
>> $ ./expl
>> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
>>
>> $ sysctl kern.version
>> kern.version=OpenBSD 5.5-beta (GENERIC) #106: Mon Feb  3 01:47:15 MST
>> 2014
>>     [hidden email]:/usr/src/sys/arch/loongson/compile/GENERIC
>
> A fix has been committed, but there's still a problem on loongson with
> libm updated:
>
> $ ls -l /usr/lib/libm.so.*
> -r--r--r--  1 root  bin  926033 Feb 12 12:17 /usr/lib/libm.so.9.0
> $ cc -o expl expl.c -O2 -pipe -lm
> $ for in in 1 2 3 4 5 6 ; do ./expl ; done
> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
> theta == 1, pr == 0.606531, pr2 == 0.606531
> theta == 1, pr == 0.606531, pr2 == 0.606531
> theta == 1, pr == 0.606531, pr2 == 0.606531
> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
> theta == 1, pr == 0.606531, pr2 == 0.606531

This isn't related to exp/expl.  Looks like a bug in either compiler
or libc/printf/gdtoa.  I don't have the hardware so I couldn't tell
much more.

As usual, there are a few ways to fix it:
1. Debug it and provide a diff,
2. Donate hardware.  Oh, BTW I have a request for loongson in
http://www.openbsd.org/want.html sitting for over a year...

[Miod Vallat]

> A fix has been committed, but there's still a problem on loongson with
> libm updated:
>
> $ ls -l /usr/lib/libm.so.*
> -r--r--r--  1 root  bin  926033 Feb 12 12:17 /usr/lib/libm.so.9.0
> $ cc -o expl expl.c -O2 -pipe -lm
> $ for in in 1 2 3 4 5 6 ; do ./expl ; done
> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
> theta == 1, pr == 0.606531, pr2 == 0.606531
> theta == 1, pr == 0.606531, pr2 == 0.606531
> theta == 1, pr == 0.606531, pr2 == 0.606531
> theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
> theta == 1, pr == 0.606531, pr2 == 0.606531
>

This is definitely a problem in gdtoa. The values are always computed
the same, but do not get printed correctly, as shown by this modified
test program. Both lines should print 0.606531 for the long doubles, but
don't always; this smells like an uninitialized value somewhere.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int
main(void)
{
        double theta = 1;
        long double lambda, pr, pr2;
        static const uint64_t zpr[2] = {
#ifdef __MIPSEL__
        0xa000000000000000ULL, 0x3ffe368b2fc6f960ULL
#else
        0x3ffe368b2fc6f960ULL, 0xa000000000000000ULL
#endif
        }, zpr2[2] = {
#ifdef __MIPSEL__
        0x9fe7aceb46aa619cULL, 0x3ffe368b2fc6f960ULL
#else
        0x3ffe368b2fc6f960ULL, 0x9fe7aceb46aa619cULL
#endif
        };

        lambda = (0.5*theta);
        pr = exp(-lambda);
        pr2 = expl(-lambda);

        printf("theta == %g, pr == %Lg, pr2 == %Lg\n", theta, pr, pr2);
        printf("theta == %g, pr == %Lg, pr2 == %Lg\n", theta, zpr[0], zpr[1],
            zpr2[0], zpr2[1]);
        exit(0);
}

[Miod Vallat]

> > A fix has been committed, but there's still a problem on loongson with
> > libm updated:

[...]

> This is definitely a problem in gdtoa. The values are always computed
> the same, but do not get printed correctly, as shown by this modified
> test program. Both lines should print 0.606531 for the long doubles, but
> don't always; this smells like an uninitialized value somewhere.

Well, it turned out this was caused by the program stack being aligned
to an 8 byte boundary, instead of a 16 byte boundary, causing long
double variables on the stack to be misaligned 50% of the time.

This has been fixed now.

Miod