libgcc/config/m68k/fpgnulib.c

   1 /* This is a stripped down version of floatlib.c.  It supplies only those
   2    functions which exist in libgcc, but for which there is not assembly
   3    language versions in m68k/lb1sf68.S.
   4
   5    It also includes simplistic support for extended floats (by working in
   6    double precision).  You must compile this file again with -DEXTFLOAT
   7    to get this support.  */
   8
   9 /*
  10 ** gnulib support for software floating point.
  11 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
  12 ** Permission is granted to do *anything* you want with this file,
  13 ** commercial or otherwise, provided this message remains intact.  So there!
  14 ** I would appreciate receiving any updates/patches/changes that anyone
  15 ** makes, and am willing to be the repository for said changes (am I
  16 ** making a big mistake?).
  17 **
  18 ** Pat Wood
  19 ** Pipeline Associates, Inc.
  20 ** pipeline!phw@motown.com or
  21 ** sun!pipeline!phw or
  22 ** uunet!motown!pipeline!phw
  23 **
  24 ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
  25 ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
  26 **                  -- fixed problems with adding and subtracting zero
  27 **                  -- fixed rounding in truncdfsf2
  28 **                  -- fixed SWAP define and tested on 386
  29 */
  30
  31 /*
  32 ** The following are routines that replace the gnulib soft floating point
  33 ** routines that are called automatically when -msoft-float is selected.
  34 ** The support single and double precision IEEE format, with provisions
  35 ** for byte-swapped machines (tested on 386).  Some of the double-precision
  36 ** routines work at full precision, but most of the hard ones simply punt
  37 ** and call the single precision routines, producing a loss of accuracy.
  38 ** long long support is not assumed or included.
  39 ** Overall accuracy is close to IEEE (actually 68882) for single-precision
  40 ** arithmetic.  I think there may still be a 1 in 1000 chance of a bit
  41 ** being rounded the wrong way during a multiply.  I'm not fussy enough to
  42 ** bother with it, but if anyone is, knock yourself out.
  43 **
  44 ** Efficiency has only been addressed where it was obvious that something
  45 ** would make a big difference.  Anyone who wants to do this right for
  46 ** best speed should go in and rewrite in assembler.
  47 **
  48 ** I have tested this only on a 68030 workstation and 386/ix integrated
  49 ** in with -msoft-float.
  50 */
  51
  52 /* the following deal with IEEE single-precision numbers */
  53 #define EXCESS          126L
  54 #define SIGNBIT         0x80000000L
  55 #define HIDDEN          (1L << 23L)
  56 #define SIGN(fp)        ((fp) & SIGNBIT)
  57 #define EXP(fp)         (((fp) >> 23L) & 0xFF)
  58 #define MANT(fp)        (((fp) & 0x7FFFFFL) | HIDDEN)
  59 #define PACK(s,e,m)     ((s) | ((e) << 23L) | (m))
  60
  61 /* the following deal with IEEE double-precision numbers */
  62 #define EXCESSD         1022L
  63 #define HIDDEND         (1L << 20L)
  64 #define EXPDBITS        11
  65 #define EXPDMASK        0x7FFL
  66 #define EXPD(fp)        (((fp.l.upper) >> 20L) & 0x7FFL)
  67 #define SIGND(fp)       ((fp.l.upper) & SIGNBIT)
  68 #define MANTD(fp)       (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
  69                                 (fp.l.lower >> 22))
  70 #define MANTDMASK       0xFFFFFL /* mask of upper part */
  71
  72 /* the following deal with IEEE extended-precision numbers */
  73 #define EXCESSX         16382L
  74 #define HIDDENX         (1L << 31L)
  75 #define EXPXBITS        15
  76 #define EXPXMASK        0x7FFF
  77 #define EXPX(fp)        (((fp.l.upper) >> 16) & EXPXMASK)
  78 #define SIGNX(fp)       ((fp.l.upper) & SIGNBIT)
  79 #define MANTXMASK       0x7FFFFFFFL /* mask of upper part */
  80
  81 union double_long
  82 {
  83   double d;
  84   struct {
  85       long upper;
  86       unsigned long lower;
  87     } l;
  88 };
  89
  90 union float_long {
  91   float f;
  92   long l;
  93 };
  94
  95 union long_double_long
  96 {
  97   long double ld;
  98   struct
  99     {
 100       long upper;
 101       unsigned long middle;
 102       unsigned long lower;
 103     } l;
 104 };
 105 \f
 106 #ifndef EXTFLOAT
 107
 108 int
 109 __unordsf2(float a, float b)
 110 {
 111   union float_long fl;
 112
 113   fl.f = a;
 114   if (EXP(fl.l) == EXP(~0u) && (MANT(fl.l) & ~HIDDEN) != 0)
 115     return 1;
 116   fl.f = b;
 117   if (EXP(fl.l) == EXP(~0u) && (MANT(fl.l) & ~HIDDEN) != 0)
 118     return 1;
 119   return 0;
 120 }
 121
 122 int
 123 __unorddf2(double a, double b)
 124 {
 125   union double_long dl;
 126
 127   dl.d = a;
 128   if (EXPD(dl) == EXPDMASK
 129       && ((dl.l.upper & MANTDMASK) != 0 || dl.l.lower != 0))
 130     return 1;
 131   dl.d = b;
 132   if (EXPD(dl) == EXPDMASK
 133       && ((dl.l.upper & MANTDMASK) != 0 || dl.l.lower != 0))
 134     return 1;
 135   return 0;
 136 }
 137
 138 /* convert unsigned int to double */
 139 double
 140 __floatunsidf (unsigned long a1)
 141 {
 142   long exp = 32 + EXCESSD;
 143   union double_long dl;
 144
 145   if (!a1)
 146     {
 147       dl.l.upper = dl.l.lower = 0;
 148       return dl.d;
 149     }
 150
 151   while (a1 < 0x2000000L)
 152     {
 153       a1 <<= 4;
 154       exp -= 4;
 155     }
 156
 157   while (a1 < 0x80000000L)
 158     {
 159       a1 <<= 1;
 160       exp--;
 161     }
 162
 163   /* pack up and go home */
 164   dl.l.upper = exp << 20L;
 165   dl.l.upper |= (a1 >> 11L) & ~HIDDEND;
 166   dl.l.lower = a1 << 21L;
 167
 168   return dl.d;
 169 }
 170
 171 /* convert int to double */
 172 double
 173 __floatsidf (long a1)
 174 {
 175   long sign = 0, exp = 31 + EXCESSD;
 176   union double_long dl;
 177
 178   if (!a1)
 179     {
 180       dl.l.upper = dl.l.lower = 0;
 181       return dl.d;
 182     }
 183
 184   if (a1 < 0)
 185     {
 186       sign = SIGNBIT;
 187       a1 = (long)-(unsigned long)a1;
 188       if (a1 < 0)
 189         {
 190           dl.l.upper = SIGNBIT | ((32 + EXCESSD) << 20L);
 191           dl.l.lower = 0;
 192           return dl.d;
 193         }
 194     }
 195
 196   while (a1 < 0x1000000L)
 197     {
 198       a1 <<= 4;
 199       exp -= 4;
 200     }
 201
 202   while (a1 < 0x40000000L)
 203     {
 204       a1 <<= 1;
 205       exp--;
 206     }
 207
 208   /* pack up and go home */
 209   dl.l.upper = sign;
 210   dl.l.upper |= exp << 20L;
 211   dl.l.upper |= (a1 >> 10L) & ~HIDDEND;
 212   dl.l.lower = a1 << 22L;
 213
 214   return dl.d;
 215 }
 216
 217 /* convert unsigned int to float */
 218 float
 219 __floatunsisf (unsigned long l)
 220 {
 221   double foo = __floatunsidf (l);
 222   return foo;
 223 }
 224
 225 /* convert int to float */
 226 float
 227 __floatsisf (long l)
 228 {
 229   double foo = __floatsidf (l);
 230   return foo;
 231 }
 232
 233 /* convert float to double */
 234 double
 235 __extendsfdf2 (float a1)
 236 {
 237   register union float_long fl1;
 238   register union double_long dl;
 239   register long exp;
 240   register long mant;
 241
 242   fl1.f = a1;
 243
 244   dl.l.upper = SIGN (fl1.l);
 245   if ((fl1.l & ~SIGNBIT) == 0)
 246     {
 247       dl.l.lower = 0;
 248       return dl.d;
 249     }
 250
 251   exp = EXP(fl1.l);
 252   mant = MANT (fl1.l) & ~HIDDEN;
 253   if (exp == 0)
 254     {
 255       /* Denormal.  */
 256       exp = 1;
 257       while (!(mant & HIDDEN))
 258         {
 259           mant <<= 1;
 260           exp--;
 261         }
 262       mant &= ~HIDDEN;
 263     }
 264   exp = exp - EXCESS + EXCESSD;
 265   dl.l.upper |= exp << 20;
 266   dl.l.upper |= mant >> 3;
 267   dl.l.lower = mant << 29;
 268
 269   return dl.d;
 270 }
 271
 272 /* convert double to float */
 273 float
 274 __truncdfsf2 (double a1)
 275 {
 276   register long exp;
 277   register long mant;
 278   register union float_long fl;
 279   register union double_long dl1;
 280   int sticky;
 281   int shift;
 282
 283   dl1.d = a1;
 284
 285   if ((dl1.l.upper & ~SIGNBIT) == 0 && !dl1.l.lower)
 286     {
 287       fl.l = SIGND(dl1);
 288       return fl.f;
 289     }
 290
 291   exp = EXPD (dl1) - EXCESSD + EXCESS;
 292
 293   sticky = dl1.l.lower & ((1 << 22) - 1);
 294   mant = MANTD (dl1);
 295   /* shift double mantissa 6 bits so we can round */
 296   sticky |= mant & ((1 << 6) - 1);
 297   mant >>= 6;
 298
 299   /* Check for underflow and denormals.  */
 300   if (exp <= 0)
 301     {
 302       if (exp < -24)
 303         {
 304           sticky |= mant;
 305           mant = 0;
 306         }
 307       else
 308         {
 309           sticky |= mant & ((1 << (1 - exp)) - 1);
 310           mant >>= 1 - exp;
 311         }
 312       exp = 0;
 313     }
 314
 315   /* now round */
 316   shift = 1;
 317   if ((mant & 1) && (sticky || (mant & 2)))
 318     {
 319       int rounding = exp ? 2 : 1;
 320
 321       mant += 1;
 322
 323       /* did the round overflow? */
 324       if (mant >= (HIDDEN << rounding))
 325         {
 326           exp++;
 327           shift = rounding;
 328         }
 329     }
 330   /* shift down */
 331   mant >>= shift;
 332
 333   mant &= ~HIDDEN;
 334
 335   /* pack up and go home */
 336   fl.l = PACK (SIGND (dl1), exp, mant);
 337   return (fl.f);
 338 }
 339
 340 /* convert double to int */
 341 long
 342 __fixdfsi (double a1)
 343 {
 344   register union double_long dl1;
 345   register long exp;
 346   register long l;
 347
 348   dl1.d = a1;
 349
 350   if (!dl1.l.upper && !dl1.l.lower)
 351     return 0;
 352
 353   exp = EXPD (dl1) - EXCESSD - 31;
 354   l = MANTD (dl1);
 355
 356   if (exp > 0)
 357     {
 358       /* Return largest integer.  */
 359       return SIGND (dl1) ? 0x80000000L : 0x7fffffffL;
 360     }
 361
 362   if (exp <= -32)
 363     return 0;
 364
 365   /* shift down until exp = 0 */
 366   if (exp < 0)
 367     l >>= -exp;
 368
 369   return (SIGND (dl1) ? -l : l);
 370 }
 371
 372 /* convert float to int */
 373 long
 374 __fixsfsi (float a1)
 375 {
 376   double foo = a1;
 377   return __fixdfsi (foo);
 378 }
 379 \f
 380 #else /* EXTFLOAT */
 381
 382 /* We do not need these routines for coldfire, as it has no extended
 383    float format. */
 384 #if !defined (__mcoldfire__)
 385
 386 /* Primitive extended precision floating point support.
 387
 388    We assume all numbers are normalized, don't do any rounding, etc.  */
 389
 390 /* Prototypes for the above in case we use them.  */
 391 double __floatunsidf (unsigned long);
 392 double __floatsidf (long);
 393 float __floatsisf (long);
 394 double __extendsfdf2 (float);
 395 float __truncdfsf2 (double);
 396 long __fixdfsi (double);
 397 long __fixsfsi (float);
 398
 399 int
 400 __unordxf2(long double a, long double b)
 401 {
 402   union long_double_long ldl;
 403
 404   ldl.ld = a;
 405   if (EXPX(ldl) == EXPXMASK
 406       && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0))
 407     return 1;
 408   ldl.ld = b;
 409   if (EXPX(ldl) == EXPXMASK
 410       && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0))
 411     return 1;
 412   return 0;
 413 }
 414
 415 /* convert double to long double */
 416 long double
 417 __extenddfxf2 (double d)
 418 {
 419   register union double_long dl;
 420   register union long_double_long ldl;
 421   register long exp;
 422
 423   dl.d = d;
 424   /*printf ("dfxf in: %g\n", d);*/
 425
 426   ldl.l.upper = SIGND (dl);
 427   if ((dl.l.upper & ~SIGNBIT) == 0 && !dl.l.lower)
 428     {
 429       ldl.l.middle = 0;
 430       ldl.l.lower = 0;
 431       return ldl.ld;
 432     }
 433
 434   exp = EXPD (dl) - EXCESSD + EXCESSX;
 435   ldl.l.upper |= exp << 16;
 436   ldl.l.middle = HIDDENX;
 437   /* 31-20: # mantissa bits in ldl.l.middle - # mantissa bits in dl.l.upper */
 438   ldl.l.middle |= (dl.l.upper & MANTDMASK) << (31 - 20);
 439   /* 1+20: explicit-integer-bit + # mantissa bits in dl.l.upper */
 440   ldl.l.middle |= dl.l.lower >> (1 + 20);
 441   /* 32 - 21: # bits of dl.l.lower in ldl.l.middle */
 442   ldl.l.lower = dl.l.lower << (32 - 21);
 443
 444   /*printf ("dfxf out: %s\n", dumpxf (ldl.ld));*/
 445   return ldl.ld;
 446 }
 447
 448 /* convert long double to double */
 449 double
 450 __truncxfdf2 (long double ld)
 451 {
 452   register long exp;
 453   register union double_long dl;
 454   register union long_double_long ldl;
 455
 456   ldl.ld = ld;
 457   /*printf ("xfdf in: %s\n", dumpxf (ld));*/
 458
 459   dl.l.upper = SIGNX (ldl);
 460   if ((ldl.l.upper & ~SIGNBIT) == 0 && !ldl.l.middle && !ldl.l.lower)
 461     {
 462       dl.l.lower = 0;
 463       return dl.d;
 464     }
 465
 466   exp = EXPX (ldl) - EXCESSX + EXCESSD;
 467   /* ??? quick and dirty: keep `exp' sane */
 468   if (exp >= EXPDMASK)
 469     exp = EXPDMASK - 1;
 470   dl.l.upper |= exp << (32 - (EXPDBITS + 1));
 471   /* +1-1: add one for sign bit, but take one off for explicit-integer-bit */
 472   dl.l.upper |= (ldl.l.middle & MANTXMASK) >> (EXPDBITS + 1 - 1);
 473   dl.l.lower = (ldl.l.middle & MANTXMASK) << (32 - (EXPDBITS + 1 - 1));
 474   dl.l.lower |= ldl.l.lower >> (EXPDBITS + 1 - 1);
 475
 476   /*printf ("xfdf out: %g\n", dl.d);*/
 477   return dl.d;
 478 }
 479
 480 /* convert a float to a long double */
 481 long double
 482 __extendsfxf2 (float f)
 483 {
 484   long double foo = __extenddfxf2 (__extendsfdf2 (f));
 485   return foo;
 486 }
 487
 488 /* convert a long double to a float */
 489 float
 490 __truncxfsf2 (long double ld)
 491 {
 492   float foo = __truncdfsf2 (__truncxfdf2 (ld));
 493   return foo;
 494 }
 495
 496 /* convert an int to a long double */
 497 long double
 498 __floatsixf (long l)
 499 {
 500   double foo = __floatsidf (l);
 501   return foo;
 502 }
 503
 504 /* convert an unsigned int to a long double */
 505 long double
 506 __floatunsixf (unsigned long l)
 507 {
 508   double foo = __floatunsidf (l);
 509   return foo;
 510 }
 511
 512 /* convert a long double to an int */
 513 long
 514 __fixxfsi (long double ld)
 515 {
 516   long foo = __fixdfsi ((double) ld);
 517   return foo;
 518 }
 519
 520 /* The remaining provide crude math support by working in double precision.  */
 521
 522 long double
 523 __addxf3 (long double x1, long double x2)
 524 {
 525   return (double) x1 + (double) x2;
 526 }
 527
 528 long double
 529 __subxf3 (long double x1, long double x2)
 530 {
 531   return (double) x1 - (double) x2;
 532 }
 533
 534 long double
 535 __mulxf3 (long double x1, long double x2)
 536 {
 537   return (double) x1 * (double) x2;
 538 }
 539
 540 long double
 541 __divxf3 (long double x1, long double x2)
 542 {
 543   return (double) x1 / (double) x2;
 544 }
 545
 546 long double
 547 __negxf2 (long double x1)
 548 {
 549   return - (double) x1;
 550 }
 551
 552 long
 553 __cmpxf2 (long double x1, long double x2)
 554 {
 555   return __cmpdf2 ((double) x1, (double) x2);
 556 }
 557
 558 long
 559 __eqxf2 (long double x1, long double x2)
 560 {
 561   return __cmpdf2 ((double) x1, (double) x2);
 562 }
 563
 564 long
 565 __nexf2 (long double x1, long double x2)
 566 {
 567   return __cmpdf2 ((double) x1, (double) x2);
 568 }
 569
 570 long
 571 __ltxf2 (long double x1, long double x2)
 572 {
 573   return __cmpdf2 ((double) x1, (double) x2);
 574 }
 575
 576 long
 577 __lexf2 (long double x1, long double x2)
 578 {
 579   return __cmpdf2 ((double) x1, (double) x2);
 580 }
 581
 582 long
 583 __gtxf2 (long double x1, long double x2)
 584 {
 585   return __cmpdf2 ((double) x1, (double) x2);
 586 }
 587
 588 long
 589 __gexf2 (long double x1, long double x2)
 590 {
 591   return __cmpdf2 ((double) x1, (double) x2);
 592 }
 593
 594 #endif /* !__mcoldfire__ */
 595 #endif /* EXTFLOAT */