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