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 81union double_long 82{ 83 double d; 84 struct { 85 long upper; 86 unsigned long lower; 87 } l; 88}; 89 90union float_long { 91 float f; 92 long l; 93}; 94 95union 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 108int 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 122int 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 */ 139double 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 */ 172double 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 */ 218float 219__floatunsisf (unsigned long l) 220{ 221 double foo = __floatunsidf (l); 222 return foo; 223} 224 225/* convert int to float */ 226float 227__floatsisf (long l) 228{ 229 double foo = __floatsidf (l); 230 return foo; 231} 232 233/* convert float to double */ 234double 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 */ 273float 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 */ 341long 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 */ 373long 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. */ 391double __floatunsidf (unsigned long); 392double __floatsidf (long); 393float __floatsisf (long); 394double __extendsfdf2 (float); 395float __truncdfsf2 (double); 396long __fixdfsi (double); 397long __fixsfsi (float); 398 399int 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 */ 416long 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 */ 449double 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 */ 481long 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 */ 489float 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 */ 497long 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 */ 505long 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 */ 513long 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 522long double 523__addxf3 (long double x1, long double x2) 524{ 525 return (double) x1 + (double) x2; 526} 527 528long double 529__subxf3 (long double x1, long double x2) 530{ 531 return (double) x1 - (double) x2; 532} 533 534long double 535__mulxf3 (long double x1, long double x2) 536{ 537 return (double) x1 * (double) x2; 538} 539 540long double 541__divxf3 (long double x1, long double x2) 542{ 543 return (double) x1 / (double) x2; 544} 545 546long double 547__negxf2 (long double x1) 548{ 549 return - (double) x1; 550} 551 552long 553__cmpxf2 (long double x1, long double x2) 554{ 555 return __cmpdf2 ((double) x1, (double) x2); 556} 557 558long 559__eqxf2 (long double x1, long double x2) 560{ 561 return __cmpdf2 ((double) x1, (double) x2); 562} 563 564long 565__nexf2 (long double x1, long double x2) 566{ 567 return __cmpdf2 ((double) x1, (double) x2); 568} 569 570long 571__ltxf2 (long double x1, long double x2) 572{ 573 return __cmpdf2 ((double) x1, (double) x2); 574} 575 576long 577__lexf2 (long double x1, long double x2) 578{ 579 return __cmpdf2 ((double) x1, (double) x2); 580} 581 582long 583__gtxf2 (long double x1, long double x2) 584{ 585 return __cmpdf2 ((double) x1, (double) x2); 586} 587 588long 589__gexf2 (long double x1, long double x2) 590{ 591 return __cmpdf2 ((double) x1, (double) x2); 592} 593 594#endif /* !__mcoldfire__ */ 595#endif /* EXTFLOAT */ 596