OpenCores
URL https://opencores.org/ocsvn/or1k/or1k/trunk

Subversion Repositories or1k

[/] [or1k/] [branches/] [oc/] [gdb-5.0/] [sim/] [ppc/] [dp-bit.c] - Blame information for rev 1771

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 106 markom
/* This is a software floating point library which can be used instead of
2
   the floating point routines in libgcc1.c for targets without hardware
3
   floating point.  */
4
 
5
/* Copyright (C) 1994 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
/* The intended way to use this file is to make two copies, add `#define FLOAT'
44
   to one copy, then compile both copies and add them to libgcc.a.  */
45
 
46
/* The following macros can be defined to change the behaviour of this file:
47
   FLOAT: Implement a `float', aka SFmode, fp library.  If this is not
48
     defined, then this file implements a `double', aka DFmode, fp library.
49
   FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
50
     don't include float->double conversion which requires the double library.
51
     This is useful only for machines which can't support doubles, e.g. some
52
     8-bit processors.
53
   CMPtype: Specify the type that floating point compares should return.
54
     This defaults to SItype, aka int.
55
   US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
56
     US Software goFast library.  If this is not defined, the entry points use
57
     the same names as libgcc1.c.
58
   _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
59
     two integers to the FLO_union_type.
60
   NO_NANS: Disable nan and infinity handling
61
   SMALL_MACHINE: Useful when operations on QIs and HIs are faster
62
     than on an SI */
63
 
64
#ifndef SFtype
65
typedef SFtype __attribute__ ((mode (SF)));
66
#endif
67
#ifndef DFtype
68
typedef DFtype __attribute__ ((mode (DF)));
69
#endif
70
 
71
#ifndef HItype
72
typedef int HItype __attribute__ ((mode (HI)));
73
#endif
74
#ifndef SItype
75
typedef int SItype __attribute__ ((mode (SI)));
76
#endif
77
#ifndef DItype
78
typedef int DItype __attribute__ ((mode (DI)));
79
#endif
80
 
81
/* The type of the result of a fp compare */
82
#ifndef CMPtype
83
#define CMPtype SItype
84
#endif
85
 
86
#ifndef UHItype
87
typedef unsigned int UHItype __attribute__ ((mode (HI)));
88
#endif
89
#ifndef USItype
90
typedef unsigned int USItype __attribute__ ((mode (SI)));
91
#endif
92
#ifndef UDItype
93
typedef unsigned int UDItype __attribute__ ((mode (DI)));
94
#endif
95
 
96
#define MAX_SI_INT   ((SItype) ((unsigned) (~0)>>1))
97
#define MAX_USI_INT  ((USItype) ~0)
98
 
99
 
100
#ifdef FLOAT_ONLY
101
#define NO_DI_MODE
102
#endif
103
 
104
#ifdef FLOAT
105
#       define NGARDS    7L
106
#       define GARDROUND 0x3f
107
#       define GARDMASK  0x7f
108
#       define GARDMSB   0x40
109
#       define EXPBITS 8
110
#       define EXPBIAS 127
111
#       define FRACBITS 23
112
#       define EXPMAX (0xff)
113
#       define QUIET_NAN 0x100000L
114
#       define FRAC_NBITS 32
115
#       define FRACHIGH  0x80000000L
116
#       define FRACHIGH2 0xc0000000L
117
        typedef USItype fractype;
118
        typedef UHItype halffractype;
119
        typedef SFtype FLO_type;
120
        typedef SItype intfrac;
121
 
122
#else
123
#       define PREFIXFPDP dp
124
#       define PREFIXSFDF df
125
#       define NGARDS 8L
126
#       define GARDROUND 0x7f
127
#       define GARDMASK  0xff
128
#       define GARDMSB   0x80
129
#       define EXPBITS 11
130
#       define EXPBIAS 1023
131
#       define FRACBITS 52
132
#       define EXPMAX (0x7ff)
133
#       define QUIET_NAN 0x8000000000000LL
134
#       define FRAC_NBITS 64
135
#       define FRACHIGH  0x8000000000000000LL
136
#       define FRACHIGH2 0xc000000000000000LL
137
        typedef UDItype fractype;
138
        typedef USItype halffractype;
139
        typedef DFtype FLO_type;
140
        typedef DItype intfrac;
141
#endif
142
 
143
#ifdef US_SOFTWARE_GOFAST
144
#       ifdef FLOAT
145
#               define add              fpadd
146
#               define sub              fpsub
147
#               define multiply         fpmul
148
#               define divide           fpdiv
149
#               define compare          fpcmp
150
#               define si_to_float      sitofp
151
#               define float_to_si      fptosi
152
#               define float_to_usi     fptoui
153
#               define negate           __negsf2
154
#               define sf_to_df         fptodp
155
#               define dptofp           dptofp
156
#else
157
#               define add              dpadd
158
#               define sub              dpsub
159
#               define multiply         dpmul
160
#               define divide           dpdiv
161
#               define compare          dpcmp
162
#               define si_to_float      litodp
163
#               define float_to_si      dptoli
164
#               define float_to_usi     dptoul
165
#               define negate           __negdf2
166
#               define df_to_sf         dptofp
167
#endif
168
#else
169
#       ifdef FLOAT
170
#               define add              __addsf3
171
#               define sub              __subsf3
172
#               define multiply         __mulsf3
173
#               define divide           __divsf3
174
#               define compare          __cmpsf2
175
#               define _eq_f2           __eqsf2
176
#               define _ne_f2           __nesf2
177
#               define _gt_f2           __gtsf2
178
#               define _ge_f2           __gesf2
179
#               define _lt_f2           __ltsf2
180
#               define _le_f2           __lesf2
181
#               define si_to_float      __floatsisf
182
#               define float_to_si      __fixsfsi
183
#               define float_to_usi     __fixunssfsi
184
#               define negate           __negsf2
185
#               define sf_to_df         __extendsfdf2
186
#else
187
#               define add              __adddf3
188
#               define sub              __subdf3
189
#               define multiply         __muldf3
190
#               define divide           __divdf3
191
#               define compare          __cmpdf2
192
#               define _eq_f2           __eqdf2
193
#               define _ne_f2           __nedf2
194
#               define _gt_f2           __gtdf2
195
#               define _ge_f2           __gedf2
196
#               define _lt_f2           __ltdf2
197
#               define _le_f2           __ledf2
198
#               define si_to_float      __floatsidf
199
#               define float_to_si      __fixdfsi
200
#               define float_to_usi     __fixunsdfsi
201
#               define negate           __negdf2
202
#               define df_to_sf         __truncdfsf2
203
#       endif
204
#endif
205
 
206
 
207
#ifndef INLINE
208
#define INLINE __inline__
209
#endif
210
 
211
/* Preserve the sticky-bit when shifting fractions to the right.  */
212
#define LSHIFT(a) { a = (a & 1) | (a >> 1); }
213
 
214
/* numeric parameters */
215
/* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
216
   of a float and of a double. Assumes there are only two float types.
217
   (double::FRAC_BITS+double::NGARGS-(float::FRAC_BITS-float::NGARDS))
218
 */
219
#define F_D_BITOFF (52+8-(23+7))
220
 
221
 
222
#define NORMAL_EXPMIN (-(EXPBIAS)+1)
223
#define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
224
#define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
225
 
226
/* common types */
227
 
228
typedef enum
229
{
230
  CLASS_SNAN,
231
  CLASS_QNAN,
232
  CLASS_ZERO,
233
  CLASS_NUMBER,
234
  CLASS_INFINITY
235
} fp_class_type;
236
 
237
typedef struct
238
{
239
#ifdef SMALL_MACHINE
240
  char class;
241
  unsigned char sign;
242
  short normal_exp;
243
#else
244
  fp_class_type class;
245
  unsigned int sign;
246
  int normal_exp;
247
#endif
248
 
249
  union
250
    {
251
      fractype ll;
252
      halffractype l[2];
253
    } fraction;
254
} fp_number_type;
255
 
256
typedef union
257
{
258
  FLO_type value;
259
#ifdef _DEBUG_BITFLOAT
260
  int l[2];
261
#endif
262
  struct
263
    {
264
#ifndef FLOAT_BIT_ORDER_MISMATCH
265
      unsigned int sign:1 __attribute__ ((packed));
266
      unsigned int exp:EXPBITS __attribute__ ((packed));
267
      fractype fraction:FRACBITS __attribute__ ((packed));
268
#else
269
      fractype fraction:FRACBITS __attribute__ ((packed));
270
      unsigned int exp:EXPBITS __attribute__ ((packed));
271
      unsigned int sign:1 __attribute__ ((packed));
272
#endif
273
    }
274
  bits;
275
}
276
FLO_union_type;
277
 
278
 
279
/* end of header */
280
 
281
/* IEEE "special" number predicates */
282
 
283
#ifdef NO_NANS
284
 
285
#define nan() 0
286
#define isnan(x) 0
287
#define isinf(x) 0
288
#else
289
 
290
INLINE
291
static fp_number_type *
292
nan ()
293
{
294
  static fp_number_type thenan;
295
 
296
  return &thenan;
297
}
298
 
299
INLINE
300
static int
301
isnan ( fp_number_type *  x)
302
{
303
  return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
304
}
305
 
306
INLINE
307
static int
308
isinf ( fp_number_type *  x)
309
{
310
  return x->class == CLASS_INFINITY;
311
}
312
 
313
#endif
314
 
315
INLINE
316
static int
317
iszero ( fp_number_type *  x)
318
{
319
  return x->class == CLASS_ZERO;
320
}
321
 
322
INLINE
323
static void
324
flip_sign ( fp_number_type *  x)
325
{
326
  x->sign = !x->sign;
327
}
328
 
329
static FLO_type
330
pack_d ( fp_number_type *  src)
331
{
332
  FLO_union_type dst;
333
  fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
334
 
335
  dst.bits.sign = src->sign;
336
 
337
  if (isnan (src))
338
    {
339
      dst.bits.exp = EXPMAX;
340
      dst.bits.fraction = src->fraction.ll;
341
      if (src->class == CLASS_QNAN || 1)
342
        {
343
          dst.bits.fraction |= QUIET_NAN;
344
        }
345
    }
346
  else if (isinf (src))
347
    {
348
      dst.bits.exp = EXPMAX;
349
      dst.bits.fraction = 0;
350
    }
351
  else if (iszero (src))
352
    {
353
      dst.bits.exp = 0;
354
      dst.bits.fraction = 0;
355
    }
356
  else if (fraction == 0)
357
    {
358
      dst.value = 0;
359
    }
360
  else
361
    {
362
      if (src->normal_exp < NORMAL_EXPMIN)
363
        {
364
          /* This number's exponent is too low to fit into the bits
365
             available in the number, so we'll store 0 in the exponent and
366
             shift the fraction to the right to make up for it.  */
367
 
368
          int shift = NORMAL_EXPMIN - src->normal_exp;
369
 
370
          dst.bits.exp = 0;
371
 
372
          if (shift > FRAC_NBITS - NGARDS)
373
            {
374
              /* No point shifting, since it's more that 64 out.  */
375
              fraction = 0;
376
            }
377
          else
378
            {
379
              /* Shift by the value */
380
              fraction >>= shift;
381
            }
382
          fraction >>= NGARDS;
383
          dst.bits.fraction = fraction;
384
        }
385
      else if (src->normal_exp > EXPBIAS)
386
        {
387
          dst.bits.exp = EXPMAX;
388
          dst.bits.fraction = 0;
389
        }
390
      else
391
        {
392
          dst.bits.exp = src->normal_exp + EXPBIAS;
393
          /* IF the gard bits are the all zero, but the first, then we're
394
             half way between two numbers, choose the one which makes the
395
             lsb of the answer 0.  */
396
          if ((fraction & GARDMASK) == GARDMSB)
397
            {
398
              if (fraction & (1 << NGARDS))
399
                fraction += GARDROUND + 1;
400
            }
401
          else
402
            {
403
              /* Add a one to the guards to round up */
404
              fraction += GARDROUND;
405
            }
406
          if (fraction >= IMPLICIT_2)
407
            {
408
              fraction >>= 1;
409
              dst.bits.exp += 1;
410
            }
411
          fraction >>= NGARDS;
412
          dst.bits.fraction = fraction;
413
        }
414
    }
415
  return dst.value;
416
}
417
 
418
static void
419
unpack_d (FLO_union_type * src, fp_number_type * dst)
420
{
421
  fractype fraction = src->bits.fraction;
422
 
423
  dst->sign = src->bits.sign;
424
  if (src->bits.exp == 0)
425
    {
426
      /* Hmm.  Looks like 0 */
427
      if (fraction == 0)
428
        {
429
          /* tastes like zero */
430
          dst->class = CLASS_ZERO;
431
        }
432
      else
433
        {
434
          /* Zero exponent with non zero fraction - it's denormalized,
435
             so there isn't a leading implicit one - we'll shift it so
436
             it gets one.  */
437
          dst->normal_exp = src->bits.exp - EXPBIAS + 1;
438
          fraction <<= NGARDS;
439
 
440
          dst->class = CLASS_NUMBER;
441
#if 1
442
          while (fraction < IMPLICIT_1)
443
            {
444
              fraction <<= 1;
445
              dst->normal_exp--;
446
            }
447
#endif
448
          dst->fraction.ll = fraction;
449
        }
450
    }
451
  else if (src->bits.exp == EXPMAX)
452
    {
453
      /* Huge exponent*/
454
      if (fraction == 0)
455
        {
456
          /* Attached to a zero fraction - means infinity */
457
          dst->class = CLASS_INFINITY;
458
        }
459
      else
460
        {
461
          /* Non zero fraction, means nan */
462
          if (dst->sign)
463
            {
464
              dst->class = CLASS_SNAN;
465
            }
466
          else
467
            {
468
              dst->class = CLASS_QNAN;
469
            }
470
          /* Keep the fraction part as the nan number */
471
          dst->fraction.ll = fraction;
472
        }
473
    }
474
  else
475
    {
476
      /* Nothing strange about this number */
477
      dst->normal_exp = src->bits.exp - EXPBIAS;
478
      dst->class = CLASS_NUMBER;
479
      dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
480
    }
481
}
482
 
483
static fp_number_type *
484
_fpadd_parts (fp_number_type * a,
485
              fp_number_type * b,
486
              fp_number_type * tmp)
487
{
488
  intfrac tfraction;
489
 
490
  /* Put commonly used fields in local variables.  */
491
  int a_normal_exp;
492
  int b_normal_exp;
493
  fractype a_fraction;
494
  fractype b_fraction;
495
 
496
  if (isnan (a))
497
    {
498
      return a;
499
    }
500
  if (isnan (b))
501
    {
502
      return b;
503
    }
504
  if (isinf (a))
505
    {
506
      /* Adding infinities with opposite signs yields a NaN.  */
507
      if (isinf (b) && a->sign != b->sign)
508
        return nan ();
509
      return a;
510
    }
511
  if (isinf (b))
512
    {
513
      return b;
514
    }
515
  if (iszero (b))
516
    {
517
      return a;
518
    }
519
  if (iszero (a))
520
    {
521
      return b;
522
    }
523
 
524
  /* Got two numbers. shift the smaller and increment the exponent till
525
     they're the same */
526
  {
527
    int diff;
528
 
529
    a_normal_exp = a->normal_exp;
530
    b_normal_exp = b->normal_exp;
531
    a_fraction = a->fraction.ll;
532
    b_fraction = b->fraction.ll;
533
 
534
    diff = a_normal_exp - b_normal_exp;
535
 
536
    if (diff < 0)
537
      diff = -diff;
538
    if (diff < FRAC_NBITS)
539
      {
540
        /* ??? This does shifts one bit at a time.  Optimize.  */
541
        while (a_normal_exp > b_normal_exp)
542
          {
543
            b_normal_exp++;
544
            LSHIFT (b_fraction);
545
          }
546
        while (b_normal_exp > a_normal_exp)
547
          {
548
            a_normal_exp++;
549
            LSHIFT (a_fraction);
550
          }
551
      }
552
    else
553
      {
554
        /* Somethings's up.. choose the biggest */
555
        if (a_normal_exp > b_normal_exp)
556
          {
557
            b_normal_exp = a_normal_exp;
558
            b_fraction = 0;
559
          }
560
        else
561
          {
562
            a_normal_exp = b_normal_exp;
563
            a_fraction = 0;
564
          }
565
      }
566
  }
567
 
568
  if (a->sign != b->sign)
569
    {
570
      if (a->sign)
571
        {
572
          tfraction = -a_fraction + b_fraction;
573
        }
574
      else
575
        {
576
          tfraction = a_fraction - b_fraction;
577
        }
578
      if (tfraction > 0)
579
        {
580
          tmp->sign = 0;
581
          tmp->normal_exp = a_normal_exp;
582
          tmp->fraction.ll = tfraction;
583
        }
584
      else
585
        {
586
          tmp->sign = 1;
587
          tmp->normal_exp = a_normal_exp;
588
          tmp->fraction.ll = -tfraction;
589
        }
590
      /* and renormalize it */
591
 
592
      while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
593
        {
594
          tmp->fraction.ll <<= 1;
595
          tmp->normal_exp--;
596
        }
597
    }
598
  else
599
    {
600
      tmp->sign = a->sign;
601
      tmp->normal_exp = a_normal_exp;
602
      tmp->fraction.ll = a_fraction + b_fraction;
603
    }
604
  tmp->class = CLASS_NUMBER;
605
  /* Now the fraction is added, we have to shift down to renormalize the
606
     number */
607
 
608
  if (tmp->fraction.ll >= IMPLICIT_2)
609
    {
610
      LSHIFT (tmp->fraction.ll);
611
      tmp->normal_exp++;
612
    }
613
  return tmp;
614
 
615
}
616
 
617
FLO_type
618
add (FLO_type arg_a, FLO_type arg_b)
619
{
620
  fp_number_type a;
621
  fp_number_type b;
622
  fp_number_type tmp;
623
  fp_number_type *res;
624
 
625
  unpack_d ((FLO_union_type *) & arg_a, &a);
626
  unpack_d ((FLO_union_type *) & arg_b, &b);
627
 
628
  res = _fpadd_parts (&a, &b, &tmp);
629
 
630
  return pack_d (res);
631
}
632
 
633
FLO_type
634
sub (FLO_type arg_a, FLO_type arg_b)
635
{
636
  fp_number_type a;
637
  fp_number_type b;
638
  fp_number_type tmp;
639
  fp_number_type *res;
640
 
641
  unpack_d ((FLO_union_type *) & arg_a, &a);
642
  unpack_d ((FLO_union_type *) & arg_b, &b);
643
 
644
  b.sign ^= 1;
645
 
646
  res = _fpadd_parts (&a, &b, &tmp);
647
 
648
  return pack_d (res);
649
}
650
 
651
static fp_number_type *
652
_fpmul_parts ( fp_number_type *  a,
653
               fp_number_type *  b,
654
               fp_number_type * tmp)
655
{
656
  fractype low = 0;
657
  fractype high = 0;
658
 
659
  if (isnan (a))
660
    {
661
      a->sign = a->sign != b->sign;
662
      return a;
663
    }
664
  if (isnan (b))
665
    {
666
      b->sign = a->sign != b->sign;
667
      return b;
668
    }
669
  if (isinf (a))
670
    {
671
      if (iszero (b))
672
        return nan ();
673
      a->sign = a->sign != b->sign;
674
      return a;
675
    }
676
  if (isinf (b))
677
    {
678
      if (iszero (a))
679
        {
680
          return nan ();
681
        }
682
      b->sign = a->sign != b->sign;
683
      return b;
684
    }
685
  if (iszero (a))
686
    {
687
      a->sign = a->sign != b->sign;
688
      return a;
689
    }
690
  if (iszero (b))
691
    {
692
      b->sign = a->sign != b->sign;
693
      return b;
694
    }
695
 
696
  /* Calculate the mantissa by multiplying both 64bit numbers to get a
697
     128 bit number */
698
  {
699
    fractype x = a->fraction.ll;
700
    fractype ylow = b->fraction.ll;
701
    fractype yhigh = 0;
702
    int bit;
703
 
704
#if defined(NO_DI_MODE)
705
    {
706
      /* ??? This does multiplies one bit at a time.  Optimize.  */
707
      for (bit = 0; bit < FRAC_NBITS; bit++)
708
        {
709
          int carry;
710
 
711
          if (x & 1)
712
            {
713
              carry = (low += ylow) < ylow;
714
              high += yhigh + carry;
715
            }
716
          yhigh <<= 1;
717
          if (ylow & FRACHIGH)
718
            {
719
              yhigh |= 1;
720
            }
721
          ylow <<= 1;
722
          x >>= 1;
723
        }
724
    }
725
#elif defined(FLOAT) 
726
    {
727
      /* Multiplying two 32 bit numbers to get a 64 bit number  on
728
        a machine with DI, so we're safe */
729
 
730
      DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
731
 
732
      high = answer >> 32;
733
      low = answer;
734
    }
735
#else
736
    /* Doing a 64*64 to 128 */
737
    {
738
      UDItype nl = a->fraction.ll & 0xffffffff;
739
      UDItype nh = a->fraction.ll >> 32;
740
      UDItype ml = b->fraction.ll & 0xffffffff;
741
      UDItype mh = b->fraction.ll >>32;
742
      UDItype pp_ll = ml * nl;
743
      UDItype pp_hl = mh * nl;
744
      UDItype pp_lh = ml * nh;
745
      UDItype pp_hh = mh * nh;
746
      UDItype res2 = 0;
747
      UDItype res0 = 0;
748
      UDItype ps_hh__ = pp_hl + pp_lh;
749
      if (ps_hh__ < pp_hl)
750
        res2 += 0x100000000LL;
751
      pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
752
      res0 = pp_ll + pp_hl;
753
      if (res0 < pp_ll)
754
        res2++;
755
      res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
756
      high = res2;
757
      low = res0;
758
    }
759
#endif
760
  }
761
 
762
  tmp->normal_exp = a->normal_exp + b->normal_exp;
763
  tmp->sign = a->sign != b->sign;
764
#ifdef FLOAT
765
  tmp->normal_exp += 2;         /* ??????????????? */
766
#else
767
  tmp->normal_exp += 4;         /* ??????????????? */
768
#endif
769
  while (high >= IMPLICIT_2)
770
    {
771
      tmp->normal_exp++;
772
      if (high & 1)
773
        {
774
          low >>= 1;
775
          low |= FRACHIGH;
776
        }
777
      high >>= 1;
778
    }
779
  while (high < IMPLICIT_1)
780
    {
781
      tmp->normal_exp--;
782
 
783
      high <<= 1;
784
      if (low & FRACHIGH)
785
        high |= 1;
786
      low <<= 1;
787
    }
788
  /* rounding is tricky. if we only round if it won't make us round later. */
789
#if 0
790
  if (low & FRACHIGH2)
791
    {
792
      if (((high & GARDMASK) != GARDMSB)
793
          && (((high + 1) & GARDMASK) == GARDMSB))
794
        {
795
          /* don't round, it gets done again later. */
796
        }
797
      else
798
        {
799
          high++;
800
        }
801
    }
802
#endif
803
  if ((high & GARDMASK) == GARDMSB)
804
    {
805
      if (high & (1 << NGARDS))
806
        {
807
          /* half way, so round to even */
808
          high += GARDROUND + 1;
809
        }
810
      else if (low)
811
        {
812
          /* but we really weren't half way */
813
          high += GARDROUND + 1;
814
        }
815
    }
816
  tmp->fraction.ll = high;
817
  tmp->class = CLASS_NUMBER;
818
  return tmp;
819
}
820
 
821
FLO_type
822
multiply (FLO_type arg_a, FLO_type arg_b)
823
{
824
  fp_number_type a;
825
  fp_number_type b;
826
  fp_number_type tmp;
827
  fp_number_type *res;
828
 
829
  unpack_d ((FLO_union_type *) & arg_a, &a);
830
  unpack_d ((FLO_union_type *) & arg_b, &b);
831
 
832
  res = _fpmul_parts (&a, &b, &tmp);
833
 
834
  return pack_d (res);
835
}
836
 
837
static fp_number_type *
838
_fpdiv_parts (fp_number_type * a,
839
              fp_number_type * b,
840
              fp_number_type * tmp)
841
{
842
  fractype low = 0;
843
  fractype high = 0;
844
  fractype r0, r1, y0, y1, bit;
845
  fractype q;
846
  fractype numerator;
847
  fractype denominator;
848
  fractype quotient;
849
  fractype remainder;
850
 
851
  if (isnan (a))
852
    {
853
      return a;
854
    }
855
  if (isnan (b))
856
    {
857
      return b;
858
    }
859
  if (isinf (a) || iszero (a))
860
    {
861
      if (a->class == b->class)
862
        return nan ();
863
      return a;
864
    }
865
  a->sign = a->sign ^ b->sign;
866
 
867
  if (isinf (b))
868
    {
869
      a->fraction.ll = 0;
870
      a->normal_exp = 0;
871
      return a;
872
    }
873
  if (iszero (b))
874
    {
875
      a->class = CLASS_INFINITY;
876
      return b;
877
    }
878
 
879
  /* Calculate the mantissa by multiplying both 64bit numbers to get a
880
     128 bit number */
881
  {
882
    int carry;
883
    intfrac d0, d1;             /* weren't unsigned before ??? */
884
 
885
    /* quotient =
886
       ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
887
     */
888
 
889
    a->normal_exp = a->normal_exp - b->normal_exp;
890
    numerator = a->fraction.ll;
891
    denominator = b->fraction.ll;
892
 
893
    if (numerator < denominator)
894
      {
895
        /* Fraction will be less than 1.0 */
896
        numerator *= 2;
897
        a->normal_exp--;
898
      }
899
    bit = IMPLICIT_1;
900
    quotient = 0;
901
    /* ??? Does divide one bit at a time.  Optimize.  */
902
    while (bit)
903
      {
904
        if (numerator >= denominator)
905
          {
906
            quotient |= bit;
907
            numerator -= denominator;
908
          }
909
        bit >>= 1;
910
        numerator *= 2;
911
      }
912
 
913
    if ((quotient & GARDMASK) == GARDMSB)
914
      {
915
        if (quotient & (1 << NGARDS))
916
          {
917
            /* half way, so round to even */
918
            quotient += GARDROUND + 1;
919
          }
920
        else if (numerator)
921
          {
922
            /* but we really weren't half way, more bits exist */
923
            quotient += GARDROUND + 1;
924
          }
925
      }
926
 
927
    a->fraction.ll = quotient;
928
    return (a);
929
  }
930
}
931
 
932
FLO_type
933
divide (FLO_type arg_a, FLO_type arg_b)
934
{
935
  fp_number_type a;
936
  fp_number_type b;
937
  fp_number_type tmp;
938
  fp_number_type *res;
939
 
940
  unpack_d ((FLO_union_type *) & arg_a, &a);
941
  unpack_d ((FLO_union_type *) & arg_b, &b);
942
 
943
  res = _fpdiv_parts (&a, &b, &tmp);
944
 
945
  return pack_d (res);
946
}
947
 
948
/* according to the demo, fpcmp returns a comparison with 0... thus
949
   a<b -> -1
950
   a==b -> 0
951
   a>b -> +1
952
 */
953
 
954
static int
955
_fpcmp_parts (fp_number_type * a, fp_number_type * b)
956
{
957
#if 0
958
  /* either nan -> unordered. Must be checked outside of this routine. */
959
  if (isnan (a) && isnan (b))
960
    {
961
      return 1;                 /* still unordered! */
962
    }
963
#endif
964
 
965
  if (isnan (a) || isnan (b))
966
    {
967
      return 1;                 /* how to indicate unordered compare? */
968
    }
969
  if (isinf (a) && isinf (b))
970
    {
971
      /* +inf > -inf, but +inf != +inf */
972
      /* b    \a| +inf(0)| -inf(1)
973
       ______\+--------+--------
974
       +inf(0)| a==b(0)| a<b(-1)
975
       -------+--------+--------
976
       -inf(1)| a>b(1) | a==b(0)
977
       -------+--------+--------
978
       So since unordered must be non zero, just line up the columns...
979
       */
980
      return b->sign - a->sign;
981
    }
982
  /* but not both... */
983
  if (isinf (a))
984
    {
985
      return a->sign ? -1 : 1;
986
    }
987
  if (isinf (b))
988
    {
989
      return b->sign ? 1 : -1;
990
    }
991
  if (iszero (a) && iszero (b))
992
    {
993
      return 0;
994
    }
995
  if (iszero (a))
996
    {
997
      return b->sign ? 1 : -1;
998
    }
999
  if (iszero (b))
1000
    {
1001
      return a->sign ? -1 : 1;
1002
    }
1003
  /* now both are "normal". */
1004
  if (a->sign != b->sign)
1005
    {
1006
      /* opposite signs */
1007
      return a->sign ? -1 : 1;
1008
    }
1009
  /* same sign; exponents? */
1010
  if (a->normal_exp > b->normal_exp)
1011
    {
1012
      return a->sign ? -1 : 1;
1013
    }
1014
  if (a->normal_exp < b->normal_exp)
1015
    {
1016
      return a->sign ? 1 : -1;
1017
    }
1018
  /* same exponents; check size. */
1019
  if (a->fraction.ll > b->fraction.ll)
1020
    {
1021
      return a->sign ? -1 : 1;
1022
    }
1023
  if (a->fraction.ll < b->fraction.ll)
1024
    {
1025
      return a->sign ? 1 : -1;
1026
    }
1027
  /* after all that, they're equal. */
1028
  return 0;
1029
}
1030
 
1031
CMPtype
1032
compare (FLO_type arg_a, FLO_type arg_b)
1033
{
1034
  fp_number_type a;
1035
  fp_number_type b;
1036
 
1037
  unpack_d ((FLO_union_type *) & arg_a, &a);
1038
  unpack_d ((FLO_union_type *) & arg_b, &b);
1039
 
1040
  return _fpcmp_parts (&a, &b);
1041
}
1042
 
1043
#ifndef US_SOFTWARE_GOFAST
1044
 
1045
/* These should be optimized for their specific tasks someday.  */
1046
 
1047
CMPtype
1048
_eq_f2 (FLO_type arg_a, FLO_type arg_b)
1049
{
1050
  fp_number_type a;
1051
  fp_number_type b;
1052
 
1053
  unpack_d ((FLO_union_type *) & arg_a, &a);
1054
  unpack_d ((FLO_union_type *) & arg_b, &b);
1055
 
1056
  if (isnan (&a) || isnan (&b))
1057
    return 1;                   /* false, truth == 0 */
1058
 
1059
  return _fpcmp_parts (&a, &b) ;
1060
}
1061
 
1062
CMPtype
1063
_ne_f2 (FLO_type arg_a, FLO_type arg_b)
1064
{
1065
  fp_number_type a;
1066
  fp_number_type b;
1067
 
1068
  unpack_d ((FLO_union_type *) & arg_a, &a);
1069
  unpack_d ((FLO_union_type *) & arg_b, &b);
1070
 
1071
  if (isnan (&a) || isnan (&b))
1072
    return 1;                   /* true, truth != 0 */
1073
 
1074
  return  _fpcmp_parts (&a, &b) ;
1075
}
1076
 
1077
CMPtype
1078
_gt_f2 (FLO_type arg_a, FLO_type arg_b)
1079
{
1080
  fp_number_type a;
1081
  fp_number_type b;
1082
 
1083
  unpack_d ((FLO_union_type *) & arg_a, &a);
1084
  unpack_d ((FLO_union_type *) & arg_b, &b);
1085
 
1086
  if (isnan (&a) || isnan (&b))
1087
    return -1;                  /* false, truth > 0 */
1088
 
1089
  return _fpcmp_parts (&a, &b);
1090
}
1091
 
1092
CMPtype
1093
_ge_f2 (FLO_type arg_a, FLO_type arg_b)
1094
{
1095
  fp_number_type a;
1096
  fp_number_type b;
1097
 
1098
  unpack_d ((FLO_union_type *) & arg_a, &a);
1099
  unpack_d ((FLO_union_type *) & arg_b, &b);
1100
 
1101
  if (isnan (&a) || isnan (&b))
1102
    return -1;                  /* false, truth >= 0 */
1103
  return _fpcmp_parts (&a, &b) ;
1104
}
1105
 
1106
CMPtype
1107
_lt_f2 (FLO_type arg_a, FLO_type arg_b)
1108
{
1109
  fp_number_type a;
1110
  fp_number_type b;
1111
 
1112
  unpack_d ((FLO_union_type *) & arg_a, &a);
1113
  unpack_d ((FLO_union_type *) & arg_b, &b);
1114
 
1115
  if (isnan (&a) || isnan (&b))
1116
    return 1;                   /* false, truth < 0 */
1117
 
1118
  return _fpcmp_parts (&a, &b);
1119
}
1120
 
1121
CMPtype
1122
_le_f2 (FLO_type arg_a, FLO_type arg_b)
1123
{
1124
  fp_number_type a;
1125
  fp_number_type b;
1126
 
1127
  unpack_d ((FLO_union_type *) & arg_a, &a);
1128
  unpack_d ((FLO_union_type *) & arg_b, &b);
1129
 
1130
  if (isnan (&a) || isnan (&b))
1131
    return 1;                   /* false, truth <= 0 */
1132
 
1133
  return _fpcmp_parts (&a, &b) ;
1134
}
1135
 
1136
#endif /* ! US_SOFTWARE_GOFAST */
1137
 
1138
FLO_type
1139
si_to_float (SItype arg_a)
1140
{
1141
  fp_number_type in;
1142
 
1143
  in.class = CLASS_NUMBER;
1144
  in.sign = arg_a < 0;
1145
  if (!arg_a)
1146
    {
1147
      in.class = CLASS_ZERO;
1148
    }
1149
  else
1150
    {
1151
      in.normal_exp = FRACBITS + NGARDS;
1152
      if (in.sign)
1153
        {
1154
          /* Special case for minint, since there is no +ve integer
1155
             representation for it */
1156
          if (arg_a == 0x80000000)
1157
            {
1158
              return -2147483648.0;
1159
            }
1160
          in.fraction.ll = (-arg_a);
1161
        }
1162
      else
1163
        in.fraction.ll = arg_a;
1164
 
1165
      while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1166
        {
1167
          in.fraction.ll <<= 1;
1168
          in.normal_exp -= 1;
1169
        }
1170
    }
1171
  return pack_d (&in);
1172
}
1173
 
1174
SItype
1175
float_to_si (FLO_type arg_a)
1176
{
1177
  fp_number_type a;
1178
  SItype tmp;
1179
 
1180
  unpack_d ((FLO_union_type *) & arg_a, &a);
1181
  if (iszero (&a))
1182
    return 0;
1183
  if (isnan (&a))
1184
    return 0;
1185
  /* get reasonable MAX_SI_INT... */
1186
  if (isinf (&a))
1187
    return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
1188
  /* it is a number, but a small one */
1189
  if (a.normal_exp < 0)
1190
    return 0;
1191
  if (a.normal_exp > 30)
1192
    return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1193
  tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1194
  return a.sign ? (-tmp) : (tmp);
1195
}
1196
 
1197
#ifdef US_SOFTWARE_GOFAST
1198
/* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1199
   we also define them for GOFAST because the ones in libgcc2.c have the
1200
   wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1201
   out of libgcc2.c.  We can't define these here if not GOFAST because then
1202
   there'd be duplicate copies.  */
1203
 
1204
USItype
1205
float_to_usi (FLO_type arg_a)
1206
{
1207
  fp_number_type a;
1208
 
1209
  unpack_d ((FLO_union_type *) & arg_a, &a);
1210
  if (iszero (&a))
1211
    return 0;
1212
  if (isnan (&a))
1213
    return 0;
1214
  /* get reasonable MAX_USI_INT... */
1215
  if (isinf (&a))
1216
    return a.sign ? MAX_USI_INT : 0;
1217
  /* it is a negative number */
1218
  if (a.sign)
1219
    return 0;
1220
  /* it is a number, but a small one */
1221
  if (a.normal_exp < 0)
1222
    return 0;
1223
  if (a.normal_exp > 31)
1224
    return MAX_USI_INT;
1225
  else if (a.normal_exp > (FRACBITS + NGARDS))
1226
    return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
1227
  else
1228
    return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1229
}
1230
#endif
1231
 
1232
FLO_type
1233
negate (FLO_type arg_a)
1234
{
1235
  fp_number_type a;
1236
 
1237
  unpack_d ((FLO_union_type *) & arg_a, &a);
1238
  flip_sign (&a);
1239
  return pack_d (&a);
1240
}
1241
 
1242
#ifdef FLOAT
1243
 
1244
SFtype
1245
__make_fp(fp_class_type class,
1246
             unsigned int sign,
1247
             int exp,
1248
             USItype frac)
1249
{
1250
  fp_number_type in;
1251
 
1252
  in.class = class;
1253
  in.sign = sign;
1254
  in.normal_exp = exp;
1255
  in.fraction.ll = frac;
1256
  return pack_d (&in);
1257
}
1258
 
1259
#ifndef FLOAT_ONLY
1260
 
1261
/* This enables one to build an fp library that supports float but not double.
1262
   Otherwise, we would get an undefined reference to __make_dp.
1263
   This is needed for some 8-bit ports that can't handle well values that
1264
   are 8-bytes in size, so we just don't support double for them at all.  */
1265
 
1266
extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1267
 
1268
DFtype
1269
sf_to_df (SFtype arg_a)
1270
{
1271
  fp_number_type in;
1272
 
1273
  unpack_d ((FLO_union_type *) & arg_a, &in);
1274
  return __make_dp (in.class, in.sign, in.normal_exp,
1275
                    ((UDItype) in.fraction.ll) << F_D_BITOFF);
1276
}
1277
 
1278
#endif
1279
#endif
1280
 
1281
#ifndef FLOAT
1282
 
1283
extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1284
 
1285
DFtype
1286
__make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1287
{
1288
  fp_number_type in;
1289
 
1290
  in.class = class;
1291
  in.sign = sign;
1292
  in.normal_exp = exp;
1293
  in.fraction.ll = frac;
1294
  return pack_d (&in);
1295
}
1296
 
1297
SFtype
1298
df_to_sf (DFtype arg_a)
1299
{
1300
  fp_number_type in;
1301
 
1302
  unpack_d ((FLO_union_type *) & arg_a, &in);
1303
  return __make_fp (in.class, in.sign, in.normal_exp,
1304
                    in.fraction.ll >> F_D_BITOFF);
1305
}
1306
 
1307
#endif

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.