This commit was generated by cvs2svn to track changes on a CVS vendor
[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 >>= 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   val.i = pack_fpu (s, 1);
2094   return val.d;
2095 }
2096
2097
2098 #if 0
2099 INLINE_SIM_FPU (void)
2100 sim_fpu_f2 (sim_fpu *f,
2101             float s)
2102 {
2103   sim_fpu_map val;
2104   val.d = s;
2105   unpack_fpu (f, val.i, 1);
2106 }
2107 #endif
2108
2109
2110 INLINE_SIM_FPU (void)
2111 sim_fpu_d2 (sim_fpu *f,
2112             double d)
2113 {
2114   sim_fpu_map val;
2115   val.d = d;
2116   unpack_fpu (f, val.i, 1);
2117 }
2118
2119
2120 /* General */
2121
2122 INLINE_SIM_FPU (int)
2123 sim_fpu_is_nan (const sim_fpu *d)
2124 {
2125   switch (d->class)
2126     {
2127     case sim_fpu_class_qnan:
2128     case sim_fpu_class_snan:
2129       return 1;
2130     default:
2131       return 0;
2132     }
2133 }
2134
2135 INLINE_SIM_FPU (int)
2136 sim_fpu_is_qnan (const sim_fpu *d)
2137 {
2138   switch (d->class)
2139     {
2140     case sim_fpu_class_qnan:
2141       return 1;
2142     default:
2143       return 0;
2144     }
2145 }
2146
2147 INLINE_SIM_FPU (int)
2148 sim_fpu_is_snan (const sim_fpu *d)
2149 {
2150   switch (d->class)
2151     {
2152     case sim_fpu_class_snan:
2153       return 1;
2154     default:
2155       return 0;
2156     }
2157 }
2158
2159 INLINE_SIM_FPU (int)
2160 sim_fpu_is_zero (const sim_fpu *d)
2161 {
2162   switch (d->class)
2163     {
2164     case sim_fpu_class_zero:
2165       return 1;
2166     default:
2167       return 0;
2168     }
2169 }
2170
2171 INLINE_SIM_FPU (int)
2172 sim_fpu_is_infinity (const sim_fpu *d)
2173 {
2174   switch (d->class)
2175     {
2176     case sim_fpu_class_infinity:
2177       return 1;
2178     default:
2179       return 0;
2180     }
2181 }
2182
2183 INLINE_SIM_FPU (int)
2184 sim_fpu_is_number (const sim_fpu *d)
2185 {
2186   switch (d->class)
2187     {
2188     case sim_fpu_class_denorm:
2189     case sim_fpu_class_number:
2190       return 1;
2191     default:
2192       return 0;
2193     }
2194 }
2195
2196 INLINE_SIM_FPU (int)
2197 sim_fpu_is_denorm (const sim_fpu *d)
2198 {
2199   switch (d->class)
2200     {
2201     case sim_fpu_class_denorm:
2202       return 1;
2203     default:
2204       return 0;
2205     }
2206 }
2207
2208
2209 INLINE_SIM_FPU (int)
2210 sim_fpu_sign (const sim_fpu *d)
2211 {
2212   return d->sign;
2213 }
2214
2215
2216 INLINE_SIM_FPU (int)
2217 sim_fpu_exp (const sim_fpu *d)
2218 {
2219   return d->normal_exp;
2220 }
2221
2222
2223
2224 INLINE_SIM_FPU (int)
2225 sim_fpu_is (const sim_fpu *d)
2226 {
2227   switch (d->class)
2228     {
2229     case sim_fpu_class_qnan:
2230       return SIM_FPU_IS_QNAN;
2231     case sim_fpu_class_snan:
2232       return SIM_FPU_IS_SNAN;
2233     case sim_fpu_class_infinity:
2234       if (d->sign)
2235         return SIM_FPU_IS_NINF;
2236       else
2237         return SIM_FPU_IS_PINF;
2238     case sim_fpu_class_number:
2239       if (d->sign)
2240         return SIM_FPU_IS_NNUMBER;
2241       else
2242         return SIM_FPU_IS_PNUMBER;
2243     case sim_fpu_class_denorm:
2244       if (d->sign)
2245         return SIM_FPU_IS_NDENORM;
2246       else
2247         return SIM_FPU_IS_PDENORM;
2248     case sim_fpu_class_zero:
2249       if (d->sign)
2250         return SIM_FPU_IS_NZERO;
2251       else
2252         return SIM_FPU_IS_PZERO;
2253     default:
2254       return -1;
2255       abort ();
2256     }
2257 }
2258
2259 INLINE_SIM_FPU (int)
2260 sim_fpu_cmp (const sim_fpu *l, const sim_fpu *r)
2261 {
2262   sim_fpu res;
2263   sim_fpu_sub (&res, l, r);
2264   return sim_fpu_is (&res);
2265 }
2266
2267 INLINE_SIM_FPU (int)
2268 sim_fpu_is_lt (const sim_fpu *l, const sim_fpu *r)
2269 {
2270   int status;
2271   sim_fpu_lt (&status, l, r);
2272   return status;
2273 }
2274
2275 INLINE_SIM_FPU (int)
2276 sim_fpu_is_le (const sim_fpu *l, const sim_fpu *r)
2277 {
2278   int is;
2279   sim_fpu_le (&is, l, r);
2280   return is;
2281 }
2282
2283 INLINE_SIM_FPU (int)
2284 sim_fpu_is_eq (const sim_fpu *l, const sim_fpu *r)
2285 {
2286   int is;
2287   sim_fpu_eq (&is, l, r);
2288   return is;
2289 }
2290
2291 INLINE_SIM_FPU (int)
2292 sim_fpu_is_ne (const sim_fpu *l, const sim_fpu *r)
2293 {
2294   int is;
2295   sim_fpu_ne (&is, l, r);
2296   return is;
2297 }
2298
2299 INLINE_SIM_FPU (int)
2300 sim_fpu_is_ge (const sim_fpu *l, const sim_fpu *r)
2301 {
2302   int is;
2303   sim_fpu_ge (&is, l, r);
2304   return is;
2305 }
2306
2307 INLINE_SIM_FPU (int)
2308 sim_fpu_is_gt (const sim_fpu *l, const sim_fpu *r)
2309 {
2310   int is;
2311   sim_fpu_gt (&is, l, r);
2312   return is;
2313 }
2314
2315
2316 /* Compare operators */
2317
2318 INLINE_SIM_FPU (int)
2319 sim_fpu_lt (int *is,
2320             const sim_fpu *l,
2321             const sim_fpu *r)
2322 {
2323   if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2324     {
2325       sim_fpu_map lval;
2326       sim_fpu_map rval;
2327       lval.i = pack_fpu (l, 1);
2328       rval.i = pack_fpu (r, 1);
2329       (*is) = (lval.d < rval.d);
2330       return 0;
2331     }
2332   else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2333     {
2334       *is = 0;
2335       return sim_fpu_status_invalid_snan;
2336     }
2337   else
2338     {
2339       *is = 0;
2340       return sim_fpu_status_invalid_qnan;
2341     }
2342 }
2343
2344 INLINE_SIM_FPU (int)
2345 sim_fpu_le (int *is,
2346             const sim_fpu *l,
2347             const sim_fpu *r)
2348 {
2349   if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2350     {
2351       sim_fpu_map lval;
2352       sim_fpu_map rval;
2353       lval.i = pack_fpu (l, 1);
2354       rval.i = pack_fpu (r, 1);
2355       *is = (lval.d <= rval.d);
2356       return 0;
2357     }
2358   else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2359     {
2360       *is = 0;
2361       return sim_fpu_status_invalid_snan;
2362     }
2363   else
2364     {
2365       *is = 0;
2366       return sim_fpu_status_invalid_qnan;
2367     }
2368 }
2369
2370 INLINE_SIM_FPU (int)
2371 sim_fpu_eq (int *is,
2372             const sim_fpu *l,
2373             const sim_fpu *r)
2374 {
2375   if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2376     {
2377       sim_fpu_map lval;
2378       sim_fpu_map rval;
2379       lval.i = pack_fpu (l, 1);
2380       rval.i = pack_fpu (r, 1);
2381       (*is) = (lval.d == rval.d);
2382       return 0;
2383     }
2384   else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2385     {
2386       *is = 0;
2387       return sim_fpu_status_invalid_snan;
2388     }
2389   else
2390     {
2391       *is = 0;
2392       return sim_fpu_status_invalid_qnan;
2393     }
2394 }
2395
2396 INLINE_SIM_FPU (int)
2397 sim_fpu_ne (int *is,
2398             const sim_fpu *l,
2399             const sim_fpu *r)
2400 {
2401   if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r))
2402     {
2403       sim_fpu_map lval;
2404       sim_fpu_map rval;
2405       lval.i = pack_fpu (l, 1);
2406       rval.i = pack_fpu (r, 1);
2407       (*is) = (lval.d != rval.d);
2408       return 0;
2409     }
2410   else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r))
2411     {
2412       *is = 0;
2413       return sim_fpu_status_invalid_snan;
2414     }
2415   else
2416     {
2417       *is = 0;
2418       return sim_fpu_status_invalid_qnan;
2419     }
2420 }
2421
2422 INLINE_SIM_FPU (int)
2423 sim_fpu_ge (int *is,
2424             const sim_fpu *l,
2425             const sim_fpu *r)
2426 {
2427   return sim_fpu_le (is, r, l);
2428 }
2429
2430 INLINE_SIM_FPU (int)
2431 sim_fpu_gt (int *is,
2432             const sim_fpu *l,
2433             const sim_fpu *r)
2434 {
2435   return sim_fpu_lt (is, r, l);
2436 }
2437
2438
2439 /* A number of useful constants */
2440
2441 EXTERN_SIM_FPU (const sim_fpu) sim_fpu_zero = {
2442   sim_fpu_class_zero,
2443 };
2444 EXTERN_SIM_FPU (const sim_fpu) sim_fpu_qnan = {
2445   sim_fpu_class_qnan,
2446 };
2447 EXTERN_SIM_FPU (const sim_fpu) sim_fpu_one = {
2448   sim_fpu_class_number, 0, IMPLICIT_1, 1
2449 };
2450 EXTERN_SIM_FPU (const sim_fpu) sim_fpu_two = {
2451   sim_fpu_class_number, 0, IMPLICIT_1, 2
2452 };
2453 EXTERN_SIM_FPU (const sim_fpu) sim_fpu_max32 = {
2454   sim_fpu_class_number, 0, LSMASK64 (NR_FRAC_GUARD, NR_GUARDS32), NORMAL_EXPMAX32
2455 };
2456 EXTERN_SIM_FPU (const sim_fpu) sim_fpu_max64 = {
2457   sim_fpu_class_number, 0, LSMASK64 (NR_FRAC_GUARD, NR_GUARDS64), NORMAL_EXPMAX64
2458 };
2459
2460
2461 /* For debugging */
2462
2463 INLINE_SIM_FPU (void)
2464 sim_fpu_print_fpu (const sim_fpu *f,
2465                    sim_fpu_print_func *print,
2466                    void *arg)
2467 {
2468   print (arg, "%s", f->sign ? "-" : "+");
2469   switch (f->class)
2470     {
2471     case sim_fpu_class_qnan:
2472       print (arg, "0.");
2473       print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg);
2474       print (arg, "*QuietNaN");
2475       break;
2476     case sim_fpu_class_snan:
2477       print (arg, "0.");
2478       print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg);
2479       print (arg, "*SignalNaN");
2480       break;
2481     case sim_fpu_class_zero:
2482       print (arg, "0.0");
2483       break;
2484     case sim_fpu_class_infinity:
2485       print (arg, "INF");
2486       break;
2487     case sim_fpu_class_number:
2488     case sim_fpu_class_denorm:
2489       print (arg, "1.");
2490       print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg);
2491       print (arg, "*2^%+-5d", f->normal_exp);
2492       ASSERT (f->fraction >= IMPLICIT_1);
2493       ASSERT (f->fraction < IMPLICIT_2);
2494     }
2495 }
2496
2497
2498 INLINE_SIM_FPU (void)
2499 sim_fpu_print_status (int status,
2500                       sim_fpu_print_func *print,
2501                       void *arg)
2502 {
2503   int i = 1;
2504   char *prefix = "";
2505   while (status >= i)
2506     {
2507       switch ((sim_fpu_status) (status & i))
2508         {
2509         case sim_fpu_status_denorm:
2510           print (arg, "%sD", prefix);
2511           break;
2512         case sim_fpu_status_invalid_snan:
2513           print (arg, "%sSNaN", prefix);
2514           break;
2515         case sim_fpu_status_invalid_qnan:
2516           print (arg, "%sQNaN", prefix);
2517           break;
2518         case sim_fpu_status_invalid_isi:
2519           print (arg, "%sISI", prefix);
2520           break;
2521         case sim_fpu_status_invalid_idi:
2522           print (arg, "%sIDI", prefix);
2523           break;
2524         case sim_fpu_status_invalid_zdz:
2525           print (arg, "%sZDZ", prefix);
2526           break;
2527         case sim_fpu_status_invalid_imz:
2528           print (arg, "%sIMZ", prefix);
2529           break;
2530         case sim_fpu_status_invalid_cvi:
2531           print (arg, "%sCVI", prefix);
2532           break;
2533         case sim_fpu_status_invalid_cmp:
2534           print (arg, "%sCMP", prefix);
2535           break;
2536         case sim_fpu_status_invalid_sqrt:
2537           print (arg, "%sSQRT", prefix);
2538           break;
2539           break;
2540         case sim_fpu_status_inexact:
2541           print (arg, "%sX", prefix);
2542           break;
2543           break;
2544         case sim_fpu_status_overflow:
2545           print (arg, "%sO", prefix);
2546           break;
2547           break;
2548         case sim_fpu_status_underflow:
2549           print (arg, "%sU", prefix);
2550           break;
2551           break;
2552         case sim_fpu_status_invalid_div0:
2553           print (arg, "%s/", prefix);
2554           break;
2555           break;
2556         case sim_fpu_status_rounded:
2557           print (arg, "%sR", prefix);
2558           break;
2559           break;
2560         }
2561       i <<= 1;
2562       prefix = ",";
2563     }
2564 }
2565
2566 #endif