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