Merkwürdige Effekte durch __restrict
-
Hallo allerseits!
Ich bin auf ein merkwürdiges Verhalten meines Programmes durch die Benutzung von restrict (bzw. seiner GCC-Erweiterung für C++) gestoßen. Folgender Code führt eine Bereichnung durch, wobei alle 5 beteiligten Arrays garantiert disjunkt sind:
template <typename T, typename Idx> inline void apply(const Idx& idx, const T* const __restrict H, const T* const __restrict in_re, const T* const __restrict in_im, T* const __restrict out_re, T* const __restrict out_im) { for(int j = 1; j <= idx.cols(); ++j) { double re_tmp = 0., im_tmp = 0.; for(int i = idx.lower_idx(j)+1; i <= idx.upper_idx(j); ++i) { const double Hij = H[idx(i, j)]; out_re[i-1] += Hij * in_im[j-1]; out_im[i-1] -= Hij * in_re[j-1]; re_tmp += Hij * in_im[i-1]; im_tmp -= Hij * in_re[i-1]; } out_re[j-1] += re_tmp; out_im[j-1] += im_tmp; const double Hjj = H[idx(j, j)]; out_re[j-1] += Hjj * in_im[j-1]; out_im[j-1] -= Hjj * in_re[j-1]; } }
Der Code ist absolut performance-kritisch, weshalb ich es für eine gute Idee hielt, __restrict zu verwenden. Tatsächlich führt das dazu, dass andere Zahlen aus der Rechnung resultieren (die sich zwar nur knapp oberhalb vom Maschinenepsilon unterscheiden), als wenn ich kein __restrict nähme. Leider macht das aber im weiteren Verlauf mein Programm gut 100% langsamer.
Kompiliert wird das ganze bei mir mit -O3 -march=native (insb. mit SSE/AVX2 Instructions). Eigentlich auch noch mit -ffast-math, aber das würde hier nur für Verwirrung sorgen.UB kann ich soweit ausschließen, dass es weder Warnings auf dem Level -Wall -Wextra -Wpedantic gibt, noch, dass -fsanitize=address irgendetwas beanstandet hätte. Der Effekt tritt auf, solange auch nur irgendein Funktionsargument mit __restrict kommt, egal welches. Sobald alle __restricts entfernt sind, wird der Code doppelt so schnell, zeigt minimal andere Zahlen und auch der ASM-Output ist bedeutent kürzer:
mov eax, DWORD PTR 4[rdx] test eax, eax jle .L224 mov r14d, DWORD PTR 12[rdx] mov r15d, DWORD PTR 8[rdx] dec eax lea ebx, 1[r14+r15] mov r13d, DWORD PTR [rdx] add rax, 2 movsx rdx, ebx mov DWORD PTR -72[rbp], ebx mov QWORD PTR -80[rbp], rax lea rbx, 0[0+rdx*8] mov rdx, QWORD PTR -64[rbp] mov QWORD PTR -88[rbp], rbx movsx rbx, r14d lea r11, [rdx+rbx*8] mov QWORD PTR -96[rbp], rbx mov r12, -8 xor ebx, ebx mov edx, 1 .p2align 4,,10 .p2align 3 .L223: mov eax, edx sub eax, r14d test eax, eax mov r10d, 1 cmovle eax, r10d lea r9d, [r15+rdx] inc eax cmp r9d, r13d cmovg r9d, r13d cmp eax, r9d jg .L267 movsx r10, ebx add r10, QWORD PTR -96[rbp] lea r10, [r12+r10*8] cdqe add r10, QWORD PTR -64[rbp] vmovapd xmm1, xmm4 vmovapd xmm2, xmm4 .p2align 4,,10 .p2align 3 .L222: vmovsd xmm0, QWORD PTR [r10+rax*8] vmovsd xmm3, QWORD PTR -8[r8+rdx*8] vfmadd213sd xmm3, xmm0, QWORD PTR -8[rsi+rax*8] vmovsd QWORD PTR -8[rsi+rax*8], xmm3 vmovsd xmm3, QWORD PTR -8[rdi+rdx*8] vfnmadd213sd xmm3, xmm0, QWORD PTR -8[rcx+rax*8] vmovsd QWORD PTR -8[rcx+rax*8], xmm3 vfmadd231sd xmm2, xmm0, QWORD PTR -8[r8+rax*8] vfnmadd231sd xmm1, xmm0, QWORD PTR -8[rdi+rax*8] inc rax cmp r9d, eax jge .L222 .L221: vaddsd xmm2, xmm2, QWORD PTR -8[rsi+rdx*8] add ebx, DWORD PTR -72[rbp] sub r12, 8 vmovsd QWORD PTR -8[rsi+rdx*8], xmm2 vaddsd xmm1, xmm1, QWORD PTR -8[rcx+rdx*8] vmovsd QWORD PTR -8[rcx+rdx*8], xmm1 vmovsd xmm0, QWORD PTR [r11] vmovsd xmm1, QWORD PTR -8[r8+rdx*8] add r11, QWORD PTR -88[rbp] vfmadd213sd xmm1, xmm0, QWORD PTR -8[rsi+rdx*8] vmovsd QWORD PTR -8[rsi+rdx*8], xmm1 vmovsd xmm5, QWORD PTR -8[rcx+rdx*8] vfnmadd132sd xmm0, xmm5, QWORD PTR -8[rdi+rdx*8] vmovsd QWORD PTR -8[rcx+rdx*8], xmm0 inc rdx cmp QWORD PTR -80[rbp], rdx jne .L223
gegen mehr als dreimal so viele Instructions mit __restrict:
mov eax, DWORD PTR 4[rdx] test eax, eax jle .L228 mov ebx, DWORD PTR [rdx] mov edi, DWORD PTR 12[rdx] mov DWORD PTR -240[rbp], ebx mov ebx, DWORD PTR 8[rdx] mov DWORD PTR -204[rbp], edi mov DWORD PTR -248[rbp], ebx lea ebx, 1[rdi+rbx] movsx rdx, ebx mov DWORD PTR -208[rbp], ebx lea rbx, 0[0+rdx*8] mov QWORD PTR -256[rbp], rbx dec eax movsx rbx, edi mov rdi, QWORD PTR -176[rbp] add rax, 2 mov QWORD PTR -224[rbp], rbx mov QWORD PTR -264[rbp], rax lea rbx, [rdi+rbx*8] mov rax, QWORD PTR -336[rbp] mov rdi, QWORD PTR -232[rbp] mov QWORD PTR -168[rbp], rbx add rax, rdi mov QWORD PTR -280[rbp], rax mov rax, QWORD PTR -344[rbp] mov DWORD PTR -132[rbp], 0 add rax, rdi mov QWORD PTR -288[rbp], rax mov QWORD PTR -64[rbp], r11 mov r12, r11 mov r11, QWORD PTR -312[rbp] mov r14, r15 mov r13d, 1 mov r15, r9 .p2align 4,,10 .p2align 3 .L227: mov eax, r13d sub eax, DWORD PTR -204[rbp] test eax, eax mov edi, 1 cmovle eax, edi mov edi, DWORD PTR -248[rbp] mov ecx, DWORD PTR -240[rbp] lea edx, [rdi+r13] cmp edx, ecx mov edi, edx cmovg edi, ecx mov rdx, QWORD PTR -152[rbp] lea ebx, 1[rax] mov DWORD PTR -216[rbp], r13d mov DWORD PTR -136[rbp], ebx mov DWORD PTR -160[rbp], edi vmovsd xmm4, QWORD PTR -8[rdx+r13*8] cmp ebx, edi jg .L271 sub edi, eax mov rax, QWORD PTR -280[rbp] movsx rcx, ebx lea rsi, [rax+rcx] lea r10, 0[0+rsi*8] lea rbx, -8[r10] mov rax, QWORD PTR -232[rbp] lea rdx, [r11+rbx] mov QWORD PTR -184[rbp], rbx add r10, 24 mov rbx, QWORD PTR -64[rbp] mov QWORD PTR -112[rbp], rcx mov r9d, edi add rcx, rax lea r8, [r11+r10] mov QWORD PTR -88[rbp], r9 mov QWORD PTR -96[rbp], r8 add rbx, 8 mov r8, QWORD PTR -120[rbp] lea r9, 0[0+rcx*8] mov QWORD PTR -80[rbp], rbx lea rbx, -8[r9] add r9, 24 mov DWORD PTR -72[rbp], edi lea rdi, [r8+rbx] add r8, r9 mov QWORD PTR -104[rbp], r8 movsx r8, DWORD PTR -132[rbp] mov QWORD PTR -192[rbp], rsi mov rsi, QWORD PTR -224[rbp] mov QWORD PTR -272[rbp], r8 sub r8, r13 mov rax, QWORD PTR -176[rbp] add rsi, r8 add rsi, QWORD PTR -112[rbp] sal rsi, 3 mov r8, rax add r8, rsi lea rax, 32[rax+rsi] cmp rbx, r10 lea rsi, [r11+r9] mov QWORD PTR -128[rbp], rsi setge sil cmp r9, QWORD PTR -184[rbp] setle r9b or r9d, esi cmp QWORD PTR -96[rbp], rdi setbe sil cmp rdx, QWORD PTR -104[rbp] setnb r10b mov QWORD PTR -200[rbp], rcx lea rcx, [r11+rbx] mov ebx, DWORD PTR -72[rbp] or esi, r10d and r9d, esi lea esi, -1[rbx] cmp esi, 2 seta sil and r9d, esi cmp rdi, QWORD PTR -128[rbp] setnb bl cmp rcx, QWORD PTR -104[rbp] setnb sil or ebx, esi mov rsi, QWORD PTR -192[rbp] mov r10, QWORD PTR -80[rbp] add rsi, QWORD PTR -88[rbp] lea rsi, -8[r11+rsi*8] and r9d, ebx cmp QWORD PTR -64[rbp], rsi setnb sil cmp rdx, r10 setnb bl or esi, ebx and esi, r9d mov r9, QWORD PTR -88[rbp] add r9, QWORD PTR -200[rbp] lea r9, -8[r11+r9*8] cmp QWORD PTR -64[rbp], r9 setnb r9b cmp rcx, r10 setnb bl or r9d, ebx and esi, r9d cmp QWORD PTR -96[rbp], r8 setbe r9b cmp rdx, rax setnb r10b or r9d, r10d test sil, r9b je .L222 cmp r8, QWORD PTR -128[rbp] setnb sil cmp rcx, rax setnb r9b or sil, r9b je .L222 mov rax, QWORD PTR -112[rbp] mov rsi, QWORD PTR -120[rbp] add rax, QWORD PTR -288[rbp] lea r9, -8[rsi+rax*8] mov esi, DWORD PTR -72[rbp] vbroadcastsd ymm5, QWORD PTR -8[r12+r13*8] shr esi, 2 vbroadcastsd ymm6, xmm4 sal rsi, 5 xor eax, eax vmovapd xmm1, xmm7 vmovapd xmm3, xmm7 .p2align 4,,10 .p2align 3 .L223: vmovupd ymm0, YMMWORD PTR [r8+rax] vmovapd ymm2, ymm0 vfmadd213pd ymm2, ymm6, YMMWORD PTR [rcx+rax] vmovupd YMMWORD PTR [rcx+rax], ymm2 vmovapd ymm2, ymm0 vfnmadd213pd ymm2, ymm5, YMMWORD PTR [rdx+rax] vmovupd YMMWORD PTR [rdx+rax], ymm2 vmulpd ymm2, ymm0, YMMWORD PTR [r9+rax] vmulpd ymm0, ymm0, YMMWORD PTR [rdi+rax] add rax, 32 vaddsd xmm3, xmm3, xmm2 vsubsd xmm1, xmm1, xmm0 vunpckhpd xmm9, xmm2, xmm2 vunpckhpd xmm8, xmm0, xmm0 vaddsd xmm3, xmm3, xmm9 vsubsd xmm1, xmm1, xmm8 vextractf128 xmm2, ymm2, 0x1 vextractf128 xmm0, ymm0, 0x1 vaddsd xmm3, xmm3, xmm2 vsubsd xmm1, xmm1, xmm0 vunpckhpd xmm2, xmm2, xmm2 vunpckhpd xmm0, xmm0, xmm0 vaddsd xmm3, xmm3, xmm2 vsubsd xmm1, xmm1, xmm0 cmp rax, rsi jne .L223 mov edi, DWORD PTR -72[rbp] mov eax, DWORD PTR -136[rbp] mov edx, edi and edx, -4 add eax, edx cmp edi, edx je .L221 mov ebx, DWORD PTR -216[rbp] mov r8d, DWORD PTR -204[rbp] mov edx, eax mov r10d, DWORD PTR -132[rbp] sub edx, ebx add edx, r8d mov rdi, QWORD PTR -176[rbp] add edx, r10d movsx rdx, edx vmovsd xmm0, QWORD PTR [rdi+rdx*8] movsx rdx, eax sal rdx, 3 lea rcx, -8[rdx] lea rsi, [r15+rcx] vmovapd xmm2, xmm0 vfmadd213sd xmm2, xmm4, QWORD PTR [rsi] add rcx, r14 vmovsd QWORD PTR [rsi], xmm2 vmovsd xmm2, QWORD PTR -8[r12+r13*8] vfnmadd213sd xmm2, xmm0, QWORD PTR [rcx] vmovsd QWORD PTR [rcx], xmm2 mov rcx, QWORD PTR -152[rbp] vfnmadd231sd xmm1, xmm0, QWORD PTR -8[r12+rdx] vfmadd231sd xmm3, xmm0, QWORD PTR -8[rcx+rdx] lea edx, 1[rax] cmp DWORD PTR -160[rbp], edx jl .L221 mov r9d, edx sub r9d, ebx mov ecx, r9d add ecx, r8d add ecx, r10d movsx rdx, edx lea rsi, 0[0+rdx*8] movsx rcx, ecx vmovsd xmm0, QWORD PTR [rdi+rcx*8] lea rcx, -8[rsi] mov r9, rdi vmovapd xmm2, xmm4 lea rdi, [r15+rcx] vfmadd213sd xmm2, xmm0, QWORD PTR [rdi] add rcx, r14 add eax, 2 vmovsd QWORD PTR [rdi], xmm2 vmovsd xmm2, QWORD PTR -8[r12+r13*8] mov rdi, QWORD PTR -152[rbp] vfnmadd213sd xmm2, xmm0, QWORD PTR [rcx] vfmadd231sd xmm3, xmm0, QWORD PTR -8[rdi+rsi] vmovsd QWORD PTR [rcx], xmm2 vfnmadd231sd xmm1, xmm0, QWORD PTR -8[r12+rsi] cmp DWORD PTR -160[rbp], eax jl .L221 sub eax, ebx add eax, r8d add eax, r10d cdqe vmovsd xmm0, QWORD PTR [r9+rax*8] lea rcx, [r15+rsi] vmovapd xmm2, xmm4 vfmadd213sd xmm2, xmm0, QWORD PTR [rcx] mov rax, rsi add rax, r14 vfmadd231sd xmm3, xmm0, QWORD PTR [rdi+rdx*8] vmovsd QWORD PTR [rcx], xmm2 vmovsd xmm2, QWORD PTR -8[r12+r13*8] vfnmadd213sd xmm2, xmm0, QWORD PTR [rax] vmovsd QWORD PTR [rax], xmm2 vfnmadd231sd xmm1, xmm0, QWORD PTR [r12+rdx*8] .L221: vaddsd xmm3, xmm3, QWORD PTR -8[r15+r13*8] mov rax, QWORD PTR -168[rbp] mov rdi, QWORD PTR -256[rbp] vmovsd QWORD PTR -8[r15+r13*8], xmm3 vaddsd xmm0, xmm1, QWORD PTR -8[r14+r13*8] mov ebx, DWORD PTR -208[rbp] vmovsd QWORD PTR -8[r14+r13*8], xmm0 vmovsd xmm0, QWORD PTR [rax] add rax, rdi vfmadd213sd xmm4, xmm0, QWORD PTR -8[r15+r13*8] mov QWORD PTR -168[rbp], rax mov rax, QWORD PTR -80[rbp] add DWORD PTR -132[rbp], ebx mov QWORD PTR -64[rbp], rax vmovsd QWORD PTR -8[r15+r13*8], xmm4 vmovsd xmm4, QWORD PTR -8[r14+r13*8] vfnmadd132sd xmm0, xmm4, QWORD PTR -8[r12+r13*8] vmovsd QWORD PTR -8[r14+r13*8], xmm0 inc r13 cmp r13, QWORD PTR -264[rbp] jne .L227
Die nahelegende Lösung wäre sicherlich, eben auf __restrict zu pfeifen. Aber ich möchte gerne verstehen, was hier schiefgeht bzw. wie ich das weiter untersuchen könnte. Der ASM-Code sagt mir im Detail jetzt nicht furchtbar viel. Vielen Dank!
p.s. Ich verwende GCC 9.3.
-
@Jodocus sagte in Merkwürdige Effekte durch __restrict:
Tatsächlich führt das dazu, dass andere Zahlen aus der Rechnung resultieren (die sich zwar nur knapp oberhalb vom Maschinenepsilon unterscheiden), als wenn ich kein __restrict nähme. Leider macht das aber im weiteren Verlauf mein Programm gut 100% langsamer.
Hm. Denormals? Die sollten mit
-ffast-math
eigentlich abgedreht sein. Das funktioniert aber soweit ich weiss nur wenn du das Executable selbst mit-ffast-math
baust, und es kann auf jeden Fall passieren dass Denormals nicht deaktiviert sind wenn dein Code als Library in einem anderen Executable läuft.Probier vielleicht einfach mal
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON)
und_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON)
vor deinen Berechnungen zu machen. Siehe https://scc.ustc.edu.cn/zlsc/chinagrid/intel/compiler_c/main_cls/GUID-1659EAE1-583E-44EE-BDEA-7C68C46061C7.htmWenns daran liegt, musst du dir aber 'was überlegen. Denn einfach so FTZ und DAZ zu setzen (und nie wieder zu restoren) ist auf jeden Fall "pfui!".
-
Mal eine andere Herangehensweise:
Wieso übergibst du nicht ein struct und returnst ein struct (dank RTO kein Problem)? Dann dürftest du das Aliasing-Problem nicht haben und der Compiler dürfte besser optimieren können.
-
@Steffo Dir ist schon aufgefallen dass das alles Arras sind, oder?
-
@hustbaer Habe ich mal probiert, aber macht keinen Unterschied. Zur Info: Mein Code läuft nicht als Library, also ich habe volle Kontrolle über alle Linker/Compiler Flags vom Code bis zur Executable.
Hintergrundinfo: Die Routine wird während einer Adams-Bashforth-Routine zum Lösen einer ODE aufgerufen, mit adaptiver Schrittweite. Je weniger Schritte (= Aufrufe dieser Funktion), desto besser.
Meine obigen Beobachtungen stehen nach wie vor. Mit restrict sind die Zahlen leicht anders und am Ende werden ca. 100.000 Aufrufe auf die Funktion gemacht. Nehme ich restrict raus, sind es nur noch 50.000 Aufrufe. Die Ergebnisse sind in beiden Fällen bis auf die eingestellte Genauigkeit der Konvergenz des ODE-Solvers korrekt.Zu allem Überfluss kommt jetzt hinzu: wenn ich am Code leicht etwas ändere, z.B. die Zeilen mit Hjj am Anfang statt am Ende der äußeren Loop setze (was mathematisch zwar äquivalent ist), also
for(int j = 1; j <= idx.cols(); ++j) { const double Hjj = H[idx(j, j)]; // * out_re[j-1] += Hjj * in_im[j-1]; // * out_im[j-1] -= Hjj * in_re[j-1]; // * double re_tmp = 0., im_tmp = 0.; for(int i = idx.lower_idx(j)+1; i <= idx.upper_idx(j); ++i) { const double Hij = H[idx(i, j)]; out_re[i-1] += Hij * in_im[j-1]; out_im[i-1] -= Hij * in_re[j-1]; re_tmp += Hij * in_im[i-1]; im_tmp -= Hij * in_re[i-1]; } out_re[j-1] += re_tmp; out_im[j-1] += im_tmp; }
Dann ist das Ergebnis von den __restricts völlig unabhängig (wie man es ja auch erwarten würde)! Mir ist schon klar, dass für die Fließpunktkommazahlen Addition nicht assoziativ/kommutativ sein könnte (die Reihenfolge also technisch, wenn auch nicht mathematisch eine Rolle spielt), aber ich habe explizit -ffast-math ausgeschaltet (und mehrfach überprüft und neugebaut), sodass Optimierungen, die der Compiler durch __restrict machen darf, nicht z.B. an dieser Reihenfolge rütteln dürfte (da sonst nicht mehr IEEE-konform). Warum das in einem Fall doch passiert, riecht für mich irgendwie schon arg nach Compiler-Bug.
Ganz abgesehen davon, dass die ODE-Routine numerisch alles andere als stabil ist, wenn solche Effekte zu 100% Laufzeitsteigerung führen können...
-
Hm, ich hatte deinen Code mal in Godbolt gekippt und habe diese großen Assember-Unterschiede nicht gesehen. Jeweils 95 Zeilen. Allerdings musste ich dazu natürlich das template und das inline rausnehmen und mir eine Fake-Index-Klasse machen. Vielleicht tritt der Effekt, den du siehst, nur bei bestimmten Klassen für Idx und T auf - vielleicht trifft der Compiler irgendwelche Optimierungen? Oder das Problem kommt von irgendwo anders her? Finde ich schon komisch, dass die Änderung im Epsilon-Bereich gleich doppelt so viele Calls zur Folge haben soll. Dann würde ich denken, dass der Solver nicht korrekt arbeitet - oder dass du irgendeinen Spezialfall triffst, wo das tatsächlich große Relevanz hat.
-
@Jodocus sagte in Merkwürdige Effekte durch __restrict:
Zu allem Überfluss kommt jetzt hinzu: wenn ich am Code leicht etwas ändere, z.B. die Zeilen mit Hjj am Anfang statt am Ende der äußeren Loop setze (was mathematisch zwar äquivalent ist), also
Wie groß sind die Felder? Der Loop über i liest vier Felder auf einmal, wenn diese groß sind, ist das schlecht für die Cache Hits.
-
@john-0 sagte in Merkwürdige Effekte durch __restrict:
@Jodocus sagte in Merkwürdige Effekte durch __restrict:
Zu allem Überfluss kommt jetzt hinzu: wenn ich am Code leicht etwas ändere, z.B. die Zeilen mit Hjj am Anfang statt am Ende der äußeren Loop setze (was mathematisch zwar äquivalent ist), also
Wie groß sind die Felder? Der Loop über i liest vier Felder auf einmal, wenn diese groß sind, ist das schlecht für die Cache Hits.
Das variiert, ich kann da keine Rücksicht auf die Cachegröße nehmen.
Ich glaube, ich mache hier Schluss. Es ist nicht ausgeschlossen, dass ich noch das eine oder andere an dem Code ändern muss, und da der Effekt so empfindlich ist (und evtl. am Ende gar nicht mehr auftritt), lohnt sich die Nachforschung an der Stelle nicht. Es bleibt nur ein komischer Nachgeschmack. Trotzdem danke an alle für die Empfehlungen!
-
@Jodocus sagte in Merkwürdige Effekte durch __restrict:
Das variiert, ich kann da keine Rücksicht auf die Cachegröße nehmen.
Das solltest Du aber, weil das zum Teil drastische Auswirkungen hat. Die Optimierung auf Cache Hits ist bei modernen Computer bei HPC die wichtigste Optimierung, die man machen kann. Es nützt Dir nichts einige wenige Taktzyklen einzusparen, wenn Du am Ende tausende Zyklen auf das Laden aus dem RAM warten musst.
-
@Jodocus
Im ersten Beitrag hast du geschrieben du baust mit-ffast-math
und nun schreibst du du baust mit-ffast-math
ausgeschaltet. Jetzt bin ich verwirrt
Meinst du du hast beide Varianten (-ffast-math
und-fno-fast-math
) probiert?Mit restrict sind die Zahlen leicht anders und am Ende werden ca. 100.000 Aufrufe auf die Funktion gemacht. Nehme ich restrict raus, sind es nur noch 50.000 Aufrufe.
OK. Das heisst die Funktion läuft gar nicht unbedingt langsamer, sie erzeugt nur leicht andere Ergebnisse so dass die Werte langsamer konvergieren und du mehr Iterationen brauchst. Das kann mit
-ffast-math
natürlich passieren.Ich bin mir nicht ganz sicher, aber ich denke dass es auch mit
-fno-fast-math
Unterschiede geben kann. Und zwar ist bei sowas wiea = x * y + z
glaube ich nicht definiert ob das Zwischenergebnisx * y
erstmal gerundet wird vor der Addition oder nicht. Der Compiler generiert für deinen Code (mit und ohne restrict) fused-multiply-add Instruktionen (die ganzenvfmaddXXX
undvfnmaddXXX
Dinger), und diese runden das Zwischenergebnis nicht. Da kann ich mir schon vorstellen das was anderes rauskommt wenn z.B. mehr oder weniger fused-multiply-add generiert werden, abhängig von den Annahmen die der Compiler durch restrict oder nicht-restrict treffen kann.Was vielleicht auch noch einen Versucht wert wäre: du könntest
std::fma
in deinem Code verwenden: https://en.cppreference.com/w/cpp/numeric/math/fma
Dann müsste der Compiler immer sicherstellen dass die Zwischenergebnisse nicht gerundet werden, dastd::fma
das garantiert.
Dann bliebe nur noch die Frage obstd::fma
zu mehr oder weniger Iterationen führt, verglichen mit dem ursprünglichen Code ohne restrict.