1/* fpu.c --- FPU emulator for stand-alone RX simulator.
2
3Copyright (C) 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by Red Hat, Inc.
5
6This file is part of the GNU simulators.
7
8This program is free software; you can redistribute it and/or modify
9it under the terms of the GNU General Public License as published by
10the Free Software Foundation; either version 3 of the License, or
11(at your option) any later version.
12
13This program is distributed in the hope that it will be useful,
14but WITHOUT ANY WARRANTY; without even the implied warranty of
15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16GNU General Public License for more details.
17
18You should have received a copy of the GNU General Public License
19along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
20
21#include "config.h"
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)
62static void
63fp_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.  */
79typedef 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
91static 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
102typedef 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
110static void
111fp_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
150static fp_t
151fp_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
196typedef union {
197  unsigned long long ll;
198  double d;
199} U_d_ll;
200
201static int checked_format = 0;
202
203/* We assume a double format like this:
204   S[1] E[11] M[52]
205*/
206
207static double
208fp_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
236static void
237double_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
332typedef 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
349static 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.  */
357int
358check_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;
436      else if  (a->type == FP_SNAN)
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
464static 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
473fp_t
474rxfp_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
493static 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
502fp_t
503rxfp_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
523static 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
532fp_t
533rxfp_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
553static 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
562fp_t
563rxfp_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
583static 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
592void
593rxfp_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
644long
645rxfp_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
728fp_t
729rxfp_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
793