emulated s390 mpfr problem

[Greg Smith]

Hello,

I'm back!  I am still having trouble getting the s390 emulator to pass
mpfr checks that pass on a real s390.  I composed a testcase,
mpfr_add.c, that is the first addendum below.  I get one result on a
real machine and another on the emulator.

But one thing is confusing me.  Running the testcase on intel (pentium
4) I am getting a weird result.  I don't see how the particular test
passes on intel.  So I hacked mpfr/tests/tadd.c to run the same exact
test and it passes!  I am thinking if I understand what is going on, I
can fix the emulator.

Per your previous suggestions, I am working with gmp-4.1.4.  The
emulator does not emulate floating point but rather uses the host's
floating point implementation (in this case pentium 4).  For example,
the s390 instruction that sets the rounding mode, SRNM, calls
fesetround.  (This was my original problem, host rounding mode was not
being set).

I fully understand that the host architecture may give different results
in floating point arithmetic than s390 due to the number of guard bits,
etc.  

The testcase is a special instance where the precision is px == py == pz
== 53.  Also the bad assumption is made that sizeof(double) == sizeof
(unsigned long long) == 8.

Below are the results of the testcase and the hacked tadd.c:

Real s390
------------------------------------------
~/gmp-4.1.4/mpfr/tests> ./mpfr_add
 x=89e9b6bd4a167678  y=07fdb24440062c08
z1=89e9b6bd49db11ef z2=89e9b6bd49db11ef

~/gmp-4.1.4/mpfr/tests> ./tadd
rounding mode set to GMP_RNDN
 x=89e9b6bd4a167678  y=07fdb24440062c08
z1=89e9b6bd49db11ef z2=89e9b6bd49db11ef


Intel pentium 4
------------------------------------------
$ ./mpfr_add
 x=89e9b6bd4a167678  y=07fdb24440062c08
z1=89e9b6bd49db11f0 z2=89e9b6bd49db11ef

$ ./tadd
rounding mode set to GMP_RNDN
 x=89e9b6bd4a167678  y=07fdb24440062c08
z1=89e9b6bd49db11ef z2=89e9b6bd49db11ef


Emulated s390
------------------------------------------
~/gmp-4.1.4/mpfr/tests> ./mpfr_add
 x=89e9b6bd4a167678  y=07fdb24440062c08
z1=89e9b6bd49db11f0 z2=89e9b6bd49db11ef

~/gmp-4.1.4/mpfr/tests> ./tadd
rounding mode set to GMP_RNDN
 x=89e9b6bd4a167678  y=07fdb24440062c08
z1=89e9b6bd49db11f0 z2=89e9b6bd49db11ef

z1 is the double from adding doubles x + y; z2 is the value extracted
from the result of mpfr_add.

What I can't understand is why the pentium 4 z1 values differ between
the two tests.  I am thinking that if I can understand this then I have
a chance for fixing the emulator.

Finally, please let me know if I should rather attach the addenda below
as a gzipped file.  Oh, and my mailer is probably going to wrap some of
the lines in the tadd.c hack.

Thanks,
Greg Smith



[1] mpfr_add.c
/*
 *  rnd_mode=GMP_RNDN  px=53 py=53 pz=53
 *  x=89e9b6bd4a167678  y=07fdb24440062c08
 * z1=89e9b6bd49db11f0 z2=89e9b6bd49db11ef
 */
#include <stdio.h>
#include <fenv.h>
#include "gmp.h"
#include "mpfr.h"

typedef union {
        double f;
        unsigned long long i;
} DOUBLE;

int main(int argc, char *argv[])
{
 DOUBLE x,y,z1,z2;
 mpfr_t xx,yy,zz;
 mp_rnd_t rnd_mode = GMP_RNDN;

 x.i = 0x89e9b6bd4a167678ull;
 y.i = 0x07fdb24440062c08ull;

 mpfr_init2(xx,53);
 mpfr_init2(yy,53);
 mpfr_init2(zz,53);
 mpfr_set_d(xx, x.f, rnd_mode);
 mpfr_set_d(yy, y.f, rnd_mode);
 mpfr_add(zz, xx, yy, rnd_mode);

 mpfr_set_machine_rnd_mode(rnd_mode);

 z1.f = x.f + y.f;
 z2.f = mpfr_get_d1(zz);

 printf(" x=%16.16llx  y=%16.16llx\n",x.i,y.i);
 printf("z1=%16.16llx z2=%16.16llx\n",z1.i,z2.i);
}


[2] tadd.pat
--- tadd.c.orig 2002-11-21 16:35:37.000000000 -0500
+++ tadd.c      2005-05-31 13:34:01.000000000 -0400
@@ -29,6 +29,11 @@
 #include "mpfr-impl.h"
 #include "mpfr-test.h"
 
+typedef union {
+        double f;
+        unsigned long long i;
+} DOUBLE;
+
 void checknan _PROTO((double, double, mp_rnd_t, unsigned int, unsigned
int,
unsigned int)); 
 void check3 _PROTO((double, double, mp_rnd_t));
 void check4 _PROTO((double, double, mp_rnd_t));
@@ -55,28 +60,35 @@
 /* checks that x+y gives the same results in double
    and with mpfr with 53 bits of precision */
 void
-_check (double x, double y, double z1, mp_rnd_t rnd_mode, unsigned int
px, 
+_check (double _x, double _y, double _z1, mp_rnd_t rnd_mode, unsigned
int px, 
         unsigned int py, unsigned int pz)
 {
-  double z2; mpfr_t xx,yy,zz; int cert=0;
+  DOUBLE x,y,z1,z2; mpfr_t xx,yy,zz; int cert=0;
+
+  x.f = _x; y.f = _y; z1.f = _z1;
 
   mpfr_init2(xx, px);
   mpfr_init2(yy, py);
   mpfr_init2(zz, pz);
-  mpfr_set_d(xx, x, rnd_mode);
-  mpfr_set_d(yy, y, rnd_mode);
+  mpfr_set_d(xx, x.f, rnd_mode);
+  mpfr_set_d(yy, y.f, rnd_mode);
   mpfr_add(zz, xx, yy, rnd_mode);
 #ifdef MPFR_HAVE_FESETROUND
   mpfr_set_machine_rnd_mode(rnd_mode);
   if (px==53 && py==53 && pz==53) cert=1;
+printf("rounding mode set to %s\n",mpfr_print_rnd_mode(rnd_mode));
 #endif
-  if (z1==0.0) z1=x+y; else cert=1;
-  z2 = mpfr_get_d1 (zz);
-  mpfr_set_d (yy, z2, GMP_RNDN);
-  if (!mpfr_cmp (zz, yy) && cert && z1!=z2 && !(isnan(z1) && isnan
(z2))) {
-    printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
+  if (z1.f==0.0) z1.f=x.f+y.f; else cert=1;
+  z2.f = mpfr_get_d1 (zz);
+
+printf(" x=%16.16llx  y=%16.16llx\n",x.i,y.i);
+printf("z1=%16.16llx z2=%16.16llx\n",z1.i,z2.i);
+
+  mpfr_set_d (yy, z2.f, GMP_RNDN);
+  if (!mpfr_cmp (zz, yy) && cert && z1.f!=z2.f && !(isnan(z1.f) &&
isnan(z2.f))) {
+    printf("expected sum is %1.20e, got %1.20e\n",z1.f,z2.f);
     printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
-          x, y, mpfr_print_rnd_mode(rnd_mode));
+          x.f, y.f, mpfr_print_rnd_mode(rnd_mode));
     exit(1);
   }
   mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
@@ -658,7 +670,15 @@
   double x;
   int i;
 
+  DOUBLE xx,yy;
+  xx.i = 0x89e9b6bd4a167678ull;
+  yy.i = 0x07fdb24440062c08ull;
+
   mpfr_test_init ();
+
+  check(xx.f, yy.f, GMP_RNDN, 53, 53, 53, 0.0);
+  exit(0);
+
   check_inexact ();
   check_case_1b ();
   check_case_2 ();