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.asm. 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 281 dl1.d = a1; 282 283 if ((dl1.l.upper & ~SIGNBIT) == 0 && !dl1.l.lower) 284 { 285 fl.l = SIGND(dl1); 286 return fl.f; 287 } 288 289 exp = EXPD (dl1) - EXCESSD + EXCESS; 290 291 /* shift double mantissa 6 bits so we can round */ 292 mant = MANTD (dl1) >> 6; 293 294 /* Check for underflow and denormals. */ 295 if (exp <= 0) 296 { 297 if (exp < -24) 298 mant = 0; 299 else 300 mant >>= 1 - exp; 301 exp = 0; 302 } 303 304 /* now round and shift down */ 305 mant += 1; 306 mant >>= 1; 307 308 /* did the round overflow? */ 309 if (mant & 0xFF000000L) 310 { 311 mant >>= 1; 312 exp++; 313 } 314 315 mant &= ~HIDDEN; 316 317 /* pack up and go home */ 318 fl.l = PACK (SIGND (dl1), exp, mant); 319 return (fl.f); 320 } 321 322 /* convert double to int */ 323 long 324 __fixdfsi (double a1) 325 { 326 register union double_long dl1; 327 register long exp; 328 register long l; 329 330 dl1.d = a1; 331 332 if (!dl1.l.upper && !dl1.l.lower) 333 return 0; 334 335 exp = EXPD (dl1) - EXCESSD - 31; 336 l = MANTD (dl1); 337 338 if (exp > 0) 339 { 340 /* Return largest integer. */ 341 return SIGND (dl1) ? 0x80000000L : 0x7fffffffL; 342 } 343 344 if (exp <= -32) 345 return 0; 346 347 /* shift down until exp = 0 */ 348 if (exp < 0) 349 l >>= -exp; 350 351 return (SIGND (dl1) ? -l : l); 352 } 353 354 /* convert float to int */ 355 long 356 __fixsfsi (float a1) 357 { 358 double foo = a1; 359 return __fixdfsi (foo); 360 } 361 362 #else /* EXTFLOAT */ 363 364 /* Primitive extended precision floating point support. 365 366 We assume all numbers are normalized, don't do any rounding, etc. */ 367 368 /* Prototypes for the above in case we use them. */ 369 double __floatunsidf (unsigned long); 370 double __floatsidf (long); 371 float __floatsisf (long); 372 double __extendsfdf2 (float); 373 float __truncdfsf2 (double); 374 long __fixdfsi (double); 375 long __fixsfsi (float); 376 377 int 378 __unordxf2(long double a, long double b) 379 { 380 union long_double_long ldl; 381 382 ldl.ld = a; 383 if (EXPX(ldl) == EXPXMASK 384 && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0)) 385 return 1; 386 ldl.ld = b; 387 if (EXPX(ldl) == EXPXMASK 388 && ((ldl.l.middle & MANTXMASK) != 0 || ldl.l.lower != 0)) 389 return 1; 390 return 0; 391 } 392 393 /* convert double to long double */ 394 long double 395 __extenddfxf2 (double d) 396 { 397 register union double_long dl; 398 register union long_double_long ldl; 399 register long exp; 400 401 dl.d = d; 402 /*printf ("dfxf in: %g\n", d);*/ 403 404 ldl.l.upper = SIGND (dl); 405 if ((dl.l.upper & ~SIGNBIT) == 0 && !dl.l.lower) 406 { 407 ldl.l.middle = 0; 408 ldl.l.lower = 0; 409 return ldl.ld; 410 } 411 412 exp = EXPD (dl) - EXCESSD + EXCESSX; 413 ldl.l.upper |= exp << 16; 414 ldl.l.middle = HIDDENX; 415 /* 31-20: # mantissa bits in ldl.l.middle - # mantissa bits in dl.l.upper */ 416 ldl.l.middle |= (dl.l.upper & MANTDMASK) << (31 - 20); 417 /* 1+20: explicit-integer-bit + # mantissa bits in dl.l.upper */ 418 ldl.l.middle |= dl.l.lower >> (1 + 20); 419 /* 32 - 21: # bits of dl.l.lower in ldl.l.middle */ 420 ldl.l.lower = dl.l.lower << (32 - 21); 421 422 /*printf ("dfxf out: %s\n", dumpxf (ldl.ld));*/ 423 return ldl.ld; 424 } 425 426 /* convert long double to double */ 427 double 428 __truncxfdf2 (long double ld) 429 { 430 register long exp; 431 register union double_long dl; 432 register union long_double_long ldl; 433 434 ldl.ld = ld; 435 /*printf ("xfdf in: %s\n", dumpxf (ld));*/ 436 437 dl.l.upper = SIGNX (ldl); 438 if ((ldl.l.upper & ~SIGNBIT) == 0 && !ldl.l.middle && !ldl.l.lower) 439 { 440 dl.l.lower = 0; 441 return dl.d; 442 } 443 444 exp = EXPX (ldl) - EXCESSX + EXCESSD; 445 /* ??? quick and dirty: keep `exp' sane */ 446 if (exp >= EXPDMASK) 447 exp = EXPDMASK - 1; 448 dl.l.upper |= exp << (32 - (EXPDBITS + 1)); 449 /* +1-1: add one for sign bit, but take one off for explicit-integer-bit */ 450 dl.l.upper |= (ldl.l.middle & MANTXMASK) >> (EXPDBITS + 1 - 1); 451 dl.l.lower = (ldl.l.middle & MANTXMASK) << (32 - (EXPDBITS + 1 - 1)); 452 dl.l.lower |= ldl.l.lower >> (EXPDBITS + 1 - 1); 453 454 /*printf ("xfdf out: %g\n", dl.d);*/ 455 return dl.d; 456 } 457 458 /* convert a float to a long double */ 459 long double 460 __extendsfxf2 (float f) 461 { 462 long double foo = __extenddfxf2 (__extendsfdf2 (f)); 463 return foo; 464 } 465 466 /* convert a long double to a float */ 467 float 468 __truncxfsf2 (long double ld) 469 { 470 float foo = __truncdfsf2 (__truncxfdf2 (ld)); 471 return foo; 472 } 473 474 /* convert an int to a long double */ 475 long double 476 __floatsixf (long l) 477 { 478 double foo = __floatsidf (l); 479 return foo; 480 } 481 482 /* convert an unsigned int to a long double */ 483 long double 484 __floatunsixf (unsigned long l) 485 { 486 double foo = __floatunsidf (l); 487 return foo; 488 } 489 490 /* convert a long double to an int */ 491 long 492 __fixxfsi (long double ld) 493 { 494 long foo = __fixdfsi ((double) ld); 495 return foo; 496 } 497 498 /* The remaining provide crude math support by working in double precision. */ 499 500 long double 501 __addxf3 (long double x1, long double x2) 502 { 503 return (double) x1 + (double) x2; 504 } 505 506 long double 507 __subxf3 (long double x1, long double x2) 508 { 509 return (double) x1 - (double) x2; 510 } 511 512 long double 513 __mulxf3 (long double x1, long double x2) 514 { 515 return (double) x1 * (double) x2; 516 } 517 518 long double 519 __divxf3 (long double x1, long double x2) 520 { 521 return (double) x1 / (double) x2; 522 } 523 524 long double 525 __negxf2 (long double x1) 526 { 527 return - (double) x1; 528 } 529 530 long 531 __cmpxf2 (long double x1, long double x2) 532 { 533 return __cmpdf2 ((double) x1, (double) x2); 534 } 535 536 long 537 __eqxf2 (long double x1, long double x2) 538 { 539 return __cmpdf2 ((double) x1, (double) x2); 540 } 541 542 long 543 __nexf2 (long double x1, long double x2) 544 { 545 return __cmpdf2 ((double) x1, (double) x2); 546 } 547 548 long 549 __ltxf2 (long double x1, long double x2) 550 { 551 return __cmpdf2 ((double) x1, (double) x2); 552 } 553 554 long 555 __lexf2 (long double x1, long double x2) 556 { 557 return __cmpdf2 ((double) x1, (double) x2); 558 } 559 560 long 561 __gtxf2 (long double x1, long double x2) 562 { 563 return __cmpdf2 ((double) x1, (double) x2); 564 } 565 566 long 567 __gexf2 (long double x1, long double x2) 568 { 569 return __cmpdf2 ((double) x1, (double) x2); 570 } 571 572 #endif /* EXTFLOAT */ 573