A fmod() that is faster than MSVC's, glibc's fmod() in C++20

[Bonita Montero]

double myFmod( double divisor, double dividend )
{
auto bin = []( double d ) -> uint64_t { return bit_cast<uint64_t>( d ); };
auto dbl = []( uint64_t u ) -> double { return bit_cast<double>( u ); };
uint64_t
bDivisor = bin( divisor ),
bDividend = bin( dividend );
constexpr uint64_t
SIGN_BIT = (uint64_t)1 << 63,
EXP_MASK = (uint64_t)0x7FF << 52,
IMPLCIT_BIT = (uint64_t)1 << 52,
MANT_MASK = IMPLCIT_BIT - 1,
QNAN_BIT = IMPLCIT_BIT >> 1;
uint64_t sign = (bDivisor ^ bDividend) & SIGN_BIT;
if( (bDividend & EXP_MASK) == EXP_MASK ) [[unlikely]]
// ... % +/-[Inf|QNan|SNaN]
if( !(bDividend & MANT_MASK) ) [[likely]]
// ... % +/-Inf
if( (bDivisor & EXP_MASK) != EXP_MASK ) [[unlikely]]
// +/-x % +/-Inf = -/+x
return dbl( sign | bDivisor & ~SIGN_BIT );
else if( !(bDivisor & MANT_MASK) ) [[likely]]
// +/-Inf % +/-Inf = -QNaN
return dbl( SIGN_BIT | EXP_MASK | QNAN_BIT );
else
{
// +/-[QNaN|SNaN] % +/-Inf = -[QNaN|SNaN]
if( !(bDivisor & QNAN_BIT) )
feraiseexcept( FE_INVALID );
return dbl( SIGN_BIT | bDivisor );
}
else if( bDividend & QNAN_BIT ) [[likely]]
// ... % +/-QNaN
if( (bDivisor & EXP_MASK) != EXP_MASK ) [[likely]]
// +/-x % +/-QNaN = -QNaN
return dbl( SIGN_BIT | bDivisor | bDividend );
else if( !(bDivisor & MANT_MASK) ) [[likely]]
// +/-Inf % +/-QNaN = -QNaN
return dbl( SIGN_BIT | bDivisor | bDividend );
else
{
// +/-[QNaN|SNaN] % +/-QNaN = -[QNaN|SNaN]
if( ~bDivisor & QNAN_BIT ) [[unlikely]]
feraiseexcept( FE_INVALID );
return dbl( SIGN_BIT | (bDivisor | bDividend) & ~(~bDivisor &
QNAN_BIT) );
}
else
{
// ... % +/-SNaN = -SNaN
feraiseexcept( FE_INVALID );
return dbl( (SIGN_BIT | bDivisor | bDividend) & ~QNAN_BIT );
}
if( (bDivisor & EXP_MASK) == EXP_MASK ) [[unlikely]]
// +/-[Inf|QNan|SNaN] / +/-x
if( !(bDivisor & MANT_MASK) ) [[likely]]
// +/-Inf % +/-x = -QNaN
return dbl( SIGN_BIT | EXP_MASK | QNAN_BIT );
else if( bDivisor & QNAN_BIT ) [[likely]]
// +/-QNaN % +/-x = -QNaN
return dbl( SIGN_BIT | bDivisor | bDividend );
else
{
// +/-SNaN % +/-x = -SNaN
feraiseexcept( FE_INVALID );
return dbl( SIGN_BIT | (bDivisor | bDividend) & ~QNAN_BIT );
}
int
divisorExp = (bDivisor & EXP_MASK) >> 52,
dividendExp = (bDividend & EXP_MASK) >> 52;
uint64_t
divisorMant = IMPLCIT_BIT | bDivisor & MANT_MASK,
dividendMant = IMPLCIT_BIT | bDividend & MANT_MASK;
auto normalize = []( uint64_t &mant, int &exp )
{
mant &= MANT_MASK;
unsigned shift = countl_zero( mant ) - 11;
mant <<= shift;
exp -= shift - 1;
};
if( !divisorExp ) [[unlikely]]
if( !(bDivisor & MANT_MASK) ) [[likely]]
// +/-0.0 % +/-x = -/+0.0
return dbl( sign );
else
// divisor is +/-denorm -> normalize
normalize( divisorMant, divisorExp );
if( !dividendExp ) [[unlikely]]
if( !(bDividend & MANT_MASK) ) [[likely]]
// +/-x % +/-0.0 = -/+QNaN
return dbl( sign | EXP_MASK | QNAN_BIT );
else
// dividend is +/-denorm -> normalize
normalize( dividendMant, dividendExp );
int exp = divisorExp;
uint64_t remainderMant = divisorMant;
for( ; ; )
{
int below = remainderMant < dividendMant;
if( exp - below < dividendExp ) [[unlikely]]
break;
exp -= below;
remainderMant <<= below;
if( !(remainderMant -= dividendMant) ) [[unlikely]]
return dbl( sign );
unsigned shift = (unsigned)countl_zero( remainderMant ) - 11;
remainderMant <<= shift;
exp -= shift;
}
if( exp <= 0 ) [[unlikely]]
// denormal result
remainderMant >>= -exp + 1,
exp = 0;
return dbl( sign | (uint64_t)exp << 52 | remainderMant & MANT_MASK );
}