Annotation of sys/lib/libkern/softfloat.c, Revision 1.1
1.1 ! nbrk 1: /* $OpenBSD: softfloat.c,v 1.1 2002/04/28 20:55:14 pvalchev Exp $ */
! 2: /* $NetBSD: softfloat.c,v 1.1 2001/04/26 03:10:47 ross Exp $ */
! 3:
! 4: /*
! 5: * This version hacked for use with gcc -msoft-float by bjh21.
! 6: * (Mostly a case of #ifdefing out things GCC doesn't need or provides
! 7: * itself).
! 8: */
! 9:
! 10: /*
! 11: * Things you may want to define:
! 12: *
! 13: * SOFTFLOAT_FOR_GCC - build only those functions necessary for GCC (with
! 14: * -msoft-float) to work. Include "softfloat-for-gcc.h" to get them
! 15: * properly renamed.
! 16: */
! 17:
! 18: /*
! 19: ===============================================================================
! 20:
! 21: This C source file is part of the SoftFloat IEC/IEEE Floating-point
! 22: Arithmetic Package, Release 2a.
! 23:
! 24: Written by John R. Hauser. This work was made possible in part by the
! 25: International Computer Science Institute, located at Suite 600, 1947 Center
! 26: Street, Berkeley, California 94704. Funding was partially provided by the
! 27: National Science Foundation under grant MIP-9311980. The original version
! 28: of this code was written as part of a project to build a fixed-point vector
! 29: processor in collaboration with the University of California at Berkeley,
! 30: overseen by Profs. Nelson Morgan and John Wawrzynek. More information
! 31: is available through the Web page `http://HTTP.CS.Berkeley.EDU/~jhauser/
! 32: arithmetic/SoftFloat.html'.
! 33:
! 34: THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable
! 35: effort has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT
! 36: WILL AT TIMES RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS
! 37: RESTRICTED TO PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL
! 38: RESPONSIBILITY FOR ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM
! 39: THEIR OWN USE OF THE SOFTWARE, AND WHO ALSO EFFECTIVELY INDEMNIFY
! 40: (possibly via similar legal warning) JOHN HAUSER AND THE INTERNATIONAL
! 41: COMPUTER SCIENCE INSTITUTE AGAINST ALL LOSSES, COSTS, OR OTHER PROBLEMS
! 42: ARISING FROM THE USE OF THE SOFTWARE BY THEIR CUSTOMERS AND CLIENTS.
! 43:
! 44: Derivative works are acceptable, even for commercial purposes, so long as
! 45: (1) they include prominent notice that the work is derivative, and (2) they
! 46: include prominent notice akin to these four paragraphs for those parts of
! 47: this code that are retained.
! 48:
! 49: ===============================================================================
! 50: */
! 51:
! 52: #ifndef NO_IEEE
! 53:
! 54: #include <sys/cdefs.h>
! 55: #if defined(LIBC_SCCS) && !defined(lint)
! 56: __RCSID("$NetBSD: softfloat.c,v 1.1 2001/04/26 03:10:47 ross Exp $");
! 57: #endif /* LIBC_SCCS and not lint */
! 58:
! 59: #ifdef SOFTFLOAT_FOR_GCC
! 60: #include "softfloat-for-gcc.h"
! 61: #endif
! 62:
! 63: #include "milieu.h"
! 64: #include "softfloat.h"
! 65:
! 66: /*
! 67: * Conversions between floats as stored in memory and floats as
! 68: * SoftFloat uses them
! 69: */
! 70: #ifndef FLOAT64_DEMANGLE
! 71: #define FLOAT64_DEMANGLE(a) (a)
! 72: #endif
! 73: #ifndef FLOAT64_MANGLE
! 74: #define FLOAT64_MANGLE(a) (a)
! 75: #endif
! 76:
! 77: /*
! 78: -------------------------------------------------------------------------------
! 79: Floating-point rounding mode, extended double-precision rounding precision,
! 80: and exception flags.
! 81: -------------------------------------------------------------------------------
! 82: */
! 83:
! 84: /*
! 85: * XXX: This may cause options-MULTIPROCESSOR or thread problems someday.
! 86: * Right now, it does not. I've removed all other dynamic global
! 87: * variables. [ross]
! 88: */
! 89: #ifdef FLOATX80
! 90: int8 floatx80_rounding_precision = 80;
! 91: #endif
! 92:
! 93: /*
! 94: -------------------------------------------------------------------------------
! 95: Primitive arithmetic functions, including multi-word arithmetic, and
! 96: division and square root approximations. (Can be specialized to target if
! 97: desired.)
! 98: -------------------------------------------------------------------------------
! 99: */
! 100: #include "softfloat-macros.h"
! 101:
! 102: /*
! 103: -------------------------------------------------------------------------------
! 104: Functions and definitions to determine: (1) whether tininess for underflow
! 105: is detected before or after rounding by default, (2) what (if anything)
! 106: happens when exceptions are raised, (3) how signaling NaNs are distinguished
! 107: from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
! 108: are propagated from function inputs to output. These details are target-
! 109: specific.
! 110: -------------------------------------------------------------------------------
! 111: */
! 112: #include "softfloat-specialize.h"
! 113:
! 114: #ifndef SOFTFLOAT_FOR_GCC /* Not used */
! 115: /*
! 116: -------------------------------------------------------------------------------
! 117: Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
! 118: and 7, and returns the properly rounded 32-bit integer corresponding to the
! 119: input. If `zSign' is 1, the input is negated before being converted to an
! 120: integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
! 121: is simply rounded to an integer, with the inexact exception raised if the
! 122: input cannot be represented exactly as an integer. However, if the fixed-
! 123: point input is too large, the invalid exception is raised and the largest
! 124: positive or negative integer is returned.
! 125: -------------------------------------------------------------------------------
! 126: */
! 127: static int32 roundAndPackInt32( flag zSign, bits64 absZ )
! 128: {
! 129: int8 roundingMode;
! 130: flag roundNearestEven;
! 131: int8 roundIncrement, roundBits;
! 132: int32 z;
! 133:
! 134: roundingMode = float_rounding_mode();
! 135: roundNearestEven = ( roundingMode == float_round_nearest_even );
! 136: roundIncrement = 0x40;
! 137: if ( ! roundNearestEven ) {
! 138: if ( roundingMode == float_round_to_zero ) {
! 139: roundIncrement = 0;
! 140: }
! 141: else {
! 142: roundIncrement = 0x7F;
! 143: if ( zSign ) {
! 144: if ( roundingMode == float_round_up ) roundIncrement = 0;
! 145: }
! 146: else {
! 147: if ( roundingMode == float_round_down ) roundIncrement = 0;
! 148: }
! 149: }
! 150: }
! 151: roundBits = absZ & 0x7F;
! 152: absZ = ( absZ + roundIncrement )>>7;
! 153: absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
! 154: z = absZ;
! 155: if ( zSign ) z = - z;
! 156: if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
! 157: float_raise( float_flag_invalid );
! 158: return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
! 159: }
! 160: if ( roundBits ) float_set_inexact();
! 161: return z;
! 162:
! 163: }
! 164:
! 165: /*
! 166: -------------------------------------------------------------------------------
! 167: Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
! 168: `absZ1', with binary point between bits 63 and 64 (between the input words),
! 169: and returns the properly rounded 64-bit integer corresponding to the input.
! 170: If `zSign' is 1, the input is negated before being converted to an integer.
! 171: Ordinarily, the fixed-point input is simply rounded to an integer, with
! 172: the inexact exception raised if the input cannot be represented exactly as
! 173: an integer. However, if the fixed-point input is too large, the invalid
! 174: exception is raised and the largest positive or negative integer is
! 175: returned.
! 176: -------------------------------------------------------------------------------
! 177: */
! 178: static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 )
! 179: {
! 180: int8 roundingMode;
! 181: flag roundNearestEven, increment;
! 182: int64 z;
! 183:
! 184: roundingMode = float_rounding_mode();
! 185: roundNearestEven = ( roundingMode == float_round_nearest_even );
! 186: increment = ( (sbits64) absZ1 < 0 );
! 187: if ( ! roundNearestEven ) {
! 188: if ( roundingMode == float_round_to_zero ) {
! 189: increment = 0;
! 190: }
! 191: else {
! 192: if ( zSign ) {
! 193: increment = ( roundingMode == float_round_down ) && absZ1;
! 194: }
! 195: else {
! 196: increment = ( roundingMode == float_round_up ) && absZ1;
! 197: }
! 198: }
! 199: }
! 200: if ( increment ) {
! 201: ++absZ0;
! 202: if ( absZ0 == 0 ) goto overflow;
! 203: absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );
! 204: }
! 205: z = absZ0;
! 206: if ( zSign ) z = - z;
! 207: if ( z && ( ( z < 0 ) ^ zSign ) ) {
! 208: overflow:
! 209: float_raise( float_flag_invalid );
! 210: return
! 211: zSign ? (sbits64) LIT64( 0x8000000000000000 )
! 212: : LIT64( 0x7FFFFFFFFFFFFFFF );
! 213: }
! 214: if ( absZ1 ) float_set_inexact();
! 215: return z;
! 216:
! 217: }
! 218: #endif
! 219:
! 220: /*
! 221: -------------------------------------------------------------------------------
! 222: Returns the fraction bits of the single-precision floating-point value `a'.
! 223: -------------------------------------------------------------------------------
! 224: */
! 225: INLINE bits32 extractFloat32Frac( float32 a )
! 226: {
! 227:
! 228: return a & 0x007FFFFF;
! 229:
! 230: }
! 231:
! 232: /*
! 233: -------------------------------------------------------------------------------
! 234: Returns the exponent bits of the single-precision floating-point value `a'.
! 235: -------------------------------------------------------------------------------
! 236: */
! 237: INLINE int16 extractFloat32Exp( float32 a )
! 238: {
! 239:
! 240: return ( a>>23 ) & 0xFF;
! 241:
! 242: }
! 243:
! 244: /*
! 245: -------------------------------------------------------------------------------
! 246: Returns the sign bit of the single-precision floating-point value `a'.
! 247: -------------------------------------------------------------------------------
! 248: */
! 249: INLINE flag extractFloat32Sign( float32 a )
! 250: {
! 251:
! 252: return a>>31;
! 253:
! 254: }
! 255:
! 256: /*
! 257: -------------------------------------------------------------------------------
! 258: Normalizes the subnormal single-precision floating-point value represented
! 259: by the denormalized significand `aSig'. The normalized exponent and
! 260: significand are stored at the locations pointed to by `zExpPtr' and
! 261: `zSigPtr', respectively.
! 262: -------------------------------------------------------------------------------
! 263: */
! 264: static void
! 265: normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
! 266: {
! 267: int8 shiftCount;
! 268:
! 269: shiftCount = countLeadingZeros32( aSig ) - 8;
! 270: *zSigPtr = aSig<<shiftCount;
! 271: *zExpPtr = 1 - shiftCount;
! 272:
! 273: }
! 274:
! 275: /*
! 276: -------------------------------------------------------------------------------
! 277: Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
! 278: single-precision floating-point value, returning the result. After being
! 279: shifted into the proper positions, the three fields are simply added
! 280: together to form the result. This means that any integer portion of `zSig'
! 281: will be added into the exponent. Since a properly normalized significand
! 282: will have an integer portion equal to 1, the `zExp' input should be 1 less
! 283: than the desired result exponent whenever `zSig' is a complete, normalized
! 284: significand.
! 285: -------------------------------------------------------------------------------
! 286: */
! 287: INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
! 288: {
! 289:
! 290: return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
! 291:
! 292: }
! 293:
! 294: /*
! 295: -------------------------------------------------------------------------------
! 296: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 297: and significand `zSig', and returns the proper single-precision floating-
! 298: point value corresponding to the abstract input. Ordinarily, the abstract
! 299: value is simply rounded and packed into the single-precision format, with
! 300: the inexact exception raised if the abstract input cannot be represented
! 301: exactly. However, if the abstract value is too large, the overflow and
! 302: inexact exceptions are raised and an infinity or maximal finite value is
! 303: returned. If the abstract value is too small, the input value is rounded to
! 304: a subnormal number, and the underflow and inexact exceptions are raised if
! 305: the abstract input cannot be represented exactly as a subnormal single-
! 306: precision floating-point number.
! 307: The input significand `zSig' has its binary point between bits 30
! 308: and 29, which is 7 bits to the left of the usual location. This shifted
! 309: significand must be normalized or smaller. If `zSig' is not normalized,
! 310: `zExp' must be 0; in that case, the result returned is a subnormal number,
! 311: and it must not require rounding. In the usual case that `zSig' is
! 312: normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
! 313: The handling of underflow and overflow follows the IEC/IEEE Standard for
! 314: Binary Floating-Point Arithmetic.
! 315: -------------------------------------------------------------------------------
! 316: */
! 317: static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig )
! 318: {
! 319: int8 roundingMode;
! 320: flag roundNearestEven;
! 321: int8 roundIncrement, roundBits;
! 322: flag isTiny;
! 323:
! 324: roundingMode = float_rounding_mode();
! 325: roundNearestEven = ( roundingMode == float_round_nearest_even );
! 326: roundIncrement = 0x40;
! 327: if ( ! roundNearestEven ) {
! 328: if ( roundingMode == float_round_to_zero ) {
! 329: roundIncrement = 0;
! 330: }
! 331: else {
! 332: roundIncrement = 0x7F;
! 333: if ( zSign ) {
! 334: if ( roundingMode == float_round_up ) roundIncrement = 0;
! 335: }
! 336: else {
! 337: if ( roundingMode == float_round_down ) roundIncrement = 0;
! 338: }
! 339: }
! 340: }
! 341: roundBits = zSig & 0x7F;
! 342: if ( 0xFD <= (bits16) zExp ) {
! 343: if ( ( 0xFD < zExp )
! 344: || ( ( zExp == 0xFD )
! 345: && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
! 346: ) {
! 347: float_raise( float_flag_overflow | float_flag_inexact );
! 348: return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
! 349: }
! 350: if ( zExp < 0 ) {
! 351: isTiny =
! 352: ( float_detect_tininess == float_tininess_before_rounding )
! 353: || ( zExp < -1 )
! 354: || ( zSig + roundIncrement < 0x80000000 );
! 355: shift32RightJamming( zSig, - zExp, &zSig );
! 356: zExp = 0;
! 357: roundBits = zSig & 0x7F;
! 358: if ( isTiny && roundBits ) float_raise( float_flag_underflow );
! 359: }
! 360: }
! 361: if ( roundBits ) float_set_inexact();
! 362: zSig = ( zSig + roundIncrement )>>7;
! 363: zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
! 364: if ( zSig == 0 ) zExp = 0;
! 365: return packFloat32( zSign, zExp, zSig );
! 366:
! 367: }
! 368:
! 369: /*
! 370: -------------------------------------------------------------------------------
! 371: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 372: and significand `zSig', and returns the proper single-precision floating-
! 373: point value corresponding to the abstract input. This routine is just like
! 374: `roundAndPackFloat32' except that `zSig' does not have to be normalized.
! 375: Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
! 376: floating-point exponent.
! 377: -------------------------------------------------------------------------------
! 378: */
! 379: static float32
! 380: normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig )
! 381: {
! 382: int8 shiftCount;
! 383:
! 384: shiftCount = countLeadingZeros32( zSig ) - 1;
! 385: return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount );
! 386:
! 387: }
! 388:
! 389: /*
! 390: -------------------------------------------------------------------------------
! 391: Returns the fraction bits of the double-precision floating-point value `a'.
! 392: -------------------------------------------------------------------------------
! 393: */
! 394: INLINE bits64 extractFloat64Frac( float64 a )
! 395: {
! 396:
! 397: return FLOAT64_DEMANGLE(a) & LIT64( 0x000FFFFFFFFFFFFF );
! 398:
! 399: }
! 400:
! 401: /*
! 402: -------------------------------------------------------------------------------
! 403: Returns the exponent bits of the double-precision floating-point value `a'.
! 404: -------------------------------------------------------------------------------
! 405: */
! 406: INLINE int16 extractFloat64Exp( float64 a )
! 407: {
! 408:
! 409: return ( FLOAT64_DEMANGLE(a)>>52 ) & 0x7FF;
! 410:
! 411: }
! 412:
! 413: /*
! 414: -------------------------------------------------------------------------------
! 415: Returns the sign bit of the double-precision floating-point value `a'.
! 416: -------------------------------------------------------------------------------
! 417: */
! 418: INLINE flag extractFloat64Sign( float64 a )
! 419: {
! 420:
! 421: return FLOAT64_DEMANGLE(a)>>63;
! 422:
! 423: }
! 424:
! 425: /*
! 426: -------------------------------------------------------------------------------
! 427: Normalizes the subnormal double-precision floating-point value represented
! 428: by the denormalized significand `aSig'. The normalized exponent and
! 429: significand are stored at the locations pointed to by `zExpPtr' and
! 430: `zSigPtr', respectively.
! 431: -------------------------------------------------------------------------------
! 432: */
! 433: static void
! 434: normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
! 435: {
! 436: int8 shiftCount;
! 437:
! 438: shiftCount = countLeadingZeros64( aSig ) - 11;
! 439: *zSigPtr = aSig<<shiftCount;
! 440: *zExpPtr = 1 - shiftCount;
! 441:
! 442: }
! 443:
! 444: /*
! 445: -------------------------------------------------------------------------------
! 446: Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
! 447: double-precision floating-point value, returning the result. After being
! 448: shifted into the proper positions, the three fields are simply added
! 449: together to form the result. This means that any integer portion of `zSig'
! 450: will be added into the exponent. Since a properly normalized significand
! 451: will have an integer portion equal to 1, the `zExp' input should be 1 less
! 452: than the desired result exponent whenever `zSig' is a complete, normalized
! 453: significand.
! 454: -------------------------------------------------------------------------------
! 455: */
! 456: INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
! 457: {
! 458:
! 459: return FLOAT64_MANGLE( ( ( (bits64) zSign )<<63 ) +
! 460: ( ( (bits64) zExp )<<52 ) + zSig );
! 461:
! 462: }
! 463:
! 464: /*
! 465: -------------------------------------------------------------------------------
! 466: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 467: and significand `zSig', and returns the proper double-precision floating-
! 468: point value corresponding to the abstract input. Ordinarily, the abstract
! 469: value is simply rounded and packed into the double-precision format, with
! 470: the inexact exception raised if the abstract input cannot be represented
! 471: exactly. However, if the abstract value is too large, the overflow and
! 472: inexact exceptions are raised and an infinity or maximal finite value is
! 473: returned. If the abstract value is too small, the input value is rounded to
! 474: a subnormal number, and the underflow and inexact exceptions are raised if
! 475: the abstract input cannot be represented exactly as a subnormal double-
! 476: precision floating-point number.
! 477: The input significand `zSig' has its binary point between bits 62
! 478: and 61, which is 10 bits to the left of the usual location. This shifted
! 479: significand must be normalized or smaller. If `zSig' is not normalized,
! 480: `zExp' must be 0; in that case, the result returned is a subnormal number,
! 481: and it must not require rounding. In the usual case that `zSig' is
! 482: normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
! 483: The handling of underflow and overflow follows the IEC/IEEE Standard for
! 484: Binary Floating-Point Arithmetic.
! 485: -------------------------------------------------------------------------------
! 486: */
! 487: static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig )
! 488: {
! 489: int8 roundingMode;
! 490: flag roundNearestEven;
! 491: int16 roundIncrement, roundBits;
! 492: flag isTiny;
! 493:
! 494: roundingMode = float_rounding_mode();
! 495: roundNearestEven = ( roundingMode == float_round_nearest_even );
! 496: roundIncrement = 0x200;
! 497: if ( ! roundNearestEven ) {
! 498: if ( roundingMode == float_round_to_zero ) {
! 499: roundIncrement = 0;
! 500: }
! 501: else {
! 502: roundIncrement = 0x3FF;
! 503: if ( zSign ) {
! 504: if ( roundingMode == float_round_up ) roundIncrement = 0;
! 505: }
! 506: else {
! 507: if ( roundingMode == float_round_down ) roundIncrement = 0;
! 508: }
! 509: }
! 510: }
! 511: roundBits = zSig & 0x3FF;
! 512: if ( 0x7FD <= (bits16) zExp ) {
! 513: if ( ( 0x7FD < zExp )
! 514: || ( ( zExp == 0x7FD )
! 515: && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
! 516: ) {
! 517: float_raise( float_flag_overflow | float_flag_inexact );
! 518: return FLOAT64_MANGLE(
! 519: FLOAT64_DEMANGLE(packFloat64( zSign, 0x7FF, 0 )) -
! 520: ( roundIncrement == 0 ));
! 521: }
! 522: if ( zExp < 0 ) {
! 523: isTiny =
! 524: ( float_detect_tininess == float_tininess_before_rounding )
! 525: || ( zExp < -1 )
! 526: || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
! 527: shift64RightJamming( zSig, - zExp, &zSig );
! 528: zExp = 0;
! 529: roundBits = zSig & 0x3FF;
! 530: if ( isTiny && roundBits ) float_raise( float_flag_underflow );
! 531: }
! 532: }
! 533: if ( roundBits ) float_set_inexact();
! 534: zSig = ( zSig + roundIncrement )>>10;
! 535: zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
! 536: if ( zSig == 0 ) zExp = 0;
! 537: return packFloat64( zSign, zExp, zSig );
! 538:
! 539: }
! 540:
! 541: /*
! 542: -------------------------------------------------------------------------------
! 543: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 544: and significand `zSig', and returns the proper double-precision floating-
! 545: point value corresponding to the abstract input. This routine is just like
! 546: `roundAndPackFloat64' except that `zSig' does not have to be normalized.
! 547: Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
! 548: floating-point exponent.
! 549: -------------------------------------------------------------------------------
! 550: */
! 551: static float64
! 552: normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig )
! 553: {
! 554: int8 shiftCount;
! 555:
! 556: shiftCount = countLeadingZeros64( zSig ) - 1;
! 557: return roundAndPackFloat64( zSign, zExp - shiftCount, zSig<<shiftCount );
! 558:
! 559: }
! 560:
! 561: #ifdef FLOATX80
! 562:
! 563: /*
! 564: -------------------------------------------------------------------------------
! 565: Returns the fraction bits of the extended double-precision floating-point
! 566: value `a'.
! 567: -------------------------------------------------------------------------------
! 568: */
! 569: INLINE bits64 extractFloatx80Frac( floatx80 a )
! 570: {
! 571:
! 572: return a.low;
! 573:
! 574: }
! 575:
! 576: /*
! 577: -------------------------------------------------------------------------------
! 578: Returns the exponent bits of the extended double-precision floating-point
! 579: value `a'.
! 580: -------------------------------------------------------------------------------
! 581: */
! 582: INLINE int32 extractFloatx80Exp( floatx80 a )
! 583: {
! 584:
! 585: return a.high & 0x7FFF;
! 586:
! 587: }
! 588:
! 589: /*
! 590: -------------------------------------------------------------------------------
! 591: Returns the sign bit of the extended double-precision floating-point value
! 592: `a'.
! 593: -------------------------------------------------------------------------------
! 594: */
! 595: INLINE flag extractFloatx80Sign( floatx80 a )
! 596: {
! 597:
! 598: return a.high>>15;
! 599:
! 600: }
! 601:
! 602: /*
! 603: -------------------------------------------------------------------------------
! 604: Normalizes the subnormal extended double-precision floating-point value
! 605: represented by the denormalized significand `aSig'. The normalized exponent
! 606: and significand are stored at the locations pointed to by `zExpPtr' and
! 607: `zSigPtr', respectively.
! 608: -------------------------------------------------------------------------------
! 609: */
! 610: static void
! 611: normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
! 612: {
! 613: int8 shiftCount;
! 614:
! 615: shiftCount = countLeadingZeros64( aSig );
! 616: *zSigPtr = aSig<<shiftCount;
! 617: *zExpPtr = 1 - shiftCount;
! 618:
! 619: }
! 620:
! 621: /*
! 622: -------------------------------------------------------------------------------
! 623: Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
! 624: extended double-precision floating-point value, returning the result.
! 625: -------------------------------------------------------------------------------
! 626: */
! 627: INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
! 628: {
! 629: floatx80 z;
! 630:
! 631: z.low = zSig;
! 632: z.high = ( ( (bits16) zSign )<<15 ) + zExp;
! 633: return z;
! 634:
! 635: }
! 636:
! 637: /*
! 638: -------------------------------------------------------------------------------
! 639: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 640: and extended significand formed by the concatenation of `zSig0' and `zSig1',
! 641: and returns the proper extended double-precision floating-point value
! 642: corresponding to the abstract input. Ordinarily, the abstract value is
! 643: rounded and packed into the extended double-precision format, with the
! 644: inexact exception raised if the abstract input cannot be represented
! 645: exactly. However, if the abstract value is too large, the overflow and
! 646: inexact exceptions are raised and an infinity or maximal finite value is
! 647: returned. If the abstract value is too small, the input value is rounded to
! 648: a subnormal number, and the underflow and inexact exceptions are raised if
! 649: the abstract input cannot be represented exactly as a subnormal extended
! 650: double-precision floating-point number.
! 651: If `roundingPrecision' is 32 or 64, the result is rounded to the same
! 652: number of bits as single or double precision, respectively. Otherwise, the
! 653: result is rounded to the full precision of the extended double-precision
! 654: format.
! 655: The input significand must be normalized or smaller. If the input
! 656: significand is not normalized, `zExp' must be 0; in that case, the result
! 657: returned is a subnormal number, and it must not require rounding. The
! 658: handling of underflow and overflow follows the IEC/IEEE Standard for Binary
! 659: Floating-Point Arithmetic.
! 660: -------------------------------------------------------------------------------
! 661: */
! 662: static floatx80
! 663: roundAndPackFloatx80(
! 664: int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
! 665: )
! 666: {
! 667: int8 roundingMode;
! 668: flag roundNearestEven, increment, isTiny;
! 669: int64 roundIncrement, roundMask, roundBits;
! 670:
! 671: roundingMode = float_rounding_mode();
! 672: roundNearestEven = ( roundingMode == float_round_nearest_even );
! 673: if ( roundingPrecision == 80 ) goto precision80;
! 674: if ( roundingPrecision == 64 ) {
! 675: roundIncrement = LIT64( 0x0000000000000400 );
! 676: roundMask = LIT64( 0x00000000000007FF );
! 677: }
! 678: else if ( roundingPrecision == 32 ) {
! 679: roundIncrement = LIT64( 0x0000008000000000 );
! 680: roundMask = LIT64( 0x000000FFFFFFFFFF );
! 681: }
! 682: else {
! 683: goto precision80;
! 684: }
! 685: zSig0 |= ( zSig1 != 0 );
! 686: if ( ! roundNearestEven ) {
! 687: if ( roundingMode == float_round_to_zero ) {
! 688: roundIncrement = 0;
! 689: }
! 690: else {
! 691: roundIncrement = roundMask;
! 692: if ( zSign ) {
! 693: if ( roundingMode == float_round_up ) roundIncrement = 0;
! 694: }
! 695: else {
! 696: if ( roundingMode == float_round_down ) roundIncrement = 0;
! 697: }
! 698: }
! 699: }
! 700: roundBits = zSig0 & roundMask;
! 701: if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
! 702: if ( ( 0x7FFE < zExp )
! 703: || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
! 704: ) {
! 705: goto overflow;
! 706: }
! 707: if ( zExp <= 0 ) {
! 708: isTiny =
! 709: ( float_detect_tininess == float_tininess_before_rounding )
! 710: || ( zExp < 0 )
! 711: || ( zSig0 <= zSig0 + roundIncrement );
! 712: shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
! 713: zExp = 0;
! 714: roundBits = zSig0 & roundMask;
! 715: if ( isTiny && roundBits ) float_raise( float_flag_underflow );
! 716: if ( roundBits ) float_set_inexact();
! 717: zSig0 += roundIncrement;
! 718: if ( (sbits64) zSig0 < 0 ) zExp = 1;
! 719: roundIncrement = roundMask + 1;
! 720: if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
! 721: roundMask |= roundIncrement;
! 722: }
! 723: zSig0 &= ~ roundMask;
! 724: return packFloatx80( zSign, zExp, zSig0 );
! 725: }
! 726: }
! 727: if ( roundBits ) float_set_inexact();
! 728: zSig0 += roundIncrement;
! 729: if ( zSig0 < roundIncrement ) {
! 730: ++zExp;
! 731: zSig0 = LIT64( 0x8000000000000000 );
! 732: }
! 733: roundIncrement = roundMask + 1;
! 734: if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
! 735: roundMask |= roundIncrement;
! 736: }
! 737: zSig0 &= ~ roundMask;
! 738: if ( zSig0 == 0 ) zExp = 0;
! 739: return packFloatx80( zSign, zExp, zSig0 );
! 740: precision80:
! 741: increment = ( (sbits64) zSig1 < 0 );
! 742: if ( ! roundNearestEven ) {
! 743: if ( roundingMode == float_round_to_zero ) {
! 744: increment = 0;
! 745: }
! 746: else {
! 747: if ( zSign ) {
! 748: increment = ( roundingMode == float_round_down ) && zSig1;
! 749: }
! 750: else {
! 751: increment = ( roundingMode == float_round_up ) && zSig1;
! 752: }
! 753: }
! 754: }
! 755: if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
! 756: if ( ( 0x7FFE < zExp )
! 757: || ( ( zExp == 0x7FFE )
! 758: && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
! 759: && increment
! 760: )
! 761: ) {
! 762: roundMask = 0;
! 763: overflow:
! 764: float_raise( float_flag_overflow | float_flag_inexact );
! 765: if ( ( roundingMode == float_round_to_zero )
! 766: || ( zSign && ( roundingMode == float_round_up ) )
! 767: || ( ! zSign && ( roundingMode == float_round_down ) )
! 768: ) {
! 769: return packFloatx80( zSign, 0x7FFE, ~ roundMask );
! 770: }
! 771: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 772: }
! 773: if ( zExp <= 0 ) {
! 774: isTiny =
! 775: ( float_detect_tininess == float_tininess_before_rounding )
! 776: || ( zExp < 0 )
! 777: || ! increment
! 778: || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
! 779: shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
! 780: zExp = 0;
! 781: if ( isTiny && zSig1 ) float_raise( float_flag_underflow );
! 782: if ( zSig1 ) float_set_inexact();
! 783: if ( roundNearestEven ) {
! 784: increment = ( (sbits64) zSig1 < 0 );
! 785: }
! 786: else {
! 787: if ( zSign ) {
! 788: increment = ( roundingMode == float_round_down ) && zSig1;
! 789: }
! 790: else {
! 791: increment = ( roundingMode == float_round_up ) && zSig1;
! 792: }
! 793: }
! 794: if ( increment ) {
! 795: ++zSig0;
! 796: zSig0 &=
! 797: ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
! 798: if ( (sbits64) zSig0 < 0 ) zExp = 1;
! 799: }
! 800: return packFloatx80( zSign, zExp, zSig0 );
! 801: }
! 802: }
! 803: if ( zSig1 ) float_set_inexact();
! 804: if ( increment ) {
! 805: ++zSig0;
! 806: if ( zSig0 == 0 ) {
! 807: ++zExp;
! 808: zSig0 = LIT64( 0x8000000000000000 );
! 809: }
! 810: else {
! 811: zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
! 812: }
! 813: }
! 814: else {
! 815: if ( zSig0 == 0 ) zExp = 0;
! 816: }
! 817: return packFloatx80( zSign, zExp, zSig0 );
! 818:
! 819: }
! 820:
! 821: /*
! 822: -------------------------------------------------------------------------------
! 823: Takes an abstract floating-point value having sign `zSign', exponent
! 824: `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
! 825: and returns the proper extended double-precision floating-point value
! 826: corresponding to the abstract input. This routine is just like
! 827: `roundAndPackFloatx80' except that the input significand does not have to be
! 828: normalized.
! 829: -------------------------------------------------------------------------------
! 830: */
! 831: static floatx80
! 832: normalizeRoundAndPackFloatx80(
! 833: int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
! 834: )
! 835: {
! 836: int8 shiftCount;
! 837:
! 838: if ( zSig0 == 0 ) {
! 839: zSig0 = zSig1;
! 840: zSig1 = 0;
! 841: zExp -= 64;
! 842: }
! 843: shiftCount = countLeadingZeros64( zSig0 );
! 844: shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
! 845: zExp -= shiftCount;
! 846: return
! 847: roundAndPackFloatx80( roundingPrecision, zSign, zExp, zSig0, zSig1 );
! 848:
! 849: }
! 850:
! 851: #endif
! 852:
! 853: #ifdef FLOAT128
! 854:
! 855: /*
! 856: -------------------------------------------------------------------------------
! 857: Returns the least-significant 64 fraction bits of the quadruple-precision
! 858: floating-point value `a'.
! 859: -------------------------------------------------------------------------------
! 860: */
! 861: INLINE bits64 extractFloat128Frac1( float128 a )
! 862: {
! 863:
! 864: return a.low;
! 865:
! 866: }
! 867:
! 868: /*
! 869: -------------------------------------------------------------------------------
! 870: Returns the most-significant 48 fraction bits of the quadruple-precision
! 871: floating-point value `a'.
! 872: -------------------------------------------------------------------------------
! 873: */
! 874: INLINE bits64 extractFloat128Frac0( float128 a )
! 875: {
! 876:
! 877: return a.high & LIT64( 0x0000FFFFFFFFFFFF );
! 878:
! 879: }
! 880:
! 881: /*
! 882: -------------------------------------------------------------------------------
! 883: Returns the exponent bits of the quadruple-precision floating-point value
! 884: `a'.
! 885: -------------------------------------------------------------------------------
! 886: */
! 887: INLINE int32 extractFloat128Exp( float128 a )
! 888: {
! 889:
! 890: return ( a.high>>48 ) & 0x7FFF;
! 891:
! 892: }
! 893:
! 894: /*
! 895: -------------------------------------------------------------------------------
! 896: Returns the sign bit of the quadruple-precision floating-point value `a'.
! 897: -------------------------------------------------------------------------------
! 898: */
! 899: INLINE flag extractFloat128Sign( float128 a )
! 900: {
! 901:
! 902: return a.high>>63;
! 903:
! 904: }
! 905:
! 906: /*
! 907: -------------------------------------------------------------------------------
! 908: Normalizes the subnormal quadruple-precision floating-point value
! 909: represented by the denormalized significand formed by the concatenation of
! 910: `aSig0' and `aSig1'. The normalized exponent is stored at the location
! 911: pointed to by `zExpPtr'. The most significant 49 bits of the normalized
! 912: significand are stored at the location pointed to by `zSig0Ptr', and the
! 913: least significant 64 bits of the normalized significand are stored at the
! 914: location pointed to by `zSig1Ptr'.
! 915: -------------------------------------------------------------------------------
! 916: */
! 917: static void
! 918: normalizeFloat128Subnormal(
! 919: bits64 aSig0,
! 920: bits64 aSig1,
! 921: int32 *zExpPtr,
! 922: bits64 *zSig0Ptr,
! 923: bits64 *zSig1Ptr
! 924: )
! 925: {
! 926: int8 shiftCount;
! 927:
! 928: if ( aSig0 == 0 ) {
! 929: shiftCount = countLeadingZeros64( aSig1 ) - 15;
! 930: if ( shiftCount < 0 ) {
! 931: *zSig0Ptr = aSig1>>( - shiftCount );
! 932: *zSig1Ptr = aSig1<<( shiftCount & 63 );
! 933: }
! 934: else {
! 935: *zSig0Ptr = aSig1<<shiftCount;
! 936: *zSig1Ptr = 0;
! 937: }
! 938: *zExpPtr = - shiftCount - 63;
! 939: }
! 940: else {
! 941: shiftCount = countLeadingZeros64( aSig0 ) - 15;
! 942: shortShift128Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );
! 943: *zExpPtr = 1 - shiftCount;
! 944: }
! 945:
! 946: }
! 947:
! 948: /*
! 949: -------------------------------------------------------------------------------
! 950: Packs the sign `zSign', the exponent `zExp', and the significand formed
! 951: by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
! 952: floating-point value, returning the result. After being shifted into the
! 953: proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
! 954: added together to form the most significant 32 bits of the result. This
! 955: means that any integer portion of `zSig0' will be added into the exponent.
! 956: Since a properly normalized significand will have an integer portion equal
! 957: to 1, the `zExp' input should be 1 less than the desired result exponent
! 958: whenever `zSig0' and `zSig1' concatenated form a complete, normalized
! 959: significand.
! 960: -------------------------------------------------------------------------------
! 961: */
! 962: INLINE float128
! 963: packFloat128( flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
! 964: {
! 965: float128 z;
! 966:
! 967: z.low = zSig1;
! 968: z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;
! 969: return z;
! 970:
! 971: }
! 972:
! 973: /*
! 974: -------------------------------------------------------------------------------
! 975: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 976: and extended significand formed by the concatenation of `zSig0', `zSig1',
! 977: and `zSig2', and returns the proper quadruple-precision floating-point value
! 978: corresponding to the abstract input. Ordinarily, the abstract value is
! 979: simply rounded and packed into the quadruple-precision format, with the
! 980: inexact exception raised if the abstract input cannot be represented
! 981: exactly. However, if the abstract value is too large, the overflow and
! 982: inexact exceptions are raised and an infinity or maximal finite value is
! 983: returned. If the abstract value is too small, the input value is rounded to
! 984: a subnormal number, and the underflow and inexact exceptions are raised if
! 985: the abstract input cannot be represented exactly as a subnormal quadruple-
! 986: precision floating-point number.
! 987: The input significand must be normalized or smaller. If the input
! 988: significand is not normalized, `zExp' must be 0; in that case, the result
! 989: returned is a subnormal number, and it must not require rounding. In the
! 990: usual case that the input significand is normalized, `zExp' must be 1 less
! 991: than the ``true'' floating-point exponent. The handling of underflow and
! 992: overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 993: -------------------------------------------------------------------------------
! 994: */
! 995: static float128
! 996: roundAndPackFloat128(
! 997: flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 )
! 998: {
! 999: int8 roundingMode;
! 1000: flag roundNearestEven, increment, isTiny;
! 1001:
! 1002: roundingMode = float_rounding_mode();
! 1003: roundNearestEven = ( roundingMode == float_round_nearest_even );
! 1004: increment = ( (sbits64) zSig2 < 0 );
! 1005: if ( ! roundNearestEven ) {
! 1006: if ( roundingMode == float_round_to_zero ) {
! 1007: increment = 0;
! 1008: }
! 1009: else {
! 1010: if ( zSign ) {
! 1011: increment = ( roundingMode == float_round_down ) && zSig2;
! 1012: }
! 1013: else {
! 1014: increment = ( roundingMode == float_round_up ) && zSig2;
! 1015: }
! 1016: }
! 1017: }
! 1018: if ( 0x7FFD <= (bits32) zExp ) {
! 1019: if ( ( 0x7FFD < zExp )
! 1020: || ( ( zExp == 0x7FFD )
! 1021: && eq128(
! 1022: LIT64( 0x0001FFFFFFFFFFFF ),
! 1023: LIT64( 0xFFFFFFFFFFFFFFFF ),
! 1024: zSig0,
! 1025: zSig1
! 1026: )
! 1027: && increment
! 1028: )
! 1029: ) {
! 1030: float_raise( float_flag_overflow | float_flag_inexact );
! 1031: if ( ( roundingMode == float_round_to_zero )
! 1032: || ( zSign && ( roundingMode == float_round_up ) )
! 1033: || ( ! zSign && ( roundingMode == float_round_down ) )
! 1034: ) {
! 1035: return
! 1036: packFloat128(
! 1037: zSign,
! 1038: 0x7FFE,
! 1039: LIT64( 0x0000FFFFFFFFFFFF ),
! 1040: LIT64( 0xFFFFFFFFFFFFFFFF )
! 1041: );
! 1042: }
! 1043: return packFloat128( zSign, 0x7FFF, 0, 0 );
! 1044: }
! 1045: if ( zExp < 0 ) {
! 1046: isTiny =
! 1047: ( float_detect_tininess == float_tininess_before_rounding )
! 1048: || ( zExp < -1 )
! 1049: || ! increment
! 1050: || lt128(
! 1051: zSig0,
! 1052: zSig1,
! 1053: LIT64( 0x0001FFFFFFFFFFFF ),
! 1054: LIT64( 0xFFFFFFFFFFFFFFFF )
! 1055: );
! 1056: shift128ExtraRightJamming(
! 1057: zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );
! 1058: zExp = 0;
! 1059: if ( isTiny && zSig2 ) float_raise( float_flag_underflow );
! 1060: if ( roundNearestEven ) {
! 1061: increment = ( (sbits64) zSig2 < 0 );
! 1062: }
! 1063: else {
! 1064: if ( zSign ) {
! 1065: increment = ( roundingMode == float_round_down ) && zSig2;
! 1066: }
! 1067: else {
! 1068: increment = ( roundingMode == float_round_up ) && zSig2;
! 1069: }
! 1070: }
! 1071: }
! 1072: }
! 1073: if ( zSig2 ) float_set_inexact();
! 1074: if ( increment ) {
! 1075: add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
! 1076: zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );
! 1077: }
! 1078: else {
! 1079: if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;
! 1080: }
! 1081: return packFloat128( zSign, zExp, zSig0, zSig1 );
! 1082:
! 1083: }
! 1084:
! 1085: /*
! 1086: -------------------------------------------------------------------------------
! 1087: Takes an abstract floating-point value having sign `zSign', exponent `zExp',
! 1088: and significand formed by the concatenation of `zSig0' and `zSig1', and
! 1089: returns the proper quadruple-precision floating-point value corresponding
! 1090: to the abstract input. This routine is just like `roundAndPackFloat128'
! 1091: except that the input significand has fewer bits and does not have to be
! 1092: normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
! 1093: point exponent.
! 1094: -------------------------------------------------------------------------------
! 1095: */
! 1096: static float128
! 1097: normalizeRoundAndPackFloat128(
! 1098: flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
! 1099: {
! 1100: int8 shiftCount;
! 1101: bits64 zSig2;
! 1102:
! 1103: if ( zSig0 == 0 ) {
! 1104: zSig0 = zSig1;
! 1105: zSig1 = 0;
! 1106: zExp -= 64;
! 1107: }
! 1108: shiftCount = countLeadingZeros64( zSig0 ) - 15;
! 1109: if ( 0 <= shiftCount ) {
! 1110: zSig2 = 0;
! 1111: shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
! 1112: }
! 1113: else {
! 1114: shift128ExtraRightJamming(
! 1115: zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );
! 1116: }
! 1117: zExp -= shiftCount;
! 1118: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
! 1119:
! 1120: }
! 1121:
! 1122: #endif
! 1123:
! 1124: /*
! 1125: -------------------------------------------------------------------------------
! 1126: Returns the result of converting the 32-bit two's complement integer `a'
! 1127: to the single-precision floating-point format. The conversion is performed
! 1128: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1129: -------------------------------------------------------------------------------
! 1130: */
! 1131: float32 int32_to_float32( int32 a )
! 1132: {
! 1133: flag zSign;
! 1134:
! 1135: if ( a == 0 ) return 0;
! 1136: if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
! 1137: zSign = ( a < 0 );
! 1138: return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a );
! 1139:
! 1140: }
! 1141:
! 1142: /*
! 1143: -------------------------------------------------------------------------------
! 1144: Returns the result of converting the 32-bit two's complement integer `a'
! 1145: to the double-precision floating-point format. The conversion is performed
! 1146: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1147: -------------------------------------------------------------------------------
! 1148: */
! 1149: float64 int32_to_float64( int32 a )
! 1150: {
! 1151: flag zSign;
! 1152: uint32 absA;
! 1153: int8 shiftCount;
! 1154: bits64 zSig;
! 1155:
! 1156: if ( a == 0 ) return 0;
! 1157: zSign = ( a < 0 );
! 1158: absA = zSign ? - a : a;
! 1159: shiftCount = countLeadingZeros32( absA ) + 21;
! 1160: zSig = absA;
! 1161: return packFloat64( zSign, 0x432 - shiftCount, zSig<<shiftCount );
! 1162:
! 1163: }
! 1164:
! 1165: #ifdef FLOATX80
! 1166:
! 1167: /*
! 1168: -------------------------------------------------------------------------------
! 1169: Returns the result of converting the 32-bit two's complement integer `a'
! 1170: to the extended double-precision floating-point format. The conversion
! 1171: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1172: Arithmetic.
! 1173: -------------------------------------------------------------------------------
! 1174: */
! 1175: floatx80 int32_to_floatx80( int32 a )
! 1176: {
! 1177: flag zSign;
! 1178: uint32 absA;
! 1179: int8 shiftCount;
! 1180: bits64 zSig;
! 1181:
! 1182: if ( a == 0 ) return packFloatx80( 0, 0, 0 );
! 1183: zSign = ( a < 0 );
! 1184: absA = zSign ? - a : a;
! 1185: shiftCount = countLeadingZeros32( absA ) + 32;
! 1186: zSig = absA;
! 1187: return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
! 1188:
! 1189: }
! 1190:
! 1191: #endif
! 1192:
! 1193: #ifdef FLOAT128
! 1194:
! 1195: /*
! 1196: -------------------------------------------------------------------------------
! 1197: Returns the result of converting the 32-bit two's complement integer `a' to
! 1198: the quadruple-precision floating-point format. The conversion is performed
! 1199: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1200: -------------------------------------------------------------------------------
! 1201: */
! 1202: float128 int32_to_float128( int32 a )
! 1203: {
! 1204: flag zSign;
! 1205: uint32 absA;
! 1206: int8 shiftCount;
! 1207: bits64 zSig0;
! 1208:
! 1209: if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
! 1210: zSign = ( a < 0 );
! 1211: absA = zSign ? - a : a;
! 1212: shiftCount = countLeadingZeros32( absA ) + 17;
! 1213: zSig0 = absA;
! 1214: return packFloat128( zSign, 0x402E - shiftCount, zSig0<<shiftCount, 0 );
! 1215:
! 1216: }
! 1217:
! 1218: #endif
! 1219:
! 1220: #ifndef SOFTFLOAT_FOR_GCC /* __floatdi?f is in libgcc2.c */
! 1221: /*
! 1222: -------------------------------------------------------------------------------
! 1223: Returns the result of converting the 64-bit two's complement integer `a'
! 1224: to the single-precision floating-point format. The conversion is performed
! 1225: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1226: -------------------------------------------------------------------------------
! 1227: */
! 1228: float32 int64_to_float32( int64 a )
! 1229: {
! 1230: flag zSign;
! 1231: uint64 absA;
! 1232: int8 shiftCount;
! 1233:
! 1234: if ( a == 0 ) return 0;
! 1235: zSign = ( a < 0 );
! 1236: absA = zSign ? - a : a;
! 1237: shiftCount = countLeadingZeros64( absA ) - 40;
! 1238: if ( 0 <= shiftCount ) {
! 1239: return packFloat32( zSign, 0x95 - shiftCount, absA<<shiftCount );
! 1240: }
! 1241: else {
! 1242: shiftCount += 7;
! 1243: if ( shiftCount < 0 ) {
! 1244: shift64RightJamming( absA, - shiftCount, &absA );
! 1245: }
! 1246: else {
! 1247: absA <<= shiftCount;
! 1248: }
! 1249: return roundAndPackFloat32( zSign, 0x9C - shiftCount, absA );
! 1250: }
! 1251:
! 1252: }
! 1253:
! 1254: /*
! 1255: -------------------------------------------------------------------------------
! 1256: Returns the result of converting the 64-bit two's complement integer `a'
! 1257: to the double-precision floating-point format. The conversion is performed
! 1258: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1259: -------------------------------------------------------------------------------
! 1260: */
! 1261: float64 int64_to_float64( int64 a )
! 1262: {
! 1263: flag zSign;
! 1264:
! 1265: if ( a == 0 ) return 0;
! 1266: if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
! 1267: return packFloat64( 1, 0x43E, 0 );
! 1268: }
! 1269: zSign = ( a < 0 );
! 1270: return normalizeRoundAndPackFloat64( zSign, 0x43C, zSign ? - a : a );
! 1271:
! 1272: }
! 1273:
! 1274: #ifdef FLOATX80
! 1275:
! 1276: /*
! 1277: -------------------------------------------------------------------------------
! 1278: Returns the result of converting the 64-bit two's complement integer `a'
! 1279: to the extended double-precision floating-point format. The conversion
! 1280: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1281: Arithmetic.
! 1282: -------------------------------------------------------------------------------
! 1283: */
! 1284: floatx80 int64_to_floatx80( int64 a )
! 1285: {
! 1286: flag zSign;
! 1287: uint64 absA;
! 1288: int8 shiftCount;
! 1289:
! 1290: if ( a == 0 ) return packFloatx80( 0, 0, 0 );
! 1291: zSign = ( a < 0 );
! 1292: absA = zSign ? - a : a;
! 1293: shiftCount = countLeadingZeros64( absA );
! 1294: return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );
! 1295:
! 1296: }
! 1297:
! 1298: #endif
! 1299:
! 1300: #ifdef FLOAT128
! 1301:
! 1302: /*
! 1303: -------------------------------------------------------------------------------
! 1304: Returns the result of converting the 64-bit two's complement integer `a' to
! 1305: the quadruple-precision floating-point format. The conversion is performed
! 1306: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1307: -------------------------------------------------------------------------------
! 1308: */
! 1309: float128 int64_to_float128( int64 a )
! 1310: {
! 1311: flag zSign;
! 1312: uint64 absA;
! 1313: int8 shiftCount;
! 1314: int32 zExp;
! 1315: bits64 zSig0, zSig1;
! 1316:
! 1317: if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
! 1318: zSign = ( a < 0 );
! 1319: absA = zSign ? - a : a;
! 1320: shiftCount = countLeadingZeros64( absA ) + 49;
! 1321: zExp = 0x406E - shiftCount;
! 1322: if ( 64 <= shiftCount ) {
! 1323: zSig1 = 0;
! 1324: zSig0 = absA;
! 1325: shiftCount -= 64;
! 1326: }
! 1327: else {
! 1328: zSig1 = absA;
! 1329: zSig0 = 0;
! 1330: }
! 1331: shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
! 1332: return packFloat128( zSign, zExp, zSig0, zSig1 );
! 1333:
! 1334: }
! 1335:
! 1336: #endif
! 1337: #endif /* !SOFTFLOAT_FOR_GCC */
! 1338:
! 1339: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 1340: /*
! 1341: -------------------------------------------------------------------------------
! 1342: Returns the result of converting the single-precision floating-point value
! 1343: `a' to the 32-bit two's complement integer format. The conversion is
! 1344: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1345: Arithmetic---which means in particular that the conversion is rounded
! 1346: according to the current rounding mode. If `a' is a NaN, the largest
! 1347: positive integer is returned. Otherwise, if the conversion overflows, the
! 1348: largest integer with the same sign as `a' is returned.
! 1349: -------------------------------------------------------------------------------
! 1350: */
! 1351: int32 float32_to_int32( float32 a )
! 1352: {
! 1353: flag aSign;
! 1354: int16 aExp, shiftCount;
! 1355: bits32 aSig;
! 1356: bits64 aSig64;
! 1357:
! 1358: aSig = extractFloat32Frac( a );
! 1359: aExp = extractFloat32Exp( a );
! 1360: aSign = extractFloat32Sign( a );
! 1361: if ( ( aExp == 0xFF ) && aSig ) aSign = 0;
! 1362: if ( aExp ) aSig |= 0x00800000;
! 1363: shiftCount = 0xAF - aExp;
! 1364: aSig64 = aSig;
! 1365: aSig64 <<= 32;
! 1366: if ( 0 < shiftCount ) shift64RightJamming( aSig64, shiftCount, &aSig64 );
! 1367: return roundAndPackInt32( aSign, aSig64 );
! 1368:
! 1369: }
! 1370: #endif /* !SOFTFLOAT_FOR_GCC */
! 1371:
! 1372: /*
! 1373: -------------------------------------------------------------------------------
! 1374: Returns the result of converting the single-precision floating-point value
! 1375: `a' to the 32-bit two's complement integer format. The conversion is
! 1376: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1377: Arithmetic, except that the conversion is always rounded toward zero.
! 1378: If `a' is a NaN, the largest positive integer is returned. Otherwise, if
! 1379: the conversion overflows, the largest integer with the same sign as `a' is
! 1380: returned.
! 1381: -------------------------------------------------------------------------------
! 1382: */
! 1383: int32 float32_to_int32_round_to_zero( float32 a )
! 1384: {
! 1385: flag aSign;
! 1386: int16 aExp, shiftCount;
! 1387: bits32 aSig;
! 1388: int32 z;
! 1389:
! 1390: aSig = extractFloat32Frac( a );
! 1391: aExp = extractFloat32Exp( a );
! 1392: aSign = extractFloat32Sign( a );
! 1393: shiftCount = aExp - 0x9E;
! 1394: if ( 0 <= shiftCount ) {
! 1395: if ( a != 0xCF000000 ) {
! 1396: float_raise( float_flag_invalid );
! 1397: if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
! 1398: }
! 1399: return (sbits32) 0x80000000;
! 1400: }
! 1401: else if ( aExp <= 0x7E ) {
! 1402: if ( aExp | aSig ) float_set_inexact();
! 1403: return 0;
! 1404: }
! 1405: aSig = ( aSig | 0x00800000 )<<8;
! 1406: z = aSig>>( - shiftCount );
! 1407: if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
! 1408: float_set_inexact();
! 1409: }
! 1410: if ( aSign ) z = - z;
! 1411: return z;
! 1412:
! 1413: }
! 1414:
! 1415: #ifndef SOFTFLOAT_FOR_GCC /* __fix?fdi provided by libgcc2.c */
! 1416: /*
! 1417: -------------------------------------------------------------------------------
! 1418: Returns the result of converting the single-precision floating-point value
! 1419: `a' to the 64-bit two's complement integer format. The conversion is
! 1420: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1421: Arithmetic---which means in particular that the conversion is rounded
! 1422: according to the current rounding mode. If `a' is a NaN, the largest
! 1423: positive integer is returned. Otherwise, if the conversion overflows, the
! 1424: largest integer with the same sign as `a' is returned.
! 1425: -------------------------------------------------------------------------------
! 1426: */
! 1427: int64 float32_to_int64( float32 a )
! 1428: {
! 1429: flag aSign;
! 1430: int16 aExp, shiftCount;
! 1431: bits32 aSig;
! 1432: bits64 aSig64, aSigExtra;
! 1433:
! 1434: aSig = extractFloat32Frac( a );
! 1435: aExp = extractFloat32Exp( a );
! 1436: aSign = extractFloat32Sign( a );
! 1437: shiftCount = 0xBE - aExp;
! 1438: if ( shiftCount < 0 ) {
! 1439: float_raise( float_flag_invalid );
! 1440: if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
! 1441: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 1442: }
! 1443: return (sbits64) LIT64( 0x8000000000000000 );
! 1444: }
! 1445: if ( aExp ) aSig |= 0x00800000;
! 1446: aSig64 = aSig;
! 1447: aSig64 <<= 40;
! 1448: shift64ExtraRightJamming( aSig64, 0, shiftCount, &aSig64, &aSigExtra );
! 1449: return roundAndPackInt64( aSign, aSig64, aSigExtra );
! 1450:
! 1451: }
! 1452:
! 1453: /*
! 1454: -------------------------------------------------------------------------------
! 1455: Returns the result of converting the single-precision floating-point value
! 1456: `a' to the 64-bit two's complement integer format. The conversion is
! 1457: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1458: Arithmetic, except that the conversion is always rounded toward zero. If
! 1459: `a' is a NaN, the largest positive integer is returned. Otherwise, if the
! 1460: conversion overflows, the largest integer with the same sign as `a' is
! 1461: returned.
! 1462: -------------------------------------------------------------------------------
! 1463: */
! 1464: int64 float32_to_int64_round_to_zero( float32 a )
! 1465: {
! 1466: flag aSign;
! 1467: int16 aExp, shiftCount;
! 1468: bits32 aSig;
! 1469: bits64 aSig64;
! 1470: int64 z;
! 1471:
! 1472: aSig = extractFloat32Frac( a );
! 1473: aExp = extractFloat32Exp( a );
! 1474: aSign = extractFloat32Sign( a );
! 1475: shiftCount = aExp - 0xBE;
! 1476: if ( 0 <= shiftCount ) {
! 1477: if ( a != 0xDF000000 ) {
! 1478: float_raise( float_flag_invalid );
! 1479: if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
! 1480: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 1481: }
! 1482: }
! 1483: return (sbits64) LIT64( 0x8000000000000000 );
! 1484: }
! 1485: else if ( aExp <= 0x7E ) {
! 1486: if ( aExp | aSig ) float_set_inexact();
! 1487: return 0;
! 1488: }
! 1489: aSig64 = aSig | 0x00800000;
! 1490: aSig64 <<= 40;
! 1491: z = aSig64>>( - shiftCount );
! 1492: if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {
! 1493: float_set_inexact();
! 1494: }
! 1495: if ( aSign ) z = - z;
! 1496: return z;
! 1497:
! 1498: }
! 1499: #endif /* !SOFTFLOAT_FOR_GCC */
! 1500:
! 1501: /*
! 1502: -------------------------------------------------------------------------------
! 1503: Returns the result of converting the single-precision floating-point value
! 1504: `a' to the double-precision floating-point format. The conversion is
! 1505: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1506: Arithmetic.
! 1507: -------------------------------------------------------------------------------
! 1508: */
! 1509: float64 float32_to_float64( float32 a )
! 1510: {
! 1511: flag aSign;
! 1512: int16 aExp;
! 1513: bits32 aSig;
! 1514:
! 1515: aSig = extractFloat32Frac( a );
! 1516: aExp = extractFloat32Exp( a );
! 1517: aSign = extractFloat32Sign( a );
! 1518: if ( aExp == 0xFF ) {
! 1519: if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) );
! 1520: return packFloat64( aSign, 0x7FF, 0 );
! 1521: }
! 1522: if ( aExp == 0 ) {
! 1523: if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
! 1524: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 1525: --aExp;
! 1526: }
! 1527: return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
! 1528:
! 1529: }
! 1530:
! 1531: #ifdef FLOATX80
! 1532:
! 1533: /*
! 1534: -------------------------------------------------------------------------------
! 1535: Returns the result of converting the single-precision floating-point value
! 1536: `a' to the extended double-precision floating-point format. The conversion
! 1537: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1538: Arithmetic.
! 1539: -------------------------------------------------------------------------------
! 1540: */
! 1541: floatx80 float32_to_floatx80( float32 a )
! 1542: {
! 1543: flag aSign;
! 1544: int16 aExp;
! 1545: bits32 aSig;
! 1546:
! 1547: aSig = extractFloat32Frac( a );
! 1548: aExp = extractFloat32Exp( a );
! 1549: aSign = extractFloat32Sign( a );
! 1550: if ( aExp == 0xFF ) {
! 1551: if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a ) );
! 1552: return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 1553: }
! 1554: if ( aExp == 0 ) {
! 1555: if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
! 1556: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 1557: }
! 1558: aSig |= 0x00800000;
! 1559: return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
! 1560:
! 1561: }
! 1562:
! 1563: #endif
! 1564:
! 1565: #ifdef FLOAT128
! 1566:
! 1567: /*
! 1568: -------------------------------------------------------------------------------
! 1569: Returns the result of converting the single-precision floating-point value
! 1570: `a' to the double-precision floating-point format. The conversion is
! 1571: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 1572: Arithmetic.
! 1573: -------------------------------------------------------------------------------
! 1574: */
! 1575: float128 float32_to_float128( float32 a )
! 1576: {
! 1577: flag aSign;
! 1578: int16 aExp;
! 1579: bits32 aSig;
! 1580:
! 1581: aSig = extractFloat32Frac( a );
! 1582: aExp = extractFloat32Exp( a );
! 1583: aSign = extractFloat32Sign( a );
! 1584: if ( aExp == 0xFF ) {
! 1585: if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a ) );
! 1586: return packFloat128( aSign, 0x7FFF, 0, 0 );
! 1587: }
! 1588: if ( aExp == 0 ) {
! 1589: if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
! 1590: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 1591: --aExp;
! 1592: }
! 1593: return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );
! 1594:
! 1595: }
! 1596:
! 1597: #endif
! 1598:
! 1599: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 1600: /*
! 1601: -------------------------------------------------------------------------------
! 1602: Rounds the single-precision floating-point value `a' to an integer, and
! 1603: returns the result as a single-precision floating-point value. The
! 1604: operation is performed according to the IEC/IEEE Standard for Binary
! 1605: Floating-Point Arithmetic.
! 1606: -------------------------------------------------------------------------------
! 1607: */
! 1608: float32 float32_round_to_int( float32 a )
! 1609: {
! 1610: flag aSign;
! 1611: int16 aExp;
! 1612: bits32 lastBitMask, roundBitsMask;
! 1613: int8 roundingMode;
! 1614: float32 z;
! 1615:
! 1616: aExp = extractFloat32Exp( a );
! 1617: if ( 0x96 <= aExp ) {
! 1618: if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
! 1619: return propagateFloat32NaN( a, a );
! 1620: }
! 1621: return a;
! 1622: }
! 1623: if ( aExp <= 0x7E ) {
! 1624: if ( (bits32) ( a<<1 ) == 0 ) return a;
! 1625: float_set_inexact();
! 1626: aSign = extractFloat32Sign( a );
! 1627: switch ( float_rounding_mode() ) {
! 1628: case float_round_nearest_even:
! 1629: if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
! 1630: return packFloat32( aSign, 0x7F, 0 );
! 1631: }
! 1632: break;
! 1633: case float_round_down:
! 1634: return aSign ? 0xBF800000 : 0;
! 1635: case float_round_up:
! 1636: return aSign ? 0x80000000 : 0x3F800000;
! 1637: }
! 1638: return packFloat32( aSign, 0, 0 );
! 1639: }
! 1640: lastBitMask = 1;
! 1641: lastBitMask <<= 0x96 - aExp;
! 1642: roundBitsMask = lastBitMask - 1;
! 1643: z = a;
! 1644: roundingMode = float_rounding_mode();
! 1645: if ( roundingMode == float_round_nearest_even ) {
! 1646: z += lastBitMask>>1;
! 1647: if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
! 1648: }
! 1649: else if ( roundingMode != float_round_to_zero ) {
! 1650: if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
! 1651: z += roundBitsMask;
! 1652: }
! 1653: }
! 1654: z &= ~ roundBitsMask;
! 1655: if ( z != a ) float_set_inexact();
! 1656: return z;
! 1657:
! 1658: }
! 1659: #endif /* !SOFTFLOAT_FOR_GCC */
! 1660:
! 1661: /*
! 1662: -------------------------------------------------------------------------------
! 1663: Returns the result of adding the absolute values of the single-precision
! 1664: floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
! 1665: before being returned. `zSign' is ignored if the result is a NaN.
! 1666: The addition is performed according to the IEC/IEEE Standard for Binary
! 1667: Floating-Point Arithmetic.
! 1668: -------------------------------------------------------------------------------
! 1669: */
! 1670: static float32 addFloat32Sigs( float32 a, float32 b, flag zSign )
! 1671: {
! 1672: int16 aExp, bExp, zExp;
! 1673: bits32 aSig, bSig, zSig;
! 1674: int16 expDiff;
! 1675:
! 1676: aSig = extractFloat32Frac( a );
! 1677: aExp = extractFloat32Exp( a );
! 1678: bSig = extractFloat32Frac( b );
! 1679: bExp = extractFloat32Exp( b );
! 1680: expDiff = aExp - bExp;
! 1681: aSig <<= 6;
! 1682: bSig <<= 6;
! 1683: if ( 0 < expDiff ) {
! 1684: if ( aExp == 0xFF ) {
! 1685: if ( aSig ) return propagateFloat32NaN( a, b );
! 1686: return a;
! 1687: }
! 1688: if ( bExp == 0 ) {
! 1689: --expDiff;
! 1690: }
! 1691: else {
! 1692: bSig |= 0x20000000;
! 1693: }
! 1694: shift32RightJamming( bSig, expDiff, &bSig );
! 1695: zExp = aExp;
! 1696: }
! 1697: else if ( expDiff < 0 ) {
! 1698: if ( bExp == 0xFF ) {
! 1699: if ( bSig ) return propagateFloat32NaN( a, b );
! 1700: return packFloat32( zSign, 0xFF, 0 );
! 1701: }
! 1702: if ( aExp == 0 ) {
! 1703: ++expDiff;
! 1704: }
! 1705: else {
! 1706: aSig |= 0x20000000;
! 1707: }
! 1708: shift32RightJamming( aSig, - expDiff, &aSig );
! 1709: zExp = bExp;
! 1710: }
! 1711: else {
! 1712: if ( aExp == 0xFF ) {
! 1713: if ( aSig | bSig ) return propagateFloat32NaN( a, b );
! 1714: return a;
! 1715: }
! 1716: if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
! 1717: zSig = 0x40000000 + aSig + bSig;
! 1718: zExp = aExp;
! 1719: goto roundAndPack;
! 1720: }
! 1721: aSig |= 0x20000000;
! 1722: zSig = ( aSig + bSig )<<1;
! 1723: --zExp;
! 1724: if ( (sbits32) zSig < 0 ) {
! 1725: zSig = aSig + bSig;
! 1726: ++zExp;
! 1727: }
! 1728: roundAndPack:
! 1729: return roundAndPackFloat32( zSign, zExp, zSig );
! 1730:
! 1731: }
! 1732:
! 1733: /*
! 1734: -------------------------------------------------------------------------------
! 1735: Returns the result of subtracting the absolute values of the single-
! 1736: precision floating-point values `a' and `b'. If `zSign' is 1, the
! 1737: difference is negated before being returned. `zSign' is ignored if the
! 1738: result is a NaN. The subtraction is performed according to the IEC/IEEE
! 1739: Standard for Binary Floating-Point Arithmetic.
! 1740: -------------------------------------------------------------------------------
! 1741: */
! 1742: static float32 subFloat32Sigs( float32 a, float32 b, flag zSign )
! 1743: {
! 1744: int16 aExp, bExp, zExp;
! 1745: bits32 aSig, bSig, zSig;
! 1746: int16 expDiff;
! 1747:
! 1748: aSig = extractFloat32Frac( a );
! 1749: aExp = extractFloat32Exp( a );
! 1750: bSig = extractFloat32Frac( b );
! 1751: bExp = extractFloat32Exp( b );
! 1752: expDiff = aExp - bExp;
! 1753: aSig <<= 7;
! 1754: bSig <<= 7;
! 1755: if ( 0 < expDiff ) goto aExpBigger;
! 1756: if ( expDiff < 0 ) goto bExpBigger;
! 1757: if ( aExp == 0xFF ) {
! 1758: if ( aSig | bSig ) return propagateFloat32NaN( a, b );
! 1759: float_raise( float_flag_invalid );
! 1760: return float32_default_nan;
! 1761: }
! 1762: if ( aExp == 0 ) {
! 1763: aExp = 1;
! 1764: bExp = 1;
! 1765: }
! 1766: if ( bSig < aSig ) goto aBigger;
! 1767: if ( aSig < bSig ) goto bBigger;
! 1768: return packFloat32( float_rounding_mode() == float_round_down, 0, 0 );
! 1769: bExpBigger:
! 1770: if ( bExp == 0xFF ) {
! 1771: if ( bSig ) return propagateFloat32NaN( a, b );
! 1772: return packFloat32( zSign ^ 1, 0xFF, 0 );
! 1773: }
! 1774: if ( aExp == 0 ) {
! 1775: ++expDiff;
! 1776: }
! 1777: else {
! 1778: aSig |= 0x40000000;
! 1779: }
! 1780: shift32RightJamming( aSig, - expDiff, &aSig );
! 1781: bSig |= 0x40000000;
! 1782: bBigger:
! 1783: zSig = bSig - aSig;
! 1784: zExp = bExp;
! 1785: zSign ^= 1;
! 1786: goto normalizeRoundAndPack;
! 1787: aExpBigger:
! 1788: if ( aExp == 0xFF ) {
! 1789: if ( aSig ) return propagateFloat32NaN( a, b );
! 1790: return a;
! 1791: }
! 1792: if ( bExp == 0 ) {
! 1793: --expDiff;
! 1794: }
! 1795: else {
! 1796: bSig |= 0x40000000;
! 1797: }
! 1798: shift32RightJamming( bSig, expDiff, &bSig );
! 1799: aSig |= 0x40000000;
! 1800: aBigger:
! 1801: zSig = aSig - bSig;
! 1802: zExp = aExp;
! 1803: normalizeRoundAndPack:
! 1804: --zExp;
! 1805: return normalizeRoundAndPackFloat32( zSign, zExp, zSig );
! 1806:
! 1807: }
! 1808:
! 1809: /*
! 1810: -------------------------------------------------------------------------------
! 1811: Returns the result of adding the single-precision floating-point values `a'
! 1812: and `b'. The operation is performed according to the IEC/IEEE Standard for
! 1813: Binary Floating-Point Arithmetic.
! 1814: -------------------------------------------------------------------------------
! 1815: */
! 1816: float32 float32_add( float32 a, float32 b )
! 1817: {
! 1818: flag aSign, bSign;
! 1819:
! 1820: aSign = extractFloat32Sign( a );
! 1821: bSign = extractFloat32Sign( b );
! 1822: if ( aSign == bSign ) {
! 1823: return addFloat32Sigs( a, b, aSign );
! 1824: }
! 1825: else {
! 1826: return subFloat32Sigs( a, b, aSign );
! 1827: }
! 1828:
! 1829: }
! 1830:
! 1831: /*
! 1832: -------------------------------------------------------------------------------
! 1833: Returns the result of subtracting the single-precision floating-point values
! 1834: `a' and `b'. The operation is performed according to the IEC/IEEE Standard
! 1835: for Binary Floating-Point Arithmetic.
! 1836: -------------------------------------------------------------------------------
! 1837: */
! 1838: float32 float32_sub( float32 a, float32 b )
! 1839: {
! 1840: flag aSign, bSign;
! 1841:
! 1842: aSign = extractFloat32Sign( a );
! 1843: bSign = extractFloat32Sign( b );
! 1844: if ( aSign == bSign ) {
! 1845: return subFloat32Sigs( a, b, aSign );
! 1846: }
! 1847: else {
! 1848: return addFloat32Sigs( a, b, aSign );
! 1849: }
! 1850:
! 1851: }
! 1852:
! 1853: /*
! 1854: -------------------------------------------------------------------------------
! 1855: Returns the result of multiplying the single-precision floating-point values
! 1856: `a' and `b'. The operation is performed according to the IEC/IEEE Standard
! 1857: for Binary Floating-Point Arithmetic.
! 1858: -------------------------------------------------------------------------------
! 1859: */
! 1860: float32 float32_mul( float32 a, float32 b )
! 1861: {
! 1862: flag aSign, bSign, zSign;
! 1863: int16 aExp, bExp, zExp;
! 1864: bits32 aSig, bSig;
! 1865: bits64 zSig64;
! 1866: bits32 zSig;
! 1867:
! 1868: aSig = extractFloat32Frac( a );
! 1869: aExp = extractFloat32Exp( a );
! 1870: aSign = extractFloat32Sign( a );
! 1871: bSig = extractFloat32Frac( b );
! 1872: bExp = extractFloat32Exp( b );
! 1873: bSign = extractFloat32Sign( b );
! 1874: zSign = aSign ^ bSign;
! 1875: if ( aExp == 0xFF ) {
! 1876: if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
! 1877: return propagateFloat32NaN( a, b );
! 1878: }
! 1879: if ( ( bExp | bSig ) == 0 ) {
! 1880: float_raise( float_flag_invalid );
! 1881: return float32_default_nan;
! 1882: }
! 1883: return packFloat32( zSign, 0xFF, 0 );
! 1884: }
! 1885: if ( bExp == 0xFF ) {
! 1886: if ( bSig ) return propagateFloat32NaN( a, b );
! 1887: if ( ( aExp | aSig ) == 0 ) {
! 1888: float_raise( float_flag_invalid );
! 1889: return float32_default_nan;
! 1890: }
! 1891: return packFloat32( zSign, 0xFF, 0 );
! 1892: }
! 1893: if ( aExp == 0 ) {
! 1894: if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
! 1895: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 1896: }
! 1897: if ( bExp == 0 ) {
! 1898: if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
! 1899: normalizeFloat32Subnormal( bSig, &bExp, &bSig );
! 1900: }
! 1901: zExp = aExp + bExp - 0x7F;
! 1902: aSig = ( aSig | 0x00800000 )<<7;
! 1903: bSig = ( bSig | 0x00800000 )<<8;
! 1904: shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
! 1905: zSig = zSig64;
! 1906: if ( 0 <= (sbits32) ( zSig<<1 ) ) {
! 1907: zSig <<= 1;
! 1908: --zExp;
! 1909: }
! 1910: return roundAndPackFloat32( zSign, zExp, zSig );
! 1911:
! 1912: }
! 1913:
! 1914: /*
! 1915: -------------------------------------------------------------------------------
! 1916: Returns the result of dividing the single-precision floating-point value `a'
! 1917: by the corresponding value `b'. The operation is performed according to the
! 1918: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1919: -------------------------------------------------------------------------------
! 1920: */
! 1921: float32 float32_div( float32 a, float32 b )
! 1922: {
! 1923: flag aSign, bSign, zSign;
! 1924: int16 aExp, bExp, zExp;
! 1925: bits32 aSig, bSig, zSig;
! 1926:
! 1927: aSig = extractFloat32Frac( a );
! 1928: aExp = extractFloat32Exp( a );
! 1929: aSign = extractFloat32Sign( a );
! 1930: bSig = extractFloat32Frac( b );
! 1931: bExp = extractFloat32Exp( b );
! 1932: bSign = extractFloat32Sign( b );
! 1933: zSign = aSign ^ bSign;
! 1934: if ( aExp == 0xFF ) {
! 1935: if ( aSig ) return propagateFloat32NaN( a, b );
! 1936: if ( bExp == 0xFF ) {
! 1937: if ( bSig ) return propagateFloat32NaN( a, b );
! 1938: float_raise( float_flag_invalid );
! 1939: return float32_default_nan;
! 1940: }
! 1941: return packFloat32( zSign, 0xFF, 0 );
! 1942: }
! 1943: if ( bExp == 0xFF ) {
! 1944: if ( bSig ) return propagateFloat32NaN( a, b );
! 1945: return packFloat32( zSign, 0, 0 );
! 1946: }
! 1947: if ( bExp == 0 ) {
! 1948: if ( bSig == 0 ) {
! 1949: if ( ( aExp | aSig ) == 0 ) {
! 1950: float_raise( float_flag_invalid );
! 1951: return float32_default_nan;
! 1952: }
! 1953: float_raise( float_flag_divbyzero );
! 1954: return packFloat32( zSign, 0xFF, 0 );
! 1955: }
! 1956: normalizeFloat32Subnormal( bSig, &bExp, &bSig );
! 1957: }
! 1958: if ( aExp == 0 ) {
! 1959: if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
! 1960: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 1961: }
! 1962: zExp = aExp - bExp + 0x7D;
! 1963: aSig = ( aSig | 0x00800000 )<<7;
! 1964: bSig = ( bSig | 0x00800000 )<<8;
! 1965: if ( bSig <= ( aSig + aSig ) ) {
! 1966: aSig >>= 1;
! 1967: ++zExp;
! 1968: }
! 1969: zSig = ( ( (bits64) aSig )<<32 ) / bSig;
! 1970: if ( ( zSig & 0x3F ) == 0 ) {
! 1971: zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );
! 1972: }
! 1973: return roundAndPackFloat32( zSign, zExp, zSig );
! 1974:
! 1975: }
! 1976:
! 1977: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 1978: /*
! 1979: -------------------------------------------------------------------------------
! 1980: Returns the remainder of the single-precision floating-point value `a'
! 1981: with respect to the corresponding value `b'. The operation is performed
! 1982: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 1983: -------------------------------------------------------------------------------
! 1984: */
! 1985: float32 float32_rem( float32 a, float32 b )
! 1986: {
! 1987: flag aSign, bSign, zSign;
! 1988: int16 aExp, bExp, expDiff;
! 1989: bits32 aSig, bSig;
! 1990: bits32 q;
! 1991: bits64 aSig64, bSig64, q64;
! 1992: bits32 alternateASig;
! 1993: sbits32 sigMean;
! 1994:
! 1995: aSig = extractFloat32Frac( a );
! 1996: aExp = extractFloat32Exp( a );
! 1997: aSign = extractFloat32Sign( a );
! 1998: bSig = extractFloat32Frac( b );
! 1999: bExp = extractFloat32Exp( b );
! 2000: bSign = extractFloat32Sign( b );
! 2001: if ( aExp == 0xFF ) {
! 2002: if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
! 2003: return propagateFloat32NaN( a, b );
! 2004: }
! 2005: float_raise( float_flag_invalid );
! 2006: return float32_default_nan;
! 2007: }
! 2008: if ( bExp == 0xFF ) {
! 2009: if ( bSig ) return propagateFloat32NaN( a, b );
! 2010: return a;
! 2011: }
! 2012: if ( bExp == 0 ) {
! 2013: if ( bSig == 0 ) {
! 2014: float_raise( float_flag_invalid );
! 2015: return float32_default_nan;
! 2016: }
! 2017: normalizeFloat32Subnormal( bSig, &bExp, &bSig );
! 2018: }
! 2019: if ( aExp == 0 ) {
! 2020: if ( aSig == 0 ) return a;
! 2021: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 2022: }
! 2023: expDiff = aExp - bExp;
! 2024: aSig |= 0x00800000;
! 2025: bSig |= 0x00800000;
! 2026: if ( expDiff < 32 ) {
! 2027: aSig <<= 8;
! 2028: bSig <<= 8;
! 2029: if ( expDiff < 0 ) {
! 2030: if ( expDiff < -1 ) return a;
! 2031: aSig >>= 1;
! 2032: }
! 2033: q = ( bSig <= aSig );
! 2034: if ( q ) aSig -= bSig;
! 2035: if ( 0 < expDiff ) {
! 2036: q = ( ( (bits64) aSig )<<32 ) / bSig;
! 2037: q >>= 32 - expDiff;
! 2038: bSig >>= 2;
! 2039: aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
! 2040: }
! 2041: else {
! 2042: aSig >>= 2;
! 2043: bSig >>= 2;
! 2044: }
! 2045: }
! 2046: else {
! 2047: if ( bSig <= aSig ) aSig -= bSig;
! 2048: aSig64 = ( (bits64) aSig )<<40;
! 2049: bSig64 = ( (bits64) bSig )<<40;
! 2050: expDiff -= 64;
! 2051: while ( 0 < expDiff ) {
! 2052: q64 = estimateDiv128To64( aSig64, 0, bSig64 );
! 2053: q64 = ( 2 < q64 ) ? q64 - 2 : 0;
! 2054: aSig64 = - ( ( bSig * q64 )<<38 );
! 2055: expDiff -= 62;
! 2056: }
! 2057: expDiff += 64;
! 2058: q64 = estimateDiv128To64( aSig64, 0, bSig64 );
! 2059: q64 = ( 2 < q64 ) ? q64 - 2 : 0;
! 2060: q = q64>>( 64 - expDiff );
! 2061: bSig <<= 6;
! 2062: aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;
! 2063: }
! 2064: do {
! 2065: alternateASig = aSig;
! 2066: ++q;
! 2067: aSig -= bSig;
! 2068: } while ( 0 <= (sbits32) aSig );
! 2069: sigMean = aSig + alternateASig;
! 2070: if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
! 2071: aSig = alternateASig;
! 2072: }
! 2073: zSign = ( (sbits32) aSig < 0 );
! 2074: if ( zSign ) aSig = - aSig;
! 2075: return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig );
! 2076:
! 2077: }
! 2078: #endif /* !SOFTFLOAT_FOR_GCC */
! 2079:
! 2080: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 2081: /*
! 2082: -------------------------------------------------------------------------------
! 2083: Returns the square root of the single-precision floating-point value `a'.
! 2084: The operation is performed according to the IEC/IEEE Standard for Binary
! 2085: Floating-Point Arithmetic.
! 2086: -------------------------------------------------------------------------------
! 2087: */
! 2088: float32 float32_sqrt( float32 a )
! 2089: {
! 2090: flag aSign;
! 2091: int16 aExp, zExp;
! 2092: bits32 aSig, zSig;
! 2093: bits64 rem, term;
! 2094:
! 2095: aSig = extractFloat32Frac( a );
! 2096: aExp = extractFloat32Exp( a );
! 2097: aSign = extractFloat32Sign( a );
! 2098: if ( aExp == 0xFF ) {
! 2099: if ( aSig ) return propagateFloat32NaN( a, 0 );
! 2100: if ( ! aSign ) return a;
! 2101: float_raise( float_flag_invalid );
! 2102: return float32_default_nan;
! 2103: }
! 2104: if ( aSign ) {
! 2105: if ( ( aExp | aSig ) == 0 ) return a;
! 2106: float_raise( float_flag_invalid );
! 2107: return float32_default_nan;
! 2108: }
! 2109: if ( aExp == 0 ) {
! 2110: if ( aSig == 0 ) return 0;
! 2111: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
! 2112: }
! 2113: zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
! 2114: aSig = ( aSig | 0x00800000 )<<8;
! 2115: zSig = estimateSqrt32( aExp, aSig ) + 2;
! 2116: if ( ( zSig & 0x7F ) <= 5 ) {
! 2117: if ( zSig < 2 ) {
! 2118: zSig = 0x7FFFFFFF;
! 2119: goto roundAndPack;
! 2120: }
! 2121: aSig >>= aExp & 1;
! 2122: term = ( (bits64) zSig ) * zSig;
! 2123: rem = ( ( (bits64) aSig )<<32 ) - term;
! 2124: while ( (sbits64) rem < 0 ) {
! 2125: --zSig;
! 2126: rem += ( ( (bits64) zSig )<<1 ) | 1;
! 2127: }
! 2128: zSig |= ( rem != 0 );
! 2129: }
! 2130: shift32RightJamming( zSig, 1, &zSig );
! 2131: roundAndPack:
! 2132: return roundAndPackFloat32( 0, zExp, zSig );
! 2133:
! 2134: }
! 2135: #endif /* !SOFTFLOAT_FOR_GCC */
! 2136:
! 2137: /*
! 2138: -------------------------------------------------------------------------------
! 2139: Returns 1 if the single-precision floating-point value `a' is equal to
! 2140: the corresponding value `b', and 0 otherwise. The comparison is performed
! 2141: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 2142: -------------------------------------------------------------------------------
! 2143: */
! 2144: flag float32_eq( float32 a, float32 b )
! 2145: {
! 2146:
! 2147: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
! 2148: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
! 2149: ) {
! 2150: if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
! 2151: float_raise( float_flag_invalid );
! 2152: }
! 2153: return 0;
! 2154: }
! 2155: return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
! 2156:
! 2157: }
! 2158:
! 2159: /*
! 2160: -------------------------------------------------------------------------------
! 2161: Returns 1 if the single-precision floating-point value `a' is less than
! 2162: or equal to the corresponding value `b', and 0 otherwise. The comparison
! 2163: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2164: Arithmetic.
! 2165: -------------------------------------------------------------------------------
! 2166: */
! 2167: flag float32_le( float32 a, float32 b )
! 2168: {
! 2169: flag aSign, bSign;
! 2170:
! 2171: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
! 2172: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
! 2173: ) {
! 2174: float_raise( float_flag_invalid );
! 2175: return 0;
! 2176: }
! 2177: aSign = extractFloat32Sign( a );
! 2178: bSign = extractFloat32Sign( b );
! 2179: if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
! 2180: return ( a == b ) || ( aSign ^ ( a < b ) );
! 2181:
! 2182: }
! 2183:
! 2184: /*
! 2185: -------------------------------------------------------------------------------
! 2186: Returns 1 if the single-precision floating-point value `a' is less than
! 2187: the corresponding value `b', and 0 otherwise. The comparison is performed
! 2188: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 2189: -------------------------------------------------------------------------------
! 2190: */
! 2191: flag float32_lt( float32 a, float32 b )
! 2192: {
! 2193: flag aSign, bSign;
! 2194:
! 2195: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
! 2196: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
! 2197: ) {
! 2198: float_raise( float_flag_invalid );
! 2199: return 0;
! 2200: }
! 2201: aSign = extractFloat32Sign( a );
! 2202: bSign = extractFloat32Sign( b );
! 2203: if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
! 2204: return ( a != b ) && ( aSign ^ ( a < b ) );
! 2205:
! 2206: }
! 2207:
! 2208: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 2209: /*
! 2210: -------------------------------------------------------------------------------
! 2211: Returns 1 if the single-precision floating-point value `a' is equal to
! 2212: the corresponding value `b', and 0 otherwise. The invalid exception is
! 2213: raised if either operand is a NaN. Otherwise, the comparison is performed
! 2214: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 2215: -------------------------------------------------------------------------------
! 2216: */
! 2217: flag float32_eq_signaling( float32 a, float32 b )
! 2218: {
! 2219:
! 2220: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
! 2221: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
! 2222: ) {
! 2223: float_raise( float_flag_invalid );
! 2224: return 0;
! 2225: }
! 2226: return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
! 2227:
! 2228: }
! 2229:
! 2230: /*
! 2231: -------------------------------------------------------------------------------
! 2232: Returns 1 if the single-precision floating-point value `a' is less than or
! 2233: equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
! 2234: cause an exception. Otherwise, the comparison is performed according to the
! 2235: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 2236: -------------------------------------------------------------------------------
! 2237: */
! 2238: flag float32_le_quiet( float32 a, float32 b )
! 2239: {
! 2240: flag aSign, bSign;
! 2241:
! 2242: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
! 2243: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
! 2244: ) {
! 2245: if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
! 2246: float_raise( float_flag_invalid );
! 2247: }
! 2248: return 0;
! 2249: }
! 2250: aSign = extractFloat32Sign( a );
! 2251: bSign = extractFloat32Sign( b );
! 2252: if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
! 2253: return ( a == b ) || ( aSign ^ ( a < b ) );
! 2254:
! 2255: }
! 2256:
! 2257: /*
! 2258: -------------------------------------------------------------------------------
! 2259: Returns 1 if the single-precision floating-point value `a' is less than
! 2260: the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
! 2261: exception. Otherwise, the comparison is performed according to the IEC/IEEE
! 2262: Standard for Binary Floating-Point Arithmetic.
! 2263: -------------------------------------------------------------------------------
! 2264: */
! 2265: flag float32_lt_quiet( float32 a, float32 b )
! 2266: {
! 2267: flag aSign, bSign;
! 2268:
! 2269: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
! 2270: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
! 2271: ) {
! 2272: if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
! 2273: float_raise( float_flag_invalid );
! 2274: }
! 2275: return 0;
! 2276: }
! 2277: aSign = extractFloat32Sign( a );
! 2278: bSign = extractFloat32Sign( b );
! 2279: if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
! 2280: return ( a != b ) && ( aSign ^ ( a < b ) );
! 2281:
! 2282: }
! 2283: #endif /* !SOFTFLOAT_FOR_GCC */
! 2284:
! 2285: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 2286: /*
! 2287: -------------------------------------------------------------------------------
! 2288: Returns the result of converting the double-precision floating-point value
! 2289: `a' to the 32-bit two's complement integer format. The conversion is
! 2290: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2291: Arithmetic---which means in particular that the conversion is rounded
! 2292: according to the current rounding mode. If `a' is a NaN, the largest
! 2293: positive integer is returned. Otherwise, if the conversion overflows, the
! 2294: largest integer with the same sign as `a' is returned.
! 2295: -------------------------------------------------------------------------------
! 2296: */
! 2297: int32 float64_to_int32( float64 a )
! 2298: {
! 2299: flag aSign;
! 2300: int16 aExp, shiftCount;
! 2301: bits64 aSig;
! 2302:
! 2303: aSig = extractFloat64Frac( a );
! 2304: aExp = extractFloat64Exp( a );
! 2305: aSign = extractFloat64Sign( a );
! 2306: if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
! 2307: if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
! 2308: shiftCount = 0x42C - aExp;
! 2309: if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
! 2310: return roundAndPackInt32( aSign, aSig );
! 2311:
! 2312: }
! 2313: #endif /* !SOFTFLOAT_FOR_GCC */
! 2314:
! 2315: /*
! 2316: -------------------------------------------------------------------------------
! 2317: Returns the result of converting the double-precision floating-point value
! 2318: `a' to the 32-bit two's complement integer format. The conversion is
! 2319: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2320: Arithmetic, except that the conversion is always rounded toward zero.
! 2321: If `a' is a NaN, the largest positive integer is returned. Otherwise, if
! 2322: the conversion overflows, the largest integer with the same sign as `a' is
! 2323: returned.
! 2324: -------------------------------------------------------------------------------
! 2325: */
! 2326: int32 float64_to_int32_round_to_zero( float64 a )
! 2327: {
! 2328: flag aSign;
! 2329: int16 aExp, shiftCount;
! 2330: bits64 aSig, savedASig;
! 2331: int32 z;
! 2332:
! 2333: aSig = extractFloat64Frac( a );
! 2334: aExp = extractFloat64Exp( a );
! 2335: aSign = extractFloat64Sign( a );
! 2336: if ( 0x41E < aExp ) {
! 2337: if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
! 2338: goto invalid;
! 2339: }
! 2340: else if ( aExp < 0x3FF ) {
! 2341: if ( aExp || aSig ) float_set_inexact();
! 2342: return 0;
! 2343: }
! 2344: aSig |= LIT64( 0x0010000000000000 );
! 2345: shiftCount = 0x433 - aExp;
! 2346: savedASig = aSig;
! 2347: aSig >>= shiftCount;
! 2348: z = aSig;
! 2349: if ( aSign ) z = - z;
! 2350: if ( ( z < 0 ) ^ aSign ) {
! 2351: invalid:
! 2352: float_raise( float_flag_invalid );
! 2353: return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
! 2354: }
! 2355: if ( ( aSig<<shiftCount ) != savedASig ) {
! 2356: float_set_inexact();
! 2357: }
! 2358: return z;
! 2359:
! 2360: }
! 2361:
! 2362: #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
! 2363: /*
! 2364: -------------------------------------------------------------------------------
! 2365: Returns the result of converting the double-precision floating-point value
! 2366: `a' to the 64-bit two's complement integer format. The conversion is
! 2367: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2368: Arithmetic---which means in particular that the conversion is rounded
! 2369: according to the current rounding mode. If `a' is a NaN, the largest
! 2370: positive integer is returned. Otherwise, if the conversion overflows, the
! 2371: largest integer with the same sign as `a' is returned.
! 2372: -------------------------------------------------------------------------------
! 2373: */
! 2374: int64 float64_to_int64( float64 a )
! 2375: {
! 2376: flag aSign;
! 2377: int16 aExp, shiftCount;
! 2378: bits64 aSig, aSigExtra;
! 2379:
! 2380: aSig = extractFloat64Frac( a );
! 2381: aExp = extractFloat64Exp( a );
! 2382: aSign = extractFloat64Sign( a );
! 2383: if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
! 2384: shiftCount = 0x433 - aExp;
! 2385: if ( shiftCount <= 0 ) {
! 2386: if ( 0x43E < aExp ) {
! 2387: float_raise( float_flag_invalid );
! 2388: if ( ! aSign
! 2389: || ( ( aExp == 0x7FF )
! 2390: && ( aSig != LIT64( 0x0010000000000000 ) ) )
! 2391: ) {
! 2392: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 2393: }
! 2394: return (sbits64) LIT64( 0x8000000000000000 );
! 2395: }
! 2396: aSigExtra = 0;
! 2397: aSig <<= - shiftCount;
! 2398: }
! 2399: else {
! 2400: shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
! 2401: }
! 2402: return roundAndPackInt64( aSign, aSig, aSigExtra );
! 2403:
! 2404: }
! 2405:
! 2406: /*
! 2407: -------------------------------------------------------------------------------
! 2408: Returns the result of converting the double-precision floating-point value
! 2409: `a' to the 64-bit two's complement integer format. The conversion is
! 2410: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2411: Arithmetic, except that the conversion is always rounded toward zero.
! 2412: If `a' is a NaN, the largest positive integer is returned. Otherwise, if
! 2413: the conversion overflows, the largest integer with the same sign as `a' is
! 2414: returned.
! 2415: -------------------------------------------------------------------------------
! 2416: */
! 2417: int64 float64_to_int64_round_to_zero( float64 a )
! 2418: {
! 2419: flag aSign;
! 2420: int16 aExp, shiftCount;
! 2421: bits64 aSig;
! 2422: int64 z;
! 2423:
! 2424: aSig = extractFloat64Frac( a );
! 2425: aExp = extractFloat64Exp( a );
! 2426: aSign = extractFloat64Sign( a );
! 2427: if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
! 2428: shiftCount = aExp - 0x433;
! 2429: if ( 0 <= shiftCount ) {
! 2430: if ( 0x43E <= aExp ) {
! 2431: if ( a != LIT64( 0xC3E0000000000000 ) ) {
! 2432: float_raise( float_flag_invalid );
! 2433: if ( ! aSign
! 2434: || ( ( aExp == 0x7FF )
! 2435: && ( aSig != LIT64( 0x0010000000000000 ) ) )
! 2436: ) {
! 2437: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 2438: }
! 2439: }
! 2440: return (sbits64) LIT64( 0x8000000000000000 );
! 2441: }
! 2442: z = aSig<<shiftCount;
! 2443: }
! 2444: else {
! 2445: if ( aExp < 0x3FE ) {
! 2446: if ( aExp | aSig ) float_set_inexact();
! 2447: return 0;
! 2448: }
! 2449: z = aSig>>( - shiftCount );
! 2450: if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
! 2451: float_set_inexact();
! 2452: }
! 2453: }
! 2454: if ( aSign ) z = - z;
! 2455: return z;
! 2456:
! 2457: }
! 2458: #endif /* !SOFTFLOAT_FOR_GCC */
! 2459:
! 2460: /*
! 2461: -------------------------------------------------------------------------------
! 2462: Returns the result of converting the double-precision floating-point value
! 2463: `a' to the single-precision floating-point format. The conversion is
! 2464: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2465: Arithmetic.
! 2466: -------------------------------------------------------------------------------
! 2467: */
! 2468: float32 float64_to_float32( float64 a )
! 2469: {
! 2470: flag aSign;
! 2471: int16 aExp;
! 2472: bits64 aSig;
! 2473: bits32 zSig;
! 2474:
! 2475: aSig = extractFloat64Frac( a );
! 2476: aExp = extractFloat64Exp( a );
! 2477: aSign = extractFloat64Sign( a );
! 2478: if ( aExp == 0x7FF ) {
! 2479: if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a ) );
! 2480: return packFloat32( aSign, 0xFF, 0 );
! 2481: }
! 2482: shift64RightJamming( aSig, 22, &aSig );
! 2483: zSig = aSig;
! 2484: if ( aExp || zSig ) {
! 2485: zSig |= 0x40000000;
! 2486: aExp -= 0x381;
! 2487: }
! 2488: return roundAndPackFloat32( aSign, aExp, zSig );
! 2489:
! 2490: }
! 2491:
! 2492: #ifdef FLOATX80
! 2493:
! 2494: /*
! 2495: -------------------------------------------------------------------------------
! 2496: Returns the result of converting the double-precision floating-point value
! 2497: `a' to the extended double-precision floating-point format. The conversion
! 2498: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2499: Arithmetic.
! 2500: -------------------------------------------------------------------------------
! 2501: */
! 2502: floatx80 float64_to_floatx80( float64 a )
! 2503: {
! 2504: flag aSign;
! 2505: int16 aExp;
! 2506: bits64 aSig;
! 2507:
! 2508: aSig = extractFloat64Frac( a );
! 2509: aExp = extractFloat64Exp( a );
! 2510: aSign = extractFloat64Sign( a );
! 2511: if ( aExp == 0x7FF ) {
! 2512: if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a ) );
! 2513: return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 2514: }
! 2515: if ( aExp == 0 ) {
! 2516: if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
! 2517: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
! 2518: }
! 2519: return
! 2520: packFloatx80(
! 2521: aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
! 2522:
! 2523: }
! 2524:
! 2525: #endif
! 2526:
! 2527: #ifdef FLOAT128
! 2528:
! 2529: /*
! 2530: -------------------------------------------------------------------------------
! 2531: Returns the result of converting the double-precision floating-point value
! 2532: `a' to the quadruple-precision floating-point format. The conversion is
! 2533: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 2534: Arithmetic.
! 2535: -------------------------------------------------------------------------------
! 2536: */
! 2537: float128 float64_to_float128( float64 a )
! 2538: {
! 2539: flag aSign;
! 2540: int16 aExp;
! 2541: bits64 aSig, zSig0, zSig1;
! 2542:
! 2543: aSig = extractFloat64Frac( a );
! 2544: aExp = extractFloat64Exp( a );
! 2545: aSign = extractFloat64Sign( a );
! 2546: if ( aExp == 0x7FF ) {
! 2547: if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a ) );
! 2548: return packFloat128( aSign, 0x7FFF, 0, 0 );
! 2549: }
! 2550: if ( aExp == 0 ) {
! 2551: if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
! 2552: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
! 2553: --aExp;
! 2554: }
! 2555: shift128Right( aSig, 0, 4, &zSig0, &zSig1 );
! 2556: return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );
! 2557:
! 2558: }
! 2559:
! 2560: #endif
! 2561:
! 2562: #ifndef SOFTFLOAT_FOR_GCC
! 2563: /*
! 2564: -------------------------------------------------------------------------------
! 2565: Rounds the double-precision floating-point value `a' to an integer, and
! 2566: returns the result as a double-precision floating-point value. The
! 2567: operation is performed according to the IEC/IEEE Standard for Binary
! 2568: Floating-Point Arithmetic.
! 2569: -------------------------------------------------------------------------------
! 2570: */
! 2571: float64 float64_round_to_int( float64 a )
! 2572: {
! 2573: flag aSign;
! 2574: int16 aExp;
! 2575: bits64 lastBitMask, roundBitsMask;
! 2576: int8 roundingMode;
! 2577: float64 z;
! 2578:
! 2579: aExp = extractFloat64Exp( a );
! 2580: if ( 0x433 <= aExp ) {
! 2581: if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {
! 2582: return propagateFloat64NaN( a, a );
! 2583: }
! 2584: return a;
! 2585: }
! 2586: if ( aExp < 0x3FF ) {
! 2587: if ( (bits64) ( a<<1 ) == 0 ) return a;
! 2588: float_set_inexact();
! 2589: aSign = extractFloat64Sign( a );
! 2590: switch ( float_rounding_mode() ) {
! 2591: case float_round_nearest_even:
! 2592: if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {
! 2593: return packFloat64( aSign, 0x3FF, 0 );
! 2594: }
! 2595: break;
! 2596: case float_round_down:
! 2597: return aSign ? LIT64( 0xBFF0000000000000 ) : 0;
! 2598: case float_round_up:
! 2599: return
! 2600: aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );
! 2601: }
! 2602: return packFloat64( aSign, 0, 0 );
! 2603: }
! 2604: lastBitMask = 1;
! 2605: lastBitMask <<= 0x433 - aExp;
! 2606: roundBitsMask = lastBitMask - 1;
! 2607: z = a;
! 2608: roundingMode = float_rounding_mode();
! 2609: if ( roundingMode == float_round_nearest_even ) {
! 2610: z += lastBitMask>>1;
! 2611: if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
! 2612: }
! 2613: else if ( roundingMode != float_round_to_zero ) {
! 2614: if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {
! 2615: z += roundBitsMask;
! 2616: }
! 2617: }
! 2618: z &= ~ roundBitsMask;
! 2619: if ( z != a ) float_set_inexact();
! 2620: return z;
! 2621:
! 2622: }
! 2623: #endif
! 2624:
! 2625: /*
! 2626: -------------------------------------------------------------------------------
! 2627: Returns the result of adding the absolute values of the double-precision
! 2628: floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
! 2629: before being returned. `zSign' is ignored if the result is a NaN.
! 2630: The addition is performed according to the IEC/IEEE Standard for Binary
! 2631: Floating-Point Arithmetic.
! 2632: -------------------------------------------------------------------------------
! 2633: */
! 2634: static float64 addFloat64Sigs( float64 a, float64 b, flag zSign )
! 2635: {
! 2636: int16 aExp, bExp, zExp;
! 2637: bits64 aSig, bSig, zSig;
! 2638: int16 expDiff;
! 2639:
! 2640: aSig = extractFloat64Frac( a );
! 2641: aExp = extractFloat64Exp( a );
! 2642: bSig = extractFloat64Frac( b );
! 2643: bExp = extractFloat64Exp( b );
! 2644: expDiff = aExp - bExp;
! 2645: aSig <<= 9;
! 2646: bSig <<= 9;
! 2647: if ( 0 < expDiff ) {
! 2648: if ( aExp == 0x7FF ) {
! 2649: if ( aSig ) return propagateFloat64NaN( a, b );
! 2650: return a;
! 2651: }
! 2652: if ( bExp == 0 ) {
! 2653: --expDiff;
! 2654: }
! 2655: else {
! 2656: bSig |= LIT64( 0x2000000000000000 );
! 2657: }
! 2658: shift64RightJamming( bSig, expDiff, &bSig );
! 2659: zExp = aExp;
! 2660: }
! 2661: else if ( expDiff < 0 ) {
! 2662: if ( bExp == 0x7FF ) {
! 2663: if ( bSig ) return propagateFloat64NaN( a, b );
! 2664: return packFloat64( zSign, 0x7FF, 0 );
! 2665: }
! 2666: if ( aExp == 0 ) {
! 2667: ++expDiff;
! 2668: }
! 2669: else {
! 2670: aSig |= LIT64( 0x2000000000000000 );
! 2671: }
! 2672: shift64RightJamming( aSig, - expDiff, &aSig );
! 2673: zExp = bExp;
! 2674: }
! 2675: else {
! 2676: if ( aExp == 0x7FF ) {
! 2677: if ( aSig | bSig ) return propagateFloat64NaN( a, b );
! 2678: return a;
! 2679: }
! 2680: if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
! 2681: zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
! 2682: zExp = aExp;
! 2683: goto roundAndPack;
! 2684: }
! 2685: aSig |= LIT64( 0x2000000000000000 );
! 2686: zSig = ( aSig + bSig )<<1;
! 2687: --zExp;
! 2688: if ( (sbits64) zSig < 0 ) {
! 2689: zSig = aSig + bSig;
! 2690: ++zExp;
! 2691: }
! 2692: roundAndPack:
! 2693: return roundAndPackFloat64( zSign, zExp, zSig );
! 2694:
! 2695: }
! 2696:
! 2697: /*
! 2698: -------------------------------------------------------------------------------
! 2699: Returns the result of subtracting the absolute values of the double-
! 2700: precision floating-point values `a' and `b'. If `zSign' is 1, the
! 2701: difference is negated before being returned. `zSign' is ignored if the
! 2702: result is a NaN. The subtraction is performed according to the IEC/IEEE
! 2703: Standard for Binary Floating-Point Arithmetic.
! 2704: -------------------------------------------------------------------------------
! 2705: */
! 2706: static float64 subFloat64Sigs( float64 a, float64 b, flag zSign )
! 2707: {
! 2708: int16 aExp, bExp, zExp;
! 2709: bits64 aSig, bSig, zSig;
! 2710: int16 expDiff;
! 2711:
! 2712: aSig = extractFloat64Frac( a );
! 2713: aExp = extractFloat64Exp( a );
! 2714: bSig = extractFloat64Frac( b );
! 2715: bExp = extractFloat64Exp( b );
! 2716: expDiff = aExp - bExp;
! 2717: aSig <<= 10;
! 2718: bSig <<= 10;
! 2719: if ( 0 < expDiff ) goto aExpBigger;
! 2720: if ( expDiff < 0 ) goto bExpBigger;
! 2721: if ( aExp == 0x7FF ) {
! 2722: if ( aSig | bSig ) return propagateFloat64NaN( a, b );
! 2723: float_raise( float_flag_invalid );
! 2724: return float64_default_nan;
! 2725: }
! 2726: if ( aExp == 0 ) {
! 2727: aExp = 1;
! 2728: bExp = 1;
! 2729: }
! 2730: if ( bSig < aSig ) goto aBigger;
! 2731: if ( aSig < bSig ) goto bBigger;
! 2732: return packFloat64( float_rounding_mode() == float_round_down, 0, 0 );
! 2733: bExpBigger:
! 2734: if ( bExp == 0x7FF ) {
! 2735: if ( bSig ) return propagateFloat64NaN( a, b );
! 2736: return packFloat64( zSign ^ 1, 0x7FF, 0 );
! 2737: }
! 2738: if ( aExp == 0 ) {
! 2739: ++expDiff;
! 2740: }
! 2741: else {
! 2742: aSig |= LIT64( 0x4000000000000000 );
! 2743: }
! 2744: shift64RightJamming( aSig, - expDiff, &aSig );
! 2745: bSig |= LIT64( 0x4000000000000000 );
! 2746: bBigger:
! 2747: zSig = bSig - aSig;
! 2748: zExp = bExp;
! 2749: zSign ^= 1;
! 2750: goto normalizeRoundAndPack;
! 2751: aExpBigger:
! 2752: if ( aExp == 0x7FF ) {
! 2753: if ( aSig ) return propagateFloat64NaN( a, b );
! 2754: return a;
! 2755: }
! 2756: if ( bExp == 0 ) {
! 2757: --expDiff;
! 2758: }
! 2759: else {
! 2760: bSig |= LIT64( 0x4000000000000000 );
! 2761: }
! 2762: shift64RightJamming( bSig, expDiff, &bSig );
! 2763: aSig |= LIT64( 0x4000000000000000 );
! 2764: aBigger:
! 2765: zSig = aSig - bSig;
! 2766: zExp = aExp;
! 2767: normalizeRoundAndPack:
! 2768: --zExp;
! 2769: return normalizeRoundAndPackFloat64( zSign, zExp, zSig );
! 2770:
! 2771: }
! 2772:
! 2773: /*
! 2774: -------------------------------------------------------------------------------
! 2775: Returns the result of adding the double-precision floating-point values `a'
! 2776: and `b'. The operation is performed according to the IEC/IEEE Standard for
! 2777: Binary Floating-Point Arithmetic.
! 2778: -------------------------------------------------------------------------------
! 2779: */
! 2780: float64 float64_add( float64 a, float64 b )
! 2781: {
! 2782: flag aSign, bSign;
! 2783:
! 2784: aSign = extractFloat64Sign( a );
! 2785: bSign = extractFloat64Sign( b );
! 2786: if ( aSign == bSign ) {
! 2787: return addFloat64Sigs( a, b, aSign );
! 2788: }
! 2789: else {
! 2790: return subFloat64Sigs( a, b, aSign );
! 2791: }
! 2792:
! 2793: }
! 2794:
! 2795: /*
! 2796: -------------------------------------------------------------------------------
! 2797: Returns the result of subtracting the double-precision floating-point values
! 2798: `a' and `b'. The operation is performed according to the IEC/IEEE Standard
! 2799: for Binary Floating-Point Arithmetic.
! 2800: -------------------------------------------------------------------------------
! 2801: */
! 2802: float64 float64_sub( float64 a, float64 b )
! 2803: {
! 2804: flag aSign, bSign;
! 2805:
! 2806: aSign = extractFloat64Sign( a );
! 2807: bSign = extractFloat64Sign( b );
! 2808: if ( aSign == bSign ) {
! 2809: return subFloat64Sigs( a, b, aSign );
! 2810: }
! 2811: else {
! 2812: return addFloat64Sigs( a, b, aSign );
! 2813: }
! 2814:
! 2815: }
! 2816:
! 2817: /*
! 2818: -------------------------------------------------------------------------------
! 2819: Returns the result of multiplying the double-precision floating-point values
! 2820: `a' and `b'. The operation is performed according to the IEC/IEEE Standard
! 2821: for Binary Floating-Point Arithmetic.
! 2822: -------------------------------------------------------------------------------
! 2823: */
! 2824: float64 float64_mul( float64 a, float64 b )
! 2825: {
! 2826: flag aSign, bSign, zSign;
! 2827: int16 aExp, bExp, zExp;
! 2828: bits64 aSig, bSig, zSig0, zSig1;
! 2829:
! 2830: aSig = extractFloat64Frac( a );
! 2831: aExp = extractFloat64Exp( a );
! 2832: aSign = extractFloat64Sign( a );
! 2833: bSig = extractFloat64Frac( b );
! 2834: bExp = extractFloat64Exp( b );
! 2835: bSign = extractFloat64Sign( b );
! 2836: zSign = aSign ^ bSign;
! 2837: if ( aExp == 0x7FF ) {
! 2838: if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
! 2839: return propagateFloat64NaN( a, b );
! 2840: }
! 2841: if ( ( bExp | bSig ) == 0 ) {
! 2842: float_raise( float_flag_invalid );
! 2843: return float64_default_nan;
! 2844: }
! 2845: return packFloat64( zSign, 0x7FF, 0 );
! 2846: }
! 2847: if ( bExp == 0x7FF ) {
! 2848: if ( bSig ) return propagateFloat64NaN( a, b );
! 2849: if ( ( aExp | aSig ) == 0 ) {
! 2850: float_raise( float_flag_invalid );
! 2851: return float64_default_nan;
! 2852: }
! 2853: return packFloat64( zSign, 0x7FF, 0 );
! 2854: }
! 2855: if ( aExp == 0 ) {
! 2856: if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
! 2857: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
! 2858: }
! 2859: if ( bExp == 0 ) {
! 2860: if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );
! 2861: normalizeFloat64Subnormal( bSig, &bExp, &bSig );
! 2862: }
! 2863: zExp = aExp + bExp - 0x3FF;
! 2864: aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
! 2865: bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
! 2866: mul64To128( aSig, bSig, &zSig0, &zSig1 );
! 2867: zSig0 |= ( zSig1 != 0 );
! 2868: if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
! 2869: zSig0 <<= 1;
! 2870: --zExp;
! 2871: }
! 2872: return roundAndPackFloat64( zSign, zExp, zSig0 );
! 2873:
! 2874: }
! 2875:
! 2876: /*
! 2877: -------------------------------------------------------------------------------
! 2878: Returns the result of dividing the double-precision floating-point value `a'
! 2879: by the corresponding value `b'. The operation is performed according to
! 2880: the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 2881: -------------------------------------------------------------------------------
! 2882: */
! 2883: float64 float64_div( float64 a, float64 b )
! 2884: {
! 2885: flag aSign, bSign, zSign;
! 2886: int16 aExp, bExp, zExp;
! 2887: bits64 aSig, bSig, zSig;
! 2888: bits64 rem0, rem1;
! 2889: bits64 term0, term1;
! 2890:
! 2891: aSig = extractFloat64Frac( a );
! 2892: aExp = extractFloat64Exp( a );
! 2893: aSign = extractFloat64Sign( a );
! 2894: bSig = extractFloat64Frac( b );
! 2895: bExp = extractFloat64Exp( b );
! 2896: bSign = extractFloat64Sign( b );
! 2897: zSign = aSign ^ bSign;
! 2898: if ( aExp == 0x7FF ) {
! 2899: if ( aSig ) return propagateFloat64NaN( a, b );
! 2900: if ( bExp == 0x7FF ) {
! 2901: if ( bSig ) return propagateFloat64NaN( a, b );
! 2902: float_raise( float_flag_invalid );
! 2903: return float64_default_nan;
! 2904: }
! 2905: return packFloat64( zSign, 0x7FF, 0 );
! 2906: }
! 2907: if ( bExp == 0x7FF ) {
! 2908: if ( bSig ) return propagateFloat64NaN( a, b );
! 2909: return packFloat64( zSign, 0, 0 );
! 2910: }
! 2911: if ( bExp == 0 ) {
! 2912: if ( bSig == 0 ) {
! 2913: if ( ( aExp | aSig ) == 0 ) {
! 2914: float_raise( float_flag_invalid );
! 2915: return float64_default_nan;
! 2916: }
! 2917: float_raise( float_flag_divbyzero );
! 2918: return packFloat64( zSign, 0x7FF, 0 );
! 2919: }
! 2920: normalizeFloat64Subnormal( bSig, &bExp, &bSig );
! 2921: }
! 2922: if ( aExp == 0 ) {
! 2923: if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
! 2924: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
! 2925: }
! 2926: zExp = aExp - bExp + 0x3FD;
! 2927: aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
! 2928: bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
! 2929: if ( bSig <= ( aSig + aSig ) ) {
! 2930: aSig >>= 1;
! 2931: ++zExp;
! 2932: }
! 2933: zSig = estimateDiv128To64( aSig, 0, bSig );
! 2934: if ( ( zSig & 0x1FF ) <= 2 ) {
! 2935: mul64To128( bSig, zSig, &term0, &term1 );
! 2936: sub128( aSig, 0, term0, term1, &rem0, &rem1 );
! 2937: while ( (sbits64) rem0 < 0 ) {
! 2938: --zSig;
! 2939: add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
! 2940: }
! 2941: zSig |= ( rem1 != 0 );
! 2942: }
! 2943: return roundAndPackFloat64( zSign, zExp, zSig );
! 2944:
! 2945: }
! 2946:
! 2947: #ifndef SOFTFLOAT_FOR_GCC
! 2948: /*
! 2949: -------------------------------------------------------------------------------
! 2950: Returns the remainder of the double-precision floating-point value `a'
! 2951: with respect to the corresponding value `b'. The operation is performed
! 2952: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 2953: -------------------------------------------------------------------------------
! 2954: */
! 2955: float64 float64_rem( float64 a, float64 b )
! 2956: {
! 2957: flag aSign, bSign, zSign;
! 2958: int16 aExp, bExp, expDiff;
! 2959: bits64 aSig, bSig;
! 2960: bits64 q, alternateASig;
! 2961: sbits64 sigMean;
! 2962:
! 2963: aSig = extractFloat64Frac( a );
! 2964: aExp = extractFloat64Exp( a );
! 2965: aSign = extractFloat64Sign( a );
! 2966: bSig = extractFloat64Frac( b );
! 2967: bExp = extractFloat64Exp( b );
! 2968: bSign = extractFloat64Sign( b );
! 2969: if ( aExp == 0x7FF ) {
! 2970: if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
! 2971: return propagateFloat64NaN( a, b );
! 2972: }
! 2973: float_raise( float_flag_invalid );
! 2974: return float64_default_nan;
! 2975: }
! 2976: if ( bExp == 0x7FF ) {
! 2977: if ( bSig ) return propagateFloat64NaN( a, b );
! 2978: return a;
! 2979: }
! 2980: if ( bExp == 0 ) {
! 2981: if ( bSig == 0 ) {
! 2982: float_raise( float_flag_invalid );
! 2983: return float64_default_nan;
! 2984: }
! 2985: normalizeFloat64Subnormal( bSig, &bExp, &bSig );
! 2986: }
! 2987: if ( aExp == 0 ) {
! 2988: if ( aSig == 0 ) return a;
! 2989: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
! 2990: }
! 2991: expDiff = aExp - bExp;
! 2992: aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;
! 2993: bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
! 2994: if ( expDiff < 0 ) {
! 2995: if ( expDiff < -1 ) return a;
! 2996: aSig >>= 1;
! 2997: }
! 2998: q = ( bSig <= aSig );
! 2999: if ( q ) aSig -= bSig;
! 3000: expDiff -= 64;
! 3001: while ( 0 < expDiff ) {
! 3002: q = estimateDiv128To64( aSig, 0, bSig );
! 3003: q = ( 2 < q ) ? q - 2 : 0;
! 3004: aSig = - ( ( bSig>>2 ) * q );
! 3005: expDiff -= 62;
! 3006: }
! 3007: expDiff += 64;
! 3008: if ( 0 < expDiff ) {
! 3009: q = estimateDiv128To64( aSig, 0, bSig );
! 3010: q = ( 2 < q ) ? q - 2 : 0;
! 3011: q >>= 64 - expDiff;
! 3012: bSig >>= 2;
! 3013: aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
! 3014: }
! 3015: else {
! 3016: aSig >>= 2;
! 3017: bSig >>= 2;
! 3018: }
! 3019: do {
! 3020: alternateASig = aSig;
! 3021: ++q;
! 3022: aSig -= bSig;
! 3023: } while ( 0 <= (sbits64) aSig );
! 3024: sigMean = aSig + alternateASig;
! 3025: if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
! 3026: aSig = alternateASig;
! 3027: }
! 3028: zSign = ( (sbits64) aSig < 0 );
! 3029: if ( zSign ) aSig = - aSig;
! 3030: return normalizeRoundAndPackFloat64( aSign ^ zSign, bExp, aSig );
! 3031:
! 3032: }
! 3033:
! 3034: /*
! 3035: -------------------------------------------------------------------------------
! 3036: Returns the square root of the double-precision floating-point value `a'.
! 3037: The operation is performed according to the IEC/IEEE Standard for Binary
! 3038: Floating-Point Arithmetic.
! 3039: -------------------------------------------------------------------------------
! 3040: */
! 3041: float64 float64_sqrt( float64 a )
! 3042: {
! 3043: flag aSign;
! 3044: int16 aExp, zExp;
! 3045: bits64 aSig, zSig, doubleZSig;
! 3046: bits64 rem0, rem1, term0, term1;
! 3047:
! 3048: aSig = extractFloat64Frac( a );
! 3049: aExp = extractFloat64Exp( a );
! 3050: aSign = extractFloat64Sign( a );
! 3051: if ( aExp == 0x7FF ) {
! 3052: if ( aSig ) return propagateFloat64NaN( a, a );
! 3053: if ( ! aSign ) return a;
! 3054: float_raise( float_flag_invalid );
! 3055: return float64_default_nan;
! 3056: }
! 3057: if ( aSign ) {
! 3058: if ( ( aExp | aSig ) == 0 ) return a;
! 3059: float_raise( float_flag_invalid );
! 3060: return float64_default_nan;
! 3061: }
! 3062: if ( aExp == 0 ) {
! 3063: if ( aSig == 0 ) return 0;
! 3064: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
! 3065: }
! 3066: zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
! 3067: aSig |= LIT64( 0x0010000000000000 );
! 3068: zSig = estimateSqrt32( aExp, aSig>>21 );
! 3069: aSig <<= 9 - ( aExp & 1 );
! 3070: zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );
! 3071: if ( ( zSig & 0x1FF ) <= 5 ) {
! 3072: doubleZSig = zSig<<1;
! 3073: mul64To128( zSig, zSig, &term0, &term1 );
! 3074: sub128( aSig, 0, term0, term1, &rem0, &rem1 );
! 3075: while ( (sbits64) rem0 < 0 ) {
! 3076: --zSig;
! 3077: doubleZSig -= 2;
! 3078: add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
! 3079: }
! 3080: zSig |= ( ( rem0 | rem1 ) != 0 );
! 3081: }
! 3082: return roundAndPackFloat64( 0, zExp, zSig );
! 3083:
! 3084: }
! 3085: #endif
! 3086:
! 3087: /*
! 3088: -------------------------------------------------------------------------------
! 3089: Returns 1 if the double-precision floating-point value `a' is equal to the
! 3090: corresponding value `b', and 0 otherwise. The comparison is performed
! 3091: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3092: -------------------------------------------------------------------------------
! 3093: */
! 3094: flag float64_eq( float64 a, float64 b )
! 3095: {
! 3096:
! 3097: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
! 3098: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
! 3099: ) {
! 3100: if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
! 3101: float_raise( float_flag_invalid );
! 3102: }
! 3103: return 0;
! 3104: }
! 3105: return ( a == b ) ||
! 3106: ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) == 0 );
! 3107:
! 3108: }
! 3109:
! 3110: /*
! 3111: -------------------------------------------------------------------------------
! 3112: Returns 1 if the double-precision floating-point value `a' is less than or
! 3113: equal to the corresponding value `b', and 0 otherwise. The comparison is
! 3114: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 3115: Arithmetic.
! 3116: -------------------------------------------------------------------------------
! 3117: */
! 3118: flag float64_le( float64 a, float64 b )
! 3119: {
! 3120: flag aSign, bSign;
! 3121:
! 3122: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
! 3123: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
! 3124: ) {
! 3125: float_raise( float_flag_invalid );
! 3126: return 0;
! 3127: }
! 3128: aSign = extractFloat64Sign( a );
! 3129: bSign = extractFloat64Sign( b );
! 3130: if ( aSign != bSign )
! 3131: return aSign ||
! 3132: ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) ==
! 3133: 0 );
! 3134: return ( a == b ) ||
! 3135: ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) );
! 3136:
! 3137: }
! 3138:
! 3139: /*
! 3140: -------------------------------------------------------------------------------
! 3141: Returns 1 if the double-precision floating-point value `a' is less than
! 3142: the corresponding value `b', and 0 otherwise. The comparison is performed
! 3143: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3144: -------------------------------------------------------------------------------
! 3145: */
! 3146: flag float64_lt( float64 a, float64 b )
! 3147: {
! 3148: flag aSign, bSign;
! 3149:
! 3150: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
! 3151: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
! 3152: ) {
! 3153: float_raise( float_flag_invalid );
! 3154: return 0;
! 3155: }
! 3156: aSign = extractFloat64Sign( a );
! 3157: bSign = extractFloat64Sign( b );
! 3158: if ( aSign != bSign )
! 3159: return aSign &&
! 3160: ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) !=
! 3161: 0 );
! 3162: return ( a != b ) &&
! 3163: ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) );
! 3164:
! 3165: }
! 3166:
! 3167: #ifndef SOFTFLOAT_FOR_GCC
! 3168: /*
! 3169: -------------------------------------------------------------------------------
! 3170: Returns 1 if the double-precision floating-point value `a' is equal to the
! 3171: corresponding value `b', and 0 otherwise. The invalid exception is raised
! 3172: if either operand is a NaN. Otherwise, the comparison is performed
! 3173: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3174: -------------------------------------------------------------------------------
! 3175: */
! 3176: flag float64_eq_signaling( float64 a, float64 b )
! 3177: {
! 3178:
! 3179: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
! 3180: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
! 3181: ) {
! 3182: float_raise( float_flag_invalid );
! 3183: return 0;
! 3184: }
! 3185: return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
! 3186:
! 3187: }
! 3188:
! 3189: /*
! 3190: -------------------------------------------------------------------------------
! 3191: Returns 1 if the double-precision floating-point value `a' is less than or
! 3192: equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
! 3193: cause an exception. Otherwise, the comparison is performed according to the
! 3194: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3195: -------------------------------------------------------------------------------
! 3196: */
! 3197: flag float64_le_quiet( float64 a, float64 b )
! 3198: {
! 3199: flag aSign, bSign;
! 3200:
! 3201: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
! 3202: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
! 3203: ) {
! 3204: if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
! 3205: float_raise( float_flag_invalid );
! 3206: }
! 3207: return 0;
! 3208: }
! 3209: aSign = extractFloat64Sign( a );
! 3210: bSign = extractFloat64Sign( b );
! 3211: if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
! 3212: return ( a == b ) || ( aSign ^ ( a < b ) );
! 3213:
! 3214: }
! 3215:
! 3216: /*
! 3217: -------------------------------------------------------------------------------
! 3218: Returns 1 if the double-precision floating-point value `a' is less than
! 3219: the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
! 3220: exception. Otherwise, the comparison is performed according to the IEC/IEEE
! 3221: Standard for Binary Floating-Point Arithmetic.
! 3222: -------------------------------------------------------------------------------
! 3223: */
! 3224: flag float64_lt_quiet( float64 a, float64 b )
! 3225: {
! 3226: flag aSign, bSign;
! 3227:
! 3228: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
! 3229: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
! 3230: ) {
! 3231: if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
! 3232: float_raise( float_flag_invalid );
! 3233: }
! 3234: return 0;
! 3235: }
! 3236: aSign = extractFloat64Sign( a );
! 3237: bSign = extractFloat64Sign( b );
! 3238: if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
! 3239: return ( a != b ) && ( aSign ^ ( a < b ) );
! 3240:
! 3241: }
! 3242: #endif
! 3243:
! 3244: #ifdef FLOATX80
! 3245:
! 3246: /*
! 3247: -------------------------------------------------------------------------------
! 3248: Returns the result of converting the extended double-precision floating-
! 3249: point value `a' to the 32-bit two's complement integer format. The
! 3250: conversion is performed according to the IEC/IEEE Standard for Binary
! 3251: Floating-Point Arithmetic---which means in particular that the conversion
! 3252: is rounded according to the current rounding mode. If `a' is a NaN, the
! 3253: largest positive integer is returned. Otherwise, if the conversion
! 3254: overflows, the largest integer with the same sign as `a' is returned.
! 3255: -------------------------------------------------------------------------------
! 3256: */
! 3257: int32 floatx80_to_int32( floatx80 a )
! 3258: {
! 3259: flag aSign;
! 3260: int32 aExp, shiftCount;
! 3261: bits64 aSig;
! 3262:
! 3263: aSig = extractFloatx80Frac( a );
! 3264: aExp = extractFloatx80Exp( a );
! 3265: aSign = extractFloatx80Sign( a );
! 3266: if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
! 3267: shiftCount = 0x4037 - aExp;
! 3268: if ( shiftCount <= 0 ) shiftCount = 1;
! 3269: shift64RightJamming( aSig, shiftCount, &aSig );
! 3270: return roundAndPackInt32( aSign, aSig );
! 3271:
! 3272: }
! 3273:
! 3274: /*
! 3275: -------------------------------------------------------------------------------
! 3276: Returns the result of converting the extended double-precision floating-
! 3277: point value `a' to the 32-bit two's complement integer format. The
! 3278: conversion is performed according to the IEC/IEEE Standard for Binary
! 3279: Floating-Point Arithmetic, except that the conversion is always rounded
! 3280: toward zero. If `a' is a NaN, the largest positive integer is returned.
! 3281: Otherwise, if the conversion overflows, the largest integer with the same
! 3282: sign as `a' is returned.
! 3283: -------------------------------------------------------------------------------
! 3284: */
! 3285: int32 floatx80_to_int32_round_to_zero( floatx80 a )
! 3286: {
! 3287: flag aSign;
! 3288: int32 aExp, shiftCount;
! 3289: bits64 aSig, savedASig;
! 3290: int32 z;
! 3291:
! 3292: aSig = extractFloatx80Frac( a );
! 3293: aExp = extractFloatx80Exp( a );
! 3294: aSign = extractFloatx80Sign( a );
! 3295: if ( 0x401E < aExp ) {
! 3296: if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
! 3297: goto invalid;
! 3298: }
! 3299: else if ( aExp < 0x3FFF ) {
! 3300: if ( aExp || aSig ) float_set_inexact();
! 3301: return 0;
! 3302: }
! 3303: shiftCount = 0x403E - aExp;
! 3304: savedASig = aSig;
! 3305: aSig >>= shiftCount;
! 3306: z = aSig;
! 3307: if ( aSign ) z = - z;
! 3308: if ( ( z < 0 ) ^ aSign ) {
! 3309: invalid:
! 3310: float_raise( float_flag_invalid );
! 3311: return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
! 3312: }
! 3313: if ( ( aSig<<shiftCount ) != savedASig ) {
! 3314: float_set_inexact();
! 3315: }
! 3316: return z;
! 3317:
! 3318: }
! 3319:
! 3320: /*
! 3321: -------------------------------------------------------------------------------
! 3322: Returns the result of converting the extended double-precision floating-
! 3323: point value `a' to the 64-bit two's complement integer format. The
! 3324: conversion is performed according to the IEC/IEEE Standard for Binary
! 3325: Floating-Point Arithmetic---which means in particular that the conversion
! 3326: is rounded according to the current rounding mode. If `a' is a NaN,
! 3327: the largest positive integer is returned. Otherwise, if the conversion
! 3328: overflows, the largest integer with the same sign as `a' is returned.
! 3329: -------------------------------------------------------------------------------
! 3330: */
! 3331: int64 floatx80_to_int64( floatx80 a )
! 3332: {
! 3333: flag aSign;
! 3334: int32 aExp, shiftCount;
! 3335: bits64 aSig, aSigExtra;
! 3336:
! 3337: aSig = extractFloatx80Frac( a );
! 3338: aExp = extractFloatx80Exp( a );
! 3339: aSign = extractFloatx80Sign( a );
! 3340: shiftCount = 0x403E - aExp;
! 3341: if ( shiftCount <= 0 ) {
! 3342: if ( shiftCount ) {
! 3343: float_raise( float_flag_invalid );
! 3344: if ( ! aSign
! 3345: || ( ( aExp == 0x7FFF )
! 3346: && ( aSig != LIT64( 0x8000000000000000 ) ) )
! 3347: ) {
! 3348: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 3349: }
! 3350: return (sbits64) LIT64( 0x8000000000000000 );
! 3351: }
! 3352: aSigExtra = 0;
! 3353: }
! 3354: else {
! 3355: shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
! 3356: }
! 3357: return roundAndPackInt64( aSign, aSig, aSigExtra );
! 3358:
! 3359: }
! 3360:
! 3361: /*
! 3362: -------------------------------------------------------------------------------
! 3363: Returns the result of converting the extended double-precision floating-
! 3364: point value `a' to the 64-bit two's complement integer format. The
! 3365: conversion is performed according to the IEC/IEEE Standard for Binary
! 3366: Floating-Point Arithmetic, except that the conversion is always rounded
! 3367: toward zero. If `a' is a NaN, the largest positive integer is returned.
! 3368: Otherwise, if the conversion overflows, the largest integer with the same
! 3369: sign as `a' is returned.
! 3370: -------------------------------------------------------------------------------
! 3371: */
! 3372: int64 floatx80_to_int64_round_to_zero( floatx80 a )
! 3373: {
! 3374: flag aSign;
! 3375: int32 aExp, shiftCount;
! 3376: bits64 aSig;
! 3377: int64 z;
! 3378:
! 3379: aSig = extractFloatx80Frac( a );
! 3380: aExp = extractFloatx80Exp( a );
! 3381: aSign = extractFloatx80Sign( a );
! 3382: shiftCount = aExp - 0x403E;
! 3383: if ( 0 <= shiftCount ) {
! 3384: aSig &= LIT64( 0x7FFFFFFFFFFFFFFF );
! 3385: if ( ( a.high != 0xC03E ) || aSig ) {
! 3386: float_raise( float_flag_invalid );
! 3387: if ( ! aSign || ( ( aExp == 0x7FFF ) && aSig ) ) {
! 3388: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 3389: }
! 3390: }
! 3391: return (sbits64) LIT64( 0x8000000000000000 );
! 3392: }
! 3393: else if ( aExp < 0x3FFF ) {
! 3394: if ( aExp | aSig ) float_set_inexact();
! 3395: return 0;
! 3396: }
! 3397: z = aSig>>( - shiftCount );
! 3398: if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
! 3399: float_set_inexact();
! 3400: }
! 3401: if ( aSign ) z = - z;
! 3402: return z;
! 3403:
! 3404: }
! 3405:
! 3406: /*
! 3407: -------------------------------------------------------------------------------
! 3408: Returns the result of converting the extended double-precision floating-
! 3409: point value `a' to the single-precision floating-point format. The
! 3410: conversion is performed according to the IEC/IEEE Standard for Binary
! 3411: Floating-Point Arithmetic.
! 3412: -------------------------------------------------------------------------------
! 3413: */
! 3414: float32 floatx80_to_float32( floatx80 a )
! 3415: {
! 3416: flag aSign;
! 3417: int32 aExp;
! 3418: bits64 aSig;
! 3419:
! 3420: aSig = extractFloatx80Frac( a );
! 3421: aExp = extractFloatx80Exp( a );
! 3422: aSign = extractFloatx80Sign( a );
! 3423: if ( aExp == 0x7FFF ) {
! 3424: if ( (bits64) ( aSig<<1 ) ) {
! 3425: return commonNaNToFloat32( floatx80ToCommonNaN( a ) );
! 3426: }
! 3427: return packFloat32( aSign, 0xFF, 0 );
! 3428: }
! 3429: shift64RightJamming( aSig, 33, &aSig );
! 3430: if ( aExp || aSig ) aExp -= 0x3F81;
! 3431: return roundAndPackFloat32( aSign, aExp, aSig );
! 3432:
! 3433: }
! 3434:
! 3435: /*
! 3436: -------------------------------------------------------------------------------
! 3437: Returns the result of converting the extended double-precision floating-
! 3438: point value `a' to the double-precision floating-point format. The
! 3439: conversion is performed according to the IEC/IEEE Standard for Binary
! 3440: Floating-Point Arithmetic.
! 3441: -------------------------------------------------------------------------------
! 3442: */
! 3443: float64 floatx80_to_float64( floatx80 a )
! 3444: {
! 3445: flag aSign;
! 3446: int32 aExp;
! 3447: bits64 aSig, zSig;
! 3448:
! 3449: aSig = extractFloatx80Frac( a );
! 3450: aExp = extractFloatx80Exp( a );
! 3451: aSign = extractFloatx80Sign( a );
! 3452: if ( aExp == 0x7FFF ) {
! 3453: if ( (bits64) ( aSig<<1 ) ) {
! 3454: return commonNaNToFloat64( floatx80ToCommonNaN( a ) );
! 3455: }
! 3456: return packFloat64( aSign, 0x7FF, 0 );
! 3457: }
! 3458: shift64RightJamming( aSig, 1, &zSig );
! 3459: if ( aExp || aSig ) aExp -= 0x3C01;
! 3460: return roundAndPackFloat64( aSign, aExp, zSig );
! 3461:
! 3462: }
! 3463:
! 3464: #ifdef FLOAT128
! 3465:
! 3466: /*
! 3467: -------------------------------------------------------------------------------
! 3468: Returns the result of converting the extended double-precision floating-
! 3469: point value `a' to the quadruple-precision floating-point format. The
! 3470: conversion is performed according to the IEC/IEEE Standard for Binary
! 3471: Floating-Point Arithmetic.
! 3472: -------------------------------------------------------------------------------
! 3473: */
! 3474: float128 floatx80_to_float128( floatx80 a )
! 3475: {
! 3476: flag aSign;
! 3477: int16 aExp;
! 3478: bits64 aSig, zSig0, zSig1;
! 3479:
! 3480: aSig = extractFloatx80Frac( a );
! 3481: aExp = extractFloatx80Exp( a );
! 3482: aSign = extractFloatx80Sign( a );
! 3483: if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
! 3484: return commonNaNToFloat128( floatx80ToCommonNaN( a ) );
! 3485: }
! 3486: shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
! 3487: return packFloat128( aSign, aExp, zSig0, zSig1 );
! 3488:
! 3489: }
! 3490:
! 3491: #endif
! 3492:
! 3493: /*
! 3494: -------------------------------------------------------------------------------
! 3495: Rounds the extended double-precision floating-point value `a' to an integer,
! 3496: and returns the result as an extended quadruple-precision floating-point
! 3497: value. The operation is performed according to the IEC/IEEE Standard for
! 3498: Binary Floating-Point Arithmetic.
! 3499: -------------------------------------------------------------------------------
! 3500: */
! 3501: floatx80 floatx80_round_to_int( floatx80 a )
! 3502: {
! 3503: flag aSign;
! 3504: int32 aExp;
! 3505: bits64 lastBitMask, roundBitsMask;
! 3506: int8 roundingMode;
! 3507: floatx80 z;
! 3508:
! 3509: aExp = extractFloatx80Exp( a );
! 3510: if ( 0x403E <= aExp ) {
! 3511: if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
! 3512: return propagateFloatx80NaN( a, a );
! 3513: }
! 3514: return a;
! 3515: }
! 3516: if ( aExp < 0x3FFF ) {
! 3517: if ( ( aExp == 0 )
! 3518: && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
! 3519: return a;
! 3520: }
! 3521: float_set_inexact();
! 3522: aSign = extractFloatx80Sign( a );
! 3523: switch ( float_rounding_mode() ) {
! 3524: case float_round_nearest_even:
! 3525: if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
! 3526: ) {
! 3527: return
! 3528: packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
! 3529: }
! 3530: break;
! 3531: case float_round_down:
! 3532: return
! 3533: aSign ?
! 3534: packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
! 3535: : packFloatx80( 0, 0, 0 );
! 3536: case float_round_up:
! 3537: return
! 3538: aSign ? packFloatx80( 1, 0, 0 )
! 3539: : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
! 3540: }
! 3541: return packFloatx80( aSign, 0, 0 );
! 3542: }
! 3543: lastBitMask = 1;
! 3544: lastBitMask <<= 0x403E - aExp;
! 3545: roundBitsMask = lastBitMask - 1;
! 3546: z = a;
! 3547: roundingMode = float_rounding_mode();
! 3548: if ( roundingMode == float_round_nearest_even ) {
! 3549: z.low += lastBitMask>>1;
! 3550: if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
! 3551: }
! 3552: else if ( roundingMode != float_round_to_zero ) {
! 3553: if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) {
! 3554: z.low += roundBitsMask;
! 3555: }
! 3556: }
! 3557: z.low &= ~ roundBitsMask;
! 3558: if ( z.low == 0 ) {
! 3559: ++z.high;
! 3560: z.low = LIT64( 0x8000000000000000 );
! 3561: }
! 3562: if ( z.low != a.low ) float_set_inexact();
! 3563: return z;
! 3564:
! 3565: }
! 3566:
! 3567: /*
! 3568: -------------------------------------------------------------------------------
! 3569: Returns the result of adding the absolute values of the extended double-
! 3570: precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
! 3571: negated before being returned. `zSign' is ignored if the result is a NaN.
! 3572: The addition is performed according to the IEC/IEEE Standard for Binary
! 3573: Floating-Point Arithmetic.
! 3574: -------------------------------------------------------------------------------
! 3575: */
! 3576: static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign )
! 3577: {
! 3578: int32 aExp, bExp, zExp;
! 3579: bits64 aSig, bSig, zSig0, zSig1;
! 3580: int32 expDiff;
! 3581:
! 3582: aSig = extractFloatx80Frac( a );
! 3583: aExp = extractFloatx80Exp( a );
! 3584: bSig = extractFloatx80Frac( b );
! 3585: bExp = extractFloatx80Exp( b );
! 3586: expDiff = aExp - bExp;
! 3587: if ( 0 < expDiff ) {
! 3588: if ( aExp == 0x7FFF ) {
! 3589: if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3590: return a;
! 3591: }
! 3592: if ( bExp == 0 ) --expDiff;
! 3593: shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
! 3594: zExp = aExp;
! 3595: }
! 3596: else if ( expDiff < 0 ) {
! 3597: if ( bExp == 0x7FFF ) {
! 3598: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3599: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 3600: }
! 3601: if ( aExp == 0 ) ++expDiff;
! 3602: shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
! 3603: zExp = bExp;
! 3604: }
! 3605: else {
! 3606: if ( aExp == 0x7FFF ) {
! 3607: if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
! 3608: return propagateFloatx80NaN( a, b );
! 3609: }
! 3610: return a;
! 3611: }
! 3612: zSig1 = 0;
! 3613: zSig0 = aSig + bSig;
! 3614: if ( aExp == 0 ) {
! 3615: normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
! 3616: goto roundAndPack;
! 3617: }
! 3618: zExp = aExp;
! 3619: goto shiftRight1;
! 3620: }
! 3621: zSig0 = aSig + bSig;
! 3622: if ( (sbits64) zSig0 < 0 ) goto roundAndPack;
! 3623: shiftRight1:
! 3624: shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
! 3625: zSig0 |= LIT64( 0x8000000000000000 );
! 3626: ++zExp;
! 3627: roundAndPack:
! 3628: return
! 3629: roundAndPackFloatx80(
! 3630: floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
! 3631:
! 3632: }
! 3633:
! 3634: /*
! 3635: -------------------------------------------------------------------------------
! 3636: Returns the result of subtracting the absolute values of the extended
! 3637: double-precision floating-point values `a' and `b'. If `zSign' is 1, the
! 3638: difference is negated before being returned. `zSign' is ignored if the
! 3639: result is a NaN. The subtraction is performed according to the IEC/IEEE
! 3640: Standard for Binary Floating-Point Arithmetic.
! 3641: -------------------------------------------------------------------------------
! 3642: */
! 3643: static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign )
! 3644: {
! 3645: int32 aExp, bExp, zExp;
! 3646: bits64 aSig, bSig, zSig0, zSig1;
! 3647: int32 expDiff;
! 3648: floatx80 z;
! 3649:
! 3650: aSig = extractFloatx80Frac( a );
! 3651: aExp = extractFloatx80Exp( a );
! 3652: bSig = extractFloatx80Frac( b );
! 3653: bExp = extractFloatx80Exp( b );
! 3654: expDiff = aExp - bExp;
! 3655: if ( 0 < expDiff ) goto aExpBigger;
! 3656: if ( expDiff < 0 ) goto bExpBigger;
! 3657: if ( aExp == 0x7FFF ) {
! 3658: if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
! 3659: return propagateFloatx80NaN( a, b );
! 3660: }
! 3661: float_raise( float_flag_invalid );
! 3662: z.low = floatx80_default_nan_low;
! 3663: z.high = floatx80_default_nan_high;
! 3664: return z;
! 3665: }
! 3666: if ( aExp == 0 ) {
! 3667: aExp = 1;
! 3668: bExp = 1;
! 3669: }
! 3670: zSig1 = 0;
! 3671: if ( bSig < aSig ) goto aBigger;
! 3672: if ( aSig < bSig ) goto bBigger;
! 3673: return packFloatx80( float_rounding_mode() == float_round_down, 0, 0 );
! 3674: bExpBigger:
! 3675: if ( bExp == 0x7FFF ) {
! 3676: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3677: return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 3678: }
! 3679: if ( aExp == 0 ) ++expDiff;
! 3680: shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
! 3681: bBigger:
! 3682: sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
! 3683: zExp = bExp;
! 3684: zSign ^= 1;
! 3685: goto normalizeRoundAndPack;
! 3686: aExpBigger:
! 3687: if ( aExp == 0x7FFF ) {
! 3688: if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3689: return a;
! 3690: }
! 3691: if ( bExp == 0 ) --expDiff;
! 3692: shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
! 3693: aBigger:
! 3694: sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
! 3695: zExp = aExp;
! 3696: normalizeRoundAndPack:
! 3697: return
! 3698: normalizeRoundAndPackFloatx80(
! 3699: floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
! 3700:
! 3701: }
! 3702:
! 3703: /*
! 3704: -------------------------------------------------------------------------------
! 3705: Returns the result of adding the extended double-precision floating-point
! 3706: values `a' and `b'. The operation is performed according to the IEC/IEEE
! 3707: Standard for Binary Floating-Point Arithmetic.
! 3708: -------------------------------------------------------------------------------
! 3709: */
! 3710: floatx80 floatx80_add( floatx80 a, floatx80 b )
! 3711: {
! 3712: flag aSign, bSign;
! 3713:
! 3714: aSign = extractFloatx80Sign( a );
! 3715: bSign = extractFloatx80Sign( b );
! 3716: if ( aSign == bSign ) {
! 3717: return addFloatx80Sigs( a, b, aSign );
! 3718: }
! 3719: else {
! 3720: return subFloatx80Sigs( a, b, aSign );
! 3721: }
! 3722:
! 3723: }
! 3724:
! 3725: /*
! 3726: -------------------------------------------------------------------------------
! 3727: Returns the result of subtracting the extended double-precision floating-
! 3728: point values `a' and `b'. The operation is performed according to the
! 3729: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3730: -------------------------------------------------------------------------------
! 3731: */
! 3732: floatx80 floatx80_sub( floatx80 a, floatx80 b )
! 3733: {
! 3734: flag aSign, bSign;
! 3735:
! 3736: aSign = extractFloatx80Sign( a );
! 3737: bSign = extractFloatx80Sign( b );
! 3738: if ( aSign == bSign ) {
! 3739: return subFloatx80Sigs( a, b, aSign );
! 3740: }
! 3741: else {
! 3742: return addFloatx80Sigs( a, b, aSign );
! 3743: }
! 3744:
! 3745: }
! 3746:
! 3747: /*
! 3748: -------------------------------------------------------------------------------
! 3749: Returns the result of multiplying the extended double-precision floating-
! 3750: point values `a' and `b'. The operation is performed according to the
! 3751: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3752: -------------------------------------------------------------------------------
! 3753: */
! 3754: floatx80 floatx80_mul( floatx80 a, floatx80 b )
! 3755: {
! 3756: flag aSign, bSign, zSign;
! 3757: int32 aExp, bExp, zExp;
! 3758: bits64 aSig, bSig, zSig0, zSig1;
! 3759: floatx80 z;
! 3760:
! 3761: aSig = extractFloatx80Frac( a );
! 3762: aExp = extractFloatx80Exp( a );
! 3763: aSign = extractFloatx80Sign( a );
! 3764: bSig = extractFloatx80Frac( b );
! 3765: bExp = extractFloatx80Exp( b );
! 3766: bSign = extractFloatx80Sign( b );
! 3767: zSign = aSign ^ bSign;
! 3768: if ( aExp == 0x7FFF ) {
! 3769: if ( (bits64) ( aSig<<1 )
! 3770: || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
! 3771: return propagateFloatx80NaN( a, b );
! 3772: }
! 3773: if ( ( bExp | bSig ) == 0 ) goto invalid;
! 3774: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 3775: }
! 3776: if ( bExp == 0x7FFF ) {
! 3777: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3778: if ( ( aExp | aSig ) == 0 ) {
! 3779: invalid:
! 3780: float_raise( float_flag_invalid );
! 3781: z.low = floatx80_default_nan_low;
! 3782: z.high = floatx80_default_nan_high;
! 3783: return z;
! 3784: }
! 3785: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 3786: }
! 3787: if ( aExp == 0 ) {
! 3788: if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
! 3789: normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
! 3790: }
! 3791: if ( bExp == 0 ) {
! 3792: if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
! 3793: normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
! 3794: }
! 3795: zExp = aExp + bExp - 0x3FFE;
! 3796: mul64To128( aSig, bSig, &zSig0, &zSig1 );
! 3797: if ( 0 < (sbits64) zSig0 ) {
! 3798: shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
! 3799: --zExp;
! 3800: }
! 3801: return
! 3802: roundAndPackFloatx80(
! 3803: floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
! 3804:
! 3805: }
! 3806:
! 3807: /*
! 3808: -------------------------------------------------------------------------------
! 3809: Returns the result of dividing the extended double-precision floating-point
! 3810: value `a' by the corresponding value `b'. The operation is performed
! 3811: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3812: -------------------------------------------------------------------------------
! 3813: */
! 3814: floatx80 floatx80_div( floatx80 a, floatx80 b )
! 3815: {
! 3816: flag aSign, bSign, zSign;
! 3817: int32 aExp, bExp, zExp;
! 3818: bits64 aSig, bSig, zSig0, zSig1;
! 3819: bits64 rem0, rem1, rem2, term0, term1, term2;
! 3820: floatx80 z;
! 3821:
! 3822: aSig = extractFloatx80Frac( a );
! 3823: aExp = extractFloatx80Exp( a );
! 3824: aSign = extractFloatx80Sign( a );
! 3825: bSig = extractFloatx80Frac( b );
! 3826: bExp = extractFloatx80Exp( b );
! 3827: bSign = extractFloatx80Sign( b );
! 3828: zSign = aSign ^ bSign;
! 3829: if ( aExp == 0x7FFF ) {
! 3830: if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3831: if ( bExp == 0x7FFF ) {
! 3832: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3833: goto invalid;
! 3834: }
! 3835: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 3836: }
! 3837: if ( bExp == 0x7FFF ) {
! 3838: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3839: return packFloatx80( zSign, 0, 0 );
! 3840: }
! 3841: if ( bExp == 0 ) {
! 3842: if ( bSig == 0 ) {
! 3843: if ( ( aExp | aSig ) == 0 ) {
! 3844: invalid:
! 3845: float_raise( float_flag_invalid );
! 3846: z.low = floatx80_default_nan_low;
! 3847: z.high = floatx80_default_nan_high;
! 3848: return z;
! 3849: }
! 3850: float_raise( float_flag_divbyzero );
! 3851: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 3852: }
! 3853: normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
! 3854: }
! 3855: if ( aExp == 0 ) {
! 3856: if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
! 3857: normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
! 3858: }
! 3859: zExp = aExp - bExp + 0x3FFE;
! 3860: rem1 = 0;
! 3861: if ( bSig <= aSig ) {
! 3862: shift128Right( aSig, 0, 1, &aSig, &rem1 );
! 3863: ++zExp;
! 3864: }
! 3865: zSig0 = estimateDiv128To64( aSig, rem1, bSig );
! 3866: mul64To128( bSig, zSig0, &term0, &term1 );
! 3867: sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
! 3868: while ( (sbits64) rem0 < 0 ) {
! 3869: --zSig0;
! 3870: add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
! 3871: }
! 3872: zSig1 = estimateDiv128To64( rem1, 0, bSig );
! 3873: if ( (bits64) ( zSig1<<1 ) <= 8 ) {
! 3874: mul64To128( bSig, zSig1, &term1, &term2 );
! 3875: sub128( rem1, 0, term1, term2, &rem1, &rem2 );
! 3876: while ( (sbits64) rem1 < 0 ) {
! 3877: --zSig1;
! 3878: add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
! 3879: }
! 3880: zSig1 |= ( ( rem1 | rem2 ) != 0 );
! 3881: }
! 3882: return
! 3883: roundAndPackFloatx80(
! 3884: floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
! 3885:
! 3886: }
! 3887:
! 3888: /*
! 3889: -------------------------------------------------------------------------------
! 3890: Returns the remainder of the extended double-precision floating-point value
! 3891: `a' with respect to the corresponding value `b'. The operation is performed
! 3892: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 3893: -------------------------------------------------------------------------------
! 3894: */
! 3895: floatx80 floatx80_rem( floatx80 a, floatx80 b )
! 3896: {
! 3897: flag aSign, bSign, zSign;
! 3898: int32 aExp, bExp, expDiff;
! 3899: bits64 aSig0, aSig1, bSig;
! 3900: bits64 q, term0, term1, alternateASig0, alternateASig1;
! 3901: floatx80 z;
! 3902:
! 3903: aSig0 = extractFloatx80Frac( a );
! 3904: aExp = extractFloatx80Exp( a );
! 3905: aSign = extractFloatx80Sign( a );
! 3906: bSig = extractFloatx80Frac( b );
! 3907: bExp = extractFloatx80Exp( b );
! 3908: bSign = extractFloatx80Sign( b );
! 3909: if ( aExp == 0x7FFF ) {
! 3910: if ( (bits64) ( aSig0<<1 )
! 3911: || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
! 3912: return propagateFloatx80NaN( a, b );
! 3913: }
! 3914: goto invalid;
! 3915: }
! 3916: if ( bExp == 0x7FFF ) {
! 3917: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
! 3918: return a;
! 3919: }
! 3920: if ( bExp == 0 ) {
! 3921: if ( bSig == 0 ) {
! 3922: invalid:
! 3923: float_raise( float_flag_invalid );
! 3924: z.low = floatx80_default_nan_low;
! 3925: z.high = floatx80_default_nan_high;
! 3926: return z;
! 3927: }
! 3928: normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
! 3929: }
! 3930: if ( aExp == 0 ) {
! 3931: if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
! 3932: normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
! 3933: }
! 3934: bSig |= LIT64( 0x8000000000000000 );
! 3935: zSign = aSign;
! 3936: expDiff = aExp - bExp;
! 3937: aSig1 = 0;
! 3938: if ( expDiff < 0 ) {
! 3939: if ( expDiff < -1 ) return a;
! 3940: shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
! 3941: expDiff = 0;
! 3942: }
! 3943: q = ( bSig <= aSig0 );
! 3944: if ( q ) aSig0 -= bSig;
! 3945: expDiff -= 64;
! 3946: while ( 0 < expDiff ) {
! 3947: q = estimateDiv128To64( aSig0, aSig1, bSig );
! 3948: q = ( 2 < q ) ? q - 2 : 0;
! 3949: mul64To128( bSig, q, &term0, &term1 );
! 3950: sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
! 3951: shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
! 3952: expDiff -= 62;
! 3953: }
! 3954: expDiff += 64;
! 3955: if ( 0 < expDiff ) {
! 3956: q = estimateDiv128To64( aSig0, aSig1, bSig );
! 3957: q = ( 2 < q ) ? q - 2 : 0;
! 3958: q >>= 64 - expDiff;
! 3959: mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
! 3960: sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
! 3961: shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
! 3962: while ( le128( term0, term1, aSig0, aSig1 ) ) {
! 3963: ++q;
! 3964: sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
! 3965: }
! 3966: }
! 3967: else {
! 3968: term1 = 0;
! 3969: term0 = bSig;
! 3970: }
! 3971: sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
! 3972: if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
! 3973: || ( eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
! 3974: && ( q & 1 ) )
! 3975: ) {
! 3976: aSig0 = alternateASig0;
! 3977: aSig1 = alternateASig1;
! 3978: zSign = ! zSign;
! 3979: }
! 3980: return
! 3981: normalizeRoundAndPackFloatx80(
! 3982: 80, zSign, bExp + expDiff, aSig0, aSig1 );
! 3983:
! 3984: }
! 3985:
! 3986: /*
! 3987: -------------------------------------------------------------------------------
! 3988: Returns the square root of the extended double-precision floating-point
! 3989: value `a'. The operation is performed according to the IEC/IEEE Standard
! 3990: for Binary Floating-Point Arithmetic.
! 3991: -------------------------------------------------------------------------------
! 3992: */
! 3993: floatx80 floatx80_sqrt( floatx80 a )
! 3994: {
! 3995: flag aSign;
! 3996: int32 aExp, zExp;
! 3997: bits64 aSig0, aSig1, zSig0, zSig1, doubleZSig0;
! 3998: bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
! 3999: floatx80 z;
! 4000:
! 4001: aSig0 = extractFloatx80Frac( a );
! 4002: aExp = extractFloatx80Exp( a );
! 4003: aSign = extractFloatx80Sign( a );
! 4004: if ( aExp == 0x7FFF ) {
! 4005: if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a );
! 4006: if ( ! aSign ) return a;
! 4007: goto invalid;
! 4008: }
! 4009: if ( aSign ) {
! 4010: if ( ( aExp | aSig0 ) == 0 ) return a;
! 4011: invalid:
! 4012: float_raise( float_flag_invalid );
! 4013: z.low = floatx80_default_nan_low;
! 4014: z.high = floatx80_default_nan_high;
! 4015: return z;
! 4016: }
! 4017: if ( aExp == 0 ) {
! 4018: if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
! 4019: normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
! 4020: }
! 4021: zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
! 4022: zSig0 = estimateSqrt32( aExp, aSig0>>32 );
! 4023: shift128Right( aSig0, 0, 2 + ( aExp & 1 ), &aSig0, &aSig1 );
! 4024: zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
! 4025: doubleZSig0 = zSig0<<1;
! 4026: mul64To128( zSig0, zSig0, &term0, &term1 );
! 4027: sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
! 4028: while ( (sbits64) rem0 < 0 ) {
! 4029: --zSig0;
! 4030: doubleZSig0 -= 2;
! 4031: add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
! 4032: }
! 4033: zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
! 4034: if ( ( zSig1 & LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
! 4035: if ( zSig1 == 0 ) zSig1 = 1;
! 4036: mul64To128( doubleZSig0, zSig1, &term1, &term2 );
! 4037: sub128( rem1, 0, term1, term2, &rem1, &rem2 );
! 4038: mul64To128( zSig1, zSig1, &term2, &term3 );
! 4039: sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
! 4040: while ( (sbits64) rem1 < 0 ) {
! 4041: --zSig1;
! 4042: shortShift128Left( 0, zSig1, 1, &term2, &term3 );
! 4043: term3 |= 1;
! 4044: term2 |= doubleZSig0;
! 4045: add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
! 4046: }
! 4047: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
! 4048: }
! 4049: shortShift128Left( 0, zSig1, 1, &zSig0, &zSig1 );
! 4050: zSig0 |= doubleZSig0;
! 4051: return
! 4052: roundAndPackFloatx80(
! 4053: floatx80_rounding_precision, 0, zExp, zSig0, zSig1 );
! 4054:
! 4055: }
! 4056:
! 4057: /*
! 4058: -------------------------------------------------------------------------------
! 4059: Returns 1 if the extended double-precision floating-point value `a' is
! 4060: equal to the corresponding value `b', and 0 otherwise. The comparison is
! 4061: performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4062: Arithmetic.
! 4063: -------------------------------------------------------------------------------
! 4064: */
! 4065: flag floatx80_eq( floatx80 a, floatx80 b )
! 4066: {
! 4067:
! 4068: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
! 4069: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
! 4070: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
! 4071: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
! 4072: ) {
! 4073: if ( floatx80_is_signaling_nan( a )
! 4074: || floatx80_is_signaling_nan( b ) ) {
! 4075: float_raise( float_flag_invalid );
! 4076: }
! 4077: return 0;
! 4078: }
! 4079: return
! 4080: ( a.low == b.low )
! 4081: && ( ( a.high == b.high )
! 4082: || ( ( a.low == 0 )
! 4083: && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
! 4084: );
! 4085:
! 4086: }
! 4087:
! 4088: /*
! 4089: -------------------------------------------------------------------------------
! 4090: Returns 1 if the extended double-precision floating-point value `a' is
! 4091: less than or equal to the corresponding value `b', and 0 otherwise. The
! 4092: comparison is performed according to the IEC/IEEE Standard for Binary
! 4093: Floating-Point Arithmetic.
! 4094: -------------------------------------------------------------------------------
! 4095: */
! 4096: flag floatx80_le( floatx80 a, floatx80 b )
! 4097: {
! 4098: flag aSign, bSign;
! 4099:
! 4100: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
! 4101: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
! 4102: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
! 4103: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
! 4104: ) {
! 4105: float_raise( float_flag_invalid );
! 4106: return 0;
! 4107: }
! 4108: aSign = extractFloatx80Sign( a );
! 4109: bSign = extractFloatx80Sign( b );
! 4110: if ( aSign != bSign ) {
! 4111: return
! 4112: aSign
! 4113: || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 4114: == 0 );
! 4115: }
! 4116: return
! 4117: aSign ? le128( b.high, b.low, a.high, a.low )
! 4118: : le128( a.high, a.low, b.high, b.low );
! 4119:
! 4120: }
! 4121:
! 4122: /*
! 4123: -------------------------------------------------------------------------------
! 4124: Returns 1 if the extended double-precision floating-point value `a' is
! 4125: less than the corresponding value `b', and 0 otherwise. The comparison
! 4126: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4127: Arithmetic.
! 4128: -------------------------------------------------------------------------------
! 4129: */
! 4130: flag floatx80_lt( floatx80 a, floatx80 b )
! 4131: {
! 4132: flag aSign, bSign;
! 4133:
! 4134: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
! 4135: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
! 4136: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
! 4137: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
! 4138: ) {
! 4139: float_raise( float_flag_invalid );
! 4140: return 0;
! 4141: }
! 4142: aSign = extractFloatx80Sign( a );
! 4143: bSign = extractFloatx80Sign( b );
! 4144: if ( aSign != bSign ) {
! 4145: return
! 4146: aSign
! 4147: && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 4148: != 0 );
! 4149: }
! 4150: return
! 4151: aSign ? lt128( b.high, b.low, a.high, a.low )
! 4152: : lt128( a.high, a.low, b.high, b.low );
! 4153:
! 4154: }
! 4155:
! 4156: /*
! 4157: -------------------------------------------------------------------------------
! 4158: Returns 1 if the extended double-precision floating-point value `a' is equal
! 4159: to the corresponding value `b', and 0 otherwise. The invalid exception is
! 4160: raised if either operand is a NaN. Otherwise, the comparison is performed
! 4161: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 4162: -------------------------------------------------------------------------------
! 4163: */
! 4164: flag floatx80_eq_signaling( floatx80 a, floatx80 b )
! 4165: {
! 4166:
! 4167: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
! 4168: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
! 4169: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
! 4170: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
! 4171: ) {
! 4172: float_raise( float_flag_invalid );
! 4173: return 0;
! 4174: }
! 4175: return
! 4176: ( a.low == b.low )
! 4177: && ( ( a.high == b.high )
! 4178: || ( ( a.low == 0 )
! 4179: && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
! 4180: );
! 4181:
! 4182: }
! 4183:
! 4184: /*
! 4185: -------------------------------------------------------------------------------
! 4186: Returns 1 if the extended double-precision floating-point value `a' is less
! 4187: than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs
! 4188: do not cause an exception. Otherwise, the comparison is performed according
! 4189: to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 4190: -------------------------------------------------------------------------------
! 4191: */
! 4192: flag floatx80_le_quiet( floatx80 a, floatx80 b )
! 4193: {
! 4194: flag aSign, bSign;
! 4195:
! 4196: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
! 4197: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
! 4198: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
! 4199: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
! 4200: ) {
! 4201: if ( floatx80_is_signaling_nan( a )
! 4202: || floatx80_is_signaling_nan( b ) ) {
! 4203: float_raise( float_flag_invalid );
! 4204: }
! 4205: return 0;
! 4206: }
! 4207: aSign = extractFloatx80Sign( a );
! 4208: bSign = extractFloatx80Sign( b );
! 4209: if ( aSign != bSign ) {
! 4210: return
! 4211: aSign
! 4212: || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 4213: == 0 );
! 4214: }
! 4215: return
! 4216: aSign ? le128( b.high, b.low, a.high, a.low )
! 4217: : le128( a.high, a.low, b.high, b.low );
! 4218:
! 4219: }
! 4220:
! 4221: /*
! 4222: -------------------------------------------------------------------------------
! 4223: Returns 1 if the extended double-precision floating-point value `a' is less
! 4224: than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause
! 4225: an exception. Otherwise, the comparison is performed according to the
! 4226: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 4227: -------------------------------------------------------------------------------
! 4228: */
! 4229: flag floatx80_lt_quiet( floatx80 a, floatx80 b )
! 4230: {
! 4231: flag aSign, bSign;
! 4232:
! 4233: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
! 4234: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
! 4235: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
! 4236: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
! 4237: ) {
! 4238: if ( floatx80_is_signaling_nan( a )
! 4239: || floatx80_is_signaling_nan( b ) ) {
! 4240: float_raise( float_flag_invalid );
! 4241: }
! 4242: return 0;
! 4243: }
! 4244: aSign = extractFloatx80Sign( a );
! 4245: bSign = extractFloatx80Sign( b );
! 4246: if ( aSign != bSign ) {
! 4247: return
! 4248: aSign
! 4249: && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 4250: != 0 );
! 4251: }
! 4252: return
! 4253: aSign ? lt128( b.high, b.low, a.high, a.low )
! 4254: : lt128( a.high, a.low, b.high, b.low );
! 4255:
! 4256: }
! 4257:
! 4258: #endif
! 4259:
! 4260: #ifdef FLOAT128
! 4261:
! 4262: /*
! 4263: -------------------------------------------------------------------------------
! 4264: Returns the result of converting the quadruple-precision floating-point
! 4265: value `a' to the 32-bit two's complement integer format. The conversion
! 4266: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4267: Arithmetic---which means in particular that the conversion is rounded
! 4268: according to the current rounding mode. If `a' is a NaN, the largest
! 4269: positive integer is returned. Otherwise, if the conversion overflows, the
! 4270: largest integer with the same sign as `a' is returned.
! 4271: -------------------------------------------------------------------------------
! 4272: */
! 4273: int32 float128_to_int32( float128 a )
! 4274: {
! 4275: flag aSign;
! 4276: int32 aExp, shiftCount;
! 4277: bits64 aSig0, aSig1;
! 4278:
! 4279: aSig1 = extractFloat128Frac1( a );
! 4280: aSig0 = extractFloat128Frac0( a );
! 4281: aExp = extractFloat128Exp( a );
! 4282: aSign = extractFloat128Sign( a );
! 4283: if ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) aSign = 0;
! 4284: if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
! 4285: aSig0 |= ( aSig1 != 0 );
! 4286: shiftCount = 0x4028 - aExp;
! 4287: if ( 0 < shiftCount ) shift64RightJamming( aSig0, shiftCount, &aSig0 );
! 4288: return roundAndPackInt32( aSign, aSig0 );
! 4289:
! 4290: }
! 4291:
! 4292: /*
! 4293: -------------------------------------------------------------------------------
! 4294: Returns the result of converting the quadruple-precision floating-point
! 4295: value `a' to the 32-bit two's complement integer format. The conversion
! 4296: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4297: Arithmetic, except that the conversion is always rounded toward zero. If
! 4298: `a' is a NaN, the largest positive integer is returned. Otherwise, if the
! 4299: conversion overflows, the largest integer with the same sign as `a' is
! 4300: returned.
! 4301: -------------------------------------------------------------------------------
! 4302: */
! 4303: int32 float128_to_int32_round_to_zero( float128 a )
! 4304: {
! 4305: flag aSign;
! 4306: int32 aExp, shiftCount;
! 4307: bits64 aSig0, aSig1, savedASig;
! 4308: int32 z;
! 4309:
! 4310: aSig1 = extractFloat128Frac1( a );
! 4311: aSig0 = extractFloat128Frac0( a );
! 4312: aExp = extractFloat128Exp( a );
! 4313: aSign = extractFloat128Sign( a );
! 4314: aSig0 |= ( aSig1 != 0 );
! 4315: if ( 0x401E < aExp ) {
! 4316: if ( ( aExp == 0x7FFF ) && aSig0 ) aSign = 0;
! 4317: goto invalid;
! 4318: }
! 4319: else if ( aExp < 0x3FFF ) {
! 4320: if ( aExp || aSig0 ) float_set_inexact();
! 4321: return 0;
! 4322: }
! 4323: aSig0 |= LIT64( 0x0001000000000000 );
! 4324: shiftCount = 0x402F - aExp;
! 4325: savedASig = aSig0;
! 4326: aSig0 >>= shiftCount;
! 4327: z = aSig0;
! 4328: if ( aSign ) z = - z;
! 4329: if ( ( z < 0 ) ^ aSign ) {
! 4330: invalid:
! 4331: float_raise( float_flag_invalid );
! 4332: return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
! 4333: }
! 4334: if ( ( aSig0<<shiftCount ) != savedASig ) {
! 4335: float_set_inexact();
! 4336: }
! 4337: return z;
! 4338:
! 4339: }
! 4340:
! 4341: /*
! 4342: -------------------------------------------------------------------------------
! 4343: Returns the result of converting the quadruple-precision floating-point
! 4344: value `a' to the 64-bit two's complement integer format. The conversion
! 4345: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4346: Arithmetic---which means in particular that the conversion is rounded
! 4347: according to the current rounding mode. If `a' is a NaN, the largest
! 4348: positive integer is returned. Otherwise, if the conversion overflows, the
! 4349: largest integer with the same sign as `a' is returned.
! 4350: -------------------------------------------------------------------------------
! 4351: */
! 4352: int64 float128_to_int64( float128 a )
! 4353: {
! 4354: flag aSign;
! 4355: int32 aExp, shiftCount;
! 4356: bits64 aSig0, aSig1;
! 4357:
! 4358: aSig1 = extractFloat128Frac1( a );
! 4359: aSig0 = extractFloat128Frac0( a );
! 4360: aExp = extractFloat128Exp( a );
! 4361: aSign = extractFloat128Sign( a );
! 4362: if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
! 4363: shiftCount = 0x402F - aExp;
! 4364: if ( shiftCount <= 0 ) {
! 4365: if ( 0x403E < aExp ) {
! 4366: float_raise( float_flag_invalid );
! 4367: if ( ! aSign
! 4368: || ( ( aExp == 0x7FFF )
! 4369: && ( aSig1 || ( aSig0 != LIT64( 0x0001000000000000 ) ) )
! 4370: )
! 4371: ) {
! 4372: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 4373: }
! 4374: return (sbits64) LIT64( 0x8000000000000000 );
! 4375: }
! 4376: shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
! 4377: }
! 4378: else {
! 4379: shift64ExtraRightJamming( aSig0, aSig1, shiftCount, &aSig0, &aSig1 );
! 4380: }
! 4381: return roundAndPackInt64( aSign, aSig0, aSig1 );
! 4382:
! 4383: }
! 4384:
! 4385: /*
! 4386: -------------------------------------------------------------------------------
! 4387: Returns the result of converting the quadruple-precision floating-point
! 4388: value `a' to the 64-bit two's complement integer format. The conversion
! 4389: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4390: Arithmetic, except that the conversion is always rounded toward zero.
! 4391: If `a' is a NaN, the largest positive integer is returned. Otherwise, if
! 4392: the conversion overflows, the largest integer with the same sign as `a' is
! 4393: returned.
! 4394: -------------------------------------------------------------------------------
! 4395: */
! 4396: int64 float128_to_int64_round_to_zero( float128 a )
! 4397: {
! 4398: flag aSign;
! 4399: int32 aExp, shiftCount;
! 4400: bits64 aSig0, aSig1;
! 4401: int64 z;
! 4402:
! 4403: aSig1 = extractFloat128Frac1( a );
! 4404: aSig0 = extractFloat128Frac0( a );
! 4405: aExp = extractFloat128Exp( a );
! 4406: aSign = extractFloat128Sign( a );
! 4407: if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
! 4408: shiftCount = aExp - 0x402F;
! 4409: if ( 0 < shiftCount ) {
! 4410: if ( 0x403E <= aExp ) {
! 4411: aSig0 &= LIT64( 0x0000FFFFFFFFFFFF );
! 4412: if ( ( a.high == LIT64( 0xC03E000000000000 ) )
! 4413: && ( aSig1 < LIT64( 0x0002000000000000 ) ) ) {
! 4414: if ( aSig1 ) float_set_inexact();
! 4415: }
! 4416: else {
! 4417: float_raise( float_flag_invalid );
! 4418: if ( ! aSign || ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) ) {
! 4419: return LIT64( 0x7FFFFFFFFFFFFFFF );
! 4420: }
! 4421: }
! 4422: return (sbits64) LIT64( 0x8000000000000000 );
! 4423: }
! 4424: z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
! 4425: if ( (bits64) ( aSig1<<shiftCount ) ) {
! 4426: float_set_inexact();
! 4427: }
! 4428: }
! 4429: else {
! 4430: if ( aExp < 0x3FFF ) {
! 4431: if ( aExp | aSig0 | aSig1 ) {
! 4432: float_set_inexact();
! 4433: }
! 4434: return 0;
! 4435: }
! 4436: z = aSig0>>( - shiftCount );
! 4437: if ( aSig1
! 4438: || ( shiftCount && (bits64) ( aSig0<<( shiftCount & 63 ) ) ) ) {
! 4439: float_set_inexact();
! 4440: }
! 4441: }
! 4442: if ( aSign ) z = - z;
! 4443: return z;
! 4444:
! 4445: }
! 4446:
! 4447: /*
! 4448: -------------------------------------------------------------------------------
! 4449: Returns the result of converting the quadruple-precision floating-point
! 4450: value `a' to the single-precision floating-point format. The conversion
! 4451: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4452: Arithmetic.
! 4453: -------------------------------------------------------------------------------
! 4454: */
! 4455: float32 float128_to_float32( float128 a )
! 4456: {
! 4457: flag aSign;
! 4458: int32 aExp;
! 4459: bits64 aSig0, aSig1;
! 4460: bits32 zSig;
! 4461:
! 4462: aSig1 = extractFloat128Frac1( a );
! 4463: aSig0 = extractFloat128Frac0( a );
! 4464: aExp = extractFloat128Exp( a );
! 4465: aSign = extractFloat128Sign( a );
! 4466: if ( aExp == 0x7FFF ) {
! 4467: if ( aSig0 | aSig1 ) {
! 4468: return commonNaNToFloat32( float128ToCommonNaN( a ) );
! 4469: }
! 4470: return packFloat32( aSign, 0xFF, 0 );
! 4471: }
! 4472: aSig0 |= ( aSig1 != 0 );
! 4473: shift64RightJamming( aSig0, 18, &aSig0 );
! 4474: zSig = aSig0;
! 4475: if ( aExp || zSig ) {
! 4476: zSig |= 0x40000000;
! 4477: aExp -= 0x3F81;
! 4478: }
! 4479: return roundAndPackFloat32( aSign, aExp, zSig );
! 4480:
! 4481: }
! 4482:
! 4483: /*
! 4484: -------------------------------------------------------------------------------
! 4485: Returns the result of converting the quadruple-precision floating-point
! 4486: value `a' to the double-precision floating-point format. The conversion
! 4487: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 4488: Arithmetic.
! 4489: -------------------------------------------------------------------------------
! 4490: */
! 4491: float64 float128_to_float64( float128 a )
! 4492: {
! 4493: flag aSign;
! 4494: int32 aExp;
! 4495: bits64 aSig0, aSig1;
! 4496:
! 4497: aSig1 = extractFloat128Frac1( a );
! 4498: aSig0 = extractFloat128Frac0( a );
! 4499: aExp = extractFloat128Exp( a );
! 4500: aSign = extractFloat128Sign( a );
! 4501: if ( aExp == 0x7FFF ) {
! 4502: if ( aSig0 | aSig1 ) {
! 4503: return commonNaNToFloat64( float128ToCommonNaN( a ) );
! 4504: }
! 4505: return packFloat64( aSign, 0x7FF, 0 );
! 4506: }
! 4507: shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
! 4508: aSig0 |= ( aSig1 != 0 );
! 4509: if ( aExp || aSig0 ) {
! 4510: aSig0 |= LIT64( 0x4000000000000000 );
! 4511: aExp -= 0x3C01;
! 4512: }
! 4513: return roundAndPackFloat64( aSign, aExp, aSig0 );
! 4514:
! 4515: }
! 4516:
! 4517: #ifdef FLOATX80
! 4518:
! 4519: /*
! 4520: -------------------------------------------------------------------------------
! 4521: Returns the result of converting the quadruple-precision floating-point
! 4522: value `a' to the extended double-precision floating-point format. The
! 4523: conversion is performed according to the IEC/IEEE Standard for Binary
! 4524: Floating-Point Arithmetic.
! 4525: -------------------------------------------------------------------------------
! 4526: */
! 4527: floatx80 float128_to_floatx80( float128 a )
! 4528: {
! 4529: flag aSign;
! 4530: int32 aExp;
! 4531: bits64 aSig0, aSig1;
! 4532:
! 4533: aSig1 = extractFloat128Frac1( a );
! 4534: aSig0 = extractFloat128Frac0( a );
! 4535: aExp = extractFloat128Exp( a );
! 4536: aSign = extractFloat128Sign( a );
! 4537: if ( aExp == 0x7FFF ) {
! 4538: if ( aSig0 | aSig1 ) {
! 4539: return commonNaNToFloatx80( float128ToCommonNaN( a ) );
! 4540: }
! 4541: return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
! 4542: }
! 4543: if ( aExp == 0 ) {
! 4544: if ( ( aSig0 | aSig1 ) == 0 ) return packFloatx80( aSign, 0, 0 );
! 4545: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
! 4546: }
! 4547: else {
! 4548: aSig0 |= LIT64( 0x0001000000000000 );
! 4549: }
! 4550: shortShift128Left( aSig0, aSig1, 15, &aSig0, &aSig1 );
! 4551: return roundAndPackFloatx80( 80, aSign, aExp, aSig0, aSig1 );
! 4552:
! 4553: }
! 4554:
! 4555: #endif
! 4556:
! 4557: /*
! 4558: -------------------------------------------------------------------------------
! 4559: Rounds the quadruple-precision floating-point value `a' to an integer, and
! 4560: returns the result as a quadruple-precision floating-point value. The
! 4561: operation is performed according to the IEC/IEEE Standard for Binary
! 4562: Floating-Point Arithmetic.
! 4563: -------------------------------------------------------------------------------
! 4564: */
! 4565: float128 float128_round_to_int( float128 a )
! 4566: {
! 4567: flag aSign;
! 4568: int32 aExp;
! 4569: bits64 lastBitMask, roundBitsMask;
! 4570: int8 roundingMode;
! 4571: float128 z;
! 4572:
! 4573: aExp = extractFloat128Exp( a );
! 4574: if ( 0x402F <= aExp ) {
! 4575: if ( 0x406F <= aExp ) {
! 4576: if ( ( aExp == 0x7FFF )
! 4577: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) )
! 4578: ) {
! 4579: return propagateFloat128NaN( a, a );
! 4580: }
! 4581: return a;
! 4582: }
! 4583: lastBitMask = 1;
! 4584: lastBitMask = ( lastBitMask<<( 0x406E - aExp ) )<<1;
! 4585: roundBitsMask = lastBitMask - 1;
! 4586: z = a;
! 4587: roundingMode = float_rounding_mode();
! 4588: if ( roundingMode == float_round_nearest_even ) {
! 4589: if ( lastBitMask ) {
! 4590: add128( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
! 4591: if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
! 4592: }
! 4593: else {
! 4594: if ( (sbits64) z.low < 0 ) {
! 4595: ++z.high;
! 4596: if ( (bits64) ( z.low<<1 ) == 0 ) z.high &= ~1;
! 4597: }
! 4598: }
! 4599: }
! 4600: else if ( roundingMode != float_round_to_zero ) {
! 4601: if ( extractFloat128Sign( z )
! 4602: ^ ( roundingMode == float_round_up ) ) {
! 4603: add128( z.high, z.low, 0, roundBitsMask, &z.high, &z.low );
! 4604: }
! 4605: }
! 4606: z.low &= ~ roundBitsMask;
! 4607: }
! 4608: else {
! 4609: if ( aExp < 0x3FFF ) {
! 4610: if ( ( ( (bits64) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
! 4611: float_set_inexact();
! 4612: aSign = extractFloat128Sign( a );
! 4613: switch ( float_rounding_mode() ) {
! 4614: case float_round_nearest_even:
! 4615: if ( ( aExp == 0x3FFE )
! 4616: && ( extractFloat128Frac0( a )
! 4617: | extractFloat128Frac1( a ) )
! 4618: ) {
! 4619: return packFloat128( aSign, 0x3FFF, 0, 0 );
! 4620: }
! 4621: break;
! 4622: case float_round_down:
! 4623: return
! 4624: aSign ? packFloat128( 1, 0x3FFF, 0, 0 )
! 4625: : packFloat128( 0, 0, 0, 0 );
! 4626: case float_round_up:
! 4627: return
! 4628: aSign ? packFloat128( 1, 0, 0, 0 )
! 4629: : packFloat128( 0, 0x3FFF, 0, 0 );
! 4630: }
! 4631: return packFloat128( aSign, 0, 0, 0 );
! 4632: }
! 4633: lastBitMask = 1;
! 4634: lastBitMask <<= 0x402F - aExp;
! 4635: roundBitsMask = lastBitMask - 1;
! 4636: z.low = 0;
! 4637: z.high = a.high;
! 4638: roundingMode = float_rounding_mode();
! 4639: if ( roundingMode == float_round_nearest_even ) {
! 4640: z.high += lastBitMask>>1;
! 4641: if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
! 4642: z.high &= ~ lastBitMask;
! 4643: }
! 4644: }
! 4645: else if ( roundingMode != float_round_to_zero ) {
! 4646: if ( extractFloat128Sign( z )
! 4647: ^ ( roundingMode == float_round_up ) ) {
! 4648: z.high |= ( a.low != 0 );
! 4649: z.high += roundBitsMask;
! 4650: }
! 4651: }
! 4652: z.high &= ~ roundBitsMask;
! 4653: }
! 4654: if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
! 4655: float_set_inexact();
! 4656: }
! 4657: return z;
! 4658:
! 4659: }
! 4660:
! 4661: /*
! 4662: -------------------------------------------------------------------------------
! 4663: Returns the result of adding the absolute values of the quadruple-precision
! 4664: floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
! 4665: before being returned. `zSign' is ignored if the result is a NaN.
! 4666: The addition is performed according to the IEC/IEEE Standard for Binary
! 4667: Floating-Point Arithmetic.
! 4668: -------------------------------------------------------------------------------
! 4669: */
! 4670: static float128 addFloat128Sigs( float128 a, float128 b, flag zSign )
! 4671: {
! 4672: int32 aExp, bExp, zExp;
! 4673: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
! 4674: int32 expDiff;
! 4675:
! 4676: aSig1 = extractFloat128Frac1( a );
! 4677: aSig0 = extractFloat128Frac0( a );
! 4678: aExp = extractFloat128Exp( a );
! 4679: bSig1 = extractFloat128Frac1( b );
! 4680: bSig0 = extractFloat128Frac0( b );
! 4681: bExp = extractFloat128Exp( b );
! 4682: expDiff = aExp - bExp;
! 4683: if ( 0 < expDiff ) {
! 4684: if ( aExp == 0x7FFF ) {
! 4685: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b );
! 4686: return a;
! 4687: }
! 4688: if ( bExp == 0 ) {
! 4689: --expDiff;
! 4690: }
! 4691: else {
! 4692: bSig0 |= LIT64( 0x0001000000000000 );
! 4693: }
! 4694: shift128ExtraRightJamming(
! 4695: bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
! 4696: zExp = aExp;
! 4697: }
! 4698: else if ( expDiff < 0 ) {
! 4699: if ( bExp == 0x7FFF ) {
! 4700: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
! 4701: return packFloat128( zSign, 0x7FFF, 0, 0 );
! 4702: }
! 4703: if ( aExp == 0 ) {
! 4704: ++expDiff;
! 4705: }
! 4706: else {
! 4707: aSig0 |= LIT64( 0x0001000000000000 );
! 4708: }
! 4709: shift128ExtraRightJamming(
! 4710: aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
! 4711: zExp = bExp;
! 4712: }
! 4713: else {
! 4714: if ( aExp == 0x7FFF ) {
! 4715: if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
! 4716: return propagateFloat128NaN( a, b );
! 4717: }
! 4718: return a;
! 4719: }
! 4720: add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
! 4721: if ( aExp == 0 ) return packFloat128( zSign, 0, zSig0, zSig1 );
! 4722: zSig2 = 0;
! 4723: zSig0 |= LIT64( 0x0002000000000000 );
! 4724: zExp = aExp;
! 4725: goto shiftRight1;
! 4726: }
! 4727: aSig0 |= LIT64( 0x0001000000000000 );
! 4728: add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
! 4729: --zExp;
! 4730: if ( zSig0 < LIT64( 0x0002000000000000 ) ) goto roundAndPack;
! 4731: ++zExp;
! 4732: shiftRight1:
! 4733: shift128ExtraRightJamming(
! 4734: zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
! 4735: roundAndPack:
! 4736: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
! 4737:
! 4738: }
! 4739:
! 4740: /*
! 4741: -------------------------------------------------------------------------------
! 4742: Returns the result of subtracting the absolute values of the quadruple-
! 4743: precision floating-point values `a' and `b'. If `zSign' is 1, the
! 4744: difference is negated before being returned. `zSign' is ignored if the
! 4745: result is a NaN. The subtraction is performed according to the IEC/IEEE
! 4746: Standard for Binary Floating-Point Arithmetic.
! 4747: -------------------------------------------------------------------------------
! 4748: */
! 4749: static float128 subFloat128Sigs( float128 a, float128 b, flag zSign )
! 4750: {
! 4751: int32 aExp, bExp, zExp;
! 4752: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
! 4753: int32 expDiff;
! 4754: float128 z;
! 4755:
! 4756: aSig1 = extractFloat128Frac1( a );
! 4757: aSig0 = extractFloat128Frac0( a );
! 4758: aExp = extractFloat128Exp( a );
! 4759: bSig1 = extractFloat128Frac1( b );
! 4760: bSig0 = extractFloat128Frac0( b );
! 4761: bExp = extractFloat128Exp( b );
! 4762: expDiff = aExp - bExp;
! 4763: shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
! 4764: shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
! 4765: if ( 0 < expDiff ) goto aExpBigger;
! 4766: if ( expDiff < 0 ) goto bExpBigger;
! 4767: if ( aExp == 0x7FFF ) {
! 4768: if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
! 4769: return propagateFloat128NaN( a, b );
! 4770: }
! 4771: float_raise( float_flag_invalid );
! 4772: z.low = float128_default_nan_low;
! 4773: z.high = float128_default_nan_high;
! 4774: return z;
! 4775: }
! 4776: if ( aExp == 0 ) {
! 4777: aExp = 1;
! 4778: bExp = 1;
! 4779: }
! 4780: if ( bSig0 < aSig0 ) goto aBigger;
! 4781: if ( aSig0 < bSig0 ) goto bBigger;
! 4782: if ( bSig1 < aSig1 ) goto aBigger;
! 4783: if ( aSig1 < bSig1 ) goto bBigger;
! 4784: return packFloat128( float_rounding_mode() == float_round_down, 0, 0, 0 );
! 4785: bExpBigger:
! 4786: if ( bExp == 0x7FFF ) {
! 4787: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
! 4788: return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
! 4789: }
! 4790: if ( aExp == 0 ) {
! 4791: ++expDiff;
! 4792: }
! 4793: else {
! 4794: aSig0 |= LIT64( 0x4000000000000000 );
! 4795: }
! 4796: shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
! 4797: bSig0 |= LIT64( 0x4000000000000000 );
! 4798: bBigger:
! 4799: sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
! 4800: zExp = bExp;
! 4801: zSign ^= 1;
! 4802: goto normalizeRoundAndPack;
! 4803: aExpBigger:
! 4804: if ( aExp == 0x7FFF ) {
! 4805: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b );
! 4806: return a;
! 4807: }
! 4808: if ( bExp == 0 ) {
! 4809: --expDiff;
! 4810: }
! 4811: else {
! 4812: bSig0 |= LIT64( 0x4000000000000000 );
! 4813: }
! 4814: shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
! 4815: aSig0 |= LIT64( 0x4000000000000000 );
! 4816: aBigger:
! 4817: sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
! 4818: zExp = aExp;
! 4819: normalizeRoundAndPack:
! 4820: --zExp;
! 4821: return normalizeRoundAndPackFloat128( zSign, zExp - 14, zSig0, zSig1 );
! 4822:
! 4823: }
! 4824:
! 4825: /*
! 4826: -------------------------------------------------------------------------------
! 4827: Returns the result of adding the quadruple-precision floating-point values
! 4828: `a' and `b'. The operation is performed according to the IEC/IEEE Standard
! 4829: for Binary Floating-Point Arithmetic.
! 4830: -------------------------------------------------------------------------------
! 4831: */
! 4832: float128 float128_add( float128 a, float128 b )
! 4833: {
! 4834: flag aSign, bSign;
! 4835:
! 4836: aSign = extractFloat128Sign( a );
! 4837: bSign = extractFloat128Sign( b );
! 4838: if ( aSign == bSign ) {
! 4839: return addFloat128Sigs( a, b, aSign );
! 4840: }
! 4841: else {
! 4842: return subFloat128Sigs( a, b, aSign );
! 4843: }
! 4844:
! 4845: }
! 4846:
! 4847: /*
! 4848: -------------------------------------------------------------------------------
! 4849: Returns the result of subtracting the quadruple-precision floating-point
! 4850: values `a' and `b'. The operation is performed according to the IEC/IEEE
! 4851: Standard for Binary Floating-Point Arithmetic.
! 4852: -------------------------------------------------------------------------------
! 4853: */
! 4854: float128 float128_sub( float128 a, float128 b )
! 4855: {
! 4856: flag aSign, bSign;
! 4857:
! 4858: aSign = extractFloat128Sign( a );
! 4859: bSign = extractFloat128Sign( b );
! 4860: if ( aSign == bSign ) {
! 4861: return subFloat128Sigs( a, b, aSign );
! 4862: }
! 4863: else {
! 4864: return addFloat128Sigs( a, b, aSign );
! 4865: }
! 4866:
! 4867: }
! 4868:
! 4869: /*
! 4870: -------------------------------------------------------------------------------
! 4871: Returns the result of multiplying the quadruple-precision floating-point
! 4872: values `a' and `b'. The operation is performed according to the IEC/IEEE
! 4873: Standard for Binary Floating-Point Arithmetic.
! 4874: -------------------------------------------------------------------------------
! 4875: */
! 4876: float128 float128_mul( float128 a, float128 b )
! 4877: {
! 4878: flag aSign, bSign, zSign;
! 4879: int32 aExp, bExp, zExp;
! 4880: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
! 4881: float128 z;
! 4882:
! 4883: aSig1 = extractFloat128Frac1( a );
! 4884: aSig0 = extractFloat128Frac0( a );
! 4885: aExp = extractFloat128Exp( a );
! 4886: aSign = extractFloat128Sign( a );
! 4887: bSig1 = extractFloat128Frac1( b );
! 4888: bSig0 = extractFloat128Frac0( b );
! 4889: bExp = extractFloat128Exp( b );
! 4890: bSign = extractFloat128Sign( b );
! 4891: zSign = aSign ^ bSign;
! 4892: if ( aExp == 0x7FFF ) {
! 4893: if ( ( aSig0 | aSig1 )
! 4894: || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
! 4895: return propagateFloat128NaN( a, b );
! 4896: }
! 4897: if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
! 4898: return packFloat128( zSign, 0x7FFF, 0, 0 );
! 4899: }
! 4900: if ( bExp == 0x7FFF ) {
! 4901: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
! 4902: if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
! 4903: invalid:
! 4904: float_raise( float_flag_invalid );
! 4905: z.low = float128_default_nan_low;
! 4906: z.high = float128_default_nan_high;
! 4907: return z;
! 4908: }
! 4909: return packFloat128( zSign, 0x7FFF, 0, 0 );
! 4910: }
! 4911: if ( aExp == 0 ) {
! 4912: if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
! 4913: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
! 4914: }
! 4915: if ( bExp == 0 ) {
! 4916: if ( ( bSig0 | bSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
! 4917: normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
! 4918: }
! 4919: zExp = aExp + bExp - 0x4000;
! 4920: aSig0 |= LIT64( 0x0001000000000000 );
! 4921: shortShift128Left( bSig0, bSig1, 16, &bSig0, &bSig1 );
! 4922: mul128To256( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
! 4923: add128( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
! 4924: zSig2 |= ( zSig3 != 0 );
! 4925: if ( LIT64( 0x0002000000000000 ) <= zSig0 ) {
! 4926: shift128ExtraRightJamming(
! 4927: zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
! 4928: ++zExp;
! 4929: }
! 4930: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
! 4931:
! 4932: }
! 4933:
! 4934: /*
! 4935: -------------------------------------------------------------------------------
! 4936: Returns the result of dividing the quadruple-precision floating-point value
! 4937: `a' by the corresponding value `b'. The operation is performed according to
! 4938: the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 4939: -------------------------------------------------------------------------------
! 4940: */
! 4941: float128 float128_div( float128 a, float128 b )
! 4942: {
! 4943: flag aSign, bSign, zSign;
! 4944: int32 aExp, bExp, zExp;
! 4945: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
! 4946: bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
! 4947: float128 z;
! 4948:
! 4949: aSig1 = extractFloat128Frac1( a );
! 4950: aSig0 = extractFloat128Frac0( a );
! 4951: aExp = extractFloat128Exp( a );
! 4952: aSign = extractFloat128Sign( a );
! 4953: bSig1 = extractFloat128Frac1( b );
! 4954: bSig0 = extractFloat128Frac0( b );
! 4955: bExp = extractFloat128Exp( b );
! 4956: bSign = extractFloat128Sign( b );
! 4957: zSign = aSign ^ bSign;
! 4958: if ( aExp == 0x7FFF ) {
! 4959: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b );
! 4960: if ( bExp == 0x7FFF ) {
! 4961: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
! 4962: goto invalid;
! 4963: }
! 4964: return packFloat128( zSign, 0x7FFF, 0, 0 );
! 4965: }
! 4966: if ( bExp == 0x7FFF ) {
! 4967: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
! 4968: return packFloat128( zSign, 0, 0, 0 );
! 4969: }
! 4970: if ( bExp == 0 ) {
! 4971: if ( ( bSig0 | bSig1 ) == 0 ) {
! 4972: if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
! 4973: invalid:
! 4974: float_raise( float_flag_invalid );
! 4975: z.low = float128_default_nan_low;
! 4976: z.high = float128_default_nan_high;
! 4977: return z;
! 4978: }
! 4979: float_raise( float_flag_divbyzero );
! 4980: return packFloat128( zSign, 0x7FFF, 0, 0 );
! 4981: }
! 4982: normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
! 4983: }
! 4984: if ( aExp == 0 ) {
! 4985: if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
! 4986: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
! 4987: }
! 4988: zExp = aExp - bExp + 0x3FFD;
! 4989: shortShift128Left(
! 4990: aSig0 | LIT64( 0x0001000000000000 ), aSig1, 15, &aSig0, &aSig1 );
! 4991: shortShift128Left(
! 4992: bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
! 4993: if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
! 4994: shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
! 4995: ++zExp;
! 4996: }
! 4997: zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
! 4998: mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
! 4999: sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
! 5000: while ( (sbits64) rem0 < 0 ) {
! 5001: --zSig0;
! 5002: add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
! 5003: }
! 5004: zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
! 5005: if ( ( zSig1 & 0x3FFF ) <= 4 ) {
! 5006: mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
! 5007: sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
! 5008: while ( (sbits64) rem1 < 0 ) {
! 5009: --zSig1;
! 5010: add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
! 5011: }
! 5012: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
! 5013: }
! 5014: shift128ExtraRightJamming( zSig0, zSig1, 0, 15, &zSig0, &zSig1, &zSig2 );
! 5015: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
! 5016:
! 5017: }
! 5018:
! 5019: /*
! 5020: -------------------------------------------------------------------------------
! 5021: Returns the remainder of the quadruple-precision floating-point value `a'
! 5022: with respect to the corresponding value `b'. The operation is performed
! 5023: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 5024: -------------------------------------------------------------------------------
! 5025: */
! 5026: float128 float128_rem( float128 a, float128 b )
! 5027: {
! 5028: flag aSign, bSign, zSign;
! 5029: int32 aExp, bExp, expDiff;
! 5030: bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
! 5031: bits64 allZero, alternateASig0, alternateASig1, sigMean1;
! 5032: sbits64 sigMean0;
! 5033: float128 z;
! 5034:
! 5035: aSig1 = extractFloat128Frac1( a );
! 5036: aSig0 = extractFloat128Frac0( a );
! 5037: aExp = extractFloat128Exp( a );
! 5038: aSign = extractFloat128Sign( a );
! 5039: bSig1 = extractFloat128Frac1( b );
! 5040: bSig0 = extractFloat128Frac0( b );
! 5041: bExp = extractFloat128Exp( b );
! 5042: bSign = extractFloat128Sign( b );
! 5043: if ( aExp == 0x7FFF ) {
! 5044: if ( ( aSig0 | aSig1 )
! 5045: || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
! 5046: return propagateFloat128NaN( a, b );
! 5047: }
! 5048: goto invalid;
! 5049: }
! 5050: if ( bExp == 0x7FFF ) {
! 5051: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
! 5052: return a;
! 5053: }
! 5054: if ( bExp == 0 ) {
! 5055: if ( ( bSig0 | bSig1 ) == 0 ) {
! 5056: invalid:
! 5057: float_raise( float_flag_invalid );
! 5058: z.low = float128_default_nan_low;
! 5059: z.high = float128_default_nan_high;
! 5060: return z;
! 5061: }
! 5062: normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
! 5063: }
! 5064: if ( aExp == 0 ) {
! 5065: if ( ( aSig0 | aSig1 ) == 0 ) return a;
! 5066: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
! 5067: }
! 5068: expDiff = aExp - bExp;
! 5069: if ( expDiff < -1 ) return a;
! 5070: shortShift128Left(
! 5071: aSig0 | LIT64( 0x0001000000000000 ),
! 5072: aSig1,
! 5073: 15 - ( expDiff < 0 ),
! 5074: &aSig0,
! 5075: &aSig1
! 5076: );
! 5077: shortShift128Left(
! 5078: bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
! 5079: q = le128( bSig0, bSig1, aSig0, aSig1 );
! 5080: if ( q ) sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
! 5081: expDiff -= 64;
! 5082: while ( 0 < expDiff ) {
! 5083: q = estimateDiv128To64( aSig0, aSig1, bSig0 );
! 5084: q = ( 4 < q ) ? q - 4 : 0;
! 5085: mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
! 5086: shortShift192Left( term0, term1, term2, 61, &term1, &term2, &allZero );
! 5087: shortShift128Left( aSig0, aSig1, 61, &aSig0, &allZero );
! 5088: sub128( aSig0, 0, term1, term2, &aSig0, &aSig1 );
! 5089: expDiff -= 61;
! 5090: }
! 5091: if ( -64 < expDiff ) {
! 5092: q = estimateDiv128To64( aSig0, aSig1, bSig0 );
! 5093: q = ( 4 < q ) ? q - 4 : 0;
! 5094: q >>= - expDiff;
! 5095: shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
! 5096: expDiff += 52;
! 5097: if ( expDiff < 0 ) {
! 5098: shift128Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
! 5099: }
! 5100: else {
! 5101: shortShift128Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 );
! 5102: }
! 5103: mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
! 5104: sub128( aSig0, aSig1, term1, term2, &aSig0, &aSig1 );
! 5105: }
! 5106: else {
! 5107: shift128Right( aSig0, aSig1, 12, &aSig0, &aSig1 );
! 5108: shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
! 5109: }
! 5110: do {
! 5111: alternateASig0 = aSig0;
! 5112: alternateASig1 = aSig1;
! 5113: ++q;
! 5114: sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
! 5115: } while ( 0 <= (sbits64) aSig0 );
! 5116: add128(
! 5117: aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 );
! 5118: if ( ( sigMean0 < 0 )
! 5119: || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
! 5120: aSig0 = alternateASig0;
! 5121: aSig1 = alternateASig1;
! 5122: }
! 5123: zSign = ( (sbits64) aSig0 < 0 );
! 5124: if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
! 5125: return
! 5126: normalizeRoundAndPackFloat128( aSign ^ zSign, bExp - 4, aSig0, aSig1 );
! 5127:
! 5128: }
! 5129:
! 5130: /*
! 5131: -------------------------------------------------------------------------------
! 5132: Returns the square root of the quadruple-precision floating-point value `a'.
! 5133: The operation is performed according to the IEC/IEEE Standard for Binary
! 5134: Floating-Point Arithmetic.
! 5135: -------------------------------------------------------------------------------
! 5136: */
! 5137: float128 float128_sqrt( float128 a )
! 5138: {
! 5139: flag aSign;
! 5140: int32 aExp, zExp;
! 5141: bits64 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
! 5142: bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
! 5143: float128 z;
! 5144:
! 5145: aSig1 = extractFloat128Frac1( a );
! 5146: aSig0 = extractFloat128Frac0( a );
! 5147: aExp = extractFloat128Exp( a );
! 5148: aSign = extractFloat128Sign( a );
! 5149: if ( aExp == 0x7FFF ) {
! 5150: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, a );
! 5151: if ( ! aSign ) return a;
! 5152: goto invalid;
! 5153: }
! 5154: if ( aSign ) {
! 5155: if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a;
! 5156: invalid:
! 5157: float_raise( float_flag_invalid );
! 5158: z.low = float128_default_nan_low;
! 5159: z.high = float128_default_nan_high;
! 5160: return z;
! 5161: }
! 5162: if ( aExp == 0 ) {
! 5163: if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( 0, 0, 0, 0 );
! 5164: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
! 5165: }
! 5166: zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFE;
! 5167: aSig0 |= LIT64( 0x0001000000000000 );
! 5168: zSig0 = estimateSqrt32( aExp, aSig0>>17 );
! 5169: shortShift128Left( aSig0, aSig1, 13 - ( aExp & 1 ), &aSig0, &aSig1 );
! 5170: zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
! 5171: doubleZSig0 = zSig0<<1;
! 5172: mul64To128( zSig0, zSig0, &term0, &term1 );
! 5173: sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
! 5174: while ( (sbits64) rem0 < 0 ) {
! 5175: --zSig0;
! 5176: doubleZSig0 -= 2;
! 5177: add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
! 5178: }
! 5179: zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
! 5180: if ( ( zSig1 & 0x1FFF ) <= 5 ) {
! 5181: if ( zSig1 == 0 ) zSig1 = 1;
! 5182: mul64To128( doubleZSig0, zSig1, &term1, &term2 );
! 5183: sub128( rem1, 0, term1, term2, &rem1, &rem2 );
! 5184: mul64To128( zSig1, zSig1, &term2, &term3 );
! 5185: sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
! 5186: while ( (sbits64) rem1 < 0 ) {
! 5187: --zSig1;
! 5188: shortShift128Left( 0, zSig1, 1, &term2, &term3 );
! 5189: term3 |= 1;
! 5190: term2 |= doubleZSig0;
! 5191: add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
! 5192: }
! 5193: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
! 5194: }
! 5195: shift128ExtraRightJamming( zSig0, zSig1, 0, 14, &zSig0, &zSig1, &zSig2 );
! 5196: return roundAndPackFloat128( 0, zExp, zSig0, zSig1, zSig2 );
! 5197:
! 5198: }
! 5199:
! 5200: /*
! 5201: -------------------------------------------------------------------------------
! 5202: Returns 1 if the quadruple-precision floating-point value `a' is equal to
! 5203: the corresponding value `b', and 0 otherwise. The comparison is performed
! 5204: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 5205: -------------------------------------------------------------------------------
! 5206: */
! 5207: flag float128_eq( float128 a, float128 b )
! 5208: {
! 5209:
! 5210: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
! 5211: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
! 5212: || ( ( extractFloat128Exp( b ) == 0x7FFF )
! 5213: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
! 5214: ) {
! 5215: if ( float128_is_signaling_nan( a )
! 5216: || float128_is_signaling_nan( b ) ) {
! 5217: float_raise( float_flag_invalid );
! 5218: }
! 5219: return 0;
! 5220: }
! 5221: return
! 5222: ( a.low == b.low )
! 5223: && ( ( a.high == b.high )
! 5224: || ( ( a.low == 0 )
! 5225: && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
! 5226: );
! 5227:
! 5228: }
! 5229:
! 5230: /*
! 5231: -------------------------------------------------------------------------------
! 5232: Returns 1 if the quadruple-precision floating-point value `a' is less than
! 5233: or equal to the corresponding value `b', and 0 otherwise. The comparison
! 5234: is performed according to the IEC/IEEE Standard for Binary Floating-Point
! 5235: Arithmetic.
! 5236: -------------------------------------------------------------------------------
! 5237: */
! 5238: flag float128_le( float128 a, float128 b )
! 5239: {
! 5240: flag aSign, bSign;
! 5241:
! 5242: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
! 5243: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
! 5244: || ( ( extractFloat128Exp( b ) == 0x7FFF )
! 5245: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
! 5246: ) {
! 5247: float_raise( float_flag_invalid );
! 5248: return 0;
! 5249: }
! 5250: aSign = extractFloat128Sign( a );
! 5251: bSign = extractFloat128Sign( b );
! 5252: if ( aSign != bSign ) {
! 5253: return
! 5254: aSign
! 5255: || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 5256: == 0 );
! 5257: }
! 5258: return
! 5259: aSign ? le128( b.high, b.low, a.high, a.low )
! 5260: : le128( a.high, a.low, b.high, b.low );
! 5261:
! 5262: }
! 5263:
! 5264: /*
! 5265: -------------------------------------------------------------------------------
! 5266: Returns 1 if the quadruple-precision floating-point value `a' is less than
! 5267: the corresponding value `b', and 0 otherwise. The comparison is performed
! 5268: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 5269: -------------------------------------------------------------------------------
! 5270: */
! 5271: flag float128_lt( float128 a, float128 b )
! 5272: {
! 5273: flag aSign, bSign;
! 5274:
! 5275: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
! 5276: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
! 5277: || ( ( extractFloat128Exp( b ) == 0x7FFF )
! 5278: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
! 5279: ) {
! 5280: float_raise( float_flag_invalid );
! 5281: return 0;
! 5282: }
! 5283: aSign = extractFloat128Sign( a );
! 5284: bSign = extractFloat128Sign( b );
! 5285: if ( aSign != bSign ) {
! 5286: return
! 5287: aSign
! 5288: && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 5289: != 0 );
! 5290: }
! 5291: return
! 5292: aSign ? lt128( b.high, b.low, a.high, a.low )
! 5293: : lt128( a.high, a.low, b.high, b.low );
! 5294:
! 5295: }
! 5296:
! 5297: /*
! 5298: -------------------------------------------------------------------------------
! 5299: Returns 1 if the quadruple-precision floating-point value `a' is equal to
! 5300: the corresponding value `b', and 0 otherwise. The invalid exception is
! 5301: raised if either operand is a NaN. Otherwise, the comparison is performed
! 5302: according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 5303: -------------------------------------------------------------------------------
! 5304: */
! 5305: flag float128_eq_signaling( float128 a, float128 b )
! 5306: {
! 5307:
! 5308: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
! 5309: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
! 5310: || ( ( extractFloat128Exp( b ) == 0x7FFF )
! 5311: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
! 5312: ) {
! 5313: float_raise( float_flag_invalid );
! 5314: return 0;
! 5315: }
! 5316: return
! 5317: ( a.low == b.low )
! 5318: && ( ( a.high == b.high )
! 5319: || ( ( a.low == 0 )
! 5320: && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
! 5321: );
! 5322:
! 5323: }
! 5324:
! 5325: /*
! 5326: -------------------------------------------------------------------------------
! 5327: Returns 1 if the quadruple-precision floating-point value `a' is less than
! 5328: or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
! 5329: cause an exception. Otherwise, the comparison is performed according to the
! 5330: IEC/IEEE Standard for Binary Floating-Point Arithmetic.
! 5331: -------------------------------------------------------------------------------
! 5332: */
! 5333: flag float128_le_quiet( float128 a, float128 b )
! 5334: {
! 5335: flag aSign, bSign;
! 5336:
! 5337: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
! 5338: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
! 5339: || ( ( extractFloat128Exp( b ) == 0x7FFF )
! 5340: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
! 5341: ) {
! 5342: if ( float128_is_signaling_nan( a )
! 5343: || float128_is_signaling_nan( b ) ) {
! 5344: float_raise( float_flag_invalid );
! 5345: }
! 5346: return 0;
! 5347: }
! 5348: aSign = extractFloat128Sign( a );
! 5349: bSign = extractFloat128Sign( b );
! 5350: if ( aSign != bSign ) {
! 5351: return
! 5352: aSign
! 5353: || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 5354: == 0 );
! 5355: }
! 5356: return
! 5357: aSign ? le128( b.high, b.low, a.high, a.low )
! 5358: : le128( a.high, a.low, b.high, b.low );
! 5359:
! 5360: }
! 5361:
! 5362: /*
! 5363: -------------------------------------------------------------------------------
! 5364: Returns 1 if the quadruple-precision floating-point value `a' is less than
! 5365: the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
! 5366: exception. Otherwise, the comparison is performed according to the IEC/IEEE
! 5367: Standard for Binary Floating-Point Arithmetic.
! 5368: -------------------------------------------------------------------------------
! 5369: */
! 5370: flag float128_lt_quiet( float128 a, float128 b )
! 5371: {
! 5372: flag aSign, bSign;
! 5373:
! 5374: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
! 5375: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
! 5376: || ( ( extractFloat128Exp( b ) == 0x7FFF )
! 5377: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
! 5378: ) {
! 5379: if ( float128_is_signaling_nan( a )
! 5380: || float128_is_signaling_nan( b ) ) {
! 5381: float_raise( float_flag_invalid );
! 5382: }
! 5383: return 0;
! 5384: }
! 5385: aSign = extractFloat128Sign( a );
! 5386: bSign = extractFloat128Sign( b );
! 5387: if ( aSign != bSign ) {
! 5388: return
! 5389: aSign
! 5390: && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
! 5391: != 0 );
! 5392: }
! 5393: return
! 5394: aSign ? lt128( b.high, b.low, a.high, a.low )
! 5395: : lt128( a.high, a.low, b.high, b.low );
! 5396:
! 5397: }
! 5398:
! 5399: #endif
! 5400:
! 5401:
! 5402: #if defined(SOFTFLOAT_FOR_GCC) && defined(SOFTFLOAT_NEED_FIXUNS)
! 5403:
! 5404: /*
! 5405: * These two routines are not part of the original softfloat distribution.
! 5406: *
! 5407: * They are based on the corresponding conversions to integer but return
! 5408: * unsigned numbers instead since these functions are required by GCC.
! 5409: *
! 5410: * Added by Mark Brinicombe <mark@netbsd.org> 27/09/97
! 5411: *
! 5412: * float64 version overhauled for SoftFloat 2a [bjh21 2000-07-15]
! 5413: */
! 5414:
! 5415: /*
! 5416: -------------------------------------------------------------------------------
! 5417: Returns the result of converting the double-precision floating-point value
! 5418: `a' to the 32-bit unsigned integer format. The conversion is
! 5419: performed according to the IEC/IEEE Standard for Binary Floating-point
! 5420: Arithmetic, except that the conversion is always rounded toward zero. If
! 5421: `a' is a NaN, the largest positive integer is returned. If the conversion
! 5422: overflows, the largest integer positive is returned.
! 5423: -------------------------------------------------------------------------------
! 5424: */
! 5425: uint32 float64_to_uint32_round_to_zero( float64 a )
! 5426: {
! 5427: flag aSign;
! 5428: int16 aExp, shiftCount;
! 5429: bits64 aSig, savedASig;
! 5430: uint32 z;
! 5431:
! 5432: aSig = extractFloat64Frac( a );
! 5433: aExp = extractFloat64Exp( a );
! 5434: aSign = extractFloat64Sign( a );
! 5435:
! 5436: if (aSign) {
! 5437: float_raise( float_flag_invalid );
! 5438: return(0);
! 5439: }
! 5440:
! 5441: if ( 0x41E < aExp ) {
! 5442: float_raise( float_flag_invalid );
! 5443: return 0xffffffff;
! 5444: }
! 5445: else if ( aExp < 0x3FF ) {
! 5446: if ( aExp || aSig ) float_set_inexact();
! 5447: return 0;
! 5448: }
! 5449: aSig |= LIT64( 0x0010000000000000 );
! 5450: shiftCount = 0x433 - aExp;
! 5451: savedASig = aSig;
! 5452: aSig >>= shiftCount;
! 5453: z = aSig;
! 5454: if ( ( aSig<<shiftCount ) != savedASig ) {
! 5455: float_set_inexact();
! 5456: }
! 5457: return z;
! 5458:
! 5459: }
! 5460:
! 5461: /*
! 5462: -------------------------------------------------------------------------------
! 5463: Returns the result of converting the single-precision floating-point value
! 5464: `a' to the 32-bit unsigned integer format. The conversion is
! 5465: performed according to the IEC/IEEE Standard for Binary Floating-point
! 5466: Arithmetic, except that the conversion is always rounded toward zero. If
! 5467: `a' is a NaN, the largest positive integer is returned. If the conversion
! 5468: overflows, the largest positive integer is returned.
! 5469: -------------------------------------------------------------------------------
! 5470: */
! 5471: uint32 float32_to_uint32_round_to_zero( float32 a )
! 5472: {
! 5473: flag aSign;
! 5474: int16 aExp, shiftCount;
! 5475: bits32 aSig;
! 5476: uint32 z;
! 5477:
! 5478: aSig = extractFloat32Frac( a );
! 5479: aExp = extractFloat32Exp( a );
! 5480: aSign = extractFloat32Sign( a );
! 5481: shiftCount = aExp - 0x9E;
! 5482:
! 5483: if (aSign) {
! 5484: float_raise( float_flag_invalid );
! 5485: return(0);
! 5486: }
! 5487: if ( 0 < shiftCount ) {
! 5488: float_raise( float_flag_invalid );
! 5489: return 0xFFFFFFFF;
! 5490: }
! 5491: else if ( aExp <= 0x7E ) {
! 5492: if ( aExp | aSig ) float_set_inexact();
! 5493: return 0;
! 5494: }
! 5495: aSig = ( aSig | 0x800000 )<<8;
! 5496: z = aSig>>( - shiftCount );
! 5497: if ( aSig<<( shiftCount & 31 ) ) {
! 5498: float_set_inexact();
! 5499: }
! 5500: return z;
! 5501:
! 5502: }
! 5503:
! 5504: #endif
! 5505:
! 5506: #endif /* !NO_IEEE */
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to softfloat.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / lib / libkern