[comp.sys.mips] floating point bug?

tasayco@phoenix.Princeton.EDU (Francisco Figueirido) (12/06/90)

Some time ago I was writing my own floating point routines for an
Atari (MC68000 processor) and while checking that it gave correct
answers, I found that the MIPS processor (machine: SGI 4D/240, OS
version 4D1-3.3) gives a slightly different answer than a
SparcStation. I am including the program I used to test the operations
and also the diffs between the output of both processors. I believe
the correctly rounded result is the one given by the Sparc (I regret to
say this because I really like the MIPS). As I am not an expert on the
architercture of either processor or, for that matter, floating point
stuff, I would appreciate any comments regarding this matter.

The following is the program I used (sorry, I guess I could have
posted a smaller version). Compile it and run with

                     testdf -n 100

(if testdf is the name of the executable). This will run a loop of 100
iterations adding two numbers (see the source code for more
information).

==========================================================================

#include <stdio.h>

double atof(/* const char * */);

#define NUMCMDLIN 4

#define ADD 1
#define SUBSTRACT 2
#define MULTIPLY 3
#define DIVIDE 4

static char cmdlin[] = 
  { 'x', 'y', 'o', 'n', };

static double fx = 0.333333333333;
static double fy = 0.123456789;
static unsigned int oflag = ADD;
static unsigned int N = 1000;

void __get_fx(argcp, argvp)
int *argcp;
char ***argvp;
{
  (*argvp)++; (*argcp)--;
  if (*(*argvp) != (char *)0) fx = atof(*(*argvp));
}

void __get_fy(argcp, argvp)
int *argcp;
char ***argvp;
{
  (*argvp)++; (*argcp)--;
  if (*(*argvp) != (char *)0) fy = atof(*(*argvp));
}

void __get_operator(argcp, argvp)
int *argcp;
char ***argvp;
{
  char op;

  (*argvp)++; (*argcp)--;
  if (*(*argvp) != (char *)0) {
    op = *(*(*argvp));
    switch (op) {
    case '+' : oflag = ADD; break;
    case '-' : oflag = SUBSTRACT; break;
    case '*' : oflag = MULTIPLY; break;
    case '/' : oflag = DIVIDE; break;
    }
  }
}

void __get_iterations(argcp, argvp)
int *argcp;
char ***argvp;
{
  (*argvp)++; (*argcp)--;
  if (*(*argvp) != (char *)0) N = atoi(*(*argvp));
}

static void (* getlin[])() = {
  __get_fx,
  __get_fy,
  __get_operator,
  __get_iterations,
};

extern int main(argc, argv)
int argc;
char **argv;
{
  unsigned int i, j;
  union { double d; int i[2]; } d2i;
  union { float f; int i; } f2i;

 /* Check the command line arguments */
  argv++;
  argc--;
  while (argc > 0) {
    if ((*argv)[0] != '-') {
      /* skip this argument */
      argc--;
      argv++;
    }
    for (i = 0; i < NUMCMDLIN; i++) {
      if ( argc > 0 && (*argv)[1] == cmdlin[i]) 
	(* getlin[i])(&argc, &argv);
    }
    /* go to next argument */
    argc--;
    argv++;
  }


  printf("\n");
  d2i.d = fx;
  f2i.f = (float)fx;
  printf("fx = %08x  %08x = %08x (%lf)  \n", 
	 d2i.i[0], d2i.i[1], f2i.i, fx); 
  d2i.d = fy;
  f2i.f = (float)fy;
  printf("fy = %08x  %08x = %08x (%lf)  \n", 
	 d2i.i[0], d2i.i[1], f2i.i, fy);

  printf("\n"); 

/* Compare the numbers as doubles */
  if (fx == fy) 
    printf("The numbers are equal (as doubles).\n\n");
  else 
    if (fx > fy) 
      printf("fx is larger than fy (as doubles).\n\n");
    else 
      if (fx < fy) 
        printf("fy is larger than fx (as doubles).\n\n");
      else printf("Cannot decide for doubles!\n\n");

/* Compare the numbers as floats */
  if ((float)fx == (float)fy) 
    printf("The numbers are equal (as floats).\n\n");
  else 
    if ((float)fx > (float)fy) 
      printf("fx is larger than fy (as floats).\n\n");
    else 
      if ((float)fx < (float)fy) 
        printf("fy is larger than fx (as floats).\n\n");
      else printf("Cannot decide for floats!\n\n");

  d2i.d = fx;
  f2i.f = (float)fx;

  for (j = 0; j < N; j++) {
    printf("  >> %08x  %08x = %08x \n", 
	   d2i.i[0], d2i.i[1], f2i.i); 
    switch (oflag) {
    case ADD :       d2i.d += fy; f2i.f += (float)fy; break;
    case SUBSTRACT : d2i.d -= fy; f2i.f -= (float)fy; break;
    case MULTIPLY :  d2i.d *= fy; f2i.f *= (float)fy; break;
    case DIVIDE :    d2i.d /= fy; f2i.f /= (float)fy; break;
    }
  }

  printf("\n"); 

}

==========================================================================

This is the diff between the output for a Sparc and a MIPS:

10c10
<   >> 3fdd3c0c  a32396b9 = 3ee9e066     (Sparc)
--- 
>   >> 3fdd3c0c  a32396b9 = 3ee9e065     (MIPS)

==========================================================================

	Francisco Figueirido
	e-mail:	tasayco@phoenix.princeton.edu
		figuei@max.physics.sunysb.edu