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