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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gdb-7.2/] [sim/] [ppc/] [dp-bit.c] - Blame information for rev 841

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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