I re-implemented fmod() as a drop-in replacement for glibc.
It requires sth. like CMOVcc and BSR on x86 on other CPUs to be fast.
On my Zen2-CPU I get +150% more Throughput. Ony my Skylake Linuc PC
the code is +105% faster.

#include <stdint.h>
#include <string.h>

#define LIKELY(x) __builtin_expect((x), 1)
#define UNLIKELY(x) __builtin_expect((x), 0)

#define EXP_MAX (0x7FF)
#define SIGN_BIT ((uint64_t)1 << 63)
#define EXP_MASK ((uint64_t)EXP_MAX << 52)
#define IMPLCIT_BIT ((uint64_t)1 << 52)
#define MANT_MASK (IMPLCIT_BIT - 1)

static inline int hasMaxExp( uint64_t b )
{ return b >= EXP_MASK;
}

static inline int hasInfMant( uint64_t b )
{ return !(b & MANT_MASK);
}

static inline uint64_t bin( double d )
{ uint64_t u;
memcpy( &u, &d, sizeof d );
return u;
}

static inline double dbl( uint64_t u )
{ double d;
memcpy( &d, &u, sizeof u );
return d;
}

static inline void normalize( uint64_t *mant, int *exp )
{ unsigned bits = __builtin_clzll( *mant ) - 11;
*mant <<= bits;
*exp -= bits;
}

double __ieee754_fmod( double counter, double denom )
{ uint64_t
bCounter = bin( counter ),
bDenom = bin( denom ) & ~SIGN_BIT,
bSign = bCounter & SIGN_BIT;
bCounter &= ~SIGN_BIT;
if( UNLIKELY(!bDenom) || UNLIKELY(hasMaxExp( bCounter )) )
return (counter * denom) / (counter * denom);
if( UNLIKELY(hasMaxExp( bDenom )) )
if( LIKELY(hasInfMant( bDenom )) )
return counter;
else
return (counter * denom) / (counter * denom);
if( UNLIKELY(!bCounter) )
return counter;
int
counterExp = bCounter >> 52 & EXP_MAX,
denomExp = bDenom >> 52 & EXP_MAX;
uint64_t
counterMant = (uint64_t)(counterExp != 0) << 52 | bCounter & MANT_MASK,
denomMant = (uint64_t)(denomExp != 0) << 52 | bDenom & MANT_MASK;
if( UNLIKELY(!counterExp) )
normalize( &counterMant, &counterExp ),
++counterExp;
if( UNLIKELY(!denomExp) )
normalize( &denomMant, &denomExp ),
++denomExp;
int remExp = counterExp;
uint64_t remMant = counterMant;
for( ; ; )
{
int below = remMant < denomMant;
if( UNLIKELY(remExp - below < denomExp) )
break;
remExp -= below;
remMant <<= below;
if( UNLIKELY(!(remMant -= denomMant)) )
{
remExp = 0;
break;
}
normalize( &remMant, &remExp );
};
if( UNLIKELY(remExp <= 0) )
// denormal result
remMant >>= -remExp + 1,
remExp = 0;
return dbl( bSign | (uint64_t)remExp << 52 | remMant & MANT_MASK );
} libm_alias_finite(__ieee754_fmod, __fmod)