* sim-fpu.h (enum sim_fpu_class): Add sim_fpu_class_denorm.
[deliverable/binutils-gdb.git] / sim / common / sim-fpu.c
1 /* This is a software floating point library which can be used instead
2 of the floating point routines in libgcc1.c for targets without
3 hardware floating point. */
4
5 /* Copyright (C) 1994,1997 Free Software Foundation, Inc.
6
7 This file is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 2, or (at your option) any
10 later version.
11
12 In addition to the permissions in the GNU General Public License, the
13 Free Software Foundation gives you unlimited permission to link the
14 compiled version of this file with other programs, and to distribute
15 those programs without any restriction coming from the use of this
16 file. (The General Public License restrictions do apply in other
17 respects; for example, they cover modification of the file, and
18 distribution when not linked into another program.)
19
20 This file is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23 General Public License for more details.
24
25 You should have received a copy of the GNU General Public License
26 along with this program; see the file COPYING. If not, write to
27 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
28
29 /* As a special exception, if you link this library with other files,
30 some of which are compiled with GCC, to produce an executable,
31 this library does not by itself cause the resulting executable
32 to be covered by the GNU General Public License.
33 This exception does not however invalidate any other reasons why
34 the executable file might be covered by the GNU General Public License. */
35
36 /* This implements IEEE 754 format arithmetic, but does not provide a
37 mechanism for setting the rounding mode, or for generating or handling
38 exceptions.
39
40 The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
41 Wilson, all of Cygnus Support. */
42
43
44 #ifndef SIM_FPU_C
45 #define SIM_FPU_C
46
47 #include "sim-basics.h"
48 #include "sim-fpu.h"
49
50 #include "sim-io.h"
51 #include "sim-assert.h"
52
53
54 /* Debugging support. */
55
56 static void
57 print_bits (unsigned64 x,
58 int msbit,
59 sim_fpu_print_func print,
60 void *arg)
61 {
62 unsigned64 bit = LSBIT64 (msbit);
63 int i = 4;
64 while (bit)
65 {
66 if (i == 0)
67 print (arg, ",");
68 if ((x & bit))
69 print (arg, "1");
70 else
71 print (arg, "0");
72 bit >>= 1;
73 i = (i + 1) % 4;
74 }
75 }
76
77
78
79 /* Quick and dirty conversion between a host double and host 64bit int */
80
81 typedef union {
82 double d;
83 unsigned64 i;
84 } sim_fpu_map;
85
86
87 /* A packed IEEE floating point number.
88
89 Form is <SIGN:1><BIASEDEXP:NR_EXPBITS><FRAC:NR_FRACBITS> for both
90 32 and 64 bit numbers. This number is interpreted as:
91
92 Normalized (0 < BIASEDEXP && BIASEDEXP < EXPMAX):
93 (sign ? '-' : '+') 1.<FRAC> x 2 ^ (BIASEDEXP - EXPBIAS)
94
95 Denormalized (0 == BIASEDEXP && FRAC != 0):
96 (sign ? "-" : "+") 0.<FRAC> x 2 ^ (- EXPBIAS)
97
98 Zero (0 == BIASEDEXP && FRAC == 0):
99 (sign ? "-" : "+") 0.0
100
101 Infinity (BIASEDEXP == EXPMAX && FRAC == 0):
102 (sign ? "-" : "+") "infinity"
103
104 SignalingNaN (BIASEDEXP == EXPMAX && FRAC > 0 && FRAC < QUIET_NAN):
105 SNaN.FRAC
106
107 QuietNaN (BIASEDEXP == EXPMAX && FRAC > 0 && FRAC > QUIET_NAN):
108 QNaN.FRAC
109
110 */
111
112 #define NR_EXPBITS (is_double ? 11 : 8)
113 #define NR_FRACBITS (is_double ? 52 : 23)
114 #define SIGNBIT (is_double ? MSBIT64 (0) : MSBIT64 (32))
115
116 #define EXPMAX ((unsigned) (is_double ? 2047 : 255))
117 #define EXPBIAS (is_double ? 1023 : 127)
118
119 #define QUIET_NAN LSBIT64 (NR_FRACBITS - 1)
120
121
122
123 /* An unpacked floating point number.
124
125 When unpacked, the fraction of both a 32 and 64 bit floating point
126 number is stored using the same format:
127
128 64 bit - <IMPLICIT_1:1><FRACBITS:52><GUARDS:8><PAD:00>
129 32 bit - <IMPLICIT_1:1><FRACBITS:23><GUARDS:7><PAD:30> */
130
131 #define NR_PAD (is_double ? 0 : 30)
132 #define PADMASK (is_double ? 0 : LSMASK64 (29, 0))
133 #define NR_GUARDS ((is_double ? 8 : 7 ) + NR_PAD)
134 #define GUARDMASK LSMASK64 (NR_GUARDS - 1, 0)
135 #define GUARDMSB LSBIT64 (NR_GUARDS - 1)
136 #define GUARDLSB LSBIT64 (NR_PAD)
137 #define GUARDROUND LSMASK64 (NR_GUARDS - 2, 0)
138
139 #define NR_FRAC_GUARD (60)
140 #define IMPLICIT_1 LSBIT64 (NR_FRAC_GUARD)
141 #define IMPLICIT_2 LSBIT64 (NR_FRAC_GUARD + 1)
142 #define IMPLICIT_4 LSBIT64 (NR_FRAC_GUARD + 2)
143 #define NR_SPARE 2
144
145 #define FRAC32MASK LSMASK64 (63, NR_FRAC_GUARD - 32 + 1)
146
147 #define NORMAL_EXPMIN (-(EXPBIAS)+1)
148 #define NORMAL_EXPMAX (EXPBIAS)
149
150
151 /* Integer constants */
152
153 #define MAX_INT32 ((signed64) LSMASK64 (30, 0))
154 #define MAX_UINT32 LSMASK64 (31, 0)
155 #define MIN_INT32 ((signed64) LSMASK64 (63, 31))
156
157 #define MAX_INT64 ((signed64) LSMASK64 (62, 0))
158 #define MAX_UINT64 LSMASK64 (63, 0)
159 #define MIN_INT64 ((signed64) LSMASK64 (63, 63))
160
161 #define MAX_INT (is_64bit ? MAX_INT64 : MAX_INT32)
162 #define MIN_INT (is_64bit ? MIN_INT64 : MIN_INT32)
163 #define MAX_UINT (is_64bit ? MAX_UINT64 : MAX_UINT32)
164 #define NR_INTBITS (is_64bit ? 64 : 32)
165
166 /* Squeese an unpacked sim_fpu struct into a 32/64 bit integer */
167 STATIC_INLINE_SIM_FPU (unsigned64)
168 pack_fpu (const sim_fpu *src,
169 int is_double)
170 {
171 int sign;
172 unsigned64 exp;
173 unsigned64 fraction;
174 unsigned64 packed;
175
176 switch (src->class)
177 {
178 /* create a NaN */
179 case sim_fpu_class_qnan:
180 sign = src->sign;
181 exp = EXPMAX;
182 /* force fraction to correct class */
183 fraction = src->fraction;
184 fraction >>= NR_GUARDS;
185 fraction |= QUIET_NAN;
186 break;
187 case sim_fpu_class_snan:
188 sign = src->sign;
189 exp = EXPMAX;
190 /* force fraction to correct class */
191 fraction = src->fraction;
192 fraction >>= NR_GUARDS;
193 fraction &= ~QUIET_NAN;
194 break;
195 case sim_fpu_class_infinity:
196 sign = src->sign;
197 exp = EXPMAX;
198 fraction = 0;
199 break;
200 case sim_fpu_class_zero:
201 sign = src->sign;
202 exp = 0;
203 fraction = 0;
204 break;
205 case sim_fpu_class_number:
206 case sim_fpu_class_denorm:
207 ASSERT (src->fraction >= IMPLICIT_1);
208 ASSERT (src->fraction < IMPLICIT_2);
209 if (src->normal_exp < NORMAL_EXPMIN)
210 {
211 /* This number's exponent is too low to fit into the bits
212 available in the number We'll denormalize the number by
213 storing zero in the exponent and shift the fraction to
214 the right to make up for it. */
215 int nr_shift = NORMAL_EXPMIN - src->normal_exp;
216 if (nr_shift > NR_FRACBITS)
217 {
218 /* underflow, just make the number zero */
219 sign = src->sign;
220 exp = 0;
221 fraction = 0;
222 }
223 else
224 {
225 sign = src->sign;
226 exp = 0;
227 /* Shift by the value */
228 fraction = src->fraction;
229 fraction >>= NR_GUARDS;
230 fraction >>= nr_shift;
231 }
232 }
233 else if (src->normal_exp > NORMAL_EXPMAX)
234 {
235 /* Infinity */
236 sign = src->sign;
237 exp = EXPMAX;
238 fraction = 0;
239 }
240 else
241 {
242 exp = (src->normal_exp + EXPBIAS);
243 sign = src->sign;
244 fraction = src->fraction;
245 /* FIXME: Need to round according to WITH_SIM_FPU_ROUNDING
246 or some such */
247 /* Round to nearest: If the guard bits are the all zero, but
248 the first, then we're half way between two numbers,
249 choose the one which makes the lsb of the answer 0. */
250 if ((fraction & GUARDMASK) == GUARDMSB)
251 {
252 if ((fraction & (GUARDMSB << 1)))
253 fraction += (GUARDMSB << 1);
254 }
255 else
256 {
257 /* Add a one to the guards to force round to nearest */
258 fraction += GUARDROUND;
259 }
260 if ((fraction & IMPLICIT_2)) /* rounding resulted in carry */
261 {
262 exp += 1;
263 fraction >>= 1;
264 }
265 fraction >>= NR_GUARDS;
266 /* When exp == EXPMAX (overflow from carry) fraction must
267 have been made zero */
268 ASSERT ((exp == EXPMAX) <= ((fraction & ~IMPLICIT_1) == 0));
269 }
270 break;
271 default:
272 abort ();
273 }
274
275 packed = ((sign ? SIGNBIT : 0)
276 | (exp << NR_FRACBITS)
277 | LSMASKED64 (fraction, NR_FRACBITS - 1, 0));
278
279 /* trace operation */
280 #if 0
281 if (is_double)
282 {
283 }
284 else
285 {
286 printf ("pack_fpu: ");
287 printf ("-> %c%0lX.%06lX\n",
288 LSMASKED32 (packed, 31, 31) ? '8' : '0',
289 (long) LSEXTRACTED32 (packed, 30, 23),
290 (long) LSEXTRACTED32 (packed, 23 - 1, 0));
291 }
292 #endif
293
294 return packed;
295 }
296
297
298 /* Unpack a 32/64 bit integer into a sim_fpu structure */
299 STATIC_INLINE_SIM_FPU (void)
300 unpack_fpu (sim_fpu *dst, unsigned64 packed, int is_double)
301 {
302 unsigned64 fraction = LSMASKED64 (packed, NR_FRACBITS - 1, 0);
303 unsigned exp = LSEXTRACTED64 (packed, NR_EXPBITS + NR_FRACBITS - 1, NR_FRACBITS);
304 int sign = (packed & SIGNBIT) != 0;
305
306 if (exp == 0)
307 {
308 /* Hmm. Looks like 0 */
309 if (fraction == 0)
310 {
311 /* tastes like zero */
312 dst->class = sim_fpu_class_zero;
313 dst->sign = sign;
314 }
315 else
316 {
317 /* Zero exponent with non zero fraction - it's denormalized,
318 so there isn't a leading implicit one - we'll shift it so
319 it gets one. */
320 dst->normal_exp = exp - EXPBIAS + 1;
321 dst->class = sim_fpu_class_denorm;
322 dst->sign = sign;
323 fraction <<= NR_GUARDS;
324 while (fraction < IMPLICIT_1)
325 {
326 fraction <<= 1;
327 dst->normal_exp--;
328 }
329 dst->fraction = fraction;
330 }
331 }
332 else if (exp == EXPMAX)
333 {
334 /* Huge exponent*/
335 if (fraction == 0)
336 {
337 /* Attached to a zero fraction - means infinity */
338 dst->class = sim_fpu_class_infinity;
339 dst->sign = sign;
340 /* dst->normal_exp = EXPBIAS; */
341 /* dst->fraction = 0; */
342 }
343 else
344 {
345 /* Non zero fraction, means NaN */
346 dst->sign = sign;
347 dst->fraction = (fraction << NR_GUARDS);
348 if (fraction >= QUIET_NAN)
349 dst->class = sim_fpu_class_qnan;
350 else
351 dst->class = sim_fpu_class_snan;
352 }
353 }
354 else
355 {
356 /* Nothing strange about this number */
357 dst->class = sim_fpu_class_number;
358 dst->sign = sign;
359 dst->fraction = ((fraction << NR_GUARDS) | IMPLICIT_1);
360 dst->normal_exp = exp - EXPBIAS;
361 }
362
363 /* trace operation */
364 #if 0
365 if (is_double)
366 {
367 }
368 else
369 {
370 printf ("unpack_fpu: %c%02lX.%06lX ->\n",
371 LSMASKED32 (packed, 31, 31) ? '8' : '0',
372 (long) LSEXTRACTED32 (packed, 30, 23),
373 (long) LSEXTRACTED32 (packed, 23 - 1, 0));
374 }
375 #endif
376
377 /* sanity checks */
378 {
379 sim_fpu_map val;
380 val.i = pack_fpu (dst, 1);
381 if (is_double)
382 {
383 ASSERT (val.i == packed);
384 }
385 else
386 {
387 unsigned32 val = pack_fpu (dst, 0);
388 unsigned32 org = packed;
389 ASSERT (val == org);
390 }
391 }
392 }
393
394
395 /* Convert a floating point into an integer */
396 STATIC_INLINE_SIM_FPU (int)
397 fpu2i (signed64 *i,
398 const sim_fpu *s,
399 int is_64bit,
400 sim_fpu_round round)
401 {
402 unsigned64 tmp;
403 int shift;
404 int status = 0;
405 if (sim_fpu_is_zero (s))
406 {
407 *i = 0;
408 return 0;
409 }
410 if (sim_fpu_is_snan (s))
411 {
412 *i = MIN_INT; /* FIXME */
413 return sim_fpu_status_invalid_cvi;
414 }
415 if (sim_fpu_is_qnan (s))
416 {
417 *i = MIN_INT; /* FIXME */
418 return sim_fpu_status_invalid_cvi;
419 }
420 /* map infinity onto MAX_INT... */
421 if (sim_fpu_is_infinity (s))
422 {
423 *i = s->sign ? MIN_INT : MAX_INT;
424 return sim_fpu_status_invalid_cvi;
425 }
426 /* it is a number, but a small one */
427 if (s->normal_exp < 0)
428 {
429 *i = 0;
430 return sim_fpu_status_inexact;
431 }
432 /* Is the floating point MIN_INT or just close? */
433 if (s->sign && s->normal_exp == (NR_INTBITS - 1))
434 {
435 *i = MIN_INT;
436 ASSERT (s->fraction >= IMPLICIT_1);
437 if (s->fraction == IMPLICIT_1)
438 return 0; /* exact */
439 if (is_64bit) /* can't round */
440 return sim_fpu_status_invalid_cvi; /* must be overflow */
441 /* For a 32bit with MAX_INT, rounding is possible */
442 switch (round)
443 {
444 case sim_fpu_round_default:
445 abort ();
446 case sim_fpu_round_zero:
447 if ((s->fraction & FRAC32MASK) != IMPLICIT_1)
448 return sim_fpu_status_invalid_cvi;
449 else
450 return sim_fpu_status_inexact;
451 break;
452 case sim_fpu_round_near:
453 {
454 if ((s->fraction & FRAC32MASK) != IMPLICIT_1)
455 return sim_fpu_status_invalid_cvi;
456 else if ((s->fraction & !FRAC32MASK) >= (~FRAC32MASK >> 1))
457 return sim_fpu_status_invalid_cvi;
458 else
459 return sim_fpu_status_inexact;
460 }
461 case sim_fpu_round_up:
462 if ((s->fraction & FRAC32MASK) == IMPLICIT_1)
463 return sim_fpu_status_inexact;
464 else
465 return sim_fpu_status_invalid_cvi;
466 case sim_fpu_round_down:
467 return sim_fpu_status_invalid_cvi;
468 }
469 }
470 /* Would right shifting result in the FRAC being shifted into
471 (through) the integer's sign bit? */
472 if (s->normal_exp > (NR_INTBITS - 2))
473 {
474 *i = s->sign ? MIN_INT : MAX_INT;
475 return sim_fpu_status_invalid_cvi;
476 }
477 /* normal number shift it into place */
478 tmp = s->fraction;
479 shift = (s->normal_exp - (NR_FRAC_GUARD));
480 if (shift > 0)
481 {
482 tmp <<= shift;
483 }
484 else
485 {
486 shift = -shift;
487 if (tmp & ((SIGNED64 (1) << shift) - 1))
488 status |= sim_fpu_status_inexact;
489 tmp >>= shift;
490 }
491 *i = s->sign ? (-tmp) : (tmp);
492 return status;
493 }
494
495 /* convert an integer into a floating point */
496 STATIC_INLINE_SIM_FPU (int)
497 i2fpu (sim_fpu *f, signed64 i, int is_64bit)
498 {
499 int status = 0;
500 if (i == 0)
501 {
502 f->class = sim_fpu_class_zero;
503 f->sign = 0;
504 }
505 else
506 {
507 f->class = sim_fpu_class_number;
508 f->sign = (i < 0);
509 f->normal_exp = NR_FRAC_GUARD;
510
511 if (f->sign)
512 {
513 /* Special case for minint, since there is no corresponding
514 +ve integer representation for it */
515 if (i == MIN_INT)
516 {
517 f->fraction = IMPLICIT_1;
518 f->normal_exp = NR_INTBITS - 1;
519 }
520 else
521 f->fraction = (-i);
522 }
523 else
524 f->fraction = i;
525
526 if (f->fraction >= IMPLICIT_2)
527 {
528 do
529 {
530 f->fraction >>= 1;
531 f->normal_exp += 1;
532 }
533 while (f->fraction >= IMPLICIT_2);
534 }
535 else if (f->fraction < IMPLICIT_1)
536 {
537 do
538 {
539 f->fraction <<= 1;
540 f->normal_exp -= 1;
541 }
542 while (f->fraction < IMPLICIT_1);
543 }
544 }
545
546 /* trace operation */
547 #if 0
548 {
549 printf ("i2fpu: 0x%08lX ->\n", (long) i);
550 }
551 #endif
552
553 /* sanity check */
554 {
555 signed64 val;
556 fpu2i (&val, f, is_64bit, sim_fpu_round_zero);
557 if (i >= MIN_INT32 && i <= MAX_INT32)
558 {
559 ASSERT (val == i);
560 }
561 }
562
563 return status;
564 }
565
566
567 /* Convert a floating point into an integer */
568 STATIC_INLINE_SIM_FPU (int)
569 fpu2u (unsigned64 *u, const sim_fpu *s, int is_64bit)
570 {
571 const int is_double = 1;
572 unsigned64 tmp;
573 int shift;
574 if (sim_fpu_is_zero (s))
575 {
576 *u = 0;
577 return 0;
578 }
579 if (sim_fpu_is_nan (s))
580 {
581 *u = 0;
582 return 0;
583 }
584 /* it is a negative number */
585 if (s->sign)
586 {
587 *u = 0;
588 return 0;
589 }
590 /* get reasonable MAX_USI_INT... */
591 if (sim_fpu_is_infinity (s))
592 {
593 *u = MAX_UINT;
594 return 0;
595 }
596 /* it is a number, but a small one */
597 if (s->normal_exp < 0)
598 {
599 *u = 0;
600 return 0;
601 }
602 /* overflow */
603 if (s->normal_exp > (NR_INTBITS - 1))
604 {
605 *u = MAX_UINT;
606 return 0;
607 }
608 /* normal number */
609 tmp = (s->fraction & ~PADMASK);
610 shift = (s->normal_exp - (NR_FRACBITS + NR_GUARDS));
611 if (shift > 0)
612 {
613 tmp <<= shift;
614 }
615 else
616 {
617 shift = -shift;
618 tmp >>= shift;
619 }
620 *u = tmp;
621 return 0;
622 }
623
624 /* Convert an unsigned integer into a floating point */
625 STATIC_INLINE_SIM_FPU (int)
626 u2fpu (sim_fpu *f, unsigned64 u, int is_64bit)
627 {
628 if (u == 0)
629 {
630 f->class = sim_fpu_class_zero;
631 f->sign = 0;
632 }
633 else
634 {
635 f->class = sim_fpu_class_number;
636 f->sign = 0;
637 f->normal_exp = NR_FRAC_GUARD;
638 f->fraction = u;
639
640 while (f->fraction < IMPLICIT_1)
641 {
642 f->fraction <<= 1;
643 f->normal_exp -= 1;
644 }
645 }
646 return 0;
647 }
648
649
650 /* register <-> sim_fpu */
651
652 INLINE_SIM_FPU (void)
653 sim_fpu_32to (sim_fpu *f, unsigned32 s)
654 {
655 unpack_fpu (f, s, 0);
656 }
657
658
659 INLINE_SIM_FPU (void)
660 sim_fpu_232to (sim_fpu *f, unsigned32 h, unsigned32 l)
661 {
662 unsigned64 s = h;
663 s = (s << 32) | l;
664 unpack_fpu (f, s, 1);
665 }
666
667
668 INLINE_SIM_FPU (void)
669 sim_fpu_64to (sim_fpu *f, unsigned64 s)
670 {
671 unpack_fpu (f, s, 1);
672 }
673
674
675 INLINE_SIM_FPU (void)
676 sim_fpu_to32 (unsigned32 *s,
677 const sim_fpu *f)
678 {
679 *s = pack_fpu (f, 0);
680 }
681
682
683 INLINE_SIM_FPU (void)
684 sim_fpu_to232 (unsigned32 *h, unsigned32 *l,
685 const sim_fpu *f)
686 {
687 unsigned64 s = pack_fpu (f, 1);
688 *l = s;
689 *h = (s >> 32);
690 }
691
692
693 INLINE_SIM_FPU (void)
694 sim_fpu_to64 (unsigned64 *u,
695 const sim_fpu *f)
696 {
697 *u = pack_fpu (f, 1);
698 }
699
700
701 /* Rounding */
702
703 STATIC_INLINE_SIM_FPU (int)
704 do_normal_overflow (sim_fpu *f,
705 int is_double,
706 sim_fpu_round round)
707 {
708 switch (round)
709 {
710 case sim_fpu_round_default:
711 return 0;
712 case sim_fpu_round_near:
713 f->class = sim_fpu_class_infinity;
714 break;
715 case sim_fpu_round_up:
716 if (!f->sign)
717 f->class = sim_fpu_class_infinity;
718 break;
719 case sim_fpu_round_down:
720 if (f->sign)
721 f->class = sim_fpu_class_infinity;
722 break;
723 case sim_fpu_round_zero:
724 break;
725 }
726 f->normal_exp = NORMAL_EXPMAX;
727 f->fraction = LSMASK64 (NR_FRAC_GUARD, NR_GUARDS);
728 return (sim_fpu_status_overflow | sim_fpu_status_inexact);
729 }
730
731 STATIC_INLINE_SIM_FPU (int)
732 do_normal_underflow (sim_fpu *f,
733 int is_double,
734 sim_fpu_round round)
735 {
736 switch (round)
737 {
738 case sim_fpu_round_default:
739 return 0;
740 case sim_fpu_round_near:
741 f->class = sim_fpu_class_zero;
742 break;
743 case sim_fpu_round_up:
744 if (f->sign)
745 f->class = sim_fpu_class_zero;
746 break;
747 case sim_fpu_round_down:
748 if (!f->sign)
749 f->class = sim_fpu_class_zero;
750 break;
751 case sim_fpu_round_zero:
752 f->class = sim_fpu_class_zero;
753 break;
754 }
755 f->normal_exp = NORMAL_EXPMIN - NR_FRACBITS;
756 f->fraction = IMPLICIT_1;
757 return (sim_fpu_status_inexact | sim_fpu_status_underflow);
758 }
759
760
761
762 /* Round a number using NR_GUARDS.
763 Will return the rounded number or F->FRACTION == 0 when underflow */
764
765 STATIC_INLINE_SIM_FPU (int)
766 do_normal_round (sim_fpu *f,
767 int nr_guards,
768 sim_fpu_round round)
769 {
770 unsigned64 guardmask = LSMASK64 (nr_guards - 1, 0);
771 unsigned64 guardmsb = LSBIT64 (nr_guards - 1);
772 unsigned64 fraclsb = guardmsb << 1;
773 if ((f->fraction & guardmask))
774 {
775 int status = sim_fpu_status_inexact;
776 switch (round)
777 {
778 case sim_fpu_round_default:
779 return 0;
780 case sim_fpu_round_near:
781 if ((f->fraction & guardmsb))
782 {
783 if ((f->fraction & fraclsb))
784 {
785 status |= sim_fpu_status_rounded;
786 }
787 else if ((f->fraction & (guardmask >> 1)))
788 {
789 status |= sim_fpu_status_rounded;
790 }
791 }
792 break;
793 case sim_fpu_round_up:
794 if (!f->sign)
795 status |= sim_fpu_status_rounded;
796 break;
797 case sim_fpu_round_down:
798 if (f->sign)
799 status |= sim_fpu_status_rounded;
800 break;
801 case sim_fpu_round_zero:
802 break;
803 }
804 f->fraction &= ~guardmask;
805 /* round if needed, handle resulting overflow */
806 if ((status & sim_fpu_status_rounded))
807 {
808 f->fraction += fraclsb;
809 if ((f->fraction & IMPLICIT_2))
810 {
811 f->fraction >>= 1;
812 f->normal_exp += 1;
813 }
814 }
815 return status;
816 }
817 else
818 return 0;
819 }
820
821
822 STATIC_INLINE_SIM_FPU (int)
823 do_round (sim_fpu *f,
824 int is_double,
825 sim_fpu_round round,
826 sim_fpu_denorm denorm)
827 {
828 switch (f->class)
829 {
830 case sim_fpu_class_qnan:
831 case sim_fpu_class_zero:
832 case sim_fpu_class_infinity:
833 return 0;
834 break;
835 case sim_fpu_class_snan:
836 /* Quieten a SignalingNaN */
837 f->class = sim_fpu_class_qnan;
838 return sim_fpu_status_invalid_snan;
839 break;
840 case sim_fpu_class_number:
841 case sim_fpu_class_denorm:
842 {
843 int status;
844 ASSERT (f->fraction < IMPLICIT_2);
845 ASSERT (f->fraction >= IMPLICIT_1);
846 if (f->normal_exp < NORMAL_EXPMIN)
847 {
848 /* This number's exponent is too low to fit into the bits
849 available in the number. Round off any bits that will be
850 discarded as a result of denormalization. Edge case is
851 the implicit bit shifted to GUARD0 and then rounded
852 up. */
853 int shift = NORMAL_EXPMIN - f->normal_exp;
854 if (shift + NR_GUARDS <= NR_FRAC_GUARD + 1
855 && !(denorm & sim_fpu_denorm_zero))
856 {
857 status = do_normal_round (f, shift + NR_GUARDS, round);
858 if (f->fraction == 0) /* rounding underflowed */
859 {
860 status |= do_normal_underflow (f, is_double, round);
861 }
862 else if (f->normal_exp < NORMAL_EXPMIN) /* still underflow? */
863 {
864 status |= sim_fpu_status_denorm;
865 /* Any loss of precision when denormalizing is
866 underflow. Some processors check for underflow
867 before rounding, some after! */
868 if (status & sim_fpu_status_inexact)
869 status |= sim_fpu_status_underflow;
870 /* Flag that resultant value has been denormalized */
871 f->class = sim_fpu_class_denorm;
872 }
873 else if ((denorm & sim_fpu_denorm_underflow_inexact))
874 {
875 if ((status & sim_fpu_status_inexact))
876 status |= sim_fpu_status_underflow;
877 }
878 }
879 else
880 {
881 status = do_normal_underflow (f, is_double, round);
882 }
883 }
884 else if (f->normal_exp > NORMAL_EXPMAX)
885 {
886 /* Infinity */
887 status = do_normal_overflow (f, is_double, round);
888 }
889 else
890 {
891 status = do_normal_round (f, NR_GUARDS, round);
892 if (f->fraction == 0)
893 /* f->class = sim_fpu_class_zero; */
894 status |= do_normal_underflow (f, is_double, round);
895 else if (f->normal_exp > NORMAL_EXPMAX)
896 /* oops! rounding caused overflow */
897 status |= do_normal_overflow (f, is_double, round);
898 }
899 ASSERT ((f->class == sim_fpu_class_number
900 || f->class == sim_fpu_class_denorm)
901 <= (f->fraction < IMPLICIT_2 && f->fraction >= IMPLICIT_1));
902 return status;
903 }
904 }
905 return 0;
906 }
907
908 INLINE_SIM_FPU (int)
909 sim_fpu_round_32 (sim_fpu *f,
910 sim_fpu_round round,
911 sim_fpu_denorm denorm)
912 {
913 return do_round (f, 0, round, denorm);
914 }
915
916 INLINE_SIM_FPU (int)
917 sim_fpu_round_64 (sim_fpu *f,
918 sim_fpu_round round,
919 sim_fpu_denorm denorm)
920 {
921 return do_round (f, 1, round, denorm);
922 }
923
924
925
926 /* Arithmetic ops */
927
928 INLINE_SIM_FPU (int)
929 sim_fpu_add (sim_fpu *f,
930 const sim_fpu *l,
931 const sim_fpu *r)
932 {
933 if (sim_fpu_is_snan (l))
934 {
935 *f = *l;
936 f->class = sim_fpu_class_qnan;
937 return sim_fpu_status_invalid_snan;
938 }
939 if (sim_fpu_is_snan (r))
940 {
941 *f = *r;
942 f->class = sim_fpu_class_qnan;
943 return sim_fpu_status_invalid_snan;
944 }
945 if (sim_fpu_is_qnan (l))
946 {
947 *f = *l;
948 return 0;
949 }
950 if (sim_fpu_is_qnan (r))
951 {
952 *f = *r;
953 return 0;
954 }
955 if (sim_fpu_is_infinity (l))
956 {
957 if (sim_fpu_is_infinity (r)
958 && l->sign != r->sign)
959 {
960 *f = sim_fpu_qnan;
961 return sim_fpu_status_invalid_isi;
962 }
963 *f = *l;
964 return 0;
965 }
966 if (sim_fpu_is_infinity (r))
967 {
968 *f = *r;
969 return 0;
970 }
971 if (sim_fpu_is_zero (l))
972 {
973 if (sim_fpu_is_zero (r))
974 {
975 *f = sim_fpu_zero;
976 f->sign = l->sign & r->sign;
977 }
978 else
979 *f = *r;
980 return 0;
981 }
982 if (sim_fpu_is_zero (r))
983 {
984 *f = *l;
985 return 0;
986 }
987 {
988 int status = 0;
989 int shift = l->normal_exp - r->normal_exp;
990 unsigned64 lfraction;
991 unsigned64 rfraction;
992 /* use exp of larger */
993 if (shift >= NR_FRAC_GUARD)
994 {
995 /* left has much bigger magnitute */
996 *f = *l;
997 return sim_fpu_status_inexact;
998 }
999 if (shift <= - NR_FRAC_GUARD)
1000 {
1001 /* right has much bigger magnitute */
1002 *f = *r;
1003 return sim_fpu_status_inexact;
1004 }
1005 lfraction = l->fraction;
1006 rfraction = r->fraction;
1007 if (shift > 0)
1008 {
1009 f->normal_exp = l->normal_exp;
1010 if (rfraction & LSMASK64 (shift - 1, 0))
1011 {
1012 status |= sim_fpu_status_inexact;
1013 rfraction |= LSBIT64 (shift); /* stick LSBit */
1014 }
1015 rfraction >>= shift;
1016 }
1017 else if (shift < 0)
1018 {
1019 f->normal_exp = r->normal_exp;
1020 if (lfraction & LSMASK64 (- shift - 1, 0))
1021 {
1022 status |= sim_fpu_status_inexact;
1023 lfraction |= LSBIT64 (- shift); /* stick LSBit */
1024 }
1025 lfraction >>= -shift;
1026 }
1027 else
1028 {
1029 f->normal_exp = r->normal_exp;
1030 }
1031
1032 /* perform the addition */
1033 if (l->sign)
1034 lfraction = - lfraction;
1035 if (r->sign)
1036 rfraction = - rfraction;
1037 f->fraction = lfraction + rfraction;
1038
1039 /* zero? */
1040 if (f->fraction == 0)
1041 {
1042 *f = sim_fpu_zero;
1043 return 0;
1044 }
1045
1046 /* sign? */
1047 f->class = sim_fpu_class_number;
1048 if ((signed64) f->fraction >= 0)
1049 f->sign = 0;
1050 else
1051 {
1052 f->sign = 1;
1053 f->fraction = - f->fraction;
1054 }
1055
1056 /* normalize it */
1057 if ((f->fraction & IMPLICIT_2))
1058 {
1059 f->fraction = (f->fraction >> 1) | (f->fraction & 1);
1060 f->normal_exp ++;
1061 }
1062 else if (f->fraction < IMPLICIT_1)
1063 {
1064 do
1065 {
1066 f->fraction <<= 1;
1067 f->normal_exp --;
1068 }
1069 while (f->fraction < IMPLICIT_1);
1070 }
1071 ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2);
1072 return status;
1073 }
1074 }
1075
1076
1077 INLINE_SIM_FPU (int)
1078 sim_fpu_sub (sim_fpu *f,
1079 const sim_fpu *l,
1080 const sim_fpu *r)
1081 {
1082 if (sim_fpu_is_snan (l))
1083 {
1084 *f = *l;
1085 f->class = sim_fpu_class_qnan;
1086 return sim_fpu_status_invalid_snan;
1087 }
1088 if (sim_fpu_is_snan (r))
1089 {
1090 *f = *r;
1091 f->class = sim_fpu_class_qnan;
1092 return sim_fpu_status_invalid_snan;
1093 }
1094 if (sim_fpu_is_qnan (l))
1095 {
1096 *f = *l;
1097 return 0;
1098 }
1099 if (sim_fpu_is_qnan (r))
1100 {
1101 *f = *r;
1102 return 0;
1103 }
1104 if (sim_fpu_is_infinity (l))
1105 {
1106 if (sim_fpu_is_infinity (r)
1107 && l->sign == r->sign)
1108 {
1109 *f = sim_fpu_qnan;
1110 return sim_fpu_status_invalid_isi;
1111 }
1112 *f = *l;
1113 return 0;
1114 }
1115 if (sim_fpu_is_infinity (r))
1116 {
1117 *f = *r;
1118 f->sign = !r->sign;
1119 return 0;
1120 }
1121 if (sim_fpu_is_zero (l))
1122 {
1123 if (sim_fpu_is_zero (r))
1124 {
1125 *f = sim_fpu_zero;
1126 f->sign = l->sign & !r->sign;
1127 }
1128 else
1129 {
1130 *f = *r;
1131 f->sign = !r->sign;
1132 }
1133 return 0;
1134 }
1135 if (sim_fpu_is_zero (r))
1136 {
1137 *f = *l;
1138 return 0;
1139 }
1140 {
1141 int status = 0;
1142 int shift = l->normal_exp - r->normal_exp;
1143 unsigned64 lfraction;
1144 unsigned64 rfraction;
1145 /* use exp of larger */
1146 if (shift >= NR_FRAC_GUARD)
1147 {
1148 /* left has much bigger magnitute */
1149 *f = *l;
1150 return sim_fpu_status_inexact;
1151 }
1152 if (shift <= - NR_FRAC_GUARD)
1153 {
1154 /* right has much bigger magnitute */
1155 *f = *r;
1156 f->sign = !r->sign;
1157 return sim_fpu_status_inexact;
1158 }
1159 lfraction = l->fraction;
1160 rfraction = r->fraction;
1161 if (shift > 0)
1162 {
1163 f->normal_exp = l->normal_exp;
1164 if (rfraction & LSMASK64 (shift - 1, 0))
1165 {
1166 status |= sim_fpu_status_inexact;
1167 rfraction |= LSBIT64 (shift); /* stick LSBit */
1168 }
1169 rfraction >>= shift;
1170 }
1171 else if (shift < 0)
1172 {
1173 f->normal_exp = r->normal_exp;
1174 if (lfraction & LSMASK64 (- shift - 1, 0))
1175 {
1176 status |= sim_fpu_status_inexact;
1177 lfraction |= LSBIT64 (- shift); /* stick LSBit */
1178 }
1179 lfraction >>= -shift;
1180 }
1181 else
1182 {
1183 f->normal_exp = r->normal_exp;
1184 }
1185
1186 /* perform the subtraction */
1187 if (l->sign)
1188 lfraction = - lfraction;
1189 if (!r->sign)
1190 rfraction = - rfraction;
1191 f->fraction = lfraction + rfraction;
1192
1193 /* zero? */
1194 if (f->fraction == 0)
1195 {
1196 *f = sim_fpu_zero;
1197 return 0;
1198 }
1199
1200 /* sign? */
1201 f->class = sim_fpu_class_number;
1202 if ((signed64) f->fraction >= 0)
1203 f->sign = 0;
1204 else
1205 {
1206 f->sign = 1;
1207 f->fraction = - f->fraction;
1208 }
1209
1210 /* normalize it */
1211 if ((f->fraction & IMPLICIT_2))
1212 {
1213 f->fraction = (f->fraction >> 1) | (f->fraction & 1);
1214 f->normal_exp ++;
1215 }
1216 else if (f->fraction < IMPLICIT_1)
1217 {
1218 do
1219 {
1220 f->fraction <<= 1;
1221 f->normal_exp --;
1222 }
1223 while (f->fraction < IMPLICIT_1);
1224 }
1225 ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2);
1226 return status;
1227 }
1228 }
1229
1230
1231 INLINE_SIM_FPU (int)
1232 sim_fpu_mul (sim_fpu *f,
1233 const sim_fpu *l,
1234 const sim_fpu *r)
1235 {
1236 if (sim_fpu_is_snan (l))
1237 {
1238 *f = *l;
1239 f->class = sim_fpu_class_qnan;
1240 return sim_fpu_status_invalid_snan;
1241 }
1242 if (sim_fpu_is_snan (r))
1243 {
1244 *f = *r;
1245 f->class = sim_fpu_class_qnan;
1246 return sim_fpu_status_invalid_snan;
1247 }
1248 if (sim_fpu_is_qnan (l))
1249 {
1250 *f = *l;
1251 return 0;
1252 }
1253 if (sim_fpu_is_qnan (r))
1254 {
1255 *f = *r;
1256 return 0;
1257 }
1258 if (sim_fpu_is_infinity (l))
1259 {
1260 if (sim_fpu_is_zero (r))
1261 {
1262 *f = sim_fpu_qnan;
1263 return sim_fpu_status_invalid_imz;
1264 }
1265 *f = *l;
1266 f->sign = l->sign ^ r->sign;
1267 return 0;
1268 }
1269 if (sim_fpu_is_infinity (r))
1270 {
1271 if (sim_fpu_is_zero (l))
1272 {
1273 *f = sim_fpu_qnan;
1274 return sim_fpu_status_invalid_imz;
1275 }
1276 *f = *r;
1277 f->sign = l->sign ^ r->sign;
1278 return 0;
1279 }
1280 if (sim_fpu_is_zero (l) || sim_fpu_is_zero (r))
1281 {
1282 *f = sim_fpu_zero;
1283 f->sign = l->sign ^ r->sign;
1284 return 0;
1285 }
1286 /* Calculate the mantissa by multiplying both 64bit numbers to get a
1287 128 bit number */
1288 {
1289 unsigned64 low;
1290 unsigned64 high;
1291 unsigned64 nl = l->fraction & 0xffffffff;
1292 unsigned64 nh = l->fraction >> 32;
1293 unsigned64 ml = r->fraction & 0xffffffff;
1294 unsigned64 mh = r->fraction >>32;
1295 unsigned64 pp_ll = ml * nl;
1296 unsigned64 pp_hl = mh * nl;
1297 unsigned64 pp_lh = ml * nh;
1298 unsigned64 pp_hh = mh * nh;
1299 unsigned64 res2 = 0;
1300 unsigned64 res0 = 0;
1301 unsigned64 ps_hh__ = pp_hl + pp_lh;
1302 if (ps_hh__ < pp_hl)
1303 res2 += UNSIGNED64 (0x100000000);
1304 pp_hl = (ps_hh__ << 32) & UNSIGNED64 (0xffffffff00000000);
1305 res0 = pp_ll + pp_hl;
1306 if (res0 < pp_ll)
1307 res2++;
1308 res2 += ((ps_hh__ >> 32) & 0xffffffff) + pp_hh;
1309 high = res2;
1310 low = res0;
1311
1312 f->normal_exp = l->normal_exp + r->normal_exp;
1313 f->sign = l->sign ^ r->sign;
1314 f->class = sim_fpu_class_number;
1315
1316 /* Input is bounded by [1,2) ; [2^60,2^61)
1317 Output is bounded by [1,4) ; [2^120,2^122) */
1318
1319 /* Adjust the exponent according to where the decimal point ended
1320 up in the high 64 bit word. In the source the decimal point
1321 was at NR_FRAC_GUARD. */
1322 f->normal_exp += NR_FRAC_GUARD + 64 - (NR_FRAC_GUARD * 2);
1323
1324 /* The high word is bounded according to the above. Consequently
1325 it has never overflowed into IMPLICIT_2. */
1326 ASSERT (high < LSBIT64 (((NR_FRAC_GUARD + 1) * 2) - 64));
1327 ASSERT (high >= LSBIT64 ((NR_FRAC_GUARD * 2) - 64));
1328 ASSERT (LSBIT64 (((NR_FRAC_GUARD + 1) * 2) - 64) < IMPLICIT_1);
1329
1330 #if 0
1331 printf ("\n");
1332 print_bits (high, 63, (sim_fpu_print_func*)fprintf, stdout);
1333 printf (";");
1334 print_bits (low, 63, (sim_fpu_print_func*)fprintf, stdout);
1335 printf ("\n");
1336 #endif
1337
1338 /* normalize */
1339 do
1340 {
1341 f->normal_exp--;
1342 high <<= 1;
1343 if (low & LSBIT64 (63))
1344 high |= 1;
1345 low <<= 1;
1346 }
1347 while (high < IMPLICIT_1);
1348
1349 #if 0
1350 print_bits (high, 63, (sim_fpu_print_func*)fprintf, stdout);
1351 printf (";");
1352 print_bits (low, 63, (sim_fpu_print_func*)fprintf, stdout);
1353 printf ("\n");
1354 #endif
1355
1356 ASSERT (high >= IMPLICIT_1 && high < IMPLICIT_2);
1357 if (low != 0)
1358 {
1359 f->fraction = (high | 1); /* sticky */
1360 return sim_fpu_status_inexact;
1361 }
1362 else
1363 {
1364 f->fraction = high;
1365 return 0;
1366 }
1367 return 0;
1368 }
1369 }
1370
1371 INLINE_SIM_FPU (int)
1372 sim_fpu_div (sim_fpu *f,
1373 const sim_fpu *l,
1374 const sim_fpu *r)
1375 {
1376 if (sim_fpu_is_snan (l))
1377 {
1378 *f = *l;
1379 f->class = sim_fpu_class_qnan;
1380 return sim_fpu_status_invalid_snan;
1381 }
1382 if (sim_fpu_is_snan (r))
1383 {
1384 *f = *r;
1385 f->class = sim_fpu_class_qnan;
1386 return sim_fpu_status_invalid_snan;
1387 }
1388 if (sim_fpu_is_qnan (l))
1389 {
1390 *f = *l;
1391 f->class = sim_fpu_class_qnan;
1392 return 0;
1393 }
1394 if (sim_fpu_is_qnan (r))
1395 {
1396 *f = *r;
1397 f->class = sim_fpu_class_qnan;
1398 return 0;
1399 }
1400 if (sim_fpu_is_infinity (l))
1401 {
1402 if (sim_fpu_is_infinity (r))
1403 {
1404 *f = sim_fpu_qnan;
1405 return sim_fpu_status_invalid_idi;
1406 }
1407 else
1408 {
1409 *f = *l;
1410 f->sign = l->sign ^ r->sign;
1411 return 0;
1412 }
1413 }
1414 if (sim_fpu_is_zero (l))
1415 {
1416 if (sim_fpu_is_zero (r))
1417 {
1418 *f = sim_fpu_qnan;
1419 return sim_fpu_status_invalid_zdz;
1420 }
1421 else
1422 {
1423 *f = *l;
1424 f->sign = l->sign ^ r->sign;
1425 return 0;
1426 }
1427 }
1428 if (sim_fpu_is_infinity (r))
1429 {
1430 *f = sim_fpu_zero;
1431 f->sign = l->sign ^ r->sign;
1432 return 0;
1433 }
1434 if (sim_fpu_is_zero (r))
1435 {
1436 f->class = sim_fpu_class_infinity;
1437 f->sign = l->sign ^ r->sign;
1438 return sim_fpu_status_invalid_div0;
1439 }
1440
1441 /* Calculate the mantissa by multiplying both 64bit numbers to get a
1442 128 bit number */
1443 {
1444 /* quotient = ( ( numerator / denominator)
1445 x 2^(numerator exponent - denominator exponent)
1446 */
1447 unsigned64 numerator;
1448 unsigned64 denominator;
1449 unsigned64 quotient;
1450 unsigned64 bit;
1451
1452 f->class = sim_fpu_class_number;
1453 f->sign = l->sign ^ r->sign;
1454 f->normal_exp = l->normal_exp - r->normal_exp;
1455
1456 numerator = l->fraction;
1457 denominator = r->fraction;
1458
1459 /* Fraction will be less than 1.0 */
1460 if (numerator < denominator)
1461 {
1462 numerator <<= 1;
1463 f->normal_exp--;
1464 }
1465 ASSERT (numerator >= denominator);
1466
1467 /* Gain extra precision, already used one spare bit */
1468 numerator <<= NR_SPARE;
1469 denominator <<= NR_SPARE;
1470
1471 /* Does divide one bit at a time. Optimize??? */
1472 quotient = 0;
1473 bit = (IMPLICIT_1 << NR_SPARE);
1474 while (bit)
1475 {
1476 if (numerator >= denominator)
1477 {
1478 quotient |= bit;
1479 numerator -= denominator;
1480 }
1481 bit >>= 1;
1482 numerator <<= 1;
1483 }
1484
1485 #if 0
1486 printf ("\n");
1487 print_bits (quotient, 63, (sim_fpu_print_func*)fprintf, stdout);
1488 printf ("\n");
1489 print_bits (numerator, 63, (sim_fpu_print_func*)fprintf, stdout);
1490 printf ("\n");
1491 print_bits (denominator, 63, (sim_fpu_print_func*)fprintf, stdout);
1492 printf ("\n");
1493 #endif
1494
1495 /* discard (but save) the extra bits */
1496 if ((quotient & LSMASK64 (NR_SPARE -1, 0)))
1497 quotient = (quotient >> NR_SPARE) | 1;
1498 else
1499 quotient = (quotient >> NR_SPARE);
1500
1501 f->fraction = quotient;
1502 ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2);
1503 if (numerator != 0)
1504 {
1505 f->fraction |= 1; /* stick remaining bits */
1506 return sim_fpu_status_inexact;
1507 }
1508 else
1509 return 0;
1510 }
1511 }
1512
1513
1514 INLINE_SIM_FPU (int)
1515 sim_fpu_neg (sim_fpu *f,
1516 const sim_fpu *r)
1517 {
1518 if (sim_fpu_is_snan (r))
1519 {
1520 *f = *r;
1521 f->class = sim_fpu_class_qnan;
1522 return sim_fpu_status_invalid_snan;
1523 }
1524 if (sim_fpu_is_qnan (r))
1525 {
1526 *f = *r;
1527 return 0;
1528 }
1529 *f = *r;
1530 f->sign = !r->sign;
1531 return 0;
1532 }
1533
1534
1535 INLINE_SIM_FPU (int)
1536 sim_fpu_abs (sim_fpu *f,
1537 const sim_fpu *r)
1538 {
1539 if (sim_fpu_is_snan (r))
1540 {
1541 *f = *r;
1542 f->class = sim_fpu_class_qnan;
1543 return sim_fpu_status_invalid_snan;
1544 }
1545 if (sim_fpu_is_qnan (r))
1546 {
1547 *f = *r;
1548 return 0;
1549 }
1550 *f = *r;
1551 f->sign = 0;
1552 return 0;
1553 }
1554
1555
1556 INLINE_SIM_FPU (int)
1557 sim_fpu_inv (sim_fpu *f,
1558 const sim_fpu *r)
1559 {
1560 if (sim_fpu_is_snan (r))
1561 {
1562 *f = *r;
1563 f->class = sim_fpu_class_qnan;
1564 return sim_fpu_status_invalid_snan;
1565 }
1566 if (sim_fpu_is_qnan (r))
1567 {
1568 *f = *r;
1569 f->class = sim_fpu_class_qnan;
1570 return 0;
1571 }
1572 if (sim_fpu_is_infinity (r))
1573 {
1574 *f = sim_fpu_zero;
1575 f->sign = r->sign;
1576 return 0;
1577 }
1578 if (sim_fpu_is_zero (r))
1579 {
1580 f->class = sim_fpu_class_infinity;
1581 f->sign = r->sign;
1582 return sim_fpu_status_invalid_div0;
1583 }
1584 *f = *r;
1585 f->normal_exp = - r->normal_exp;
1586 return 0;
1587 }
1588
1589
1590 INLINE_SIM_FPU (int)
1591 sim_fpu_sqrt (sim_fpu *f,
1592 const sim_fpu *r)
1593 {
1594 if (sim_fpu_is_snan (r))
1595 {
1596 *f = sim_fpu_qnan;
1597 return sim_fpu_status_invalid_snan;
1598 }
1599 if (sim_fpu_is_qnan (r))
1600 {
1601 *f = sim_fpu_qnan;
1602 return 0;
1603 }
1604 if (sim_fpu_is_zero (r))
1605 {
1606 f->class = sim_fpu_class_zero;
1607 f->sign = r->sign;
1608 return 0;
1609 }
1610 if (sim_fpu_is_infinity (r))
1611 {
1612 if (r->sign)
1613 {
1614 *f = sim_fpu_qnan;
1615 return sim_fpu_status_invalid_sqrt;
1616 }
1617 else
1618 {
1619 f->class = sim_fpu_class_infinity;
1620 f->sign = 0;
1621 f->sign = 0;
1622 return 0;
1623 }
1624 }
1625 if (r->sign)
1626 {
1627 *f = sim_fpu_qnan;
1628 return sim_fpu_status_invalid_sqrt;
1629 }
1630
1631 /* @(#)e_sqrt.c 5.1 93/09/24 */
1632 /*
1633 * ====================================================
1634 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1635 *
1636 * Developed at SunPro, a Sun Microsystems, Inc. business.
1637 * Permission to use, copy, modify, and distribute this
1638 * software is freely granted, provided that this notice
1639 * is preserved.
1640 * ====================================================
1641 */
1642
1643 /* __ieee754_sqrt(x)
1644 * Return correctly rounded sqrt.
1645 * ------------------------------------------
1646 * | Use the hardware sqrt if you have one |
1647 * ------------------------------------------
1648 * Method:
1649 * Bit by bit method using integer arithmetic. (Slow, but portable)
1650 * 1. Normalization
1651 * Scale x to y in [1,4) with even powers of 2:
1652 * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
1653 * sqrt(x) = 2^k * sqrt(y)
1654 -
1655 - Since:
1656 - sqrt ( x*2^(2m) ) = sqrt(x).2^m ; m even
1657 - sqrt ( x*2^(2m + 1) ) = sqrt(2.x).2^m ; m odd
1658 - Define:
1659 - y = ((m even) ? x : 2.x)
1660 - Then:
1661 - y in [1, 4) ; [IMPLICIT_1,IMPLICIT_4)
1662 - And:
1663 - sqrt (y) in [1, 2) ; [IMPLICIT_1,IMPLICIT_2)
1664 -
1665 * 2. Bit by bit computation
1666 * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
1667 * i 0
1668 * i+1 2
1669 * s = 2*q , and y = 2 * ( y - q ). (1)
1670 * i i i i
1671 *
1672 * To compute q from q , one checks whether
1673 * i+1 i
1674 *
1675 * -(i+1) 2
1676 * (q + 2 ) <= y. (2)
1677 * i
1678 * -(i+1)
1679 * If (2) is false, then q = q ; otherwise q = q + 2 .
1680 * i+1 i i+1 i
1681 *
1682 * With some algebric manipulation, it is not difficult to see
1683 * that (2) is equivalent to
1684 * -(i+1)
1685 * s + 2 <= y (3)
1686 * i i
1687 *
1688 * The advantage of (3) is that s and y can be computed by
1689 * i i
1690 * the following recurrence formula:
1691 * if (3) is false
1692 *
1693 * s = s , y = y ; (4)
1694 * i+1 i i+1 i
1695 *
1696 -
1697 - NOTE: y = 2*y
1698 - i+1 i
1699 -
1700 * otherwise,
1701 * -i -(i+1)
1702 * s = s + 2 , y = y - s - 2 (5)
1703 * i+1 i i+1 i i
1704 *
1705 -
1706 - -(i+1)
1707 - NOTE: y = 2 (y - s - 2 )
1708 - i+1 i i
1709 -
1710 * One may easily use induction to prove (4) and (5).
1711 * Note. Since the left hand side of (3) contain only i+2 bits,
1712 * it does not necessary to do a full (53-bit) comparison
1713 * in (3).
1714 * 3. Final rounding
1715 * After generating the 53 bits result, we compute one more bit.
1716 * Together with the remainder, we can decide whether the
1717 * result is exact, bigger than 1/2ulp, or less than 1/2ulp
1718 * (it will never equal to 1/2ulp).
1719 * The rounding mode can be detected by checking whether
1720 * huge + tiny is equal to huge, and whether huge - tiny is
1721 * equal to huge for some floating point number "huge" and "tiny".
1722 *
1723 * Special cases:
1724 * sqrt(+-0) = +-0 ... exact
1725 * sqrt(inf) = inf
1726 * sqrt(-ve) = NaN ... with invalid signal
1727 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
1728 *
1729 * Other methods : see the appended file at the end of the program below.
1730 *---------------
1731 */
1732
1733 {
1734 /* generate sqrt(x) bit by bit */
1735 unsigned64 y;
1736 unsigned64 q;
1737 unsigned64 s;
1738 unsigned64 b;
1739
1740 f->class = sim_fpu_class_number;
1741 f->sign = 0;
1742 y = r->fraction;
1743 f->normal_exp = (r->normal_exp >> 1); /* exp = [exp/2] */
1744
1745 /* odd exp, double x to make it even */
1746 ASSERT (y >= IMPLICIT_1 && y < IMPLICIT_4);
1747 if ((r->normal_exp & 1))
1748 {
1749 y += y;
1750 }
1751 ASSERT (y >= IMPLICIT_1 && y < (IMPLICIT_2 << 1));
1752
1753 /* Let loop determine first value of s (either 1 or 2) */
1754 b = IMPLICIT_1;
1755 q = 0;
1756 s = 0;
1757
1758 while (b)
1759 {
1760 unsigned64 t = s + b;
1761 if (t <= y)
1762 {
1763 s |= (b << 1);
1764 y -= t;
1765 q |= b;
1766 }
1767 y <<= 1;
1768 b >>= 1;
1769 }
1770
1771 ASSERT (q >= IMPLICIT_1 && q < IMPLICIT_2);
1772 f->fraction = q;
1773 if (y != 0)
1774 {
1775 f->fraction |= 1; /* stick remaining bits */
1776 return sim_fpu_status_inexact;
1777 }
1778 else
1779 return 0;
1780 }
1781 }
1782
1783
1784 /* int/long <-> sim_fpu */
1785
1786 INLINE_SIM_FPU (int)
1787 sim_fpu_i32to (sim_fpu *f,
1788 signed32 i,
1789 sim_fpu_round round)
1790 {
1791 i2fpu (f, i, 0);
1792 return 0;
1793 }
1794
1795 INLINE_SIM_FPU (int)
1796 sim_fpu_u32to (sim_fpu *f,
1797 unsigned32 u,
1798 sim_fpu_round round)
1799 {
1800 u2fpu (f, u, 0);
1801 return 0;
1802 }
1803
1804 INLINE_SIM_FPU (int)
1805 sim_fpu_i64to (sim_fpu *f,
1806 signed64 i,
1807 sim_fpu_round round)
1808 {
1809 i2fpu (f, i, 1);
1810 return 0;
1811 }
1812
1813 INLINE_SIM_FPU (int)
1814 sim_fpu_u64to (sim_fpu *f,
1815 unsigned64 u,
1816 sim_fpu_round round)
1817 {
1818 u2fpu (f, u, 1);
1819 return 0;
1820 }
1821
1822
1823 INLINE_SIM_FPU (int)
1824 sim_fpu_to32i (signed32 *i,
1825 const sim_fpu *f,
1826 sim_fpu_round round)
1827 {
1828 signed64 i64;
1829 int status = fpu2i (&i64, f, 0, round);
1830 *i = i64;
1831 return status;
1832 }
1833
1834 INLINE_SIM_FPU (int)
1835 sim_fpu_to32u (unsigned32 *u,
1836 const sim_fpu *f,
1837 sim_fpu_round round)
1838 {
1839 unsigned64 u64;
1840 int status = fpu2u (&u64, f, 0);
1841 *u = u64;
1842 return status;
1843 }
1844
1845 INLINE_SIM_FPU (int)
1846 sim_fpu_to64i (signed64 *i,
1847 const sim_fpu *f,
1848 sim_fpu_round round)
1849 {
1850 return fpu2i (i, f, 1, round);
1851 }
1852
1853
1854 INLINE_SIM_FPU (int)
1855 sim_fpu_to64u (unsigned64 *u,
1856 const sim_fpu *f,
1857 sim_fpu_round round)
1858 {
1859 return fpu2u (u, f, 1);
1860 }
1861
1862
1863
1864 /* sim_fpu -> host format */
1865
1866 #if 0
1867 INLINE_SIM_FPU (float)
1868 sim_fpu_2f (const sim_fpu *f)
1869 {
1870 return fval.d;
1871 }
1872 #endif
1873
1874
1875 INLINE_SIM_FPU (double)
1876 sim_fpu_2d (const sim_fpu *s)
1877 {
1878 sim_fpu_map val;
1879 val.i = pack_fpu (s, 1);
1880 return val.d;
1881 }
1882
1883
1884 #if 0
1885 INLINE_SIM_FPU (void)
1886 sim_fpu_f2 (sim_fpu *f,
1887 float s)
1888 {
1889 sim_fpu_map val;
1890 val.d = s;
1891 unpack_fpu (f, val.i, 1);
1892 }
1893 #endif
1894
1895
1896 INLINE_SIM_FPU (void)
1897 sim_fpu_d2 (sim_fpu *f,
1898 double d)
1899 {
1900 sim_fpu_map val;
1901 val.d = d;
1902 unpack_fpu (f, val.i, 1);
1903 }
1904
1905
1906 /* General */
1907
1908 INLINE_SIM_FPU (int)
1909 sim_fpu_is_nan (const sim_fpu *d)
1910 {
1911 switch (d->class)
1912 {
1913 case sim_fpu_class_qnan:
1914 case sim_fpu_class_snan:
1915 return 1;
1916 default:
1917 return 0;
1918 }
1919 }
1920
1921 INLINE_SIM_FPU (int)
1922 sim_fpu_is_qnan (const sim_fpu *d)
1923 {
1924 switch (d->class)
1925 {
1926 case sim_fpu_class_qnan:
1927 return 1;
1928 default:
1929 return 0;
1930 }
1931 }
1932
1933 INLINE_SIM_FPU (int)
1934 sim_fpu_is_snan (const sim_fpu *d)
1935 {
1936 switch (d->class)
1937 {
1938 case sim_fpu_class_snan:
1939 return 1;
1940 default:
1941 return 0;
1942 }
1943 }
1944
1945 INLINE_SIM_FPU (int)
1946 sim_fpu_is_zero (const sim_fpu *d)
1947 {
1948 switch (d->class)
1949 {
1950 case sim_fpu_class_zero:
1951 return 1;
1952 default:
1953 return 0;
1954 }
1955 }
1956
1957 INLINE_SIM_FPU (int)
1958 sim_fpu_is_infinity (const sim_fpu *d)
1959 {
1960 switch (d->class)
1961 {
1962 case sim_fpu_class_infinity:
1963 return 1;
1964 default:
1965 return 0;
1966 }
1967 }
1968
1969 INLINE_SIM_FPU (int)
1970 sim_fpu_is_number (const sim_fpu *d)
1971 {
1972 switch (d->class)
1973 {
1974 case sim_fpu_class_denorm:
1975 case sim_fpu_class_number:
1976 return 1;
1977 default:
1978 return 0;
1979 }
1980 }
1981
1982 INLINE_SIM_FPU (int)
1983 sim_fpu_is_denorm (const sim_fpu *d)
1984 {
1985 switch (d->class)
1986 {
1987 case sim_fpu_class_denorm:
1988 return 1;
1989 default:
1990 return 0;
1991 }
1992 }
1993
1994 INLINE_SIM_FPU (int)
1995 sim_fpu_is (const sim_fpu *d)
1996 {
1997 switch (d->class)
1998 {
1999 case sim_fpu_class_qnan:
2000 return SIM_FPU_IS_QNAN;
2001 case sim_fpu_class_snan:
2002 return SIM_FPU_IS_SNAN;
2003 case sim_fpu_class_infinity:
2004 return SIM_FPU_IS_NINF;
2005 return SIM_FPU_IS_PINF;
2006 case sim_fpu_class_number:
2007 if (d->sign)
2008 return SIM_FPU_IS_NNUMBER;
2009 else
2010 return SIM_FPU_IS_PNUMBER;
2011 case sim_fpu_class_denorm:
2012 if (d->sign)
2013 return SIM_FPU_IS_NDENORM;
2014 else
2015 return SIM_FPU_IS_PDENORM;
2016 case sim_fpu_class_zero:
2017 if (d->sign)
2018 return SIM_FPU_IS_NZERO;
2019 else
2020 return SIM_FPU_IS_PZERO;
2021 default:
2022 return -1;
2023 abort ();
2024 }
2025 }
2026
2027 INLINE_SIM_FPU (int)
2028 sim_fpu_cmp (const sim_fpu *l, const sim_fpu *r)
2029 {
2030 sim_fpu res;
2031 sim_fpu_sub (&res, l, r);
2032 return sim_fpu_is (&res);
2033 }
2034
2035 INLINE_SIM_FPU (int)
2036 sim_fpu_is_lt (const sim_fpu *l, const sim_fpu *r)
2037 {
2038 int status;
2039 sim_fpu_lt (&status, l, r);
2040 return status;
2041 }
2042
2043 INLINE_SIM_FPU (int)
2044 sim_fpu_is_le (const sim_fpu *l, const sim_fpu *r)
2045 {
2046 int is;
2047 sim_fpu_le (&is, l, r);
2048 return is;
2049 }
2050
2051 INLINE_SIM_FPU (int)
2052 sim_fpu_is_eq (const sim_fpu *l, const sim_fpu *r)
2053 {
2054 int is;
2055 sim_fpu_eq (&is, l, r);
2056 return is;
2057 }
2058
2059 INLINE_SIM_FPU (int)
2060 sim_fpu_is_ne (const sim_fpu *l, const sim_fpu *r)
2061 {
2062 int is;
2063 sim_fpu_ne (&is, l, r);
2064 return is;
2065 }
2066
2067 INLINE_SIM_FPU (int)
2068 sim_fpu_is_ge (const sim_fpu *l, const sim_fpu *r)
2069 {
2070 int is;
2071 sim_fpu_ge (&is, l, r);
2072 return is;
2073 }
2074
2075 INLINE_SIM_FPU (int)
2076 sim_fpu_is_gt (const sim_fpu *l, const sim_fpu *r)
2077 {
2078 int is;
2079 sim_fpu_gt (&is, l, r);
2080 return is;
2081 }
2082
2083
2084 /* Compare operators */
2085
2086 INLINE_SIM_FPU (int)
2087 sim_fpu_lt (int *is,
2088 const sim_fpu *l,
2089 const sim_fpu *r)
2090 {
2091 if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2092 {
2093 sim_fpu_map lval;
2094 sim_fpu_map rval;
2095 lval.i = pack_fpu (l, 1);
2096 rval.i = pack_fpu (r, 1);
2097 (*is) = (lval.d < rval.d);
2098 return 0;
2099 }
2100 else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2101 {
2102 *is = 0;
2103 return sim_fpu_status_invalid_snan;
2104 }
2105 else
2106 {
2107 *is = 0;
2108 return sim_fpu_status_invalid_qnan;
2109 }
2110 }
2111
2112 INLINE_SIM_FPU (int)
2113 sim_fpu_le (int *is,
2114 const sim_fpu *l,
2115 const sim_fpu *r)
2116 {
2117 if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2118 {
2119 sim_fpu_map lval;
2120 sim_fpu_map rval;
2121 lval.i = pack_fpu (l, 1);
2122 rval.i = pack_fpu (r, 1);
2123 *is = (lval.d <= rval.d);
2124 return 0;
2125 }
2126 else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2127 {
2128 *is = 0;
2129 return sim_fpu_status_invalid_snan;
2130 }
2131 else
2132 {
2133 *is = 0;
2134 return sim_fpu_status_invalid_qnan;
2135 }
2136 }
2137
2138 INLINE_SIM_FPU (int)
2139 sim_fpu_eq (int *is,
2140 const sim_fpu *l,
2141 const sim_fpu *r)
2142 {
2143 if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2144 {
2145 sim_fpu_map lval;
2146 sim_fpu_map rval;
2147 lval.i = pack_fpu (l, 1);
2148 rval.i = pack_fpu (r, 1);
2149 (*is) = (lval.d == rval.d);
2150 return 0;
2151 }
2152 else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2153 {
2154 *is = 0;
2155 return sim_fpu_status_invalid_snan;
2156 }
2157 else
2158 {
2159 *is = 0;
2160 return sim_fpu_status_invalid_qnan;
2161 }
2162 }
2163
2164 INLINE_SIM_FPU (int)
2165 sim_fpu_ne (int *is,
2166 const sim_fpu *l,
2167 const sim_fpu *r)
2168 {
2169 if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2170 {
2171 sim_fpu_map lval;
2172 sim_fpu_map rval;
2173 lval.i = pack_fpu (l, 1);
2174 rval.i = pack_fpu (r, 1);
2175 (*is) = (lval.d != rval.d);
2176 return 0;
2177 }
2178 else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2179 {
2180 *is = 0;
2181 return sim_fpu_status_invalid_snan;
2182 }
2183 else
2184 {
2185 *is = 0;
2186 return sim_fpu_status_invalid_qnan;
2187 }
2188 }
2189
2190 INLINE_SIM_FPU (int)
2191 sim_fpu_ge (int *is,
2192 const sim_fpu *l,
2193 const sim_fpu *r)
2194 {
2195 return sim_fpu_le (is, r, l);
2196 }
2197
2198 INLINE_SIM_FPU (int)
2199 sim_fpu_gt (int *is,
2200 const sim_fpu *l,
2201 const sim_fpu *r)
2202 {
2203 return sim_fpu_lt (is, r, l);
2204 }
2205
2206
2207 /* A number of useful constants */
2208
2209 const sim_fpu sim_fpu_zero = { sim_fpu_class_zero, };
2210 const sim_fpu sim_fpu_qnan = { sim_fpu_class_qnan, };
2211
2212
2213 /* For debugging */
2214
2215 INLINE_SIM_FPU (void)
2216 sim_fpu_print_fpu (const sim_fpu *f,
2217 sim_fpu_print_func *print,
2218 void *arg)
2219 {
2220 print (arg, "%s", f->sign ? "-" : "+");
2221 switch (f->class)
2222 {
2223 case sim_fpu_class_qnan:
2224 print (arg, "0.");
2225 print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg);
2226 print (arg, "*QuietNaN");
2227 break;
2228 case sim_fpu_class_snan:
2229 print (arg, "0.");
2230 print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg);
2231 print (arg, "*SignalNaN");
2232 break;
2233 case sim_fpu_class_zero:
2234 print (arg, "0.0");
2235 break;
2236 case sim_fpu_class_infinity:
2237 print (arg, "INF");
2238 break;
2239 case sim_fpu_class_number:
2240 case sim_fpu_class_denorm:
2241 print (arg, "1.");
2242 print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg);
2243 print (arg, "*2^%+-5d", f->normal_exp);
2244 ASSERT (f->fraction >= IMPLICIT_1);
2245 ASSERT (f->fraction < IMPLICIT_2);
2246 }
2247 }
2248
2249
2250 INLINE_SIM_FPU (void)
2251 sim_fpu_print_status (int status,
2252 sim_fpu_print_func *print,
2253 void *arg)
2254 {
2255 int i = 1;
2256 char *prefix = "";
2257 while (status >= i)
2258 {
2259 switch ((sim_fpu_status) (status & i))
2260 {
2261 case sim_fpu_status_denorm:
2262 print (arg, "%sD", prefix);
2263 break;
2264 case sim_fpu_status_invalid_snan:
2265 print (arg, "%sSNaN", prefix);
2266 break;
2267 case sim_fpu_status_invalid_qnan:
2268 print (arg, "%sQNaN", prefix);
2269 break;
2270 case sim_fpu_status_invalid_isi:
2271 print (arg, "%sISI", prefix);
2272 break;
2273 case sim_fpu_status_invalid_idi:
2274 print (arg, "%sIDI", prefix);
2275 break;
2276 case sim_fpu_status_invalid_zdz:
2277 print (arg, "%sZDZ", prefix);
2278 break;
2279 case sim_fpu_status_invalid_imz:
2280 print (arg, "%sIMZ", prefix);
2281 break;
2282 case sim_fpu_status_invalid_cvi:
2283 print (arg, "%sCVI", prefix);
2284 break;
2285 case sim_fpu_status_invalid_cmp:
2286 print (arg, "%sCMP", prefix);
2287 break;
2288 case sim_fpu_status_invalid_sqrt:
2289 print (arg, "%sSQRT", prefix);
2290 break;
2291 break;
2292 case sim_fpu_status_inexact:
2293 print (arg, "%sX", prefix);
2294 break;
2295 break;
2296 case sim_fpu_status_overflow:
2297 print (arg, "%sO", prefix);
2298 break;
2299 break;
2300 case sim_fpu_status_underflow:
2301 print (arg, "%sU", prefix);
2302 break;
2303 break;
2304 case sim_fpu_status_invalid_div0:
2305 print (arg, "%s/", prefix);
2306 break;
2307 break;
2308 case sim_fpu_status_rounded:
2309 print (arg, "%sR", prefix);
2310 break;
2311 break;
2312 }
2313 i <<= 1;
2314 prefix = ",";
2315 }
2316 }
2317
2318 #endif
This page took 0.089505 seconds and 4 git commands to generate.