a subnormal number, and the underflow and inexact exceptions are raised ifthe abstract input cannot be represented exactly as a subnormal double-precision floating-point number.    The input significand `zSig' has its binary point between bits 62and 61, which is 10 bits to the left of the usual location.  This shiftedsignificand must be normalized or smaller.  If `zSig' is not normalized,`zExp' must be 0; in that case, the result returned is a subnormal number,and it must not require rounding.  In the usual case that `zSig' isnormalized, `zExp' must be 1 less than the ``true'' floating-point exponent.The handling of underflow and overflow follows the IEC/IEEE Standard forBinary Floating-point Arithmetic.-------------------------------------------------------------------------------*/static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig ){    int8 roundingMode;    flag roundNearestEven;    int16 roundIncrement, roundBits;    flag isTiny;    roundingMode = float_rounding_mode;    roundNearestEven = ( roundingMode == float_round_nearest_even );    roundIncrement = 0x200;    if ( ! roundNearestEven ) {        if ( roundingMode == float_round_to_zero ) {            roundIncrement = 0;        }        else {            roundIncrement = 0x3FF;            if ( zSign ) {                if ( roundingMode == float_round_up ) roundIncrement = 0;            }            else {                if ( roundingMode == float_round_down ) roundIncrement = 0;            }        }    }    roundBits = zSig & 0x3FF;    if ( 0x7FD <= (bits16) zExp ) {        if (    ( 0x7FD < zExp )             || (    ( zExp == 0x7FD )                  && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )           ) {            //register int lr = __builtin_return_address(0);            //printk("roundAndPackFloat64 called from 0x%08x\n",lr);            float_raise( float_flag_overflow | float_flag_inexact );            return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );        }        if ( zExp < 0 ) {            isTiny =                   ( float_detect_tininess == float_tininess_before_rounding )                || ( zExp < -1 )                || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );            shift64RightJamming( zSig, - zExp, &zSig );            zExp = 0;            roundBits = zSig & 0x3FF;            if ( isTiny && roundBits ) float_raise( float_flag_underflow );        }    }    if ( roundBits ) float_exception_flags |= float_flag_inexact;    zSig = ( zSig + roundIncrement )>>10;    zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );    if ( zSig == 0 ) zExp = 0;    return packFloat64( zSign, zExp, zSig );}/*-------------------------------------------------------------------------------Takes an abstract floating-point value having sign `zSign', exponent `zExp',and significand `zSig', and returns the proper double-precision floating-point value corresponding to the abstract input.  This routine is just like`roundAndPackFloat64' except that `zSig' does not have to be normalized inany way.  In all cases, `zExp' must be 1 less than the ``true'' floating-point exponent.-------------------------------------------------------------------------------*/static float64 normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig ){    int8 shiftCount;    shiftCount = countLeadingZeros64( zSig ) - 1;    return roundAndPackFloat64( zSign, zExp - shiftCount, zSig<<shiftCount );}#ifdef FLOATX80/*-------------------------------------------------------------------------------Returns the fraction bits of the extended double-precision floating-pointvalue `a'.-------------------------------------------------------------------------------*/INLINE bits64 extractFloatx80Frac( floatx80 a ){    return a.low;}/*-------------------------------------------------------------------------------Returns the exponent bits of the extended double-precision floating-pointvalue `a'.-------------------------------------------------------------------------------*/INLINE int32 extractFloatx80Exp( floatx80 a ){    return a.high & 0x7FFF;}/*-------------------------------------------------------------------------------Returns the sign bit of the extended double-precision floating-point value`a'.-------------------------------------------------------------------------------*/INLINE flag extractFloatx80Sign( floatx80 a ){    return a.high>>15;}/*-------------------------------------------------------------------------------Normalizes the subnormal extended double-precision floating-point valuerepresented by the denormalized significand `aSig'.  The normalized exponentand significand are stored at the locations pointed to by `zExpPtr' and`zSigPtr', respectively.-------------------------------------------------------------------------------*/static void normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr ){    int8 shiftCount;    shiftCount = countLeadingZeros64( aSig );    *zSigPtr = aSig<<shiftCount;    *zExpPtr = 1 - shiftCount;}/*-------------------------------------------------------------------------------Packs the sign `zSign', exponent `zExp', and significand `zSig' into anextended double-precision floating-point value, returning the result.-------------------------------------------------------------------------------*/INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig ){    floatx80 z;    z.low = zSig;    z.high = ( ( (bits16) zSign )<<15 ) + zExp;    return z;}/*-------------------------------------------------------------------------------Takes an abstract floating-point value having sign `zSign', exponent `zExp',and extended significand formed by the concatenation of `zSig0' and `zSig1',and returns the proper extended double-precision floating-point valuecorresponding to the abstract input.  Ordinarily, the abstract value isrounded and packed into the extended double-precision format, with theinexact exception raised if the abstract input cannot be representedexactly.  If the abstract value is too large, however, the overflow andinexact exceptions are raised and an infinity or maximal finite value isreturned.  If the abstract value is too small, the input value is rounded toa subnormal number, and the underflow and inexact exceptions are raised ifthe abstract input cannot be represented exactly as a subnormal extendeddouble-precision floating-point number.    If `roundingPrecision' is 32 or 64, the result is rounded to the samenumber of bits as single or double precision, respectively.  Otherwise, theresult is rounded to the full precision of the extended double-precisionformat.    The input significand must be normalized or smaller.  If the inputsignificand is not normalized, `zExp' must be 0; in that case, the resultreturned is a subnormal number, and it must not require rounding.  Thehandling of underflow and overflow follows the IEC/IEEE Standard for BinaryFloating-point Arithmetic.-------------------------------------------------------------------------------*/static floatx80 roundAndPackFloatx80(     int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 ){    int8 roundingMode;    flag roundNearestEven, increment, isTiny;    int64 roundIncrement, roundMask, roundBits;    roundingMode = float_rounding_mode;    roundNearestEven = ( roundingMode == float_round_nearest_even );    if ( roundingPrecision == 80 ) goto precision80;    if ( roundingPrecision == 64 ) {        roundIncrement = LIT64( 0x0000000000000400 );        roundMask = LIT64( 0x00000000000007FF );    }    else if ( roundingPrecision == 32 ) {        roundIncrement = LIT64( 0x0000008000000000 );        roundMask = LIT64( 0x000000FFFFFFFFFF );    }    else {        goto precision80;    }    zSig0 |= ( zSig1 != 0 );    if ( ! roundNearestEven ) {        if ( roundingMode == float_round_to_zero ) {            roundIncrement = 0;        }        else {            roundIncrement = roundMask;            if ( zSign ) {                if ( roundingMode == float_round_up ) roundIncrement = 0;            }            else {                if ( roundingMode == float_round_down ) roundIncrement = 0;            }        }    }    roundBits = zSig0 & roundMask;    if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {        if (    ( 0x7FFE < zExp )             || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )           ) {            goto overflow;        }        if ( zExp <= 0 ) {            isTiny =                   ( float_detect_tininess == float_tininess_before_rounding )                || ( zExp < 0 )                || ( zSig0 <= zSig0 + roundIncrement );            shift64RightJamming( zSig0, 1 - zExp, &zSig0 );            zExp = 0;            roundBits = zSig0 & roundMask;            if ( isTiny && roundBits ) float_raise( float_flag_underflow );            if ( roundBits ) float_exception_flags |= float_flag_inexact;            zSig0 += roundIncrement;            if ( (sbits64) zSig0 < 0 ) zExp = 1;            roundIncrement = roundMask + 1;            if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {                roundMask |= roundIncrement;            }            zSig0 &= ~ roundMask;            return packFloatx80( zSign, zExp, zSig0 );        }    }    if ( roundBits ) float_exception_flags |= float_flag_inexact;    zSig0 += roundIncrement;    if ( zSig0 < roundIncrement ) {        ++zExp;        zSig0 = LIT64( 0x8000000000000000 );    }    roundIncrement = roundMask + 1;    if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {        roundMask |= roundIncrement;    }    zSig0 &= ~ roundMask;    if ( zSig0 == 0 ) zExp = 0;    return packFloatx80( zSign, zExp, zSig0 ); precision80:    increment = ( (sbits64) zSig1 < 0 );    if ( ! roundNearestEven ) {        if ( roundingMode == float_round_to_zero ) {            increment = 0;        }        else {            if ( zSign ) {                increment = ( roundingMode == float_round_down ) && zSig1;            }            else {                increment = ( roundingMode == float_round_up ) && zSig1;            }        }    }    if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {        if (    ( 0x7FFE < zExp )             || (    ( zExp == 0x7FFE )                  && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )                  && increment                )           ) {            roundMask = 0; overflow:            float_raise( float_flag_overflow | float_flag_inexact );            if (    ( roundingMode == float_round_to_zero )                 || ( zSign && ( roundingMode == float_round_up ) )                 || ( ! zSign && ( roundingMode == float_round_down ) )               ) {                return packFloatx80( zSign, 0x7FFE, ~ roundMask );            }            return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );        }        if ( zExp <= 0 ) {            isTiny =                   ( float_detect_tininess == float_tininess_before_rounding )                || ( zExp < 0 )                || ! increment                || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );            shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );            zExp = 0;            if ( isTiny && zSig1 ) float_raise( float_flag_underflow );            if ( zSig1 ) float_exception_flags |= float_flag_inexact;            if ( roundNearestEven ) {                increment = ( (sbits64) zSig1 < 0 );            }            else {                if ( zSign ) {                    increment = ( roundingMode == float_round_down ) && zSig1;                }                else {                    increment = ( roundingMode == float_round_up ) && zSig1;                }            }            if ( increment ) {                ++zSig0;                zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven );                if ( (sbits64) zSig0 < 0 ) zExp = 1;            }            return packFloatx80( zSign, zExp, zSig0 );        }    }    if ( zSig1 ) float_exception_flags |= float_flag_inexact;    if ( increment ) {        ++zSig0;        if ( zSig0 == 0 ) {            ++zExp;            zSig0 = LIT64( 0x8000000000000000 );        }        else {            zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven );        }    }    else {        if ( zSig0 == 0 ) zExp = 0;    }        return packFloatx80( zSign, zExp, zSig0 );}/*-------------------------------------------------------------------------------Takes an abstract floating-point value having sign `zSign', exponent`zExp', and significand formed by the concatenation of `zSig0' and `zSig1',and returns the proper extended double-precision floating-point valuecorresponding to the abstract input.  This routine is just like`roundAndPackFloatx80' except that the input significand does not have to benormalized.-------------------------------------------------------------------------------*/static floatx80 normalizeRoundAndPackFloatx80(     int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 ){    int8 shiftCount;    if ( zSig0 == 0 ) {        zSig0 = zSig1;        zSig1 = 0;        zExp -= 64;    }    shiftCount = countLeadingZeros64( zSig0 );    shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );    zExp -= shiftCount;    return        roundAndPackFloatx80( roundingPrecision, zSign, zExp, zSig0, zSig1 );}#endif