in nwfpe/softfloat.c [1347:1441]
float32 float32_rem( struct roundingData *roundData, 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 );
}
roundData->exception |= float_flag_invalid;
return float32_default_nan;
}
if ( bExp == 0xFF ) {
if ( bSig ) return propagateFloat32NaN( a, b );
return a;
}
if ( bExp == 0 ) {
if ( bSig == 0 ) {
roundData->exception |= 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 ) {
bits64 tmp = ( (bits64) aSig )<<32;
do_div( tmp, bSig );
q = tmp;
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( roundData, aSign ^ zSign, bExp, aSig );
}