fmod(), schneller als jede andere Implementation

[Bonita Montero]

Ich habe mal fmod( double, double ) in C++20 implementiert.
Das ganze ist auf meiner CPU ca. 2 x so schnell wie die Implementation
in MSVC und 44% schneller als die Implementation in der glibc.

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 );
}

Ich würde das gerne mal komplett in C implementieren.
Dazu bräuchte ich aber die Funktion countl_zero(), die die Anzahl
der führenden (höherwertigen) Null-Bits in einem Wert zählt. Gibt's
das ähnlich in irgendeinem C-Standard, oder muss ich da ein Intrinsic
nutzen ?
In C++20 mappt dann auf entsprechende CPU-Instruktionen moderner CPUS.
BSR gibt's seit Ewigkeiten in x86, kann aber nicht mit Null-Werten um-
gehen, also da wo BSR dann im Endeffekt bei einem 64 Bit Wert 64 ergeben
würde, weswegen der Compiler da ein Sonderfall-Handling einbaut. LZCNT
ist eine Instruktion die es seit der Core-i-Serie der Haswells gibt und
dieses Problem nicht hat.

[Bonita Montero]

Ich hab mal die alte x87-Funktion FPREM getestet, die ja auch den Divi-
sions-Rest berechnet.
Helmut hält ja von diesen x87-Funktionen einiges, bzw. der meint das sei
alles sehr viel schneller als wenn man das händisch implementiert. Ich
hab ja bereits für FSIN nachgewiesen, dass das eben nicht so ist. FPREM
ist in der Tat extrem schnell bzw. das habe ich mal mit einem kleinen
Benchmark getestet, liegt auf jeden Fall auf meiner CPU immer unter 40
Takten.
Aber FPREM ist halt nicht wie in der Intel-Dokumentation angegeben prä-
zise, bzw. das Ergebnis liegt oft außerhalb des Werte-Bereichs auf mei-
ner AMD Zen2-CPU. Wie kann man bitte so eine Scheiße bauen ?
Das C-Programm unten testet das, vielleicht kann das mal jemand für
GCC Inline Assembler anpassen.

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

void test( double (*fn)( double, double ) );
double x87( double counter, double denom );

int main()
{
test( x87 );
test( fmod );
}

uint64_t bin( double d );
double dbl( uint64_t u );

void test( double (*fn)( double, double ) )
{
double const
COUNTER = dbl( 0x7FEFFFFFFFFFFFFFu ),
DENOM = dbl( 0x0001000000000001u );
printf( "%016llx\n", bin( fn( COUNTER, DENOM ) ) );
}

double x87( double counter, double denom )
{
double result;
__asm
{
fld counter
fld denom
fprem
fstp result
fstp st(0)
}
return result;
}

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

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

Auf jeden Fall gibt dieses Programm auf meiner CPU ...
0001000000000001
0000fffbe0000001
.. aus, d.h. das Ergebnis müsste ja >= 0.0 < Dividend sein, aber
die x87-FPU, bzw. die Implementation mit meinem Zen2-Dingens, gibt
als Ergebnis den Dividenden aus, also exakt an der oberen Schranke
und nicht kleiner. Was ist das bitte für ein Müll ?
Aber grundsätzlich: so schnell wie FPREM ist kann das auch nicht
exakt sein. Man sieht ja an meinem Beispiel im vorangegangen Posting,
dass diese schriftliche "Division" die ich da mache so lange anhält
bis der Exponent des Rests kleiner ist oder der Exponent des Rests
gleich ist und die Mantisse kleiner. Theoretisch müsste man das ja
auch mit teilweise Hardware-Untestützung und teilwise im Microcode
effizient lösen können, aber für mein Gefühl gehört sowas aufwendi-
ges nicht in eine einzelne CPU-Instruktion.

[Bonita Montero]

Ich hab das Ganze hier nochmal in reinem C nachprogrammiert. Auf jeden
Fall hab ich das mal jetzt mit gleich-verteilten Exponenten für den
Zähler und den Nenner ge-benchmarkt und da ist mein Code je nach CPU
auf der ich das getestet habe zwei bis dreimal schneller (2.2 - 3.2)
als der der Gnu C Libary.
Alledings braucht das Ganze dann auch spezielle CPU-Unterstützung in
Form einer Instruktion die die führenden, also höherwertigen Nullen
einer Ganzzahl zählt. Bei x86 gab es das in Form von BSR schon immer,
aber seit der Core i Haswell Serie oder entsprechenden AMD CPUs gibt's
noch die flexiblere Variante LZCNT.

#include <stdint.h>
#include <string.h>
#include <fenv.h>
#if defined(_MSC_VER)
#include <intrin.h>
#endif

#if defined(__GNUC__) || defined(__clang__)
#define likely(x) __builtin_expect((x), 1)
#define unlikely(x) __builtin_expect((x), 0)
#else
#define likely(x) (x)
#define unlikely(x) (x)
#endif

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

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

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

inline double NaN()
{
feraiseexcept( FE_INVALID );
return dbl( SIGN_BIT | EXP_MASK | QNAN_BIT );
}

inline void normalize( uint64_t *mant, int *exp )
{
#if defined(__GNUC__) || defined(__clang__)
unsigned bits = __builtin_clz( *mant ) - 11;
*mant <<= bits;
#elif defined(_MSC_VER)
unsigned long bits;
_BitScanReverse64( &bits, *mant );
bits = (bits ^ 63) - 11;
*mant <<= bits;
#else
unsigned bits = 0;
for( ; !(*mant & IMPLCIT_BIT); *mant <<= 1, ++bits );
#endif
*exp -= bits;
}

double myFmodC( double counter, double denominator )
{
uint64_t const
bCounter = bin( counter ),
bDenom = bin( denominator );
uint64_t const sign = bCounter & SIGN_BIT;
if( unlikely((bCounter & EXP_MASK) == EXP_MASK) )
// +/-[Inf|QNaN|SNaN] % ... = -QNaN
return NaN();
if( unlikely((bDenom & EXP_MASK) == EXP_MASK) )
// +/-x % +/-[Inf|QNan|SNaN]
if( likely(!(bDenom & MANT_MASK)) )

// +/-x % +/-Inf = -/+x

return dbl( sign | bCounter & ~SIGN_BIT );
else
// +/-x % +/-[QNaN|SNaN] = -NaN
return NaN();
int
counterExp = bCounter >> 52 & MAX_EXP,
denomExp = bDenom >> 52 & MAX_EXP;
uint64_t
counterMant = (uint64_t)!!counterExp << 52 | bCounter & MANT_MASK,
denomMant = (uint64_t)!!denomExp << 52 | bDenom & MANT_MASK;
if( unlikely(!counterExp) )
// counter is denormal
if( likely(!counterMant) )
// counter == +/-0.0
if( likely(denomMant) )

// +/-0.0 % +/-x = -/+0.0
return dbl( sign );
else

// +/-0.0 % +/-0.0 = -QNaN
return NaN();
else
// normalize counter
normalize( &counterMant, &counterExp ),
++counterExp;
if( unlikely(!denomExp) )
// denominator is denormal
if( likely(!denomMant) )

// +/-x % +/-0.0 = -/+QNaN

return NaN();
else
// normalize denominator
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( sign | (uint64_t)remExp << 52 | remMant & MANT_MASK );
}