1 /* fpu.c --- FPU emulator for stand-alone RX simulator.
2
3 Copyright (C) 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
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
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)
62 static void
fp_raise(int mask)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
fp_explode(fp_t f,FP_Parts * p)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
fp_implode(FP_Parts * p)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
fp_to_double(FP_Parts * p)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
double_to_fp(double d,FP_Parts * p)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
check_exceptions(FP_Parts * a,FP_Parts * b,fp_t * c,FP_ExceptionCases ex_tab[5][5],FP_ExceptionCases * case_ret)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;
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
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
rxfp_add(fp_t fa,fp_t fb)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
rxfp_sub(fp_t fa,fp_t fb)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
rxfp_mul(fp_t fa,fp_t fb)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
rxfp_div(fp_t fa,fp_t fb)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
rxfp_cmp(fp_t fa,fp_t fb)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
rxfp_ftoi(fp_t fa,int round_mode)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
rxfp_itof(long fa,int round_mode)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
793