kernels/linux-2.4.0/arch/arm/nwfpe/softfloat.c
/*
-------------------------------------------------------------------------------
Returns the square root of the extended double-precision floating-point
value `a'.  The operation is performed according to the IEC/IEEE Standard
for Binary Floating-point Arithmetic.
-------------------------------------------------------------------------------
*/
Returns the square root of the extended double-precision floating-point
value `a'.  The operation is performed according to the IEC/IEEE Standard
for Binary Floating-point Arithmetic.
Returns the square root of the extended double - precision floating - point value ` a ' . The operation is performed according to the IEC / IEEE Standard for Binary Floating - point Arithmetic .
floatx80 floatx80_sqrt( floatx80 a )
{
    flag aSign;
    int32 aExp, zExp;
    bits64 aSig0, aSig1, zSig0, zSig1;
    bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
    bits64 shiftedRem0, shiftedRem1;
    floatx80 z;
    aSig0 = extractFloatx80Frac( a );
    aExp = extractFloatx80Exp( a );
    aSign = extractFloatx80Sign( a );
    if ( aExp == 0x7FFF ) {
        if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a );
        if ( ! aSign ) return a;
        goto invalid;
    }
    if ( aSign ) {
        if ( ( aExp | aSig0 ) == 0 ) return a;
 invalid:
        float_raise( float_flag_invalid );
        z.low = floatx80_default_nan_low;
        z.high = floatx80_default_nan_high;
        return z;
    }
    if ( aExp == 0 ) {
        if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
        normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
    }
    zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
    zSig0 = estimateSqrt32( aExp, aSig0>>32 );
    zSig0 <<= 31;
    aSig1 = 0;
    shift128Right( aSig0, 0, ( aExp & 1 ) + 2, &aSig0, &aSig1 );
    zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0 ) + zSig0 + 4;
    if ( 0 <= (sbits64) zSig0 ) zSig0 = LIT64( 0xFFFFFFFFFFFFFFFF );
    shortShift128Left( aSig0, aSig1, 2, &aSig0, &aSig1 );
    mul64To128( zSig0, zSig0, &term0, &term1 );
    sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
    while ( (sbits64) rem0 < 0 ) {
        --zSig0;
        shortShift128Left( 0, zSig0, 1, &term0, &term1 );
        term1 |= 1;
        add128( rem0, rem1, term0, term1, &rem0, &rem1 );
    }
    shortShift128Left( rem0, rem1, 63, &shiftedRem0, &shiftedRem1 );
    zSig1 = estimateDiv128To64( shiftedRem0, shiftedRem1, zSig0 );
    if ( (bits64) ( zSig1<<1 ) <= 10 ) {
        if ( zSig1 == 0 ) zSig1 = 1;
        mul64To128( zSig0, zSig1, &term1, &term2 );
        shortShift128Left( term1, term2, 1, &term1, &term2 );
        sub128( rem1, 0, term1, term2, &rem1, &rem2 );
        mul64To128( zSig1, zSig1, &term2, &term3 );
        sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
        while ( (sbits64) rem1 < 0 ) {
            --zSig1;
            shortShift192Left( 0, zSig0, zSig1, 1, &term1, &term2, &term3 );
            term3 |= 1;
            add192(
                rem1, rem2, rem3, term1, term2, term3, &rem1, &rem2, &rem3 );
        }
        zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
    }
    return
        roundAndPackFloatx80(
            floatx80_rounding_precision, 0, zExp, zSig0, zSig1 );
}
floatx80 floatx80_sqrt ( floatx80 a ) { flag aSign ; int32 aExp , zExp ; bits64 aSig0 , aSig1 , zSig0 , zSig1 ; bits64 rem0 , rem1 , rem2 , rem3 , term0 , term1 , term2 , term3 ; bits64 shiftedRem0 , shiftedRem1 ; floatx80 z ; aSig0 = extractFloatx80Frac ( a ) ; aExp = extractFloatx80Exp ( a ) ; aSign = extractFloatx80Sign ( a ) ; if ( aExp == 0x7FFF ) { if ( ( bits64 ) ( aSig0 << 1 ) ) return propagateFloatx80NaN ( a , a ) ; if ( ! aSign ) return a ; goto invalid ; } if ( aSign ) { if ( ( aExp | aSig0 ) == 0 ) return a ; invalid : float_raise ( float_flag_invalid ) ; z . low = floatx80_default_nan_low ; z . high = floatx80_default_nan_high ; return z ; } if ( aExp == 0 ) { if ( aSig0 == 0 ) return packFloatx80 ( 0 , 0 , 0 ) ; normalizeFloatx80Subnormal ( aSig0 , & aExp , & aSig0 ) ; } zExp = ( ( aExp - 0x3FFF ) >> 1 ) + 0x3FFF ; zSig0 = estimateSqrt32 ( aExp , aSig0 >> 32 ) ; zSig0 <<= 31 ; aSig1 = 0 ; shift128Right ( aSig0 , 0 , ( aExp & 1 ) + 2 , & aSig0 , & aSig1 ) ; zSig0 = estimateDiv128To64 ( aSig0 , aSig1 , zSig0 ) + zSig0 + 4 ; if ( 0 <= ( sbits64 ) zSig0 ) zSig0 = LIT64 ( 0xFFFFFFFFFFFFFFFF ) ; shortShift128Left ( aSig0 , aSig1 , 2 , & aSig0 , & aSig1 ) ; mul64To128 ( zSig0 , zSig0 , & term0 , & term1 ) ; sub128 ( aSig0 , aSig1 , term0 , term1 , & rem0 , & rem1 ) ; while ( ( sbits64 ) rem0 < 0 ) { -- zSig0 ; shortShift128Left ( 0 , zSig0 , 1 , & term0 , & term1 ) ; term1 |= 1 ; add128 ( rem0 , rem1 , term0 , term1 , & rem0 , & rem1 ) ; } shortShift128Left ( rem0 , rem1 , 63 , & shiftedRem0 , & shiftedRem1 ) ; zSig1 = estimateDiv128To64 ( shiftedRem0 , shiftedRem1 , zSig0 ) ; if ( ( bits64 ) ( zSig1 << 1 ) <= 10 ) { if ( zSig1 == 0 ) zSig1 = 1 ; mul64To128 ( zSig0 , zSig1 , & term1 , & term2 ) ; shortShift128Left ( term1 , term2 , 1 , & term1 , & term2 ) ; sub128 ( rem1 , 0 , term1 , term2 , & rem1 , & rem2 ) ; mul64To128 ( zSig1 , zSig1 , & term2 , & term3 ) ; sub192 ( rem1 , rem2 , 0 , 0 , term2 , term3 , & rem1 , & rem2 , & rem3 ) ; while ( ( sbits64 ) rem1 < 0 ) { -- zSig1 ; shortShift192Left ( 0 , zSig0 , zSig1 , 1 , & term1 , & term2 , & term3 ) ; term3 |= 1 ; add192 ( rem1 , rem2 , rem3 , term1 , term2 , term3 , & rem1 , & rem2 , & rem3 ) ; } zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 ) ; } return roundAndPackFloatx80 ( floatx80_rounding_precision , 0 , zExp , zSig0 , zSig1 ) ; }
Returns the square root of the extended double-precision floating-point
value `a'.
Returns the square root of the extended double - precision floating - point value ` a ' .