Mein kleines Meister-Stück

[Bonita Montero]

Ich wollte einfach mal die C-Funktion fmod() zur Berechnung des Divi-
sions-Rests zweier doubles möglichst effizient in C++20 nachprogram-
mieren. Ich hab das so hingekriegt, dass ca. 38% weniger CPU-Zeit
benötigt als die glibc-Funktion.
Grundsätzlich kann man nicht einfach hingehen und die Formel
counter - trunc( counter / denominator ) * denominator nehmen,
denn das ist zu un-präzise weil die Division ein Ergebnis haben kann
das auf der ersten Vorkomma-Stelle gerundet wurde, weil niedrigwertige
Vorkommastellen rechts von der Mantisse wegfallen oder weil bei er
Multiplikation mit dem Nenner rechts wiederum Bits rechts von der
Mantisse rausfallen.
Das Besondere an fmod() ist eben, dass man das durchführt wie eine
schriftliche Division, aber man merkt sich immer nur den Divisions
-Rest. Der ist über alle Rechen-Schritte hinweg 53 Bits breit, d.h.
das Ergebnis ist absolut exakt entsprechend den Parametern ! Ist
doch irgendwie schräg, dass alle Standard FPU Operationen das inexact
Bit setzen können wenn rechts von der Mantisse Einsen wegfallen, das
aber so eine vergleichsweise komplexe Operation eben exakt ist.

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 ); };
auto invalid = [&]( uint64_t b ) -> double { feraiseexcept( FE_INVALID
); return dbl( b ); };
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( (bDivisor & EXP_MASK) == EXP_MASK ) [[unlikely]]
// +/-[Inf|QNaN|SNaN] % ... = -QNaN
// follow SSE/AVX-rules, first NaN rules
return invalid( SIGN_BIT | bDivisor | QNAN_BIT );
if( (bDividend & EXP_MASK) == EXP_MASK ) [[unlikely]]
// +/-x % +/-[Inf|QNan|SNaN]
if( !(bDividend & MANT_MASK) ) [[likely]]
// +/-x % +/-Inf = -/+x
return dbl( SIGN_BIT | bDivisor & ~SIGN_BIT );
else
// +/-x % +/-[QNaN|SNaN] = -NaN
return invalid( SIGN_BIT | bDividend | QNAN_BIT );
int
divisorExp = (bDivisor & EXP_MASK) >> 52,
dividendExp = (bDividend & EXP_MASK) >> 52;
uint64_t
divisorMant = (divisorExp ? (uint64_t)1 : 0) << 52 | bDivisor & MANT_MASK,
dividendMant = (dividendExp ? (uint64_t)1 : 0) << 52 | bDividend &
MANT_MASK;
auto normalize = []( uint64_t &v, int &exp )
{
unsigned bits;
#if defined(HAS_FIND_FIRST_ONE)
bits = countl_zero( v ) - 11;
v <<= bits;
#else
for( bits = 0; !(v & IMPLCIT_BIT); v <<= 1, ++bits );
#endif
exp -= bits;
};
if( !divisorExp ) [[unlikely]]
// divisor is denormal
if( !divisorMant ) [[likely]]
// divisor == +/-0.0
if( dividendMant ) [[likely]]
// +/-0.0 % +/-x = -/+0.0
return dbl( sign );
else
// +/-0.0 % +/-0.0 = -QNaN
return invalid( SIGN_BIT | EXP_MASK | QNAN_BIT );
else
// normalize divisor
normalize( divisorMant, divisorExp ),
++divisorExp;
if( !dividendExp ) [[unlikely]]
// dividend is denormal
if( !dividendMant ) [[likely]]
// +/-x % +/-0.0 = -/+QNaN
return invalid( SIGN_BIT | EXP_MASK | QNAN_BIT );
else
// normalize dividend
normalize( dividendMant, dividendExp ),
++dividendExp;
int exp = divisorExp;
uint64_t remainderMant = divisorMant;
for( ; ; )
{
int mBelow = (remainderMant < dividendMant ? -1 : 0);
if( exp + mBelow < dividendExp ) [[unlikely]]
break;
exp += mBelow;
remainderMant <<= -mBelow;
if( !(remainderMant -= dividendMant) ) [[unlikely]]
{
exp = 0;
break;
}
normalize( remainderMant, exp );
};
if( exp <= 0 ) [[unlikely]]
// denormal result
remainderMant >>= -exp + 1,
exp = 0;
return dbl( sign | (uint64_t)exp << 52 | remainderMant & MANT_MASK );
}

Ich hab da ein wenig getrickst: das mBelow-Flag ist entweder -1 oder
Null, entsprechend dem Ergebnis des Vergleichs. CPUs mit CPU-Flags
können das recht gut abbilden, denn auf x86 z.B. ist das dann einfach
ein "sbb reg, reg", also eine Subtraktion eines Registers von sich
selbst minus dem Carry-Flag, was ja gesetzt wird wenn der vorangegan-
gene Verlgeich kleiner ausfällt.
Auf CPUs wie RISC-V gibt es keine CPU-Flags, aber ich hab mir das
auf godbolt.org mal im Disassembly angeschaut und man kann zumindest
ein Register in Abhängigkeit eines Vergleichs auf Null oder Eins set-
zen, was der dann auch macht und das Ergebnis negiert, dass man auf
so einer Krüppels-Architektur eben an der Stelle auch keinen beding-
ten Sprung hat der fehl-vorhergesagt werden könnte.

Ich glaub ich schreibt das nochmal auf C um und geb das den glibc
-Maintainern.