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

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

Copyright (C) 2008-2021 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/>.  */

/* This must come before any other includes.  */
#include "defs.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.  */
static 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 = {0};

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