if ( aExp == 0 ) {        ++expDiff;    }    else {        aSig |= 0x40000000;    }    shift32RightJamming( aSig, - expDiff, &aSig );    bSig |= 0x40000000; bBigger:    zSig = bSig - aSig;    zExp = bExp;    zSign ^= 1;    goto normalizeRoundAndPack; aExpBigger:    if ( aExp == 0xFF ) {        if ( aSig ) return propagateFloat32NaN( a, b );        return a;    }    if ( bExp == 0 ) {        --expDiff;    }    else {        bSig |= 0x40000000;    }    shift32RightJamming( bSig, expDiff, &bSig );    aSig |= 0x40000000; aBigger:    zSig = aSig - bSig;    zExp = aExp; normalizeRoundAndPack:    --zExp;    return normalizeRoundAndPackFloat32( zSign, zExp, zSig );}/*-------------------------------------------------------------------------------Returns the result of adding the single-precision floating-point values `a'and `b'.  The operation is performed according to the IEC/IEEE Standard forBinary Floating-point Arithmetic.-------------------------------------------------------------------------------*/float32 float32_add( float32 a, float32 b ){    flag aSign, bSign;    aSign = extractFloat32Sign( a );    bSign = extractFloat32Sign( b );    if ( aSign == bSign ) {        return addFloat32Sigs( a, b, aSign );    }    else {        return subFloat32Sigs( a, b, aSign );    }}/*-------------------------------------------------------------------------------Returns the result of subtracting the single-precision floating-point values`a' and `b'.  The operation is performed according to the IEC/IEEE Standardfor Binary Floating-point Arithmetic.-------------------------------------------------------------------------------*/float32 float32_sub( float32 a, float32 b ){    flag aSign, bSign;    aSign = extractFloat32Sign( a );    bSign = extractFloat32Sign( b );    if ( aSign == bSign ) {        return subFloat32Sigs( a, b, aSign );    }    else {        return addFloat32Sigs( a, b, aSign );    }}/*-------------------------------------------------------------------------------Returns the result of multiplying the single-precision floating-point values`a' and `b'.  The operation is performed according to the IEC/IEEE Standardfor Binary Floating-point Arithmetic.-------------------------------------------------------------------------------*/float32 float32_mul( float32 a, float32 b ){    flag aSign, bSign, zSign;    int16 aExp, bExp, zExp;    bits32 aSig, bSig;    bits64 zSig64;    bits32 zSig;    aSig = extractFloat32Frac( a );    aExp = extractFloat32Exp( a );    aSign = extractFloat32Sign( a );    bSig = extractFloat32Frac( b );    bExp = extractFloat32Exp( b );    bSign = extractFloat32Sign( b );    zSign = aSign ^ bSign;    if ( aExp == 0xFF ) {        if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {            return propagateFloat32NaN( a, b );        }        if ( ( bExp | bSig ) == 0 ) {            float_raise( float_flag_invalid );            return float32_default_nan;        }        return packFloat32( zSign, 0xFF, 0 );    }    if ( bExp == 0xFF ) {        if ( bSig ) return propagateFloat32NaN( a, b );        if ( ( aExp | aSig ) == 0 ) {            float_raise( float_flag_invalid );            return float32_default_nan;        }        return packFloat32( zSign, 0xFF, 0 );    }    if ( aExp == 0 ) {        if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );        normalizeFloat32Subnormal( aSig, &aExp, &aSig );    }    if ( bExp == 0 ) {        if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );        normalizeFloat32Subnormal( bSig, &bExp, &bSig );    }    zExp = aExp + bExp - 0x7F;    aSig = ( aSig | 0x00800000 )<<7;    bSig = ( bSig | 0x00800000 )<<8;    shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );    zSig = zSig64;    if ( 0 <= (sbits32) ( zSig<<1 ) ) {        zSig <<= 1;        --zExp;    }    return roundAndPackFloat32( zSign, zExp, zSig );}/*-------------------------------------------------------------------------------Returns the result of dividing the single-precision floating-point value `a'by the corresponding value `b'.  The operation is performed according to theIEC/IEEE Standard for Binary Floating-point Arithmetic.-------------------------------------------------------------------------------*/float32 float32_div( float32 a, float32 b ){    flag aSign, bSign, zSign;    int16 aExp, bExp, zExp;    bits32 aSig, bSig, zSig;    aSig = extractFloat32Frac( a );    aExp = extractFloat32Exp( a );    aSign = extractFloat32Sign( a );    bSig = extractFloat32Frac( b );    bExp = extractFloat32Exp( b );    bSign = extractFloat32Sign( b );    zSign = aSign ^ bSign;    if ( aExp == 0xFF ) {        if ( aSig ) return propagateFloat32NaN( a, b );        if ( bExp == 0xFF ) {            if ( bSig ) return propagateFloat32NaN( a, b );            float_raise( float_flag_invalid );            return float32_default_nan;        }        return packFloat32( zSign, 0xFF, 0 );    }    if ( bExp == 0xFF ) {        if ( bSig ) return propagateFloat32NaN( a, b );        return packFloat32( zSign, 0, 0 );    }    if ( bExp == 0 ) {        if ( bSig == 0 ) {            if ( ( aExp | aSig ) == 0 ) {                float_raise( float_flag_invalid );                return float32_default_nan;            }            float_raise( float_flag_divbyzero );            return packFloat32( zSign, 0xFF, 0 );        }        normalizeFloat32Subnormal( bSig, &bExp, &bSig );    }    if ( aExp == 0 ) {        if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );        normalizeFloat32Subnormal( aSig, &aExp, &aSig );    }    zExp = aExp - bExp + 0x7D;    aSig = ( aSig | 0x00800000 )<<7;    bSig = ( bSig | 0x00800000 )<<8;    if ( bSig <= ( aSig + aSig ) ) {        aSig >>= 1;        ++zExp;    }    zSig = ( ( (bits64) aSig )<<32 ) / bSig;    if ( ( zSig & 0x3F ) == 0 ) {        zSig |= ( ( (bits64) bSig ) * zSig != ( (bits64) aSig )<<32 );    }    return roundAndPackFloat32( zSign, zExp, zSig );}/*-------------------------------------------------------------------------------Returns the remainder of the single-precision floating-point value `a'with respect to the corresponding value `b'.  The operation is performedaccording to the IEC/IEEE Standard for Binary Floating-point Arithmetic.-------------------------------------------------------------------------------*/float32 float32_rem( float32 a, float32 b ){    flag aSign, bSign, zSign;    int16 aExp, bExp, expDiff;    bits32 aSig, bSig;    bits32 q;    bits64 aSig64, bSig64, q64;    bits32 alternateASig;    sbits32 sigMean;    aSig = extractFloat32Frac( a );    aExp = extractFloat32Exp( a );    aSign = extractFloat32Sign( a );    bSig = extractFloat32Frac( b );    bExp = extractFloat32Exp( b );    bSign = extractFloat32Sign( b );    if ( aExp == 0xFF ) {        if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {            return propagateFloat32NaN( a, b );        }        float_raise( float_flag_invalid );        return float32_default_nan;    }    if ( bExp == 0xFF ) {        if ( bSig ) return propagateFloat32NaN( a, b );        return a;    }    if ( bExp == 0 ) {        if ( bSig == 0 ) {            float_raise( float_flag_invalid );            return float32_default_nan;        }        normalizeFloat32Subnormal( bSig, &bExp, &bSig );    }    if ( aExp == 0 ) {        if ( aSig == 0 ) return a;        normalizeFloat32Subnormal( aSig, &aExp, &aSig );    }    expDiff = aExp - bExp;    aSig |= 0x00800000;    bSig |= 0x00800000;    if ( expDiff < 32 ) {        aSig <<= 8;        bSig <<= 8;        if ( expDiff < 0 ) {            if ( expDiff < -1 ) return a;            aSig >>= 1;        }        q = ( bSig <= aSig );        if ( q ) aSig -= bSig;        if ( 0 < expDiff ) {            q = ( ( (bits64) aSig )<<32 ) / bSig;            q >>= 32 - expDiff;            bSig >>= 2;            aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;        }        else {            aSig >>= 2;            bSig >>= 2;        }    }    else {        if ( bSig <= aSig ) aSig -= bSig;        aSig64 = ( (bits64) aSig )<<40;        bSig64 = ( (bits64) bSig )<<40;        expDiff -= 64;        while ( 0 < expDiff ) {            q64 = estimateDiv128To64( aSig64, 0, bSig64 );            q64 = ( 2 < q64 ) ? q64 - 2 : 0;            aSig64 = - ( ( bSig * q64 )<<38 );            expDiff -= 62;        }        expDiff += 64;        q64 = estimateDiv128To64( aSig64, 0, bSig64 );        q64 = ( 2 < q64 ) ? q64 - 2 : 0;        q = q64>>( 64 - expDiff );        bSig <<= 6;        aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;    }    do {        alternateASig = aSig;        ++q;        aSig -= bSig;    } while ( 0 <= (sbits32) aSig );    sigMean = aSig + alternateASig;    if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {        aSig = alternateASig;    }    zSign = ( (sbits32) aSig < 0 );    if ( zSign ) aSig = - aSig;    return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig );}/*-------------------------------------------------------------------------------Returns the square root of the single-precision floating-point value `a'.The operation is performed according to the IEC/IEEE Standard for BinaryFloating-point Arithmetic.-------------------------------------------------------------------------------*/float32 float32_sqrt( float32 a ){    flag aSign;    int16 aExp, zExp;    bits32 aSig, zSig;    bits64 rem, term;    aSig = extractFloat32Frac( a );    aExp = extractFloat32Exp( a );    aSign = extractFloat32Sign( a );    if ( aExp == 0xFF ) {        if ( aSig ) return propagateFloat32NaN( a, 0 );        if ( ! aSign ) return a;        float_raise( float_flag_invalid );        return float32_default_nan;    }    if ( aSign ) {        if ( ( aExp | aSig ) == 0 ) return a;        float_raise( float_flag_invalid );        return float32_default_nan;    }    if ( aExp == 0 ) {        if ( aSig == 0 ) return 0;        normalizeFloat32Subnormal( aSig, &aExp, &aSig );    }    zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;    aSig = ( aSig | 0x00800000 )<<8;    zSig = estimateSqrt32( aExp, aSig ) + 2;    if ( ( zSig & 0x7F ) <= 5 ) {        if ( zSig < 2 ) {            zSig = 0xFFFFFFFF;        }        else {            aSig >>= aExp & 1;            term = ( (bits64) zSig ) * zSig;            rem = ( ( (bits64) aSig )<<32 ) - term;            while ( (sbits64) rem < 0 ) {                --zSig;                rem += ( ( (bits64) zSig )<<1 ) | 1;            }            zSig |= ( rem != 0 );        }    }    shift32RightJamming( zSig, 1, &zSig );    return roundAndPackFloat32( 0, zExp, zSig );}/*-------------------------------------------------------------------------------Returns 1 if the single-precision floating-point value `a' is equal to thecorresponding value `b', and 0 otherwise.  The comparison is performedaccording to the IEC/IEEE Standard for Binary Floating-point Arithmetic.-------------------------------------------------------------------------------*/flag float32_eq( float32 a, float32 b ){    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )       ) {        if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {            float_raise( float_flag_invalid );        }        return 0;    }    return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );