contrib/gcc/floatlib.c

   1 /*
   2 ** libgcc support for software floating point.
   3 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
   4 ** Permission is granted to do *anything* you want with this file,
   5 ** commercial or otherwise, provided this message remains intact.  So there!
   6 ** I would appreciate receiving any updates/patches/changes that anyone
   7 ** makes, and am willing to be the repository for said changes (am I
   8 ** making a big mistake?).
   9
  10 Warning! Only single-precision is actually implemented.  This file
  11 won't really be much use until double-precision is supported.
  12
  13 However, once that is done, this file might eventually become a
  14 replacement for libgcc1.c.  It might also make possible
  15 cross-compilation for an IEEE target machine from a non-IEEE
  16 host such as a VAX.
  17
  18 If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu.
  19
  20
  21 **
  22 ** Pat Wood
  23 ** Pipeline Associates, Inc.
  24 ** pipeline!phw@motown.com or
  25 ** sun!pipeline!phw or
  26 ** uunet!motown!pipeline!phw
  27 **
  28 ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
  29 ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
  30 **                  -- fixed problems with adding and subtracting zero
  31 **                  -- fixed rounding in truncdfsf2
  32 **                  -- fixed SWAP define and tested on 386
  33 */
  34
  35 /*
  36 ** The following are routines that replace the libgcc soft floating point
  37 ** routines that are called automatically when -msoft-float is selected.
  38 ** The support single and double precision IEEE format, with provisions
  39 ** for byte-swapped machines (tested on 386).  Some of the double-precision
  40 ** routines work at full precision, but most of the hard ones simply punt
  41 ** and call the single precision routines, producing a loss of accuracy.
  42 ** long long support is not assumed or included.
  43 ** Overall accuracy is close to IEEE (actually 68882) for single-precision
  44 ** arithmetic.  I think there may still be a 1 in 1000 chance of a bit
  45 ** being rounded the wrong way during a multiply.  I'm not fussy enough to
  46 ** bother with it, but if anyone is, knock yourself out.
  47 **
  48 ** Efficiency has only been addressed where it was obvious that something
  49 ** would make a big difference.  Anyone who wants to do this right for
  50 ** best speed should go in and rewrite in assembler.
  51 **
  52 ** I have tested this only on a 68030 workstation and 386/ix integrated
  53 ** in with -msoft-float.
  54 */
  55
  56 /* the following deal with IEEE single-precision numbers */
  57 #define EXCESS          126
  58 #define SIGNBIT         0x80000000
  59 #define HIDDEN          (1 << 23)
  60 #define SIGN(fp)        ((fp) & SIGNBIT)
  61 #define EXP(fp)         (((fp) >> 23) & 0xFF)
  62 #define MANT(fp)        (((fp) & 0x7FFFFF) | HIDDEN)
  63 #define PACK(s,e,m)     ((s) | ((e) << 23) | (m))
  64
  65 /* the following deal with IEEE double-precision numbers */
  66 #define EXCESSD         1022
  67 #define HIDDEND         (1 << 20)
  68 #define EXPD(fp)        (((fp.l.upper) >> 20) & 0x7FF)
  69 #define SIGND(fp)       ((fp.l.upper) & SIGNBIT)
  70 #define MANTD(fp)       (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
  71                                 (fp.l.lower >> 22))
  72
  73 /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
  74 union double_long
  75   {
  76     double d;
  77 #ifdef SWAP
  78     struct {
  79       unsigned long lower;
  80       long upper;
  81     } l;
  82 #else
  83     struct {
  84       long upper;
  85       unsigned long lower;
  86     } l;
  87 #endif
  88   };
  89
  90 union float_long
  91   {
  92     float f;
  93     long l;
  94   };
  95
  96 /* add two floats */
  97 float
  98 __addsf3 (float a1, float a2)
  99 {
 100   register long mant1, mant2;
 101   register union float_long fl1, fl2;
 102   register int exp1, exp2;
 103   int sign = 0;
 104
 105   fl1.f = a1;
 106   fl2.f = a2;
 107
 108   /* check for zero args */
 109   if (!fl1.l)
 110     return (fl2.f);
 111   if (!fl2.l)
 112     return (fl1.f);
 113
 114   exp1 = EXP (fl1.l);
 115   exp2 = EXP (fl2.l);
 116
 117   if (exp1 > exp2 + 25)
 118     return (fl1.l);
 119   if (exp2 > exp1 + 25)
 120     return (fl2.l);
 121
 122   /* do everything in excess precision so's we can round later */
 123   mant1 = MANT (fl1.l) << 6;
 124   mant2 = MANT (fl2.l) << 6;
 125
 126   if (SIGN (fl1.l))
 127     mant1 = -mant1;
 128   if (SIGN (fl2.l))
 129     mant2 = -mant2;
 130
 131   if (exp1 > exp2)
 132     {
 133       mant2 >>= exp1 - exp2;
 134     }
 135   else
 136     {
 137       mant1 >>= exp2 - exp1;
 138       exp1 = exp2;
 139     }
 140   mant1 += mant2;
 141
 142   if (mant1 < 0)
 143     {
 144       mant1 = -mant1;
 145       sign = SIGNBIT;
 146     }
 147   else if (!mant1)
 148     return (0);
 149
 150   /* normalize up */
 151   while (!(mant1 & 0xE0000000))
 152     {
 153       mant1 <<= 1;
 154       exp1--;
 155     }
 156
 157   /* normalize down? */
 158   if (mant1 & (1 << 30))
 159     {
 160       mant1 >>= 1;
 161       exp1++;
 162     }
 163
 164   /* round to even */
 165   mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
 166
 167   /* normalize down? */
 168   if (mant1 & (1 << 30))
 169     {
 170       mant1 >>= 1;
 171       exp1++;
 172     }
 173
 174   /* lose extra precision */
 175   mant1 >>= 6;
 176
 177   /* turn off hidden bit */
 178   mant1 &= ~HIDDEN;
 179
 180   /* pack up and go home */
 181   fl1.l = PACK (sign, exp1, mant1);
 182   return (fl1.f);
 183 }
 184
 185 /* subtract two floats */
 186 float
 187 __subsf3 (float a1, float a2)
 188 {
 189   register union float_long fl1, fl2;
 190
 191   fl1.f = a1;
 192   fl2.f = a2;
 193
 194   /* check for zero args */
 195   if (!fl2.l)
 196     return (fl1.f);
 197   if (!fl1.l)
 198     return (-fl2.f);
 199
 200   /* twiddle sign bit and add */
 201   fl2.l ^= SIGNBIT;
 202   return __addsf3 (a1, fl2.f);
 203 }
 204
 205 /* compare two floats */
 206 long
 207 __cmpsf2 (float a1, float a2)
 208 {
 209   register union float_long fl1, fl2;
 210
 211   fl1.f = a1;
 212   fl2.f = a2;
 213
 214   if (SIGN (fl1.l) && SIGN (fl2.l))
 215     {
 216       fl1.l ^= SIGNBIT;
 217       fl2.l ^= SIGNBIT;
 218     }
 219   if (fl1.l < fl2.l)
 220     return (-1);
 221   if (fl1.l > fl2.l)
 222     return (1);
 223   return (0);
 224 }
 225
 226 /* multiply two floats */
 227 float
 228 __mulsf3 (float a1, float a2)
 229 {
 230   register union float_long fl1, fl2;
 231   register unsigned long result;
 232   register int exp;
 233   int sign;
 234
 235   fl1.f = a1;
 236   fl2.f = a2;
 237
 238   if (!fl1.l || !fl2.l)
 239     return (0);
 240
 241   /* compute sign and exponent */
 242   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 243   exp = EXP (fl1.l) - EXCESS;
 244   exp += EXP (fl2.l);
 245
 246   fl1.l = MANT (fl1.l);
 247   fl2.l = MANT (fl2.l);
 248
 249   /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
 250   result = (fl1.l >> 8) * (fl2.l >> 8);
 251   result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
 252   result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
 253
 254   if (result & 0x80000000)
 255     {
 256       /* round */
 257       result += 0x80;
 258       result >>= 8;
 259     }
 260   else
 261     {
 262       /* round */
 263       result += 0x40;
 264       result >>= 7;
 265       exp--;
 266     }
 267
 268   result &= ~HIDDEN;
 269
 270   /* pack up and go home */
 271   fl1.l = PACK (sign, exp, result);
 272   return (fl1.f);
 273 }
 274
 275 /* divide two floats */
 276 float
 277 __divsf3 (float a1, float a2)
 278 {
 279   register union float_long fl1, fl2;
 280   register int result;
 281   register int mask;
 282   register int exp, sign;
 283
 284   fl1.f = a1;
 285   fl2.f = a2;
 286
 287   /* subtract exponents */
 288   exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
 289
 290   /* compute sign */
 291   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 292
 293   /* divide by zero??? */
 294   if (!fl2.l)
 295     /* return NaN or -NaN */
 296     return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
 297
 298   /* numerator zero??? */
 299   if (!fl1.l)
 300     return (0);
 301
 302   /* now get mantissas */
 303   fl1.l = MANT (fl1.l);
 304   fl2.l = MANT (fl2.l);
 305
 306   /* this assures we have 25 bits of precision in the end */
 307   if (fl1.l < fl2.l)
 308     {
 309       fl1.l <<= 1;
 310       exp--;
 311     }
 312
 313   /* now we perform repeated subtraction of fl2.l from fl1.l */
 314   mask = 0x1000000;
 315   result = 0;
 316   while (mask)
 317     {
 318       if (fl1.l >= fl2.l)
 319         {
 320           result |= mask;
 321           fl1.l -= fl2.l;
 322         }
 323       fl1.l <<= 1;
 324       mask >>= 1;
 325     }
 326
 327   /* round */
 328   result += 1;
 329
 330   /* normalize down */
 331   exp++;
 332   result >>= 1;
 333
 334   result &= ~HIDDEN;
 335
 336   /* pack up and go home */
 337   fl1.l = PACK (sign, exp, result);
 338   return (fl1.f);
 339 }
 340
 341 /* convert int to double */
 342 double
 343 __floatsidf (register long a1)
 344 {
 345   register int sign = 0, exp = 31 + EXCESSD;
 346   union double_long dl;
 347
 348   if (!a1)
 349     {
 350       dl.l.upper = dl.l.lower = 0;
 351       return (dl.d);
 352     }
 353
 354   if (a1 < 0)
 355     {
 356       sign = SIGNBIT;
 357       a1 = -a1;
 358     }
 359
 360   while (a1 < 0x1000000)
 361     {
 362       a1 <<= 4;
 363       exp -= 4;
 364     }
 365
 366   while (a1 < 0x40000000)
 367     {
 368       a1 <<= 1;
 369       exp--;
 370     }
 371
 372   /* pack up and go home */
 373   dl.l.upper = sign;
 374   dl.l.upper |= exp << 20;
 375   dl.l.upper |= (a1 >> 10) & ~HIDDEND;
 376   dl.l.lower = a1 << 22;
 377
 378   return (dl.d);
 379 }
 380
 381 /* negate a float */
 382 float
 383 __negsf2 (float a1)
 384 {
 385   register union float_long fl1;
 386
 387   fl1.f = a1;
 388   if (!fl1.l)
 389     return (0);
 390
 391   fl1.l ^= SIGNBIT;
 392   return (fl1.f);
 393 }
 394
 395 /* negate a double */
 396 double
 397 __negdf2 (double a1)
 398 {
 399   register union double_long dl1;
 400
 401   dl1.d = a1;
 402
 403   if (!dl1.l.upper && !dl1.l.lower)
 404       return (dl1.d);
 405
 406   dl1.l.upper ^= SIGNBIT;
 407   return (dl1.d);
 408 }
 409
 410 /* convert float to double */
 411 double
 412 __extendsfdf2 (float a1)
 413 {
 414   register union float_long fl1;
 415   register union double_long dl;
 416   register int exp;
 417
 418   fl1.f = a1;
 419
 420   if (!fl1.l)
 421     {
 422       dl.l.upper = dl.l.lower = 0;
 423       return (dl.d);
 424     }
 425
 426   dl.l.upper = SIGN (fl1.l);
 427   exp = EXP (fl1.l) - EXCESS + EXCESSD;
 428   dl.l.upper |= exp << 20;
 429   dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
 430   dl.l.lower = MANT (fl1.l) << 29;
 431
 432   return (dl.d);
 433 }
 434
 435 /* convert double to float */
 436 float
 437 __truncdfsf2 (double a1)
 438 {
 439   register int exp;
 440   register long mant;
 441   register union float_long fl;
 442   register union double_long dl1;
 443
 444   dl1.d = a1;
 445
 446   if (!dl1.l.upper && !dl1.l.lower)
 447     return (0);
 448
 449   exp = EXPD (dl1) - EXCESSD + EXCESS;
 450
 451   /* shift double mantissa 6 bits so we can round */
 452   mant = MANTD (dl1) >> 6;
 453
 454   /* now round and shift down */
 455   mant += 1;
 456   mant >>= 1;
 457
 458   /* did the round overflow? */
 459   if (mant & 0xFF000000)
 460     {
 461       mant >>= 1;
 462       exp++;
 463     }
 464
 465   mant &= ~HIDDEN;
 466
 467   /* pack up and go home */
 468   fl.l = PACK (SIGND (dl1), exp, mant);
 469   return (fl.f);
 470 }
 471
 472 /* compare two doubles */
 473 long
 474 __cmpdf2 (double a1, double a2)
 475 {
 476   register union double_long dl1, dl2;
 477
 478   dl1.d = a1;
 479   dl2.d = a2;
 480
 481   if (SIGND (dl1) && SIGND (dl2))
 482     {
 483       dl1.l.upper ^= SIGNBIT;
 484       dl2.l.upper ^= SIGNBIT;
 485     }
 486   if (dl1.l.upper < dl2.l.upper)
 487     return (-1);
 488   if (dl1.l.upper > dl2.l.upper)
 489     return (1);
 490   if (dl1.l.lower < dl2.l.lower)
 491     return (-1);
 492   if (dl1.l.lower > dl2.l.lower)
 493     return (1);
 494   return (0);
 495 }
 496
 497 /* convert double to int */
 498 long
 499 __fixdfsi (double a1)
 500 {
 501   register union double_long dl1;
 502   register int exp;
 503   register long l;
 504
 505   dl1.d = a1;
 506
 507   if (!dl1.l.upper && !dl1.l.lower)
 508     return (0);
 509
 510   exp = EXPD (dl1) - EXCESSD - 31;
 511   l = MANTD (dl1);
 512
 513   if (exp > 0)
 514     return (0x7FFFFFFF | SIGND (dl1)); /* largest integer */
 515
 516   /* shift down until exp = 0 or l = 0 */
 517   if (exp < 0 && exp > -32 && l)
 518     l >>= -exp;
 519   else
 520     return (0);
 521
 522   return (SIGND (dl1) ? -l : l);
 523 }
 524
 525 /* convert double to unsigned int */
 526 unsigned
 527 long __fixunsdfsi (double a1)
 528 {
 529   register union double_long dl1;
 530   register int exp;
 531   register unsigned long l;
 532
 533   dl1.d = a1;
 534
 535   if (!dl1.l.upper && !dl1.l.lower)
 536     return (0);
 537
 538   exp = EXPD (dl1) - EXCESSD - 32;
 539   l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
 540
 541   if (exp > 0)
 542     return (0xFFFFFFFF);        /* largest integer */
 543
 544   /* shift down until exp = 0 or l = 0 */
 545   if (exp < 0 && exp > -32 && l)
 546     l >>= -exp;
 547   else
 548     return (0);
 549
 550   return (l);
 551 }
 552
 553 /* For now, the hard double-precision routines simply
 554    punt and do it in single */
 555 /* addtwo doubles */
 556 double
 557 __adddf3 (double a1, double a2)
 558 {
 559   return ((float) a1 + (float) a2);
 560 }
 561
 562 /* subtract two doubles */
 563 double
 564 __subdf3 (double a1, double a2)
 565 {
 566   return ((float) a1 - (float) a2);
 567 }
 568
 569 /* multiply two doubles */
 570 double
 571 __muldf3 (double a1, double a2)
 572 {
 573   return ((float) a1 * (float) a2);
 574 }
 575
 576 /* divide two doubles */
 577 double
 578 __divdf3 (double a1, double a2)
 579 {
 580   return ((float) a1 / (float) a2);
 581 }