sim/ppc/dp-bit.c

   1 /* 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