[thirdparty/binutils-gdb.git] / sim / rx / fpu.c

/* fpu.c --- FPU emulator for stand-alone RX simulator.

Copyright (C) 2008-2018 Free Software Foundation, Inc.
Contributed by Red Hat, Inc.

This file is part of the GNU simulators.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.  */

#include "config.h"
#include <stdio.h>
#include <stdlib.h>

#include "cpu.h"
#include "fpu.h"

/* FPU encodings are as follows:

   S EXPONENT MANTISSA
   1 12345678 12345678901234567890123

   0 00000000 00000000000000000000000	+0
   1 00000000 00000000000000000000000	-0

   X 00000000 00000000000000000000001	Denormals
   X 00000000 11111111111111111111111
 
   X 00000001 XXXXXXXXXXXXXXXXXXXXXXX	Normals
   X 11111110 XXXXXXXXXXXXXXXXXXXXXXX

   0 11111111 00000000000000000000000	+Inf
   1 11111111 00000000000000000000000	-Inf

   X 11111111 0XXXXXXXXXXXXXXXXXXXXXX	SNaN (X != 0)
   X 11111111 1XXXXXXXXXXXXXXXXXXXXXX	QNaN (X != 0)

*/

#define trace 0
#define tprintf if (trace) printf

/* Some magic numbers.  */
#define PLUS_MAX   0x7f7fffffUL
#define MINUS_MAX  0xff7fffffUL
#define PLUS_INF   0x7f800000UL
#define MINUS_INF  0xff800000UL
#define PLUS_ZERO  0x00000000UL
#define MINUS_ZERO 0x80000000UL

#define FP_RAISE(e) fp_raise(FPSWBITS_C##e)
static void
fp_raise (int mask)
{
  regs.r_fpsw |= mask;
  if (mask != FPSWBITS_CE)
    {
      if (regs.r_fpsw & (mask << FPSW_CESH))
	regs.r_fpsw |= (mask << FPSW_CFSH);
      if (regs.r_fpsw & FPSWBITS_FMASK)
	regs.r_fpsw |= FPSWBITS_FSUM;
      else
	regs.r_fpsw &= ~FPSWBITS_FSUM;
    }
}

/* We classify all numbers as one of these.  They correspond to the
   rows/colums in the exception tables.  */
typedef enum {
  FP_NORMAL,
  FP_PZERO,
  FP_NZERO,
  FP_PINFINITY,
  FP_NINFINITY,
  FP_DENORMAL,
  FP_QNAN,
  FP_SNAN
} FP_Type;

#if defined DEBUG0
static const char *fpt_names[] = {
  "Normal", "+0", "-0", "+Inf", "-Inf", "Denormal", "QNaN", "SNaN"
};
#endif

#define EXP_BIAS  127
#define EXP_ZERO -127
#define EXP_INF   128

#define MANT_BIAS 0x00080000UL

typedef struct {
  int exp;
  unsigned int mant; /* 24 bits */
  char type;
  char sign;
  fp_t orig_value;
} FP_Parts;

static void
fp_explode (fp_t f, FP_Parts *p)
{
  int exp, mant, sign;

  exp = ((f & 0x7f800000UL) >> 23);
  mant = f & 0x007fffffUL;
  sign = f & 0x80000000UL;
  /*printf("explode: %08x %x %2x %6x\n", f, sign, exp, mant);*/

  p->sign = sign ? -1 : 1;
  p->exp = exp - EXP_BIAS;
  p->orig_value = f;
  p->mant = mant | 0x00800000UL;

  if (p->exp == EXP_ZERO)
    {
      if (regs.r_fpsw & FPSWBITS_DN)
	mant = 0;
      if (mant)
	p->type = FP_DENORMAL;
      else
	{
	  p->mant = 0;
	  p->type = sign ? FP_NZERO : FP_PZERO;
	}
    }
  else if (p->exp == EXP_INF)
    {
      if (mant == 0)
	p->type = sign ? FP_NINFINITY : FP_PINFINITY;
      else if (mant & 0x00400000UL)
	p->type = FP_QNAN;
      else
	p->type = FP_SNAN;
    }
  else
    p->type = FP_NORMAL;
}

static fp_t
fp_implode (FP_Parts *p)
{
  int exp, mant;

  exp = p->exp + EXP_BIAS;
  mant = p->mant;
  /*printf("implode: exp %d mant 0x%x\n", exp, mant);*/
  if (p->type == FP_NORMAL)
    {
      while (mant
	     && exp > 0
	     && mant < 0x00800000UL)
	{
	  mant <<= 1;
	  exp --;
	}
      while (mant > 0x00ffffffUL)
	{
	  mant >>= 1;
	  exp ++;
	}
      if (exp < 0)
	{
	  /* underflow */
	  exp = 0;
	  mant = 0;
	  FP_RAISE (E);
	}
      if (exp >= 255)
	{
	  /* overflow */
	  exp = 255;
	  mant = 0;
	  FP_RAISE (O);
	}
    }
  mant &= 0x007fffffUL;
  exp &= 0xff;
  mant |= exp << 23;
  if (p->sign < 0)
    mant |= 0x80000000UL;

  return mant;
}

typedef union {
  unsigned long long ll;
  double d;
} U_d_ll;

static int checked_format = 0;

/* We assume a double format like this:
   S[1] E[11] M[52]
*/

static double
fp_to_double (FP_Parts *p)
{
  U_d_ll u;

  if (!checked_format)
    {
      u.d = 1.5;
      if (u.ll != 0x3ff8000000000000ULL)
	abort ();
      u.d = -225;
      if (u.ll != 0xc06c200000000000ULL)
	abort ();
      u.d = 10.1;
      if (u.ll != 0x4024333333333333ULL)
	abort ();
      checked_format = 1;
    }

  u.ll = 0;
  if (p->sign < 0)
    u.ll |= (1ULL << 63);
  /* Make sure a zero encoding stays a zero.  */
  if (p->exp != -EXP_BIAS)
    u.ll |= ((unsigned long long)p->exp + 1023ULL) << 52;
  u.ll |= (unsigned long long) (p->mant & 0x007fffffUL) << (52 - 23);
  return u.d;
}

static void
double_to_fp (double d, FP_Parts *p)
{
  int exp;
  U_d_ll u;
  int sign;

  u.d = d;

  sign = (u.ll & 0x8000000000000000ULL) ? 1 : 0;
  exp = u.ll >> 52;
  exp = (exp & 0x7ff);

  if (exp == 0)
    {
      /* A generated denormal should show up as an underflow, not
	 here.  */
      if (sign)
	fp_explode (MINUS_ZERO, p);
      else
	fp_explode (PLUS_ZERO, p);
      return;
    }

  exp = exp - 1023;
  if ((exp + EXP_BIAS) > 254)
    {
      FP_RAISE (O);
      switch (regs.r_fpsw & FPSWBITS_RM)
	{
	case FPRM_NEAREST:
	  if (sign)
	    fp_explode (MINUS_INF, p);
	  else
	    fp_explode (PLUS_INF, p);
	  break;
	case FPRM_ZERO:
	  if (sign)
	    fp_explode (MINUS_MAX, p);
	  else
	    fp_explode (PLUS_MAX, p);
	  break;
	case FPRM_PINF:
	  if (sign)
	    fp_explode (MINUS_MAX, p);
	  else
	    fp_explode (PLUS_INF, p);
	  break;
	case FPRM_NINF:
	  if (sign)
	    fp_explode (MINUS_INF, p);
	  else
	    fp_explode (PLUS_MAX, p);
	  break;
	}
      return;
    }
  if ((exp + EXP_BIAS) < 1)
    {
      if (sign)
	fp_explode (MINUS_ZERO, p);
      else
	fp_explode (PLUS_ZERO, p);
      FP_RAISE (U);
    }

  p->sign = sign ? -1 : 1;
  p->exp = exp;
  p->mant = u.ll >> (52-23) & 0x007fffffUL;
  p->mant |= 0x00800000UL;
  p->type = FP_NORMAL;

  if (u.ll & 0x1fffffffULL)
    {
      switch (regs.r_fpsw & FPSWBITS_RM)
	{
	case FPRM_NEAREST:
	  if (u.ll & 0x10000000ULL)
	    p->mant ++;
	  break;
	case FPRM_ZERO:
	  break;
	case FPRM_PINF:
	  if (sign == 1)
	    p->mant ++;
	  break;
	case FPRM_NINF:
	  if (sign == -1)
	    p->mant ++;
	  break;
	}
      FP_RAISE (X);
    }

}

typedef enum {
  eNR,		/* Use the normal result.  */
  ePZ, eNZ,	/* +- zero */
  eSZ,		/* signed zero - XOR signs of ops together.  */
  eRZ,		/* +- zero depending on rounding mode.  */
  ePI, eNI,	/* +- Infinity */
  eSI,		/* signed infinity - XOR signs of ops together.  */
  eQN, eSN,	/* Quiet/Signalling NANs */
  eIn,		/* Invalid.  */
  eUn,		/* Unimplemented.  */
  eDZ,		/* Divide-by-zero.  */
  eLT,		/* less than */
  eGT,		/* greater than */
  eEQ,		/* equal to */
} FP_ExceptionCases;

#if defined DEBUG0
static const char *ex_names[] = {
  "NR", "PZ", "NZ", "SZ", "RZ", "PI", "NI", "SI", "QN", "SN", "IN", "Un", "DZ", "LT", "GT", "EQ"
};
#endif

/* This checks for all exceptional cases (not all FP exceptions) and
   returns TRUE if it is providing the result in *c.  If it returns
   FALSE, the caller should do the "normal" operation.  */
int
check_exceptions (FP_Parts *a, FP_Parts *b, fp_t *c,
		  FP_ExceptionCases ex_tab[5][5], 
		  FP_ExceptionCases *case_ret)
{
  FP_ExceptionCases fpec;

  if (a->type == FP_SNAN
      || b->type == FP_SNAN)
    fpec = eIn;
  else if (a->type == FP_QNAN
	   || b->type == FP_QNAN)
    fpec = eQN;
  else if (a->type == FP_DENORMAL
	   || b->type == FP_DENORMAL)
    fpec = eUn;
  else
    fpec = ex_tab[(int)(a->type)][(int)(b->type)];

  /*printf("%s %s -> %s\n", fpt_names[(int)(a->type)], fpt_names[(int)(b->type)], ex_names[(int)(fpec)]);*/

  if (case_ret)
    *case_ret = fpec;

  switch (fpec)
    {
    case eNR:	/* Use the normal result.  */
      return 0;

    case ePZ:	/* + zero */
      *c = 0x00000000;
      return 1;

    case eNZ:	/* - zero */
      *c = 0x80000000;
      return 1;

    case eSZ:	/* signed zero */
      *c = (a->sign == b->sign) ? PLUS_ZERO : MINUS_ZERO;
      return 1;

    case eRZ:	/* +- zero depending on rounding mode.  */
      if ((regs.r_fpsw & FPSWBITS_RM) == FPRM_NINF)
	*c = 0x80000000;
      else
	*c = 0x00000000;
      return 1;

    case ePI:	/* + Infinity */
      *c = 0x7F800000;
      return 1;

    case eNI:	/* - Infinity */
      *c = 0xFF800000;
      return 1;

    case eSI:	/* sign Infinity */
      *c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF;
      return 1;

    case eQN:	/* Quiet NANs */
      if (a->type == FP_QNAN)
	*c = a->orig_value;
      else
	*c = b->orig_value;
      return 1;

    case eSN:	/* Signalling NANs */
      if (a->type == FP_SNAN)
	*c = a->orig_value;
      else
	*c = b->orig_value;
      FP_RAISE (V);
      return 1;

    case eIn:	/* Invalid.  */
      FP_RAISE (V);
      if (a->type == FP_SNAN)
	*c = a->orig_value | 0x00400000;
      else if  (b->type == FP_SNAN)
	*c = b->orig_value | 0x00400000;
      else
	*c = 0x7fc00000;
      return 1;

    case eUn:	/* Unimplemented.  */
      FP_RAISE (E);
      return 1;

    case eDZ:	/* Division-by-zero.  */
      *c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF;
      FP_RAISE (Z);
      return 1;

    default:
      return 0;
    }
}

#define CHECK_EXCEPTIONS(FPPa, FPPb, fpc, ex_tab) \
  if (check_exceptions (&FPPa, &FPPb, &fpc, ex_tab, 0))	\
    return fpc;

/* For each operation, we have two tables of how nonnormal cases are
   handled.  The DN=0 case is first, followed by the DN=1 case, with
   each table using the following layout: */

static FP_ExceptionCases ex_add_tab[5][5] = {
  /* N   +0   -0   +In  -In */
  { eNR, eNR, eNR, ePI, eNI }, /* Normal */
  { eNR, ePZ, eRZ, ePI, eNI }, /* +0   */
  { eNR, eRZ, eNZ, ePI, eNI }, /* -0   */
  { ePI, ePI, ePI, ePI, eIn }, /* +Inf */
  { eNI, eNI, eNI, eIn, eNI }, /* -Inf */
};

fp_t
rxfp_add (fp_t fa, fp_t fb)
{
  FP_Parts a, b, c;
  fp_t rv;
  double da, db;

  fp_explode (fa, &a);
  fp_explode (fb, &b);
  CHECK_EXCEPTIONS (a, b, rv, ex_add_tab);

  da = fp_to_double (&a);
  db = fp_to_double (&b);
  tprintf("%g + %g = %g\n", da, db, da+db);

  double_to_fp (da+db, &c);
  rv = fp_implode (&c);
  return rv;
}

static FP_ExceptionCases ex_sub_tab[5][5] = {
  /* N   +0   -0   +In  -In */
  { eNR, eNR, eNR, eNI, ePI }, /* Normal */
  { eNR, eRZ, ePZ, eNI, ePI }, /* +0   */
  { eNR, eNZ, eRZ, eNI, ePI }, /* -0   */
  { ePI, ePI, ePI, eIn, ePI }, /* +Inf */
  { eNI, eNI, eNI, eNI, eIn }, /* -Inf */
};

fp_t
rxfp_sub (fp_t fa, fp_t fb)
{
  FP_Parts a, b, c;
  fp_t rv;
  double da, db;

  fp_explode (fa, &a);
  fp_explode (fb, &b);
  CHECK_EXCEPTIONS (a, b, rv, ex_sub_tab);

  da = fp_to_double (&a);
  db = fp_to_double (&b);
  tprintf("%g - %g = %g\n", da, db, da-db);

  double_to_fp (da-db, &c);
  rv = fp_implode (&c);

  return rv;
}

static FP_ExceptionCases ex_mul_tab[5][5] = {
  /* N   +0   -0   +In  -In */
  { eNR, eNR, eNR, eSI, eSI }, /* Normal */
  { eNR, ePZ, eNZ, eIn, eIn }, /* +0   */
  { eNR, eNZ, ePZ, eIn, eIn }, /* -0   */
  { eSI, eIn, eIn, ePI, eNI }, /* +Inf */
  { eSI, eIn, eIn, eNI, ePI }, /* -Inf */
};

fp_t
rxfp_mul (fp_t fa, fp_t fb)
{
  FP_Parts a, b, c;
  fp_t rv;
  double da, db;

  fp_explode (fa, &a);
  fp_explode (fb, &b);
  CHECK_EXCEPTIONS (a, b, rv, ex_mul_tab);

  da = fp_to_double (&a);
  db = fp_to_double (&b);
  tprintf("%g x %g = %g\n", da, db, da*db);

  double_to_fp (da*db, &c);
  rv = fp_implode (&c);

  return rv;
}

static FP_ExceptionCases ex_div_tab[5][5] = {
  /* N   +0   -0   +In  -In */
  { eNR, eDZ, eDZ, eSZ, eSZ }, /* Normal */
  { eSZ, eIn, eIn, ePZ, eNZ }, /* +0   */
  { eSZ, eIn, eIn, eNZ, ePZ }, /* -0   */
  { eSI, ePI, eNI, eIn, eIn }, /* +Inf */
  { eSI, eNI, ePI, eIn, eIn }, /* -Inf */
};

fp_t
rxfp_div (fp_t fa, fp_t fb)
{
  FP_Parts a, b, c;
  fp_t rv;
  double da, db;

  fp_explode (fa, &a);
  fp_explode (fb, &b);
  CHECK_EXCEPTIONS (a, b, rv, ex_div_tab);

  da = fp_to_double (&a);
  db = fp_to_double (&b);
  tprintf("%g / %g = %g\n", da, db, da/db);

  double_to_fp (da/db, &c);
  rv = fp_implode (&c);

  return rv;
}

static FP_ExceptionCases ex_cmp_tab[5][5] = {
  /* N   +0   -0   +In  -In */
  { eNR, eNR, eNR, eLT, eGT }, /* Normal */
  { eNR, eEQ, eEQ, eLT, eGT }, /* +0   */
  { eNR, eEQ, eEQ, eLT, eGT }, /* -0   */
  { eGT, eGT, eGT, eEQ, eGT }, /* +Inf */
  { eLT, eLT, eLT, eLT, eEQ }, /* -Inf */
};

void
rxfp_cmp (fp_t fa, fp_t fb)
{
  FP_Parts a, b;
  fp_t c;
  FP_ExceptionCases reason;
  int flags = 0;
  double da, db;

  fp_explode (fa, &a);
  fp_explode (fb, &b);

  if (check_exceptions (&a, &b, &c, ex_cmp_tab, &reason))
    {
      if (reason == eQN)
	{
	  /* Special case - incomparable.  */
	  set_flags (FLAGBIT_Z | FLAGBIT_S | FLAGBIT_O, FLAGBIT_O);
	  return;
	}
      return;
    }

  switch (reason)
    {
    case eEQ:
      flags = FLAGBIT_Z;
      break;
    case eLT:
      flags = FLAGBIT_S;
      break;
    case eGT:
      flags = 0;
      break;
    case eNR:
      da = fp_to_double (&a);
      db = fp_to_double (&b);
      tprintf("fcmp: %g cmp %g\n", da, db);
      if (da < db)
	flags = FLAGBIT_S;
      else if (da == db)
	flags = FLAGBIT_Z;
      else
	flags = 0;
      break;
    default:
      abort();
    }

  set_flags (FLAGBIT_Z | FLAGBIT_S | FLAGBIT_O, flags);
}

long
rxfp_ftoi (fp_t fa, int round_mode)
{
  FP_Parts a;
  fp_t rv;
  int sign;
  int whole_bits, frac_bits;

  fp_explode (fa, &a);
  sign = fa & 0x80000000UL;

  switch (a.type)
    {
    case FP_NORMAL:
      break;
    case FP_PZERO:
    case FP_NZERO:
      return 0;
    case FP_PINFINITY:
      FP_RAISE (V);
      return 0x7fffffffL;
    case FP_NINFINITY:
      FP_RAISE (V);
      return 0x80000000L;
    case FP_DENORMAL:
      FP_RAISE (E);
      return 0;
    case FP_QNAN:
    case FP_SNAN:
      FP_RAISE (V);
      return sign ? 0x80000000U : 0x7fffffff;
    }

  if (a.exp >= 31)
    {
      FP_RAISE (V);
      return sign ? 0x80000000U : 0x7fffffff;
    }

  a.exp -= 23;

  if (a.exp <= -25)
    {
      /* Less than 0.49999 */
      frac_bits = a.mant;
      whole_bits = 0;
    }
  else if (a.exp < 0)
    {
      frac_bits = a.mant << (32 + a.exp);
      whole_bits = a.mant >> (-a.exp);
    }
  else
    {
      frac_bits = 0;
      whole_bits = a.mant << a.exp;
    }

  if (frac_bits)
    {
      switch (round_mode & 3)
	{
	case FPRM_NEAREST:
	  if (frac_bits & 0x80000000UL)
	    whole_bits ++;
	  break;
	case FPRM_ZERO:
	  break;
	case FPRM_PINF:
	  if (!sign)
	    whole_bits ++;
	  break;
	case FPRM_NINF:
	  if (sign)
	    whole_bits ++;
	  break;
	}
    }

  rv = sign ? -whole_bits : whole_bits;
  
  return rv;
}

fp_t
rxfp_itof (long fa, int round_mode)
{
  fp_t rv;
  int sign = 0;
  unsigned int frac_bits;
  volatile unsigned int whole_bits;
  FP_Parts a;

  if (fa == 0)
    return PLUS_ZERO;

  if (fa < 0)
    {
      fa = -fa;
      sign = 1;
      a.sign = -1;
    }
  else
    a.sign = 1;

  whole_bits = fa;
  a.exp = 31;

  while (! (whole_bits & 0x80000000UL))
    {
      a.exp --;
      whole_bits <<= 1;
    }
  frac_bits = whole_bits & 0xff;
  whole_bits = whole_bits >> 8;

  if (frac_bits)
    {
      /* We must round */
      switch (round_mode & 3)
	{
	case FPRM_NEAREST:
	  if (frac_bits & 0x80)
	    whole_bits ++;
	  break;
	case FPRM_ZERO:
	  break;
	case FPRM_PINF:
	  if (!sign)
	    whole_bits ++;
	  break;
	case FPRM_NINF:
	  if (sign)
	    whole_bits ++;
	  break;
	}
    }

  a.mant = whole_bits;
  if (whole_bits & 0xff000000UL)
    {
      a.mant >>= 1;
      a.exp ++;
    }

  rv = fp_implode (&a);
  return rv;
}
Commit	Line	Data
4f8d4a38 DD	1	/* fpu.c --- FPU emulator for stand-alone RX simulator.
4f8d4a38 DD	2
e2882c85	3	Copyright (C) 2008-2018 Free Software Foundation, Inc.
4f8d4a38 DD	4	Contributed by Red Hat, Inc.
	5
	6	This file is part of the GNU simulators.
	7
	8	This program is free software; you can redistribute it and/or modify
	9	it under the terms of the GNU General Public License as published by
	10	the Free Software Foundation; either version 3 of the License, or
	11	(at your option) any later version.
	12
	13	This program is distributed in the hope that it will be useful,
	14	but WITHOUT ANY WARRANTY; without even the implied warranty of
	15	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	16	GNU General Public License for more details.
	17
	18	You should have received a copy of the GNU General Public License
	19	along with this program. If not, see <http://www.gnu.org/licenses/>. */
	20
8d794149	21	#include "config.h"
4f8d4a38 DD	22	#include <stdio.h>
	23	#include <stdlib.h>
	24
	25	#include "cpu.h"
	26	#include "fpu.h"
	27
	28	/* FPU encodings are as follows:
	29
	30	S EXPONENT MANTISSA
	31	1 12345678 12345678901234567890123
	32
	33	0 00000000 00000000000000000000000 +0
	34	1 00000000 00000000000000000000000 -0
	35
	36	X 00000000 00000000000000000000001 Denormals
	37	X 00000000 11111111111111111111111
	38
	39	X 00000001 XXXXXXXXXXXXXXXXXXXXXXX Normals
	40	X 11111110 XXXXXXXXXXXXXXXXXXXXXXX
	41
	42	0 11111111 00000000000000000000000 +Inf
	43	1 11111111 00000000000000000000000 -Inf
	44
	45	X 11111111 0XXXXXXXXXXXXXXXXXXXXXX SNaN (X != 0)
	46	X 11111111 1XXXXXXXXXXXXXXXXXXXXXX QNaN (X != 0)
	47
	48	*/
	49
	50	#define trace 0
	51	#define tprintf if (trace) printf
	52
	53	/* Some magic numbers. */
	54	#define PLUS_MAX 0x7f7fffffUL
	55	#define MINUS_MAX 0xff7fffffUL
	56	#define PLUS_INF 0x7f800000UL
	57	#define MINUS_INF 0xff800000UL
	58	#define PLUS_ZERO 0x00000000UL
	59	#define MINUS_ZERO 0x80000000UL
	60
	61	#define FP_RAISE(e) fp_raise(FPSWBITS_C##e)
	62	static void
	63	fp_raise (int mask)
	64	{
	65	regs.r_fpsw \|= mask;
	66	if (mask != FPSWBITS_CE)
	67	{
	68	if (regs.r_fpsw & (mask << FPSW_CESH))
	69	regs.r_fpsw \|= (mask << FPSW_CFSH);
	70	if (regs.r_fpsw & FPSWBITS_FMASK)
	71	regs.r_fpsw \|= FPSWBITS_FSUM;
	72	else
	73	regs.r_fpsw &= ~FPSWBITS_FSUM;
	74	}
	75	}
	76
	77	/* We classify all numbers as one of these. They correspond to the
	78	rows/colums in the exception tables. */
	79	typedef enum {
	80	FP_NORMAL,
	81	FP_PZERO,
	82	FP_NZERO,
	83	FP_PINFINITY,
	84	FP_NINFINITY,
	85	FP_DENORMAL,
86	FP_QNAN,
87	FP_SNAN
88	} FP_Type;
89
90	#if defined DEBUG0
91	static const char *fpt_names[] = {
92	"Normal", "+0", "-0", "+Inf", "-Inf", "Denormal", "QNaN", "SNaN"
93	};
94	#endif
95
96	#define EXP_BIAS 127
97	#define EXP_ZERO -127
98	#define EXP_INF 128
99
100	#define MANT_BIAS 0x00080000UL
101
102	typedef struct {
103	int exp;
104	unsigned int mant; /* 24 bits */
105	char type;
106	char sign;
107	fp_t orig_value;
108	} FP_Parts;
109
110	static void
111	fp_explode (fp_t f, FP_Parts *p)
112	{
113	int exp, mant, sign;
114
115	exp = ((f & 0x7f800000UL) >> 23);
116	mant = f & 0x007fffffUL;
117	sign = f & 0x80000000UL;
118	/printf("explode: %08x %x %2x %6x\n", f, sign, exp, mant);/
119
120	p->sign = sign ? -1 : 1;
121	p->exp = exp - EXP_BIAS;
122	p->orig_value = f;
123	p->mant = mant \| 0x00800000UL;
124
125	if (p->exp == EXP_ZERO)
126	{
127	if (regs.r_fpsw & FPSWBITS_DN)
128	mant = 0;
129	if (mant)
130	p->type = FP_DENORMAL;
131	else
132	{
133	p->mant = 0;
134	p->type = sign ? FP_NZERO : FP_PZERO;
135	}
136	}
137	else if (p->exp == EXP_INF)
138	{
139	if (mant == 0)
140	p->type = sign ? FP_NINFINITY : FP_PINFINITY;
141	else if (mant & 0x00400000UL)
142	p->type = FP_QNAN;
143	else
144	p->type = FP_SNAN;
145	}
146	else
147	p->type = FP_NORMAL;
148	}
149
150	static fp_t
151	fp_implode (FP_Parts *p)
152	{
153	int exp, mant;
154
155	exp = p->exp + EXP_BIAS;
156	mant = p->mant;
157	/printf("implode: exp %d mant 0x%x\n", exp, mant);/
158	if (p->type == FP_NORMAL)
159	{
160	while (mant
161	&& exp > 0
162	&& mant < 0x00800000UL)
163	{
164	mant <<= 1;
165	exp --;
166	}
167	while (mant > 0x00ffffffUL)
168	{
169	mant >>= 1;
170	exp ++;
171	}
172	if (exp < 0)
173	{
174	/* underflow */
175	exp = 0;
176	mant = 0;
177	FP_RAISE (E);
178	}
179	if (exp >= 255)
180	{
181	/* overflow */
182	exp = 255;
183	mant = 0;
184	FP_RAISE (O);
185	}
186	}
187	mant &= 0x007fffffUL;
188	exp &= 0xff;
189	mant \|= exp << 23;
190	if (p->sign < 0)
191	mant \|= 0x80000000UL;
192
193	return mant;
194	}
195
196	typedef union {
197	unsigned long long ll;
198	double d;
199	} U_d_ll;
200
201	static int checked_format = 0;
202
203	/* We assume a double format like this:
204	S[1] E[11] M[52]
205	*/
206
207	static double
208	fp_to_double (FP_Parts *p)
209	{
210	U_d_ll u;
211
212	if (!checked_format)
213	{
214	u.d = 1.5;
215	if (u.ll != 0x3ff8000000000000ULL)
216	abort ();
217	u.d = -225;
218	if (u.ll != 0xc06c200000000000ULL)
219	abort ();
220	u.d = 10.1;
221	if (u.ll != 0x4024333333333333ULL)
222	abort ();
223	checked_format = 1;
224	}
225
226	u.ll = 0;
227	if (p->sign < 0)
228	u.ll \|= (1ULL << 63);
229	/* Make sure a zero encoding stays a zero. */
230	if (p->exp != -EXP_BIAS)
231	u.ll \|= ((unsigned long long)p->exp + 1023ULL) << 52;
232	u.ll \|= (unsigned long long) (p->mant & 0x007fffffUL) << (52 - 23);
233	return u.d;
234	}
235
236	static void
237	double_to_fp (double d, FP_Parts *p)
238	{
239	int exp;
240	U_d_ll u;
241	int sign;
242
243	u.d = d;
244
245	sign = (u.ll & 0x8000000000000000ULL) ? 1 : 0;
246	exp = u.ll >> 52;
247	exp = (exp & 0x7ff);
248
249	if (exp == 0)
250	{
251	/* A generated denormal should show up as an underflow, not
252	here. */
253	if (sign)
254	fp_explode (MINUS_ZERO, p);
255	else
256	fp_explode (PLUS_ZERO, p);
257	return;
258	}
259
260	exp = exp - 1023;
261	if ((exp + EXP_BIAS) > 254)
262	{
263	FP_RAISE (O);
264	switch (regs.r_fpsw & FPSWBITS_RM)
265	{
266	case FPRM_NEAREST:
267	if (sign)
268	fp_explode (MINUS_INF, p);
269	else
270	fp_explode (PLUS_INF, p);
271	break;
272	case FPRM_ZERO:
273	if (sign)
274	fp_explode (MINUS_MAX, p);
275	else
276	fp_explode (PLUS_MAX, p);
277	break;
278	case FPRM_PINF:
279	if (sign)
280	fp_explode (MINUS_MAX, p);
281	else
282	fp_explode (PLUS_INF, p);
283	break;
284	case FPRM_NINF:
285	if (sign)
286	fp_explode (MINUS_INF, p);
287	else
288	fp_explode (PLUS_MAX, p);
289	break;
290	}
291	return;
292	}
293	if ((exp + EXP_BIAS) < 1)
294	{
295	if (sign)
296	fp_explode (MINUS_ZERO, p);
297	else
298	fp_explode (PLUS_ZERO, p);
299	FP_RAISE (U);
300	}
301
302	p->sign = sign ? -1 : 1;
303	p->exp = exp;
304	p->mant = u.ll >> (52-23) & 0x007fffffUL;
305	p->mant \|= 0x00800000UL;
306	p->type = FP_NORMAL;
307
308	if (u.ll & 0x1fffffffULL)
309	{
310	switch (regs.r_fpsw & FPSWBITS_RM)
311	{
312	case FPRM_NEAREST:
313	if (u.ll & 0x10000000ULL)
314	p->mant ++;
315	break;
316	case FPRM_ZERO:
317	break;
318	case FPRM_PINF:
319	if (sign == 1)
320	p->mant ++;
321	break;
322	case FPRM_NINF:
323	if (sign == -1)
324	p->mant ++;
325	break;
326	}
327	FP_RAISE (X);
328	}
329
330	}
331
332	typedef enum {
333	eNR, /* Use the normal result. */
334	ePZ, eNZ, /* +- zero */
335	eSZ, /* signed zero - XOR signs of ops together. */
336	eRZ, /* +- zero depending on rounding mode. */
337	ePI, eNI, /* +- Infinity */
338	eSI, /* signed infinity - XOR signs of ops together. */
339	eQN, eSN, /* Quiet/Signalling NANs */
340	eIn, /* Invalid. */
341	eUn, /* Unimplemented. */
342	eDZ, /* Divide-by-zero. */
343	eLT, /* less than */
344	eGT, /* greater than */
345	eEQ, /* equal to */
346	} FP_ExceptionCases;
347
348	#if defined DEBUG0
349	static const char *ex_names[] = {
350	"NR", "PZ", "NZ", "SZ", "RZ", "PI", "NI", "SI", "QN", "SN", "IN", "Un", "DZ", "LT", "GT", "EQ"
351	};
352	#endif
353
354	/* This checks for all exceptional cases (not all FP exceptions) and
355	returns TRUE if it is providing the result in *c. If it returns
356	FALSE, the caller should do the "normal" operation. */
357	int
358	check_exceptions (FP_Parts a, FP_Parts b, fp_t *c,
359	FP_ExceptionCases ex_tab[5][5],
360	FP_ExceptionCases *case_ret)
361	{
362	FP_ExceptionCases fpec;
363
364	if (a->type == FP_SNAN
365	\|\| b->type == FP_SNAN)
366	fpec = eIn;
367	else if (a->type == FP_QNAN
368	\|\| b->type == FP_QNAN)
369	fpec = eQN;
370	else if (a->type == FP_DENORMAL
371	\|\| b->type == FP_DENORMAL)
372	fpec = eUn;
373	else
374	fpec = ex_tab[(int)(a->type)][(int)(b->type)];
375
376	/printf("%s %s -> %s\n", fpt_names[(int)(a->type)], fpt_names[(int)(b->type)], ex_names[(int)(fpec)]);/
377
378	if (case_ret)
379	*case_ret = fpec;
380
381	switch (fpec)
382	{
383	case eNR: /* Use the normal result. */
384	return 0;
385
386	case ePZ: /* + zero */
387	*c = 0x00000000;
388	return 1;
389
390	case eNZ: /* - zero */
391	*c = 0x80000000;
392	return 1;
393
394	case eSZ: /* signed zero */
395	*c = (a->sign == b->sign) ? PLUS_ZERO : MINUS_ZERO;
396	return 1;
397
398	case eRZ: /* +- zero depending on rounding mode. */
399	if ((regs.r_fpsw & FPSWBITS_RM) == FPRM_NINF)
400	*c = 0x80000000;
401	else
402	*c = 0x00000000;
403	return 1;
404
405	case ePI: /* + Infinity */
406	*c = 0x7F800000;
407	return 1;
408
409	case eNI: /* - Infinity */
410	*c = 0xFF800000;
411	return 1;
412
413	case eSI: /* sign Infinity */
414	*c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF;
415	return 1;
416
417	case eQN: /* Quiet NANs */
418	if (a->type == FP_QNAN)
419	*c = a->orig_value;
420	else
421	*c = b->orig_value;
422	return 1;
423
424	case eSN: /* Signalling NANs */
425	if (a->type == FP_SNAN)
426	*c = a->orig_value;
427	else
428	*c = b->orig_value;
429	FP_RAISE (V);
430	return 1;
431
432	case eIn: /* Invalid. */
433	FP_RAISE (V);
434	if (a->type == FP_SNAN)
435	*c = a->orig_value \| 0x00400000;
7881f69e	436	else if (b->type == FP_SNAN)
4f8d4a38 DD	437	*c = b->orig_value \| 0x00400000;
	438	else
	439	*c = 0x7fc00000;
	440	return 1;
	441
	442	case eUn: /* Unimplemented. */
	443	FP_RAISE (E);
	444	return 1;
	445
	446	case eDZ: /* Division-by-zero. */
	447	*c = (a->sign == b->sign) ? PLUS_INF : MINUS_INF;
	448	FP_RAISE (Z);
	449	return 1;
	450
	451	default:
	452	return 0;
	453	}
	454	}
	455
	456	#define CHECK_EXCEPTIONS(FPPa, FPPb, fpc, ex_tab) \
	457	if (check_exceptions (&FPPa, &FPPb, &fpc, ex_tab, 0)) \
	458	return fpc;
	459
	460	/* For each operation, we have two tables of how nonnormal cases are
	461	handled. The DN=0 case is first, followed by the DN=1 case, with
	462	each table using the following layout: */
	463
	464	static FP_ExceptionCases ex_add_tab[5][5] = {
	465	/* N +0 -0 +In -In */
	466	{ eNR, eNR, eNR, ePI, eNI }, /* Normal */
	467	{ eNR, ePZ, eRZ, ePI, eNI }, /* +0 */
	468	{ eNR, eRZ, eNZ, ePI, eNI }, /* -0 */
	469	{ ePI, ePI, ePI, ePI, eIn }, /* +Inf */
	470	{ eNI, eNI, eNI, eIn, eNI }, /* -Inf */
	471	};
	472
	473	fp_t
	474	rxfp_add (fp_t fa, fp_t fb)
	475	{
	476	FP_Parts a, b, c;
	477	fp_t rv;
	478	double da, db;
	479
	480	fp_explode (fa, &a);
	481	fp_explode (fb, &b);
	482	CHECK_EXCEPTIONS (a, b, rv, ex_add_tab);
	483
	484	da = fp_to_double (&a);
	485	db = fp_to_double (&b);
	486	tprintf("%g + %g = %g\n", da, db, da+db);
	487
	488	double_to_fp (da+db, &c);
	489	rv = fp_implode (&c);
	490	return rv;
	491	}
	492
	493	static FP_ExceptionCases ex_sub_tab[5][5] = {
	494	/* N +0 -0 +In -In */
	495	{ eNR, eNR, eNR, eNI, ePI }, /* Normal */
	496	{ eNR, eRZ, ePZ, eNI, ePI }, /* +0 */
	497	{ eNR, eNZ, eRZ, eNI, ePI }, /* -0 */
	498	{ ePI, ePI, ePI, eIn, ePI }, /* +Inf */
	499	{ eNI, eNI, eNI, eNI, eIn }, /* -Inf */
	500	};
501
502	fp_t
503	rxfp_sub (fp_t fa, fp_t fb)
504	{
505	FP_Parts a, b, c;
506	fp_t rv;
507	double da, db;
508
509	fp_explode (fa, &a);
510	fp_explode (fb, &b);
511	CHECK_EXCEPTIONS (a, b, rv, ex_sub_tab);
512
513	da = fp_to_double (&a);
514	db = fp_to_double (&b);
515	tprintf("%g - %g = %g\n", da, db, da-db);
516
517	double_to_fp (da-db, &c);
518	rv = fp_implode (&c);
519
520	return rv;
521	}
522
523	static FP_ExceptionCases ex_mul_tab[5][5] = {
524	/* N +0 -0 +In -In */
525	{ eNR, eNR, eNR, eSI, eSI }, /* Normal */
526	{ eNR, ePZ, eNZ, eIn, eIn }, /* +0 */
527	{ eNR, eNZ, ePZ, eIn, eIn }, /* -0 */
528	{ eSI, eIn, eIn, ePI, eNI }, /* +Inf */
529	{ eSI, eIn, eIn, eNI, ePI }, /* -Inf */
530	};
531
532	fp_t
533	rxfp_mul (fp_t fa, fp_t fb)
534	{
535	FP_Parts a, b, c;
536	fp_t rv;
537	double da, db;
538
539	fp_explode (fa, &a);
540	fp_explode (fb, &b);
541	CHECK_EXCEPTIONS (a, b, rv, ex_mul_tab);
542
543	da = fp_to_double (&a);
544	db = fp_to_double (&b);
545	tprintf("%g x %g = %g\n", da, db, da*db);
546
547	double_to_fp (da*db, &c);
548	rv = fp_implode (&c);
549
550	return rv;
551	}
552
553	static FP_ExceptionCases ex_div_tab[5][5] = {
554	/* N +0 -0 +In -In */
555	{ eNR, eDZ, eDZ, eSZ, eSZ }, /* Normal */
556	{ eSZ, eIn, eIn, ePZ, eNZ }, /* +0 */
557	{ eSZ, eIn, eIn, eNZ, ePZ }, /* -0 */
558	{ eSI, ePI, eNI, eIn, eIn }, /* +Inf */
559	{ eSI, eNI, ePI, eIn, eIn }, /* -Inf */
560	};
561
562	fp_t
563	rxfp_div (fp_t fa, fp_t fb)
564	{
565	FP_Parts a, b, c;
566	fp_t rv;
567	double da, db;
568
569	fp_explode (fa, &a);
570	fp_explode (fb, &b);
571	CHECK_EXCEPTIONS (a, b, rv, ex_div_tab);
572
573	da = fp_to_double (&a);
574	db = fp_to_double (&b);
575	tprintf("%g / %g = %g\n", da, db, da/db);
576
577	double_to_fp (da/db, &c);
578	rv = fp_implode (&c);
579
580	return rv;
581	}
582
583	static FP_ExceptionCases ex_cmp_tab[5][5] = {
584	/* N +0 -0 +In -In */
585	{ eNR, eNR, eNR, eLT, eGT }, /* Normal */
586	{ eNR, eEQ, eEQ, eLT, eGT }, /* +0 */
587	{ eNR, eEQ, eEQ, eLT, eGT }, /* -0 */
588	{ eGT, eGT, eGT, eEQ, eGT }, /* +Inf */
589	{ eLT, eLT, eLT, eLT, eEQ }, /* -Inf */
590	};
591
592	void
593	rxfp_cmp (fp_t fa, fp_t fb)
594	{
595	FP_Parts a, b;
596	fp_t c;
597	FP_ExceptionCases reason;
598	int flags = 0;
599	double da, db;
600
601	fp_explode (fa, &a);
602	fp_explode (fb, &b);
603
604	if (check_exceptions (&a, &b, &c, ex_cmp_tab, &reason))
605	{
606	if (reason == eQN)
607	{
608	/* Special case - incomparable. */
609	set_flags (FLAGBIT_Z \| FLAGBIT_S \| FLAGBIT_O, FLAGBIT_O);
610	return;
611	}
612	return;
613	}
614
615	switch (reason)
616	{
617	case eEQ:
618	flags = FLAGBIT_Z;
619	break;
620	case eLT:
621	flags = FLAGBIT_S;
622	break;
623	case eGT:
624	flags = 0;
625	break;
626	case eNR:
627	da = fp_to_double (&a);
628	db = fp_to_double (&b);
629	tprintf("fcmp: %g cmp %g\n", da, db);
630	if (da < db)
631	flags = FLAGBIT_S;
632	else if (da == db)
633	flags = FLAGBIT_Z;
634	else
635	flags = 0;
636	break;
637	default:
638	abort();
639	}
640
641	set_flags (FLAGBIT_Z \| FLAGBIT_S \| FLAGBIT_O, flags);
642	}
643
644	long
645	rxfp_ftoi (fp_t fa, int round_mode)
646	{
647	FP_Parts a;
648	fp_t rv;
649	int sign;
650	int whole_bits, frac_bits;
651
652	fp_explode (fa, &a);
653	sign = fa & 0x80000000UL;
654
655	switch (a.type)
656	{
657	case FP_NORMAL:
658	break;
659	case FP_PZERO:
660	case FP_NZERO:
661	return 0;
662	case FP_PINFINITY:
663	FP_RAISE (V);
664	return 0x7fffffffL;
665	case FP_NINFINITY:
666	FP_RAISE (V);
667	return 0x80000000L;
668	case FP_DENORMAL:
669	FP_RAISE (E);
670	return 0;
671	case FP_QNAN:
672	case FP_SNAN:
673	FP_RAISE (V);
674	return sign ? 0x80000000U : 0x7fffffff;
675	}
676
677	if (a.exp >= 31)
678	{
679	FP_RAISE (V);
680	return sign ? 0x80000000U : 0x7fffffff;
681	}
682
683	a.exp -= 23;
684
685	if (a.exp <= -25)
686	{
687	/* Less than 0.49999 */
688	frac_bits = a.mant;
689	whole_bits = 0;
690	}
691	else if (a.exp < 0)
692	{
693	frac_bits = a.mant << (32 + a.exp);
694	whole_bits = a.mant >> (-a.exp);
695	}
696	else
697	{
698	frac_bits = 0;
699	whole_bits = a.mant << a.exp;
700	}
701
702	if (frac_bits)
703	{
704	switch (round_mode & 3)
705	{
706	case FPRM_NEAREST:
707	if (frac_bits & 0x80000000UL)
708	whole_bits ++;
709	break;
710	case FPRM_ZERO:
711	break;
712	case FPRM_PINF:
713	if (!sign)
714	whole_bits ++;
715	break;
716	case FPRM_NINF:
717	if (sign)
718	whole_bits ++;
719	break;
720	}
721	}
722
723	rv = sign ? -whole_bits : whole_bits;
724
725	return rv;
726	}
727
728	fp_t
729	rxfp_itof (long fa, int round_mode)
730	{
731	fp_t rv;
732	int sign = 0;
733	unsigned int frac_bits;
734	volatile unsigned int whole_bits;
735	FP_Parts a;
736
737	if (fa == 0)
738	return PLUS_ZERO;
739
740	if (fa < 0)
741	{
742	fa = -fa;
743	sign = 1;
744	a.sign = -1;
745	}
746	else
747	a.sign = 1;
748
749	whole_bits = fa;
750	a.exp = 31;
751
752	while (! (whole_bits & 0x80000000UL))
753	{
754	a.exp --;
755	whole_bits <<= 1;
756	}
757	frac_bits = whole_bits & 0xff;
758	whole_bits = whole_bits >> 8;
759
760	if (frac_bits)
761	{
762	/* We must round */
763	switch (round_mode & 3)
764	{
765	case FPRM_NEAREST:
766	if (frac_bits & 0x80)
767	whole_bits ++;
768	break;
769	case FPRM_ZERO:
770	break;
771	case FPRM_PINF:
772	if (!sign)
773	whole_bits ++;
774	break;
775	case FPRM_NINF:
776	if (sign)
777	whole_bits ++;
778	break;
779	}
780	}
781
782	a.mant = whole_bits;
783	if (whole_bits & 0xff000000UL)
784	{
785	a.mant >>= 1;
786	a.exp ++;
787	}
788
789	rv = fp_implode (&a);
790	return rv;
791	}
792