Bug#153022: Note that this bug is arch-specific

[Daniel Martin <martin@snowplow.org>]

Or at least, that the various manifestations of it are arch-specific.

I noticed that this bug was unreproducible on my own (i386) machine,
so I decided to try it on other archs.  I used the i386 results as a
baseline on the assumption that i386 code is probably more heavily
tested.

First off, the i386 results:
cush:/tmp$ ./testrnd 1 2 3 4 5
N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 2.7182818284590455
N: exp(2) = 7.3890560989306504
Z: exp(2) = 7.3890560989306495
D: exp(2) = 7.3890560989306495
U: exp(2) = 7.3890560989306504
N: exp(3) = 20.085536923187668
Z: exp(3) = 20.085536923187664
D: exp(3) = 20.085536923187664
U: exp(3) = 20.085536923187668
N: exp(4) = 54.598150033144236
Z: exp(4) = 54.598150033144236
D: exp(4) = 54.598150033144236
U: exp(4) = 54.598150033144243
N: exp(5) = 148.4131591025766
Z: exp(5) = 148.4131591025766
D: exp(5) = 148.4131591025766
U: exp(5) = 148.41315910257663

On voltaire (a ppc), results are wrong as the submitter reported
(voltaire is running libc6 2.2.5-14).

On auric (sparc64), we have:
fizbin@auric:~$ ./testrnd 1 2 3 4 5
N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 1.3099120373162434e+151
N: exp(2) = 7.3890560989306504
Z: exp(2) = 4.0037745305985499
D: exp(2) = 4.0037745305985499
U: exp(2) = 7.3890560989306504
N: exp(3) = 20.085536923187668
Z: exp(3) = 20.085536923187664
D: exp(3) = 20.085536923187664
U: exp(3) = 0.015794933022372046
N: exp(4) = 54.598150033144236
Z: exp(4) = 32.045580221125
D: exp(4) = 32.045580221125
U: exp(4) = 54.598150033144243
N: exp(5) = 148.4131591025766
Z: exp(5) = 148.4131591025766
D: exp(5) = 148.4131591025766
U: exp(5) = 21.332962833859369

On caballero (ia64), results happen identically to my home i386
system.

On debussy (arm), only FE_ROUNDTONEAREST is defined (is this a bug?
bits/fenv.h seems to imply that the other rounding modes are available
on arm but accessed through a different API), and it gives the correct
results.

On faure (alpha), reasonably correct results happen, although every
rounding mode produces the same result as ROUNDTONEAREST does on i386.
(i.e. the results are not the same as my home i386 system)  (is this a
bug?)

I don't have access to an m68k machine; perhaps someone else with
access to kullervo can test it there.

On sarti (hppa), the results are the same as on alpha - every rounding
mode produces the same result as ROUNDTONEAREST.

On trex (s390), the following (different, wrong) results occur:
(base)fizbin@trex:~$ ./testrnd 1 2 3 4 5
N: exp(1) = 2.7182818284590451
Z: exp(1) = 2.7182818284590451
D: exp(1) = 2.7182818284590451
U: exp(1) = 6.2906695523485417e-26
N: exp(2) = 7.3890560989306504
Z: exp(2) = 4.0037745305985499
D: exp(2) = 4.0037745305985499
U: exp(2) = 7.3890560989306504
N: exp(3) = 20.085536923187668
Z: exp(3) = 20.085536923187664
D: exp(3) = 20.085536923187664
U: exp(3) = 57.005568249477001
N: exp(4) = 54.598150033144236
Z: exp(4) = 32.045580221125
D: exp(4) = 32.045580221125
U: exp(4) = 54.598150033144243
N: exp(5) = 148.4131591025766
Z: exp(5) = 148.4131591025766
D: exp(5) = 148.4131591025766
U: exp(5) = 0.87712424561247715

So, the summary is:
sparc, ppc, s390: each wrong in their own unique way
alpha, hppa: every rounding mode produces the same result as
   FE_ROUNDTONEAREST (possible bug?)
arm: only FE_ROUNDTONEAREST exists (possible bug?)
ia64, i386: expected results
m68k: untested

Also, I modified the test program to that it would compile on archs
missing one or more rounding modes (and then test only those modes
that exist on that architecture).  The modified version is attached.
It should be buildable anywhere with
  make LDFLAGS=-lm testrnd

It might be worth it at this point for someone with access to the IEEE
spec for how to implement these various rounding modes to implement
exp() using integer arithmetic and these various rounding modes just
as a crosscheck that the i386 results are indeed correct to begin
with.

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

static int rnd[] = { 
#ifdef FE_TONEAREST
	FE_TONEAREST,
#endif
#ifdef FE_TOWARDZERO
	FE_TOWARDZERO,
#endif
#ifdef FE_DOWNWARD
	FE_DOWNWARD,
#endif
#ifdef FE_UPWARD
	FE_UPWARD,
#endif
	0
};
static char rc[] = { 
#ifdef FE_TONEAREST
	'N',
#endif
#ifdef FE_TOWARDZERO
	'Z',
#endif
#ifdef FE_DOWNWARD
	'D',
#endif
#ifdef FE_UPWARD
	'U', 
#endif
	0
};

int main (int argc, char *argv[])
{
  int i;

  for (i = 1; i < argc; i++)
    {
      int r;
      double x;
      char *end;

      x = strtod (argv[i], &end);
      if (*end != '\0')
        exit (EXIT_FAILURE);

      for (r = 0; r < sizeof(rc) - 1; r++)
        {
          double y;

          fesetround (rnd[r]);
          y = exp (x);
          printf ("%c: exp(%.17g) = %.17g\n", rc[r], x, y);
        }
    }

  return 0;
}