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 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 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 long __cmpdf2 (double, double); 399 400 int 401 __unordxf2(long double a, long double b) 402 { 403 union long_double_long ldl; 404 405 ldl.ld = a; 406 if (EXPX(ldl) == EXPXMASK 407 && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0)) 408 return 1; 409 ldl.ld = b; 410 if (EXPX(ldl) == EXPXMASK 411 && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0)) 412 return 1; 413 return 0; 414 } 415 416 /* convert double to long double */ 417 long double 418 __extenddfxf2 (double d) 419 { 420 register union double_long dl; 421 register union long_double_long ldl; 422 register long exp; 423 424 dl.d = d; 425 /*printf ("dfxf in: %g\n", d);*/ 426 427 ldl.l.upper = SIGND (dl); 428 if ((dl.l.upper & ~SIGNBIT) == 0 && !dl.l.lower) 429 { 430 ldl.l.middle = 0; 431 ldl.l.lower = 0; 432 return ldl.ld; 433 } 434 435 exp = EXPD (dl) - EXCESSD + EXCESSX; 436 ldl.l.upper |= exp << 16; 437 ldl.l.middle = HIDDENX; 438 /* 31-20: # mantissa bits in ldl.l.middle - # mantissa bits in dl.l.upper */ 439 ldl.l.middle |= (dl.l.upper & MANTDMASK) << (31 - 20); 440 /* 1+20: explicit-integer-bit + # mantissa bits in dl.l.upper */ 441 ldl.l.middle |= dl.l.lower >> (1 + 20); 442 /* 32 - 21: # bits of dl.l.lower in ldl.l.middle */ 443 ldl.l.lower = dl.l.lower << (32 - 21); 444 445 /*printf ("dfxf out: %s\n", dumpxf (ldl.ld));*/ 446 return ldl.ld; 447 } 448 449 /* convert long double to double */ 450 double 451 __truncxfdf2 (long double ld) 452 { 453 register long exp; 454 register union double_long dl; 455 register union long_double_long ldl; 456 457 ldl.ld = ld; 458 /*printf ("xfdf in: %s\n", dumpxf (ld));*/ 459 460 dl.l.upper = SIGNX (ldl); 461 if ((ldl.l.upper & ~SIGNBIT) == 0 && !ldl.l.middle && !ldl.l.lower) 462 { 463 dl.l.lower = 0; 464 return dl.d; 465 } 466 467 exp = EXPX (ldl) - EXCESSX + EXCESSD; 468 /* ??? quick and dirty: keep `exp' sane */ 469 if (exp >= EXPDMASK) 470 exp = EXPDMASK - 1; 471 dl.l.upper |= exp << (32 - (EXPDBITS + 1)); 472 /* +1-1: add one for sign bit, but take one off for explicit-integer-bit */ 473 dl.l.upper |= (ldl.l.middle & MANTXMASK) >> (EXPDBITS + 1 - 1); 474 dl.l.lower = (ldl.l.middle & MANTXMASK) << (32 - (EXPDBITS + 1 - 1)); 475 dl.l.lower |= ldl.l.lower >> (EXPDBITS + 1 - 1); 476 477 /*printf ("xfdf out: %g\n", dl.d);*/ 478 return dl.d; 479 } 480 481 /* convert a float to a long double */ 482 long double 483 __extendsfxf2 (float f) 484 { 485 long double foo = __extenddfxf2 (__extendsfdf2 (f)); 486 return foo; 487 } 488 489 /* convert a long double to a float */ 490 float 491 __truncxfsf2 (long double ld) 492 { 493 float foo = __truncdfsf2 (__truncxfdf2 (ld)); 494 return foo; 495 } 496 497 /* convert an int to a long double */ 498 long double 499 __floatsixf (long l) 500 { 501 double foo = __floatsidf (l); 502 return foo; 503 } 504 505 /* convert an unsigned int to a long double */ 506 long double 507 __floatunsixf (unsigned long l) 508 { 509 double foo = __floatunsidf (l); 510 return foo; 511 } 512 513 /* convert a long double to an int */ 514 long 515 __fixxfsi (long double ld) 516 { 517 long foo = __fixdfsi ((double) ld); 518 return foo; 519 } 520 521 /* The remaining provide crude math support by working in double precision. */ 522 523 long double 524 __addxf3 (long double x1, long double x2) 525 { 526 return (double) x1 + (double) x2; 527 } 528 529 long double 530 __subxf3 (long double x1, long double x2) 531 { 532 return (double) x1 - (double) x2; 533 } 534 535 long double 536 __mulxf3 (long double x1, long double x2) 537 { 538 return (double) x1 * (double) x2; 539 } 540 541 long double 542 __divxf3 (long double x1, long double x2) 543 { 544 return (double) x1 / (double) x2; 545 } 546 547 long double 548 __negxf2 (long double x1) 549 { 550 return - (double) x1; 551 } 552 553 long 554 __cmpxf2 (long double x1, long double x2) 555 { 556 return __cmpdf2 ((double) x1, (double) x2); 557 } 558 559 long 560 __eqxf2 (long double x1, long double x2) 561 { 562 return __cmpdf2 ((double) x1, (double) x2); 563 } 564 565 long 566 __nexf2 (long double x1, long double x2) 567 { 568 return __cmpdf2 ((double) x1, (double) x2); 569 } 570 571 long 572 __ltxf2 (long double x1, long double x2) 573 { 574 return __cmpdf2 ((double) x1, (double) x2); 575 } 576 577 long 578 __lexf2 (long double x1, long double x2) 579 { 580 return __cmpdf2 ((double) x1, (double) x2); 581 } 582 583 long 584 __gtxf2 (long double x1, long double x2) 585 { 586 return __cmpdf2 ((double) x1, (double) x2); 587 } 588 589 long 590 __gexf2 (long double x1, long double x2) 591 { 592 return __cmpdf2 ((double) x1, (double) x2); 593 } 594 595 #endif /* !__mcoldfire__ */ 596 #endif /* EXTFLOAT */ 597