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