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