Ich hab mal in die Tabellen auf agner.org geguckt und dabei gesehen,
dass die x87-Instruktion FPREM ziemlich schnell ist, nämlich ca. 70
Takte auf meiner Zen4-CPU. Das mag auf Anhieb langsam wirken, aber
hinter der Operation steckt echt mehr Komplexität als bei einer FP-
Division, denn die bricht nach 53 Ergebnis-Bits ab, während das da-
hinterstehende Subtract and Shift bei einem Divisionsrest so lange
rechnet bis der Exponent des Rests kleiner als der des Divisors ist
oder die Mantisse kleiner als die des Divisors ist.
Wenn man mal in die glibc schaut, dann sieht man, dass der ganze
Krempel nicht über die x87-FPU läuft, sondern über Integer-Operati-
onen. Das ist relativ trickreich, aber eben bei weitem nicht so
schnell wie die x87-FPU angeblich sein soll.
Also habe ich mal mit VC++ ein C++-Programm geschrieben, das fmod()
mit einer kleinen MASM-Assembler-Routine vergleicht die FPREM um-
setzt.

#include <iostream>
#include <random>
#include <cstdint>
#include <cmath>
#include <bit>

extern double x87Fprem( double x, double y );

using namespace std;

int main()
{ mt19937_64 mt;
uniform_int_distribution<uint64_t> uid( 0, -1 );
constexpr size_t N = 10'000'000;
double accuracySum = 0;
for( size_t r = N; r; )
{
uint64_t
aB = uid( mt ),
bB = uid( mt );
auto chk = []( uint64_t b ) { return (b >> 52 & 0x7FF) != 0x7FF; };
if( !chk( aB ) || !chk( bB ) )
continue;
double
a = bit_cast<double>( aB ),
b = bit_cast<double>( bB ),
std = fmod( a, b ),
x87 = x87Fprem( a, b );
if( isnan( a ) || isnan( b ) || signbit( a ) != signbit( b ) && a != b )
continue;
std = abs( std );
x87 = abs( x87 );
accuracySum += x87 != std ? log2( x87 > std ? x87 : std ) - log2( abs(
x87 - std ) ) : 53;
--r;
}
cout << accuracySum / 1.0e7 << endl;
return 0;
}

PUBLIC ?x87Fprem@@YANNN@Z

_TEXT SEGMENT
?x87Fprem@@YANNN@Z PROC
movsd qword ptr [rsp - 16], xmm0
movsd qword ptr [rsp - 8], xmm1
fld qword ptr [rsp - 8]
fld qword ptr [rsp - 16]
fprem
fstp qword ptr [rsp -8]
fstp st
movsd xmm0, qword ptr [rsp - 8]
ret
?x87Fprem@@YANNN@Z ENDP
_TEXT ENDS

END

Der Assembler-Code mag ein wenig umständlich sein, aber ich krieg die
Parameter für x und y nicht über die x87 FPU-Register übergeben. Die
Logarithmen-Trickserei oben gibt mir an, wieivel Bits Genauigkeit die
beiden Ergebnisse teilen. Das Ergebnis der Runtime hat bei VC++ und
g++ nie Präzisionsverlust da bei den einzelnen Subtract and Shift
Zwischenergebnissen nie Mantissen-Bits fallengelassen werden müssen.
Da hat es mich schon ziemlich überrascht, dass FPREM im Schnitt nur
ca. 28 Bits Genauigkeit mit fmod() teilt. Das ist ja wirklich unter
aller Sau, was die x87-FPU da macht.