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