Helmut meinte ja mal, dass FSIN auf x87-FPUs schneller als das sin()
der Standardbibliothkek ist. Ich hab mir mal die Timings auf modernen
CPUs bei agner.org angeschaut und das kam mir alles ein wenig seltsam
vor bzw. falls das intern über den CORDIC-Algorithmus gehen sollte
sollte das eigentlich schneller gehen als mit bis zu 200 Takten auf
meiner Zen4-CPU.
Hab dann mal einen kleinen Benchmark in C++ und Assembler geschrieben,
und der sieht so aus:

#define _USE_MATH_DEFINES
#include <iostream>
#include <cmath>
#include <chrono>
#include <atomic>
#include <vector>

using namespace std;
using namespace chrono;

double x87Sin( double angle );

atomic<double> aSum;

int main()
{ vector<double> values;
for( int i = 0; i < 1000; ++i )
values.emplace_back( M_PI * i / 1000.0 );
auto bench = [&]( double (*sinFn)( double ) )
{
auto start = high_resolution_clock::now();
double sum = 0.0;
for( int r = 0; r < 10'000; ++r )
for( double d : values )
sum += sinFn( d );
double t = duration_cast<nanoseconds>( high_resolution_clock::now() -
start ).count() / 1.0e7;
aSum = sum;
return t;
};
cout << bench( x87Sin ) / bench( sin ) << endl;
unsigned unequal = 0;
for( double d : values )
unequal += x87Sin( d ) != sin( d );
cout << unequal << endl;
}

PUBLIC ?x87Sin@@YANN@Z

_TEXT SEGMENT
?x87Sin@@YANN@Z PROC
movsd qword ptr [rsp - 8], xmm0
fld qword ptr [rsp - 8]
fsin
fstp qword ptr [rsp - 8]
movsd xmm0, qword ptr [rsp - 8]
ret
?x87Sin@@YANN@Z ENDP
_TEXT ENDS
END

Das Ergebnis ist, dass das sin() der MSVC-Runtime ca. 11,3 mal
schneller ist als das x87 FSIN.
Was mich auch interssiert hatte ist, wie oft sich die Ergebnisse
für gleiche Eingaben aufgrund unterschiedlicher Rechenwege unter-
scheiden. Und da war ich echt überrascht, dass das sin() der Run-
time auf 1.000 Eingaben sich nur 27 mal unterschied; ich hätte
echt damit gerechnet, dass sich die meisten Endergebnisse irgendwo
auf den letzten Bits unterscheiden.

Am 26.04.24 um 17:01 schrieb Bonita Montero:
> Das Ergebnis ist, dass das sin() der MSVC-Runtime ca. 11,3 mal
> schneller ist als das x87 FSIN.

Was erwartest du, wenn jedes mal die Daten erst durch den Stack-Wolf drehst?

> Was mich auch interssiert hatte ist, wie oft sich die Ergebnisse
> für gleiche Eingaben aufgrund unterschiedlicher Rechenwege unter-
> scheiden. Und da war ich echt überrascht, dass das sin() der Run-
> time auf 1.000 Eingaben sich nur 27 mal unterschied; ich hätte
> echt damit gerechnet, dass sich die meisten Endergebnisse irgendwo
> auf den letzten Bits unterscheiden.

Probier mal große Zahlen, nicht [0,Pi), dann kommen die Unterschiede
zutage. fsin ist da eigen.

Marcel

Am 26.04.2024 um 20:37 schrieb Marcel Mueller:

> Was erwartest du, wenn jedes mal die Daten erst durch den Stack-Wolf
> drehst?

Was ?

> Probier mal große Zahlen, nicht [0,Pi), dann kommen die Unterschiede
> zutage. fsin ist da eigen.

Is aber unwahrscheinlich, dass jemand Zahlen jenseits von 360 Grad
nimmt.

Ich hab jetzt mal das FSIN als gcc Inline-Assembly gemacht. Der FSIN
-Teil ist mit dem g++ natürlich genauso "schnell", aber das std::sin()
ist sehr viel langsamer, dass das std::sin() im Endeffekt nur fünfmal
so schnell ist wie das FSIN.

#include <iostream>
#include <cmath>
#include <chrono>
#include <atomic>
#include <vector>

constexpr double PI = 3.14159265358979323846;

using namespace std;
using namespace chrono;

#if defined(_MSC_VER)
double x87Sin( double angle );
#else
double x87Sin( double angle )
{ double sinval;
asm( "fsin" : "=t" (sinval) : "0" (angle) );
return sinval;
} #endif

atomic<double> aSum;

int main()
{ vector<double> values;
values.reserve( 1'000 );
for( int i = 0; i < 1000; ++i )
values.emplace_back( PI * i / 1e3 );
auto bench = [&]( double (*sinFn)( double ) )
{
auto start = high_resolution_clock::now();
double sum = 0.0;
for( int r = 0; r < 10'000; ++r )
for( double d : values )
sum += sinFn( d );
double t = duration_cast<nanoseconds>( high_resolution_clock::now() -
start ).count() / 1.0e7;
aSum.store( sum, memory_order_relaxed );
return t;
};
double
x87 = bench( x87Sin ),
std = bench( sin );
cout << "x87: " << x87 << endl;
cout << "std: " << std << endl;
cout << "*: " << x87 / std << endl;
unsigned unequal = 0;
for( double d : values )
unequal += x87Sin( d ) != sin( d );
cout << unequal << endl;
}

Am 27.04.24 um 06:55 schrieb Bonita Montero:
> Am 26.04.2024 um 20:37 schrieb Marcel Mueller:
>
>> Was erwartest du, wenn jedes mal die Daten erst durch den Stack-Wolf
>> drehst?
>
> Was ?

Naja, der Assembler-Code verhindert jegliche Optimierungen für den
Datentransport. Und dass die zu testende Funktion nur einige Prozent der
ausgeführten Instruktionen ausmacht, macht die Sache keineswegs besser.

Es bleibt aber auch festzuhalten, dass die ruhmreichen Zeiten der x86
FPU gezählt sind. Nachdem der 486 schon langsamer als der 387 war, hat
sich das ganze nicht mehr sonderlich zum Guten weiterentwickelt. Die 586
haben nochmal Schwung gebracht (und sich Anfangs verrechnet), aber
danach kam nicht mehr viel.

>> Probier mal große Zahlen, nicht [0,Pi), dann kommen die Unterschiede
>> zutage. fsin ist da eigen.
>
> Is aber unwahrscheinlich, dass jemand Zahlen jenseits von 360 Grad
> nimmt.

Noch nicht viel mit trigonometrischen Funktionen gemacht, oder?
Bereits die Berechnung eines einfachen Sinustons bedingt größere Werte.

Die eingebaute sin-Funktion wird dann auch langsamer, weil fmod recht
teuer ist.

Mit 8-fach Unrolling komme ich bei kleinen Zahlen auf Faktor 3 langsamer
als sin(), aber bei großen Zahlen dreht sich das Bild. Dann wird sin()
immer langsamer. Das hängt natürlich von der jeweiligen CPU ab. Das war
jetzt ein alter AMD FX.

Aber zumindest scheint das Große-Zahlen-Problem mittlerweile halbwegs im
Griff zu sein. Früher hat fsin ziemlichen Blödsinn ausgerechnet, wenn
die Argumente zu groß wurden.

Marcel

Am 27.04.2024 um 11:41 schrieb Marcel Mueller:

> Noch nicht viel mit trigonometrischen Funktionen gemacht, oder?
> Bereits die Berechnung eines einfachen Sinustons bedingt größere Werte.

Ne, aber imm mal Computergrafik, da hat man Bounds beim Sinus die
immer >= 0 bis < 360 Graf oder >= -180 bis < 180 Grad liegen.

> Die eingebaute sin-Funktion wird dann auch langsamer, weil fmod recht
> teuer ist.

Kommt aber insgesamt sicher selten vor. Ich mein wenn Du schon nichts
sinnvolleres kennst als dieses Beispiel ...

> Mit 8-fach Unrolling komme ich bei kleinen Zahlen auf Faktor 3 langsamer
> als sin(), ...

So ein dummes Zeug, die Funktion dahinter hat eine derart große Kom-
plexität, dass Du an der Stelle mit Loop-Unrolling sicher gar nichts
erreichst. Zumal Loop-Unrolling auf heutigen CPUs sowieso nur selten
was bringt und wenn sehr wenig.

> Aber zumindest scheint das Große-Zahlen-Problem mittlerweile halbwegs
> im Griff zu sein. Früher hat fsin ziemlichen Blödsinn ausgerechnet,
> wenn die Argumente zu groß wurden.

Grundsätzlich arbeitet der CORDIC-Algorithmus mit sukzessiver Approxima-
tion. Dazu wird Pi / 2 = 90 Grad zu 1.0, dass man Bit-weise sich annä-
hrern kann. Wie groß der Wert dabei ist ist dann komplett egal da man
höherwertige Bits dann komplett ignorieren kann. Das war früher sicher
nicht anders als heute denn CORDIC ist uralt.

Am 27.04.24 um 14:04 schrieb Bonita Montero:
> Am 27.04.2024 um 11:41 schrieb Marcel Mueller:
>
>> Noch nicht viel mit trigonometrischen Funktionen gemacht, oder?
>> Bereits die Berechnung eines einfachen Sinustons bedingt größere Werte.
>
> Ne, aber imm mal Computergrafik, da hat man Bounds beim Sinus die
> immer >= 0 bis < 360 Graf oder >= -180 bis < 180 Grad liegen.
>
>> Die eingebaute sin-Funktion wird dann auch langsamer, weil fmod recht
>> teuer ist.
>
> Kommt aber insgesamt sicher selten vor. Ich mein wenn Du schon nichts
> sinnvolleres kennst als dieses Beispiel ...

Naja, Berechnung der FFT-Koeffizienten tut's auch.
Ich denke mal im Bereich Digital Signal Processing hat man es selten mit
Werten im Bereich ±π zu tun.

Gleichwohl stimme ich zu, dass wirklich große Werte selten sind. Und
genau darauf dürften die Implementierungen auch optimiert sein.

>> Mit 8-fach Unrolling komme ich bei kleinen Zahlen auf Faktor 3
>> langsamer als sin(), ...
>
> So ein dummes Zeug, die Funktion dahinter hat eine derart große Kom-
> plexität, dass Du an der Stelle mit Loop-Unrolling sicher gar nichts
> erreichst. Zumal Loop-Unrolling auf heutigen CPUs sowieso nur selten
> was bringt und wenn sehr wenig.

Naja, der Funktionsaufruf mit Übertragung über Stack ist schon einiger
Overhead.

> Grundsätzlich arbeitet der CORDIC-Algorithmus mit sukzessiver Approxima-
> tion. Dazu wird Pi / 2 = 90 Grad zu 1.0, dass man Bit-weise sich  annä-
> hrern kann. Wie groß der Wert dabei ist ist dann komplett egal da man
> höherwertige Bits dann komplett ignorieren kann. Das war früher sicher
> nicht anders als heute denn CORDIC ist uralt.

Keine Ahnung, wie das die Runtime macht. Aber die Berechnungszeit von
sin() steigt hier von etwa 10ns bei kleinen Zahlen auf 150ns und mehr
bei großen. Das ist schon eklatant.

Marcel

Am 27.04.2024 um 16:04 schrieb Marcel Mueller:

> Naja, Berechnung der FFT-Koeffizienten tut's auch.

Ja, super, aber dann hat man erst recht einen kleinen Wertebereich.

> Naja, der Funktionsaufruf mit Übertragung über Stack ist schon einiger
> Overhead.

Selbst wenn spielt das anteilig bei so ner komplexen Operation nicht die
Rolle.

> Keine Ahnung, wie das die Runtime macht. Aber die Berechnungszeit von
> sin() steigt hier von etwa 10ns bei kleinen Zahlen auf 150ns und mehr
> bei großen. Das ist schon eklatant.

Du musst echt irgendwas verkehrt messen oder eine uralte CPU haben.
Meine 7950X Zen4-CPU braucht für eine sin()-Operation ca. 7ns mit
g++ 12, das wären ca. 40 Takte bei 5.7GHz; x87 FSIN ist mit diesem
Compiler nur fünfmal langsamer.

Am 27.04.24 um 16:28 schrieb Bonita Montero:
>> Keine Ahnung, wie das die Runtime macht. Aber die Berechnungszeit von
>> sin() steigt hier von etwa 10ns bei kleinen Zahlen auf 150ns und mehr
>> bei großen. Das ist schon eklatant.
>
> Du musst echt irgendwas verkehrt messen oder eine uralte CPU haben.

Knapp 10 Jahre schätze ich.

> Meine 7950X Zen4-CPU braucht für eine sin()-Operation ca. 7ns mit
> g++ 12, das wären ca. 40 Takte bei 5.7GHz;

Aber 10ns vs. 7ns finde ich dafür ziemlich gut. Das wären ungefähr
genauso viele Takte.

> x87 FSIN ist mit diesem
> Compiler nur fünfmal langsamer.

War hier nur gut 3-mal, mit Parameterübergabe am Stack eher 4-mal.

Marcel

Am 27.04.2024 um 21:44 schrieb Marcel Mueller:

> Aber 10ns vs. 7ns finde ich dafür ziemlich gut. Das wären ungefähr
> genauso viele Takte.

Ja, aber es gibt keine großen Timings bei großen Werten, auch wenn
ich 1E10 hinzuaddiere.

> War hier nur gut 3-mal, mit Parameterübergabe am Stack eher 4-mal.

Die Parameterübergabe am Stack macht bei so einer komplexen
Operation gar nichts aus. Der Code für den g++ nutzt ja inline
-Assembly und ist genauso schnell wie der unter Windows mit dem
Umkopieren.

Am 27.04.2024 um 11:41 schrieb Marcel Mueller:

> Die eingebaute sin-Funktion wird dann auch langsamer, weil fmod recht
> teuer ist.

fmod() ist bei kleinen Exponenten-Differenzen bis einschließlich 11
bei der glibc sehr schnell, denn dann kann der Dividend einfach nach
links geschoben werden und mit Integer-Operationen dividiert werden.

Ich hab mal ein sin() gebastelt das mit dem Taylor-Algorithmus arbeitet
und den ganzen Wertebereich eines doubles handeln kann. Hier der Code
mit nem kleinen Benchmark:

#include <iostream>
#include <span>
#include <array>
#include <bit>
#include <numeric>
#include <chrono>
#include <atomic>

using namespace std;
using namespace chrono;

double taylor( double angle );

atomic<double> aSum( 0.0 );

int main()
{ constexpr double PI = 3.14159265358979323846264338327950288;
vector<double> angles;
for( int i = 0; i != 100; ++i )
angles.emplace_back( i * PI / 100 );
double sum = 0.0, l = 0.0;
auto start = high_resolution_clock::now();
for( size_t i = 1'000'000; i; --i )
for( double d : angles )
sum += l = taylor( d + l ),
l *= numeric_limits<double>::denorm_min();
cout << duration_cast<nanoseconds>( high_resolution_clock::now() -
start ).count() / 1.0e8 << endl;
aSum = sum;
}

double taylor( double angle )
{ double neg = 1.0;
if( angle < 0.0 )
angle = -angle,
neg = -1.0;
constexpr double PI = 3.14159265358979323846264338327950288;
angle = fmod( angle, 2 * PI );
if( angle > PI )
angle -= PI,
neg = -neg;
if( angle > PI / 2 )
angle = PI - angle;
static double const rawReverseFacs[] =
{
0x1.0000000000000p+0,
0x1.5555555555555p-3,
0x1.1111111111111p-7,
0x1.a01a01a01a01ap-13,
0x1.71de3a556c734p-19,
0x1.ae64567f544e4p-26,
0x1.6124613a86d09p-33,
0x1.ae7f3e733b81fp-41,
0x1.952c77030ad4ap-49,
0x1.2f49b46814157p-57,
0x1.71b8ef6dcf572p-66
};
span reverseFacs( rawReverseFacs );
array<double, 11> chain;
double poly = angle;
for( size_t i = 0; i != 11; ++i )
chain[i] = bit_cast<double>( (uint64_t)(i % 2) << 63 | 0x3FFull << 52
) * poly * reverseFacs[i],
poly *= angle * angle;
return neg * accumulate( chain.rbegin(), chain.rend(), 0.0 );
}

Auf jeden Fall ist bei mir dann die mittlere Zeit bis Zum Ergebnis ca.
30,8ns bei 5.7GHz, also ca. 175 Takte. Für wirklich große Winkel kann
das fmod() schon reinhauen, aber der Regelfall ist das aus besagtem
Grund wirklich nicht.
Ich mein in Anbetracht der Performance obigen Algorithmusses und der
Tatsache, dass FSIN ähnlich langsam ist liegt der Verdacht nahe, dass
FSIN wirklich mit einer Polynomlösung arbeitet, wie folngende Intel
-Seite meint:
https://www.intel.com/content/www/us/en/developer/articles/technical/the-difference-between-x87-instructions-and-mathematical-functions.html
Ich frag mich aber wieso die meinen, dass die Taylor-Lösung schneller
ist als die Cordic-Lösung, ist sie ja nachweislich nicht.

Total faszinierend: sin() gibt mit der glibc Bit-identische Ergebnisse
zu FSIN, braucht aber nur 40% der CPU-Zeit weils auf Cordic basiert.