Die Mathe-Redaktion [Home] - Neues auf einen Blick 30.09.2023 21:28 [Neues] [Forum] [Suche] - Registrieren/Login
 
Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Kubische Gleichung: numerische Näherung schneller als Cardanische Formeln?
Seite 1   [1 2 3 4]   4 Seiten
Autor
Kein bestimmter Bereich Kubische Gleichung: numerische Näherung schneller als Cardanische Formeln?
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Themenstart: 2021-10-13

Hallo zusammen, Anlass für diesen Thread ist die Aussage aus dem Titel dieses Threads, die man im Wikipedia-Artikel zur kubischen Gleichung findet, nämlich hier. Demnach soll die dort beschriebene Näherungslösung mit dem Halleyschen Verfahren "bei sorgfältiger Implementierung" um den Faktor 1,2 bis 10 (!!) schneller sein, als die Lösungen direkt per expliziter Cardanischer Formel zu berechnen. Ich habe daran Zweifel. In dem Wiki-Artikel wird direkt auf den C++-Code dieser "sorgfältig implementierten" Näherungslösung verlinkt, so dass man ihn sich direkt herauskopieren kann. Ich habe jedoch kein C++, um einen Vergleich anzustellen. Was ich habe, ist Python-Code, den ich zur Lösung der allgemeinen kubischen Gleichung geschrieben habe, siehe hier: \sourceon Python \numberson # ax^3+bx^2+cx+d=0 def cbrt(a,b,c,d): # Trivialfall quadratische Gleichung: if a==0: # Trivialfall lineare Gleichung: if b==0: return [-d/c] # Echte quadratische Gleichung: y1=-(c+(c**2-4*b*d)**0.5)/(2*b) return [y1,-c/b-y1] # Echte kubische Gleichung: x=y+r -> y^3+3p-2q=0 r=-b/(3*a) p=(c+b*r)/(3*a) q=-0.5*(d/a+3*p*r+r**3) if p==q==0: return [r,r,r] # Dreifach-Nullstelle s=(q**2+p**3)**0.5 q=(2*q)**(1/3) if s+q==0 else (q+s)**(1/3) y1=q-p/q y2=-0.5*(y1+(-3*y1**2-12*p)**0.5) return [y1+r,-y1-y2+r,y2+r] \sourceoff Ich übergebe der Zeile 1 entsprechend die Koeffizienten der kubischen Gleichung an die Funktion, und alle drei Lösungen werden per Liste zurückgegeben. Die Lösung funktioniert fehlerfrei, soweit ich das bisher beurteilen kann, auch für komplexe Koeffizienten und komplexe Lösungen. Wenn man für einen Geschwindigkeitstest die Trivialfälle (quadratische oder lineare Gleichung) wegließe, beginnt der eigentliche Code erst in Zeile 11, dann noch einen Trivialfall rausfischen in Zeile 14, und in insgesamt 9 Zeilen Code habe ich alle drei Lösungen berechnet und zurückgegeben, ohne irgendeine Schleife mit genauigkeitsrelevantem Abbruchkriterium. Ich kann mir nicht mal entfernt vorstellen, dass das im Wikipedia-Artikel beschriebene Verfahren schneller sein soll (Wendepunkt berechnen, Ober- und Untergrenze der Lösungen bestimmen, Startwert des Näherungsverfahrens nach bestimmten Kriterien aus drei Werten auswählen, und dann noch die eigentliche Iteration inkl. dafür notwendiger Berechnung der Funktionswerte, ersten und zweiten Ableitung). Hat jemand Lust, meinen Code mal in C++ zu übersetzen und gegen den verlinkten Code aus dem Wikipedia-Artikel von der Uni Köln zu vergleichen? Ich werde meinerseits wohl mal den C++-Code aus dem Artikel in Python übersetzen, aber in Sachen nackter Geschwindigkeit ist der Vergleich in C++ sicherlich am aussagekräftigsten. Oder hat jemand eine eigene Erfahrung zu diesem Thema? Ciao, Thomas

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.1, eingetragen 2021-10-13

Hallo MontyPythagoras, eine interessante Fragestellung. Auf mathC, der Seite mit dem C++ Programm, ist dieser Artikel verlinkt. In Table 4 werden die CPU-Zeiten verschiedener Algorithmen verglichen. Anstatt den Source-Code aus dem pdf zu kopieren, habe ich mir die gesamte mathC-Lib heruntergeladen. Das Paket beinhaltet Tests für alle Programme. Die Tests für eqn_cubic muss man mindestens 100000 Mal wiederholen, damit die wahre CPU-Zeit größer als die Systemzeit ist. \sourceon $ cat README data1: x = -1 data2: x = 0.1, 0.3, 0.6 data3: x = -1, 1, 2 data4: x = 2, 2, 10 (1 double root!) data5: x = -1, -1, -1 (1 triple root!) data6: x = 1.0, 1.000001, 1.000002 (nasty test proposed by Flocke, 2015, middle root is inflection point) data7: x = 1, 1.0e+07, 1.0e+14 (nasty test proposed by Strobach, 2001) \sourceoff \sourceon $ cat data5 100000 1 3 3 1 $ ./test_p3 < data5 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 CPU time: 5 ms \sourceoff Sonst wird immer 2 ms ausgegeben. Hier die Daten meiner CPU: \sourceon $ cat /proc/cpuinfo processor : 0 vendor_id : GenuineIntel cpu family : 6 model : 58 model name : Intel(R) Core(TM) i3-3240 CPU @ 3.40GHz stepping : 9 microcode : 0x21 cpu MHz : 1600.004 cache size : 3072 KB \sourceoff Im Makefile wird die lib mit g++ mit der Option-Ofast kompiliert. Ich bin gerne bereit weitere Tests zu machen. Dazu muss ich mich aber erst einmal einlesen. Ich habe auch noch einen Rechner mit einem Core i5.

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.2, vom Themenstarter, eingetragen 2021-10-13

Hallo AlphaSigma, danke für Deine Mühe! Ich bin dabei, den Code noch in leichten Variationen zu testen und die letzten Nanosekunden aus Python rauszuquetschen. Was der C++ Compiler später daraus macht, ist sicher noch einmal eine andere Nummer, aber mein i7-Prozessor in meinem Surface Pro schafft derzeit eine Zeit von 2,24µs. Ich lasse $2^{22}$ mal die gleiche Lösung berechnen, das dauert etwa 9 Sekunden. Ciao, Thomas

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.3, vom Themenstarter, eingetragen 2021-10-13

...so, ich habe noch einmal gute 7% rausgequetscht. Divisionen durch Multiplikationen mit gespeichertem Kehrwert ersetzt, Anzahl der arithmetischen Operationen generell reduziert. \sourceon Python \numberson def cbrt(a,b,c,d): # Trivialfall quadratische Gleichung: if a==0: # Trivialfall lineare Gleichung: if b==0: return [-d/c] # Echte quadratische Gleichung: y1=-(c+(c**2-4*b*d)**0.5)/(2*b) return [y1,-c/b-y1] # Echte kubische Gleichung: x=y+r -> y^3+3p-2q=0 s=1/(3*a) r=-b*s p=c*s-r**2 q=-1.5*s*(d+c*r)+r**3 if p==q==0: return [r,r,r] # Dreifach-Nullstelle s=(q**2+p**3)**0.5+q q=(2*q)**(1/3) if s==0 else s**(1/3) y1=q-p/q y2=-0.5*(y1+(-3*y1**2-12*p)**0.5) return [y1+r,-y1-y2+r,y2+r] \sourceoff Ciao, Thomas

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.4, eingetragen 2021-10-13

Hallo MontyPythagoras, beim mathc Paket kann man die Methode für die Funktion eqn_cubic konfigurieren. \sourceon C #define METHOD 1 // Deiters & Macias-Salinas #if METHOD == 0 // Cardano's method \sourceoff Mit METHOD 1 benötigt das Programm für 1000000000 Iterationen ca. 12101 ms, mit METHOD 0 ca. 33727 ms für data5. \sourceon Polynomkoeffizienten 1 3 3 1 Lösung -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 \sourceoff CPU wie in Beitrag No.1.

   Profil
pzktupel
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 02.09.2017
Mitteilungen: 2445
Wohnort: Thüringen
  Beitrag No.5, eingetragen 2021-10-13

Muss aber einräumen, wo jetzt der Vorteil liegt ? Die Berechnung der Nullstelle ist nun 300% schneller...aber bzgl 1 Mrd Durchläufe, nun das einzelne Ergebnis 0,000000072s eher vorliegt erscheint mir nicht schlüssig.

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.6, eingetragen 2021-10-13

Hallo pzktupel, was genau erscheint Dir nicht schlüssig? Die numerische Methode nach Deiters & Macias-Salinas benötigt für 1000000000 Iterationen ca. 12101 ms, also 12 ns für 1 Iteration. Die Berechnung mit Cardano benötigt für 1000000000 Iterationen ca. 33727 ms, also 34 ns für 1 Iteration. Der Faktor 3 passt doch zum Faktor 1,2 bis 10 aus dem Themenstart.

   Profil
schnitzel
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 26.02.2009
Mitteilungen: 227
  Beitrag No.7, eingetragen 2021-10-14

Hi, \quoteon(2021-10-13 21:23 - MontyPythagoras in Beitrag No. 3) ...so, ich habe noch einmal gute 7% rausgequetscht. Divisionen durch Multiplikationen mit gespeichertem Kehrwert ersetzt, Anzahl der arithmetischen Operationen generell reduziert. \quoteoff Es könnte sich noch lohnen Potenzen wie `x**2` oder `x**3` durch `x*x` und `x*x*x` zu ersetzen. Das ist meist etwas schneller. Gruß

   Profil
pzktupel
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 02.09.2017
Mitteilungen: 2445
Wohnort: Thüringen
  Beitrag No.8, eingetragen 2021-10-14

@AlphaSigma Hast recht. Mir war nur so, das PCs heute so schnell arbeiten und um eine Nullstelle einer kubischen Gleichung zu lösen, ist der zeitliche Aspekt nachrangig geworden. Eben weil in 12s oder 30s 1000 Mio mal das Ergebnis berechnet wurde.

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.9, vom Themenstarter, eingetragen 2021-10-14

Hallo AlphaSigma, danke schon einmal für Deine Mühe. Das Beispiel von Dir wirft aber genau die Fragen und Zweifel auf, die ich schon angedeutet hatte. Hier ist die Berechnung mit dem numerischen Algorithmus augenscheinlich schneller. Aber ist sie das wirklich? Wenn man DEREN Implementation der cardanischen Formel verwendet, war das zu erwarten. Sie haben ja schließlich veröffentlicht, dass der Näherungsalgorithmus schneller sein soll, und man kann davon ausgehen, dass sie den Algorithmus auch perfekt optimiert haben. Aber haben sie das auch für die Umsetzung der cardanischen Formeln getan? Ich denke nicht. Nachfolgende Fragen und Kritikpunkte drängen sich auf: 1. Die von Dir berechnete Lösung rechnet mit den Koeffizienten 1 3 3 1. Das ist ein Trivialfall, nämlich dreifache Nullstelle. Berechne ich den Wendepunkt, was der erste Schritt in deren Algorithmus ist, habe ich, oh Wunder, auch direkt die Lösung. Die iterative Schleife, die die eigentliche Zeit kostet, wird nicht ein einziges Mal durchlaufen. 2. Die Koeffizienten sind in dem Beispiel reell, die Lösungen auch. Funktioniert der Algorithmus überhaupt mit komplexen Koeffizienten und komplexen Lösungen? 3. Wie haben sie die cardanischen Formeln umgesetzt? Etwa wie in den Wikipediaartikeln zur kubischen Gleichung (hier) oder den Cardanischen Formeln (hier) beschrieben, mit den ganzen Fallunterscheidungen und $\sin$ und $\cos$ und $\cosh$ und so weiter? Wo man $\Delta=\left(\frac{q}2\right)^2 + \left(\frac{p}3\right)^3 = \frac{27 A^2 D^2+4 B^3D-18 A B C D+4 A C^3-B^2 C^2}{108 A^4}$ berechnet, und die drei Lösungen im komplexen Fall durch Multiplikation mit den dritten Einheitswurzeln berechnet? Geht es denn noch uneffektiver? Dann muss ich mich nicht wundern, wenn die cardanische Lösung langsam ist. 4. Wenn ich das richtig sehe, wird das relative Abbruchkriterium für den Näherungsalgorithmus bei $10^{-6}$ gesetzt. Das ist schon sehr grob, könnte aber aufgrund der schnellen, weil kubischen Konvergenz des Halley-Verfahrens ausreichen, um die Lösungen auf $10^{-18}$ genau zu treffen. Könnte! Es könnte aber auch Situationen geben, wo die Konvergenz ungünstig ist. Ich würde hier mal das Konvergenzkriterium noch weiter runterschrauben, um sicherzustellen, dass die maximale, mit "float" erreichbare Genauigkeit auch sicher erreicht ist, oder man müsste noch eine Prüfung durch Einsetzen der Nullstelle hinzufügen, ob der Funktionswert an dieser Stelle tatsächlich nahe null ist. Daher habe ich meinen Code hier mal reingeschrieben, wo ich eine wirklich effiziente Umsetzung der cardanischen Formel angestrebt habe. Keine Fallunterscheidungen, keine riesigen Formeln wie die oben gezeigte, das wäre alles dramatisch uneffektiv. 9 Zeilen Code mit einer Handvoll algebraischen Anweisungen - und das war's. Ich habe daher arge Zweifel, dass die cardanische Lösung von den Machern des Algorithmus optimiert wurde, und daher würde ich gerne meinen Code umgesetzt sehen. Ich übersetze es auch gerne selbst in C++. Wir sollten außerdem ein Beispiel wählen, bei dem mit komplexen Koeffizienten drei Lösungen berechnet werden, die ebenfalls komplex sind und keine trivialen Nullstellen enthalten. \quoteon(2021-10-14 07:26 - pzktupel in Beitrag No. 8) (...) Mir war nur so, das PCs heute so schnell arbeiten und um eine Nullstelle einer kubischen Gleichung zu lösen, ist der zeitliche Aspekt nachrangig geworden. Eben weil in 12s oder 30s 1000 Mio mal das Ergebnis berechnet wurde. \quoteoff Es war klar, dass dieser Einwand früher oder später kommt. Wenn man privat nur eine einzige Lösung braucht, ist die Rechenzeit völlig unerheblich, das stimmt. Aber wenn ein Großrechner aufwendige Berechnungen macht, zum Beispiel bei Wettervorhersagen, Flugverkehr, Navigation etc. und er Milliarden von solchen Berechnungen in möglichst kurzer Zeit machen muss, dann ist es eben doch relevant. (Dass dieser Einwand ausgerechnet von demjenigen kommt, der was weiß ich wie viel Zeit in die Frage steckt, ob es irgendwo Primzahl-X-linge gibt, ist schon ein wenig ironisch. 😛) In einem Wikipedia-Artikel sollte nichts drinstehen, was nicht unabhängig überprüft werden kann. Ciao, Thomas

   Profil
pzktupel
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 02.09.2017
Mitteilungen: 2445
Wohnort: Thüringen
  Beitrag No.10, eingetragen 2021-10-14

"Aber wenn ein Großrechner aufwendige Berechnungen macht, zum Beispiel bei Wettervorhersagen, Flugverkehr, Navigation etc. und er Milliarden von solchen Berechnungen in möglichst kurzer Zeit machen muss, dann ist es eben doch relevant." Das sehe ich auf jeden Fall ein, wenn kubische Gleichungen Anwedungen finden. 🙂

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.11, vom Themenstarter, eingetragen 2021-10-14

Hallo pzktupel, das war von mir nicht böse gemeint in Deine Richtung, ich hoffe, Du hast das nicht so verstanden. Ob in meinen Beispielen nun wirklich kubische Gleichungen gelöst werden, weiß ich nicht, aber Tatsache ist fraglos, dass heute genauso wie vor 50 Jahren die Rechenleistung möglichst effizient genutzt werden musste und muss, erst recht in Zeiten von "Big Data". Vielleicht bin ich am Ende auch überzeugt, gleichzeitig sicher auch überrascht, dass der Näherungsalgorithmus tatsächlich schneller ist, aber im Moment glaube ich das noch nicht. Natürlich ist mir auch klar, dass tief unter der Hochsprachenoberfläche, unten im Maschinendeck, der Prozessor auch nicht anders kann als den Sinus und Kosinus durch eine Reihe zu berechnen, wenn ich eine komplexe Zahl mit (1/3) potenziere. Aber auf Prozessorlevel ist das so perfekt durchoptimiert, da wird keine pico-Sekunde verloren. Wer hier unnötigerweise auf Hochsprachenlevel diese ganzen Fallunterscheidungen macht und die cardanischen Lösungen wie in den Artikeln beschrieben programmiert, hat in Sachen Effizienz schon den ersten gravierenden Fehler gemacht. Ich bin gespannt. Ciao, Thomas

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.12, vom Themenstarter, eingetragen 2021-10-14

Hallo schnitzel, \quoteon(2021-10-14 06:49 - schnitzel in Beitrag No. 7) Hi, \quoteon(2021-10-13 21:23 - MontyPythagoras in Beitrag No. 3) ...so, ich habe noch einmal gute 7% rausgequetscht. Divisionen durch Multiplikationen mit gespeichertem Kehrwert ersetzt, Anzahl der arithmetischen Operationen generell reduziert. \quoteoff Es könnte sich noch lohnen Potenzen wie `x**2` oder `x**3` durch `x*x` und `x*x*x` zu ersetzen. Das ist meist etwas schneller. Gruß \quoteoff Das stimmt, und bei Visual Basic mache ich das auch immer. Aber bei Python sagt mir meine (im Moment zugegebenermaßen noch recht überschaubare) Erfahrung mit ein paar Geschwindigkeitstests, dass der Unterschied kaum feststellbar ist. Bei C++ mag das anders aussehen, da sollte man das auf jeden Fall berücksichtigen. Divisionen brauchen auch länger als Multiplikationen, was den Unterschied zwischen dem Code im Threadstart und in Beitrag #3 erklärt. Ciao, Thomas

   Profil
schnitzel
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 26.02.2009
Mitteilungen: 227
  Beitrag No.13, eingetragen 2021-10-14

Hi, das \quoteon(2021-10-14 10:10 - MontyPythagoras in Beitrag No. 12) Aber bei Python sagt mir meine (im Moment zugegebenermaßen noch recht überschaubare) Erfahrung mit ein paar Geschwindigkeitstests, dass der Unterschied kaum feststellbar ist. \quoteoff stimmt tatsächlich nicht. Es ist sogar verblüffend, wie groß der Unterschied sein kann: \sourceon python # runs in jupyter xs = [45, 43434, 4564956845860] for x in xs: x2 = %timeit -o x**2 xx = %timeit -o x*x print(xx.average / x2.average) x3 = %timeit -o x**3 xxx = %timeit -o x*x*x print(xxx.average / x3.average) # 206 ns ± 1.97 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) # 53.3 ns ± 1.48 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) # 0.2588890796874734 # 217 ns ± 3.14 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) # 81 ns ± 1.14 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) # 0.3727411192757285 # 218 ns ± 2.75 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) # 55.5 ns ± 1.66 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) # 0.2546575625405366 # 247 ns ± 7.83 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) # 98.7 ns ± 2.79 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) # 0.3998553218182748 # 262 ns ± 12.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) # 68.5 ns ± 2.46 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) # 0.26097106665871633 # 291 ns ± 3.91 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) # 119 ns ± 1.96 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) # 0.4078796332879913 \sourceoff Bei deiner Funktion: \sourceon python # jupyter t = (749456, 356767, 69656, 565566) assert cbrt(*t) == cbrt2(*t) t1 = %timeit -r 10 -o cbrt(*t) t2 = %timeit -r 10 -o cbrt2(*t) print(t2.average/t1.average) # 1.51 µs ± 21.5 ns per loop (mean ± std. dev. of 10 runs, 1000000 loops each) # 1.26 µs ± 11.7 ns per loop (mean ± std. dev. of 10 runs, 1000000 loops each) # 0.8358631635653873 \sourceoff Gruß

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.14, eingetragen 2021-10-14

\quoteon(2021-10-14 08:10 - MontyPythagoras in Beitrag No. 9) Hallo AlphaSigma, danke schon einmal für Deine Mühe. Das Beispiel von Dir wirft aber genau die Fragen und Zweifel auf, die ich schon angedeutet hatte. Hier ist die Berechnung mit dem numerischen Algorithmus augenscheinlich schneller. Aber ist sie das wirklich? Wenn man DEREN Implementation der cardanischen Formel verwendet, war das zu erwarten. Sie haben ja schließlich veröffentlicht, dass der Näherungsalgorithmus schneller sein soll, und man kann davon ausgehen, dass sie den Algorithmus auch perfekt optimiert haben. Aber haben sie das auch für die Umsetzung der cardanischen Formeln getan? Ich denke nicht. \quoteoff Ohne Beleg würde ich erst einmal davon ausgehen, dass sie die Fragestellung objektiv nach den Regeln der Wissenschaft bearbeitet haben. \quoteon Nachfolgende Fragen und Kritikpunkte drängen sich auf: 1. Die von Dir berechnete Lösung rechnet mit den Koeffizienten 1 3 3 1. Das ist ein Trivialfall, nämlich dreifache Nullstelle. Berechne ich den Wendepunkt, was der erste Schritt in deren Algorithmus ist, habe ich, oh Wunder, auch direkt die Lösung. Die iterative Schleife, die die eigentliche Zeit kostet, wird nicht ein einziges Mal durchlaufen. \quoteoff Dieses Beispiel habe ich als erstes spontan aus den 7 Testfällen ausgewählt. Dafür können die Autoren nichts. \quoteon 2. Die Koeffizienten sind in dem Beispiel reell, die Lösungen auch. Funktioniert der Algorithmus überhaupt mit komplexen Koeffizienten und komplexen Lösungen? \quoteoff Die Autoren wollen Aufgaben aus der Thermodynamik lösen (siehe in Beitrag No. 0 verlinkten Artikel). Die Koeffizienten sind reell. Die Lösungen können komplex sein. Die Autoren sind für ihre Anwendung hpts. an Gleichungen mit rellen Lösungen interessiert. Sie haben ihren Artikel und das Programm nicht für den Wikipedia-Artikel geschrieben. \quoteon 3. Wie haben sie die cardanischen Formeln umgesetzt? Etwa wie in den Wikipediaartikeln zur kubischen Gleichung (hier) oder den Cardanischen Formeln (hier) beschrieben, mit den ganzen Fallunterscheidungen und $\sin$ und $\cos$ und $\cosh$ und so weiter? Wo man $\Delta=\left(\frac{q}2\right)^2 + \left(\frac{p}3\right)^3 = \frac{27 A^2 D^2+4 B^3D-18 A B C D+4 A C^3-B^2 C^2}{108 A^4}$ berechnet, und die drei Lösungen im komplexen Fall durch Multiplikation mit den dritten Einheitswurzeln berechnet? Geht es denn noch uneffektiver? Dann muss ich mich nicht wundern, wenn die cardanische Lösung langsam ist. \quoteoff Die Implementierung der Cardanoformel in ihrem Mathepaket (nicht zwingend die gleiche wie im Artikel, kann aber sein) kann man sich von der mathc Seite herunterladen (ist ein tar.gz, erfordert tar und gzip). Ich denke es ist ok, wenn ich das hier einstelle. \sourceon C /* Quelle: http://www.uni-koeln.de/deiters/math/ */ /* eqn_cubic: real roots of a cubic equation */ #if METHOD == 0 // Cardano's method int eqn_cubic(const Polynomial& a, Vector& x) noexcept { int i, nreal; double w, h, p, q, d, phi, s, c; const double w3 = sqrt(3.0); if (a.size() < 4 || a[3] == 0.0) return eqn_quadratic(a, x); w = a[2] / (3.0 * a[3]); // x = y - w, y^3 + b1 y + b0 = 0 p = cb(a[1] / (3.0 * a[3]) - sq(w)); // (b1/3)^3 q = -0.5 * (2.0 * cb(w) - (a[1] * w - a[0]) / a[3]); // -b0/2 d = sq(q) + p; // discriminant if (d < 0.0) { // 3 real solutions h = q / sqrt(-p); phi = acos(fmax(-1.0, fmin(1.0, h))); sincos(rec3 * phi, &s, &c); p = 2.0 * pow(-p, 1.0 / 6.0); x[0] = 0.5 * p * (-c - w3 * s) - w; x[1] = 0.5 * p * (-c + w3 * s) - w; x[2] = 0.5 * p * 2.0 * c - w; nreal = 3; } else { // only one real solution d = sqrt(d); x[0] = cbrt(q + d) + cbrt(q - d) - w; nreal = 1; } for (i = 0; i < nreal; i++) { // 1 Newton iteration step if ((h = a[1] + x[i] * (2.0 * a[2] + 3.0 * x[i] * a[3])) != 0.0) x[i] -= (a[0] + x[i] * (a[1] + x[i] * (a[2] + x[i] * a[3]))) / h; } #if 1 // should be superfluous, but sometimes good for closely spaced roots if (nreal == 3) { // sort results if (x[1] < x[0]) std::swap(x[0], x[1]); if (x[2] < x[1]) std::swap(x[1], x[2]); if (x[1] < x[0]) std::swap(x[0], x[1]); } #endif return nreal; } \sourceoff \quoteon 4. Wenn ich das richtig sehe, wird das relative Abbruchkriterium für den Näherungsalgorithmus bei $10^{-6}$ gesetzt. Das ist schon sehr grob, könnte aber aufgrund der schnellen, weil kubischen Konvergenz des Halley-Verfahrens ausreichen, um die Lösungen auf $10^{-18}$ genau zu treffen. Könnte! Es könnte aber auch Situationen geben, wo die Konvergenz ungünstig ist. Ich würde hier mal das Konvergenzkriterium noch weiter runterschrauben, um sicherzustellen, dass die maximale, mit "float" erreichbare Genauigkeit auch sicher erreicht ist, oder man müsste noch eine Prüfung durch Einsetzen der Nullstelle hinzufügen, ob der Funktionswert an dieser Stelle tatsächlich nahe null ist. Daher habe ich meinen Code hier mal reingeschrieben, wo ich eine wirklich effiziente Umsetzung der cardanischen Formel angestrebt habe. Keine Fallunterscheidungen, keine riesigen Formeln wie die oben gezeigte, das wäre alles dramatisch uneffektiv. 9 Zeilen Code mit einer Handvoll algebraischen Anweisungen - und das war's. Ich habe daher arge Zweifel, dass die cardanische Lösung von den Machern des Algorithmus optimiert wurde, und daher würde ich gerne meinen Code umgesetzt sehen. Ich übersetze es auch gerne selbst in C++. Wir sollten außerdem ein Beispiel wählen, bei dem mit komplexen Koeffizienten drei Lösungen berechnet werden, die ebenfalls komplex sind und keine trivialen Nullstellen enthalten. \quoteoff In Table 4 in ihrem Artikel vergleichen die Autoren (nicht die von Wikipedia) ihre Methode (Deiters) mit der Cardano-Formel für drei Fälle: Three Real Roots (3r+2x), OneReal Root and Two Extrema (1r+2x), or One Real Root and No Extrema (1r+0x). Ich habe inzwischen alle 7 Testbeispiele für eqn_cubic mit "method 0 (Cardano)" und "method 1 (Deiters)" für jeweils 1000000000 Wiederholungen berechnet. Am Anfang eine Zusammenfassung mit Zeiten für eine Iteration. Danach eine Liste mit den 7 Datensätzen data1 bis data7, jeweils Eingabedatei und Ausgabedatei. Links Cardano und rechts Deiters. \sourceon method 0 method 1 Cardano Halley (Deiters) data1 59 ns 38 ns data2 120 ns 64 ns data3 121 ns 62 ns data4 133 ns 29 ns data5 31 ns 12 ns data6 87 ns 13 ns data7 133 ns 42 ns ==================================================== Cardano Halley (Deiters) datai nsteps a0 a1 a2 a3 x1 p(x1) x1 p(x1) x2 p(x2) x2 p(x2) ... ... CPU time: CPU time: ===================================================== data1 1000000000 1 0 0 1 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 CPU time: 58579 ms CPU time: 38227 ms data2 1000000000 -0.018 0.27 -1 1 1.00000000e-01 -3.4694e-18 1.00000000e-01 3.4694e-18 3.00000000e-01 0.0000e+00 3.00000000e-01 0.0000e+00 6.00000000e-01 -3.4694e-18 6.00000000e-01 -3.4694e-18 CPU time: 119554 ms CPU time: 63813 ms data3 1000000000 2 -1 -2 1 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 1.00000000e+00 0.0000e+00 1.00000000e+00 0.0000e+00 2.00000000e+00 0.0000e+00 2.00000000e+00 0.0000e+00 CPU time: 121193 ms CPU time: 62451 ms data4 1000000000 -40 +44 -14 1 1.99999998e+00 0.0000e+00 2.00000000e+00 0.0000e+00 2.00000002e+00 0.0000e+00 2.00000000e+00 0.0000e+00 1.00000000e+01 0.0000e+00 1.00000000e+01 0.0000e+00 CPU time: 132939 ms CPU time: 28703 ms data5 1000000000 1 3 3 1 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 -1.00000000e+00 0.0000e+00 CPU time: 30937 ms CPU time: 11788 ms data6 1000000000 -1.000003000002 3.000006000002 -3.000003 1.0 9.99778975e-01 -1.0944e-11 1.00000000e+00 2.2204e-16 9.99888963e-01 -1.4062e-12 1.00000100e+00 -2.2204e-16 1.00011304e+00 1.4064e-12 1.00000200e+00 2.2204e-16 CPU time: 86814 ms CPU time: 12977 ms data7 1000000000 -1.0e+21 1.00000010000001e+21 -1.00000010000001e+14 1.0e+00 -5.83634661e+00 -6.8364e+21 1.00176003e+00 1.7600e+18 1.00000068e+07 -6.8363e+21 1.00000000e+07 -1.6384e+11 1.00000000e+14 -3.8528e+18 1.00000000e+14 -3.8528e+18 CPU time: 133139 ms CPU time: 42168 ms \sourceoff Verwendete CPU: \sourceon $ cat /proc/cpuinfo processor : 0 vendor_id : GenuineIntel cpu family : 6 model : 58 model name : Intel(R) Core(TM) i3-3240 CPU @ 3.40GHz stepping : 9 microcode : 0x21 cpu MHz : 1600.087 cache size : 3072 KB \sourceoff \quoteon ... ... In einem Wikipedia-Artikel sollte nichts drinstehen, was nicht unabhängig überprüft werden kann. Ciao, Thomas \quoteoff Der Artikel ist ja verlinkt und man kann ihn frei von der Seite des Autors herunterladen. Der Quellcode des mathc Paketes ist auch frei verfügbar. Die Wikipedia-Autoren können auch nur auf das verweisen, was sie im Internet finden.

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.15, vom Themenstarter, eingetragen 2021-10-14

Hallo schnitzel und AlphaSigma, vielen Dank für Eure Mühe. @schnitzel: ich konnte mir auch nicht im Ernst vorstellen, dass das Multiplizieren nicht gravierend schneller sein soll als das Potenzieren. Allerdings war meine Methodik, die Schleifendauer zu messen, nicht so profimäßig wie Deine. Da habe ich direkt was wichtiges gelernt, vielen Dank. @AlphaSigma: Ich wollte nicht unterstellen, dass die Autoren die cardanische Lösung bewusst schlecht umgesetzt haben. Was hätten sie davon? Es ist ja kein Wettrennen, wo es irgendwas zu gewinnen gäbe. Am Ende haben sie ja ein Interesse daran, ihre Berechnungen schnell umzusetzen, ganz egal, welche Kopfstände der Prozessor dafür machen muss. Ein paar Dinge fallen mir auf: 1. Wenn der Algorithmus keine kubischen Gleichungen mit komplexen Koeffizienten lösen kann, ist der Vergleich schon ein wenig schief. Die hier gezeigte Umsetzung der cardanischen Gleichungen sind mit reellen Koeffizienten programmiert, und es wird tatsächlich das Kochrezept mit Fallunterscheidung $\Delta<0$ und $\Delta>0$ usw. abgearbeitet. Möglicherweise haben sie nicht wirklich darüber nachgedacht, eine komplexe Lösung umzusetzen, weil es sie für ihren Anwendungszweck nicht interessierte. Für allgemeine Lösungen taugt der Algorithmus daher vielleicht gar nicht, oder man müsste ihn erst auf komplex umprogrammieren. 2. Die Umsetzung der cardanischen Lösung, die Du hier hineinkopiert hast, ist tatsächlich sehr ineffektiv. Vergleiche das mal mit meinem Code. Hier wird der Sinus und Cosinus berechnet, der Arkuscosinus, es wird potenziert, Fälle werden unterschieden, im Fall $d>0$ wird zweimal eine dritte Wurzel gezogen, wo auch einmal reicht, und am Ende wird sogar noch eine Newton-Schleife hinterhergeschoben, vermutlich um die Präzision der Lösung zu erhöhen und Rundungsfehler zu reduzieren. Und ganz am Ende sortieren sie noch die Ergebnisse. Es war vermutlich keine Absicht der Programmierer, aber mein Verdacht hat sich erhärtet, dass die Umsetzung der cardanischen Lösung hier wirklich schlecht ist. 3. Nochmal das Thema Abbruchkriterium. Ich denke, der Halley-Algorithmus bricht zu früh ab. Es wäre schön, noch einen Vergleich zu sehen, wie genau die Null wirklich getroffen wird. Möglicherweise muss der Algorithmus noch eine Schleife mehr drehen, um die gleiche Genauigkeit zu erreichen. Setze also mal das Abbruchkriterium auf $10^{-15}$ statt $10^{-6}$. Ich übersetze mal meine Lösung in C++, setze Du das Abbruchkriterium im Halley-Algorithmus mal auf $10^{-15}$, und dann vergleichen wir nochmal. 🙂 Ciao, Thomas

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.16, eingetragen 2021-10-14

\quoteon(2021-10-14 20:39 - MontyPythagoras in Beitrag No. 15) ... @AlphaSigma: ... ... Setze also mal das Abbruchkriterium auf $10^{-15}$ statt $10^{-6}$. Ich übersetze mal meine Lösung in C++, setze Du das Abbruchkriterium im Halley-Algorithmus mal auf $10^{-15}$, und dann vergleichen wir nochmal. 🙂 Ciao, Thomas \quoteoff Hallo Thomas, erst einmal habe ich mehr Nachkommastellen (16) ausgeben lassen. method 1 (Deiters): \sourceon data1 -1.0000000000000000e+00 0.0000000000000000e+00 CPU time: 40763 ms data2 1.0000000000000002e-01 3.4694469519536142e-18 3.0000000000000010e-01 0.0000000000000000e+00 5.9999999999999987e-01 -3.4694469519536142e-18 CPU time: 64151 ms data3 -1.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00 2.0000000000000000e+00 0.0000000000000000e+00 CPU time: 62523 ms data4 2.0000000000000000e+00 0.0000000000000000e+00 2.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+01 0.0000000000000000e+00 CPU time: 28783 ms data5 -1.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 CPU time: 11786 ms data6 9.9999999973354070e-01 2.2204460492503131e-16 1.0000009999999999e+00 -2.2204460492503131e-16 1.0000020002664596e+00 2.2204460492503131e-16 CPU time: 12951 ms data7 1.0017600255087018e+00 1.7600253323895112e+18 1.0000000000000000e+07 -1.6384000000000000e+11 1.0000000000000000e+14 -3.8528000000000000e+18 CPU time: 42329 ms \sourceoff Bei data1 bis data6 ist der Wert des Polynoms an der berechneten Nullstelle jeweils kleiner oder gleich 2.22e-16 (eps double). Die Nullstellen selber sind nicht alle auf 15 Nachkommastellen genau. Das Testprogramm gibt übrigens die Nullstelle xi und den Wert des Polynoms p(xi) aus und nicht Re(xi) und Im(xi), wie ich zuerst angenommen hatte. Wahrscheinlich rechnet es rein reell. Den Quellcode habe ich mir noch nicht genauer angesehen. In Fußnote #3 in supplement.pdf - Cubic rootfinder using Halley’s method schreiben die Autoren "Halley’s method usually has a convergence order of 3. ... If two subsequent iteration steps differ by less than $10^{−6}$(relatively), the termination error of the last step is less than $10^{−18}$."

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.17, eingetragen 2021-10-15

\quoteon(2021-10-13 21:23 - MontyPythagoras in Beitrag No. 3) ...so, ich habe noch einmal gute 7% rausgequetscht. Divisionen durch Multiplikationen mit gespeichertem Kehrwert ersetzt, Anzahl der arithmetischen Operationen generell reduziert. \sourceon Python \numberson def cbrt(a,b,c,d): # Trivialfall quadratische Gleichung: if a==0: # Trivialfall lineare Gleichung: if b==0: return [-d/c] # Echte quadratische Gleichung: y1=-(c+(c**2-4*b*d)**0.5)/(2*b) return [y1,-c/b-y1] # Echte kubische Gleichung: x=y+r -> y^3+3p-2q=0 s=1/(3*a) r=-b*s p=c*s-r**2 q=-1.5*s*(d+c*r)+r**3 if p==q==0: return [r,r,r] # Dreifach-Nullstelle s=(q**2+p**3)**0.5+q q=(2*q)**(1/3) if s==0 else s**(1/3) y1=q-p/q y2=-0.5*(y1+(-3*y1**2-12*p)**0.5) return [y1+r,-y1-y2+r,y2+r] \sourceoff Ciao, Thomas \quoteoff Hallo MontyPythagoras, ich bin gerade dabei Deinen Python Code nach C zu übersetzen. Was ist das in Zeile 16 in Deinem zweiten Programm für ein Konstrukt? Soll q gleich dem ersten Term gesetzt werden, falls s == 0 und sonst gleich s**1/3? Also analog dem ? : Operator in C? \quoteon(2021-10-13 21:03 - MontyPythagoras in Beitrag No. 2) Hallo AlphaSigma, danke für Deine Mühe! Ich bin dabei, den Code noch in leichten Variationen zu testen und die letzten Nanosekunden aus Python rauszuquetschen. Was der C++ Compiler später daraus macht, ist sicher noch einmal eine andere Nummer, aber mein i7-Prozessor in meinem Surface Pro schafft derzeit eine Zeit von 2,24µs. Ich lasse $2^{22}$ mal die gleiche Lösung berechnen, das dauert etwa 9 Sekunden. Ciao, Thomas \quoteoff Warum verwendest Du eigentlich Python mit Rechenzeiten von einigen µs für eine Iteration?

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.18, vom Themenstarter, eingetragen 2021-10-15

Hallo AlphaSigma, \quoteon(2021-10-15 18:07 - AlphaSigma in Beitrag No. 17) Warum verwendest Du eigentlich Python mit Rechenzeiten von einigen µs für eine Iteration? \quoteoff Kein besonderer Grund. Ich weiß, dass Python eine eher langsame Sprache ist. Selbst Visual Basic ist etwa 3-4 mal schneller. Aber sie ist sehr intuitiv und gut lesbar, die Lernkurve ist steil, sie kann komplexe Zahlen und es gibt reichlich fertige Lösungen wie scipy usw. mit vielen mathematischen Funktionen. Ich habe früher mal C++ programmiert, aber das ist noch im letzten Jahrtausend gewesen. Ich habe vor, eventuell auch wieder in C++ einzusteigen. Was ich damals konnte, kriege ich heute auch wieder hin. Ich habe in den 90ern auch viel Geschwindigkeitsoptimierungen programmiert, aber dann war ich mit dem Studium fertig. 😂 Die Zeile 16 hast Du komplett richtig verstanden. \quoteon(2021-10-14 22:04 - AlphaSigma in Beitrag No. 16) In Fußnote #3 in supplement.pdf - Cubic rootfinder using Halley’s method schreiben die Autoren "Halley’s method usually has a convergence order of 3. ... If two subsequent iteration steps differ by less than $10^{−6}$(relatively), the termination error of the last step is less than $10^{−18}$." \quoteoff Ich weiß, das habe ich gelesen, und es ist auch logisch. Aber kubische Konvergenz ist keine Garantie. Wenn die Nullstellen ungünstig liegen, z.B. wenn in der kubischen Gleichung zwei Lösungen sehr dicht bei einander liegen. Dann taucht der parabelförmige Teil nur sehr knapp unter die x-Achse, und dann ist die Konvergenz langsamer. Daher ist es sehr optimistisch, die Schleife bei $10^{-6}$ abzubrechen und nur zu glauben, dass die anvisierte Genauigkeit der Nullstelle auf $10^{-18}$ schon erreicht wurde. Wenn das eine sicherheitskritische Anwendung wäre, wäre das gefährlich. Und der Beweis, dass der Algorithmus tatsächlich zu nachlässig programmiert wurde, findet sich in Deinem Beitrag #16, im Datensatz 7. Die Lösung wurde dramatisch verfehlt, denn statt dass der Funktionswert null wäre, beträgt er bei zwei Lösungen über $10^{18}$. Wenn ich numerische Lösungen programmiere, dann breche ich erst ab, wenn ich sicher bin, dass ich die gewünschte Genauigkeit erreicht habe. Das würde hier aber bedeuten, dass der Halley-Algorithmus schon einmal mindestens eine Schleife mehr drehen müsste, und dann sieht der Zeitvergleich schon anders aus. Und dass die Leute der Uni Köln bei der cardanischen Lösung noch eine Newtonschleife hintendran basteln, haben sie auch großzügig unterschlagen. Ciao, Thomas

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.19, eingetragen 2021-10-15

\quoteon(2021-10-15 21:30 - MontyPythagoras in Beitrag No. 18) Und dass die Leute der Uni Köln bei der cardanischen Lösung noch eine Newtonschleife hintendran basteln, haben sie auch großzügig unterschlagen. Ciao, Thomas \quoteoff Da müssen wir vorsichtig sein. Die oben angegebene Cardano-Routine stammt aus dem mathc Paket der Autoren. Ob sie das exakt so für den Artikel berechnet haben, wissen wir nicht. Bei der Transformation Deines Python Codes in C erhalte ich für Datensatz 3 mit den Koeffizienten d = 2, c = -1, b = -2, a = 1 komplexe Zwischenwerte, obwohl die Nullstellen reell sind. Entweder habe ich noch einen Fehler oder die Vereinfachung auf rein reelle Gleichungen ist nicht so einfach.

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.20, vom Themenstarter, eingetragen 2021-10-15

Hallo AlphaSigma, nein, das ist doch genau richtig, denn das ist ja das geniale an den komplexen Zahlen: reelle Koeffizienten, reelle Lösungen, aber der elegante Weg führt über komplexe Zwischenwerte. Wenn man das nicht will, dann muss man das ganze Fallunterscheidungsbrimborium machen wie im Wikipedia-Artikel beschrieben mit Casus irreducibilis und so weiter. Und dann wird es eben langsam. Ich denke, Du bist genau auf dem richtigen Weg. 🙂 Ciao, Thomas

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.21, eingetragen 2021-10-16

Hallo MontyPythagoras, ich habe jetzt eine erste lauffähige Version von Deinem Python Code in C, nicht C++, implementiert. Sie benutzt den Header complex.h und die Funktionen csqrt und cpow. Die Variablen r, s, p, q, y1, y2 sind alle komplex. Relle Konstanten c habe ich mit CMPLX(c, 0.0) in komplexe Konstanten umgewandelt. Für die Datensätze 1 bis 5 liefert das Programm auch die korrekten Ergebnisse. Es ist aber deutlich langsamer als die cardano Routine aus mathc in Beitrag No. 14. cardano1 (mathc, siehe Beitrag No. 14) und cardano2 (erster Versuch Python nach C complex): \sourceon cardano1 cardano2 data1 59 ns 196 ns data2 120 ns 293 ns data3 121 ns 297 ns data4 133 ns 260 ns data5 31 ns 31 ns \sourceoff Bevor ich den Code hier veröffentliche, schaue ich ihn mir lieber noch einmal gründlich an und versuche den Flaschenhals zu finden.

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.22, eingetragen 2021-10-16

Hallo, hier noch der Link auf eine interaktive Internetseite zur Berechnung der Nullstellen nach Cardano mit einem vergleichbaren Ansatz wie im Beitrag No. 3 von MontyPythagoras. Formel von Cardano (mit Formular zur Demo)

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.23, vom Themenstarter, eingetragen 2021-10-16

Hallo AlphaSigma, das ist auf jeden Fall ein spannendes und überraschendes Ergebnis. Bin gespannt auf Deinen Code. Wenn sich das bewahrheitet, dann könnte der Unterschied dadurch erklärbar sein, dass bei komplexen Zahlen der Rechenaufwand für den Prozessor steigt. Eine komplexe Multiplikation bringt z.B. mehr Aufwand mit sich als eine reelle Multiplikation. Letztere ist einfach eine einzige Multiplikation, aber eine komplexe Multiplikation bedeutet ja $$(a+bi)(c+di)=(ac-bd)+(bc+ad)i$$Das sind in Wirklichkeit 4 reelle Multiplikationen und 2 reelle Additionen (die Addition zwischen den Klammern ist irrelevant). Das könnte dazu führen, dass eine rein reelle Lösung selbst mit den Fallunterscheidungen unter dem Strich doch weniger Rechenaufwand bedeutet als mit komplexen Zahlen zu rechnen. Die von mir angewandte Cardano-Lösung mit komplexen Zahlen dürfte dementsprechend ihren Vorteil nur dann ausspielen können, wenn man allgemeine Lösungen für komplexe Koeffizienten sucht. Sehr spannend! Ciao, Thomas

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.24, eingetragen 2021-10-17

Da Halley’s eigentliche Methode nur etwa 6 Stellen genau zu sein scheint, ist das höchstens float-Genauigkeit statt double. Die nachträgliche Newton-Schleife bringt erst die eigentliche double-Genauigkeit. Die eqn_cubic, die Cardano's method darstellen soll, ist laut Beitrag 14: - der (ur-)alte Algorithmus mit zig Fallunterscheidungen und vielen langsamen trigonometrischen Funktionen - wie schon richtig festgestellt wurde, verfälscht die nachträgliche Newton-Iteration Aussagen zur Genauigkeit und Geschwindigkeit - sincos scheint auch noch über den veralteten & langsamen Coprozessorbefehl fsincos zu laufen Wie MontyPythagoras und und https://www.mathematik.ch/anwendungenmath/Cardano/FormelCardano.html zeigten, gibt es längst optimierte "neue Formeln von Cardano" mit komplexen Variablen. Ich habe sie damals unter https://www.lamprechts.de/gerd/php/gleichung-6-grades.php PQRST-Formel genannt. Sie kommt ohne trigonometr. Funktionen & ohne Fallunterscheidungen aus: https://www.lamprechts.de/gerd/Bilder/QuadratischeGleichung_p-q-Formel_KubischeGleichung_PQRST-Formel.png Da hier komplexe Zwischenergebnisse vorkommen, habe ich die Algorithmen zunächst mit mathematica statt mit c++ verglichen. Statt winzige Genauigkeitswerte (die dann auch oft {data1...data7} noch glatte ganze Zahlen wie 0 und 1 sind!!) zu vergleichen, wo Compiler, Cache, Maschinenbefehl, Messgenauigkeit eine größere Rolle spielen als der Algorithmus selbst, habe ich mal 401 unterschiedliche krumme Parameter ganommen, die dann für alle 3 Lösungen {x1,x2,x3} auf je 1000 Stellen genau berechnet werden. (also 401 *3* 1000 = 1203000 Ziffern pro Algorithmus und noch Summe aller Werte zur Kontrolle) Zeitmessungen sind mit echter AbsoluteTime-Ausgabe und Aufsummierung der Ergebnissummen realistischer, als die geschummelte //Timing Funktion von Mathematica (da wird nur der interne hex-Wert gemessen, bevor die Konvertierung in Dezimalzahl erfolgt!) \sourceon 401 Gleichungen * 3 Ergebnisse je 1000 Nachkommastellen Solve NSolve Monty PQRST FindRoot Newton {0.2499413,0.3592909,14.5747157,1.1403582,1.1716010,0.3749123} Werte in Sekunden Korrektur, da Zwischenergebnisse zu genau berechnet waren: {0.2499415, 0.3592910,0.0781067,0.0937281,1.1403582,0.3749123} in s \sourceoff Zunächst dachte ich, dass FindRoot die Newton-Iteration verwendet. Da der Messwert jedoch zu langsam ist, habe ich eine eigene Newton-Iteration hinten angehangen -> die dann etwa so schnell wie NSolve (also reine numerische Berechnung) ist. Mathematicas Solve scheint für die komplexe 3. Wurzel hoch optimiert zu sein. PQRST & Monty verwenden ja eine komplexe 3. Wurzel und 2 quadratische Wurzeln. Mit Komplexen Zahlen ist so eine Berechnung (meist auch Iteration) natürlich langsamer als die reelle Newton-Iteration. Für Kubische Funktionen ohne quadratischen Anteil (Parameter b=0) gibt es auch noch hypergeometrische Algorithmen, die ich auch noch auf Geschwindigkeit testen werde...

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.25, vom Themenstarter, eingetragen 2021-10-17

Hallo HyperG, Deine PQRStTW-Formel bedarf einer höheren Anzahl an Multiplikationen oder allgemein an Rechenoperationen als meine Programmierung. Unter dem Strich tun beide ungefähr das gleiche. Trotzdem soll Deine Variante fast 13 mal schneller sein? Das ist unglaubwürdig. Worin soll der Vorteil liegen? Außerdem ist nicht interessant, was Mathematica tut. Mathematica ist eine Anwendung, keine Programmiersprache. Es geht hier aber um den Geschwindigkeitsvergleich in der standardisierten "double precision" innerhalb einer Programmiersprache, also bei Zahlen, die ungefähr 16 Dezimalstellen Genauigkeit haben, und nicht 1000. Darauf bezieht sich der Absatz im Wikipediaartikel. und nur darum geht es mir hier. Hilfreich wäre, wenn Du die PQRST-Formel in C++ umsetzt und dann einen Geschwindigkeitsvergleich machst. Es sollte einleuchten, dass nur dann überhaupt ein Näherungsalgorithmus, ganz gleich ob Newton oder Halley, eine Chance gegen eine explizite Lösung haben kann, wenn die Anzahl der verlangten Stellen relativ gering ist. Der Halley-Algorithmus hat kubische Konvergenz, was bedeutet, dass er die Zahl der korrekten Nachkommastellen mit jeder Schleife etwa verdreifacht. Während er bei "double precision" also oft mit ein oder zwei Schleifen auskommen mag, braucht er bei 1000 Nachkommastellen normalerweise mindestens 6 und bei 1 Mio. Nachkommastellen mindestens 12, wenn er mit Zahlen rechnet, die die gleiche Anzahl an Stellen aufweisen. Eine exakte Formel braucht aber immer nur genau einen Durchlauf, und wenn man mit 1000-stelligen Zahlen rechnen würde, würde die cardanische Lösung nach diesem einen Durchlauf halt sofort ein Ergebnis liefern, welches auf 1000 Stellen genau stimmt, minus ein paar Stellen Rundungsfehler. Ich trete jetzt erst einmal in einen einwöchigen Funkschatten. Bin gespannt, was in der Zwischenzeit noch rauskommt. Ciao, Thomas

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.26, eingetragen 2021-10-17

\quoteon(2021-10-17 07:39 - MontyPythagoras in Beitrag No. 25) ... Trotzdem soll Deine Variante fast 13 mal schneller sein? Das ist unglaubwürdig. Worin soll der Vorteil liegen? ... \quoteoff Zwischenergebnisse waren noch zu genau berechnet -> behoben und in Tabelle korrigiert \quoteon(2021-10-17 07:39 - MontyPythagoras in Beitrag No. 25) Eine exakte Formel braucht aber immer nur genau einen Durchlauf, und wenn man mit 1000-stelligen Zahlen rechnen würde, würde die cardanische Lösung nach diesem einen Durchlauf halt sofort ein Ergebnis liefern, ... \quoteoff Eben nicht 1 Durchlauf: selbst bei double benötigen Maschinenbefehle iterne Schleifen. Selbst MUL braucht mehr als 1 Takt. Zwar wird bei Sqrt und 3. Wurzel mit Näherungsformeln gearbeitet, aber je nach Befehlssatz können das 12 bis 300 Takte sein. Deshalb ist ja das Variieren der Parameter wichtig. Bei 1000 Stellen ist selbst eine FFT MUL schneller als "normale MUL". Die 2. & 3. Wurzel wird meist mit der Newton-Iteration berechnet, da die Genauigkeit bekanntermaßen quadratisch ansteigt. Ob nun 2 Iterationen für die explizite Variante, oder 1 Iteration x = x - (d + x (c + x (b + a x)))/(c + x (2 b + 3 a x)) für die normale Newton-Iteration schneller ist, oder sogar hypergeometrische Funktionen, das ist die jedenfalls meine Hauptfrage. Den Vorteil bei Halley sehe ich darin, dass schnell ein Startwert für diese Iterationen gefunden wird, ohne sich mit komplexen Zwischenwerten herumschlagen zu müssen. Grüße

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.27, eingetragen 2021-10-17

\quoteon(2021-10-16 23:46 - MontyPythagoras in Beitrag No. 23) Hallo AlphaSigma, das ist auf jeden Fall ein spannendes und überraschendes Ergebnis. Bin gespannt auf Deinen Code. ... Sehr spannend! Ciao, Thomas \quoteoff Hallo MontyPythagoras, ich habe noch einige Tests mit complex.h gemacht. \sourceon C double c; double complex z; c*z oder CMPLX(c, 0.0)*z cpow(z, c) oder cpow(z, CMPLX(c, 0.0)) \sourceoff ergeben scheinbar keinen großen Zeitunterschied. Ich habe das Testprogramm ggü. dem in Beitrag No. 21 verwendeten geringfügig modifiziert und erhalte diese leicht verbesserten Zeiten. cardano1 (mathc, siehe Beitrag No. 14) und cardano2 (zweiter Versuch Python nach C complex): \sourceon cardano1 cardano2 data1 59 ns 118 ns data2 120 ns 207 ns data3 121 ns 211 ns data4 133 ns 178 ns data5 31 ns 23 ns data6 87 ns 59 ns data7 133 ns 191 ns ======================================================================= data1 1000000000 1 0 0 1 5.0000000000000011e-01 8.6602540378443860e-01 i -1.1102230246251565e-16 3.8857805861880479e-16 i 4.9999999999999989e-01 -8.6602540378443882e-01 i -3.3306690738754696e-16 6.1062266354383610e-16 i -1.0000000000000000e+00 1.6653345369377348e-16 i 5.5466782398352393e-32 4.9960036108132044e-16 i CPU time: 118152 ms data2 1000000000 -0.018 0.27 -1 1 6.0000000000000009e-01 2.0816681711721685e-17 i 2.0816681711721685e-17 3.1225022567582559e-18 i 2.9999999999999982e-01 -5.5511151231257827e-17 i 1.3877787807814457e-17 3.3306690738754672e-18 i 1.0000000000000012e-01 3.1225022567582602e-17 i 1.7347234759768071e-17 3.1225022567582555e-18 i CPU time: 207440 ms data3 1000000000 2 -1 -2 1 2.0000000000000000e+00 -1.1102230246251565e-16 i -2.4651903288156619e-32 -3.3306690738754696e-16 i 9.9999999999999989e-01 2.2204460492503131e-16 i 2.2204460492503131e-16 -4.4408920985006262e-16 i -9.9999999999999989e-01 -5.5511151231257852e-17 i 6.6613381477509392e-16 -3.3306690738754706e-16 i CPU time: 211127 ms data4 1000000000 -40 +44 -14 1 9.9999999999999964e+00 3.3087224502121107e-24 i -2.2737367544323206e-13 2.1175823681357471e-22 i 2.0000000596046452e+00 -2.2204460492503131e-16 i -2.8421709430404007e-14 2.1175823602471418e-22 i 1.9999999403953561e+00 2.2204460327066996e-16 i -2.8421709430404007e-14 2.1175823365813146e-22 i CPU time: 178468 ms data5 1000000000 1 3 3 1 -1.0000000000000000e+00 0.0000000000000000e+00 i 0.0000000000000000e+00 0.0000000000000000e+00 i -1.0000000000000000e+00 0.0000000000000000e+00 i 0.0000000000000000e+00 0.0000000000000000e+00 i -1.0000000000000000e+00 0.0000000000000000e+00 i 0.0000000000000000e+00 0.0000000000000000e+00 i CPU time: 23045 ms data6 1000000000 -1.000003000002 3.000006000002 -3.000003 1.0 -nan -nan i -nan -nan i -nan -nan i -nan -nan i nan nan i nan nan i CPU time: 58557 ms data7 1000000000 -1.0e+21 1.00000010000001e+21 -1.00000010000001e+14 1.0e+00 1.0000000000000000e+14 1.1175870895385742e-08 i -3.8528000000000000e+18 1.1175869777798540e+20 i 1.0028794335937500e+07 -1.1111880186945200e-01 i -2.8877241540234207e+25 1.1175869793943296e+20 i -2.8793335937500000e+04 1.1111879051327003e-01 i -2.8877244444588056e+25 1.1175869782225963e+20 i CPU time: 190586 ms ======================================================================= data6mod 1000000000 -1.0 3.0 -3.0 1.0 1.0000000000000000e+00 -0.0000000000000000e+00 i 0.0000000000000000e+00 0.0000000000000000e+00 i 1.0000000000000000e+00 -0.0000000000000000e+00 i 0.0000000000000000e+00 0.0000000000000000e+00 i 1.0000000000000000e+00 -0.0000000000000000e+00 i 0.0000000000000000e+00 0.0000000000000000e+00 i CPU time: 22958 ms \sourceoff Die Ausgabe ist creal(xi) cimag(xi) i creal(p(xi)) cimag(p(xi)) i. Hier noch der Sourcecode. Das Programm ist die Antwort auf die Frage im Themenstart und Beitrag No. 3. Es arbeitet mit double precision, nicht mit arbitrary precision. Es ist nur für die Eingangsdaten data1 bis data7 getestet, insbesondere nur für rein reelle Polynomkoeffizienten. Zur Vergleichbarkeit mit den bisherigen Tests wurden wieder data1 bis data7 als Eingangsdaten verwendet. mp_cbrt_cmplx.c \sourceon C /*******************************************************************/ /* Source: https://www.matheplanet.de/matheplanet/nuke/html/viewtopic.php?topic=255773&post_id=1858001 */ /* Translation of cbrt(a,b,c,d) from Python to C */ /*******************************************************************/ /* Solve a*x^3 + b*x^2 + c*x + d = 0 (Python) */ /* Solve a[3]*x^3 + a[2]*x^2 + a[1]*x + a[0] = 0 (C) */ /*******************************************************************/ #include #include #include #include int mp_cbrt_cmplx(const double complex *ac, double complex *xc); int main(int argc, char **argv) { int i, n, nrep, t0, t1; double a[4]; double complex ac[4]; double complex xc[3]; scanf("%d", &nrep); for (i = 0; i < 4; ++i) { scanf("%lf", &a[i]); ac[i] = CMPLX(a[i], 0.0); } t0 = clock(); for (i = 0; i < nrep; i++) n = mp_cbrt_cmplx(ac, xc); t1 = clock(); for (i = 0; i < n; i++) { printf("%24.16e %24.16e i %24.16e %24.16e i\n", creal(xc[i]), cimag(xc[i]), creal(ac[0] + xc[i] * (ac[1] + xc[i] * (ac[2] + xc[i] * ac[3]))), cimag(ac[0] + xc[i] * (ac[1] + xc[i] * (ac[2] + xc[i] * ac[3])))); } printf("CPU time: %d ms\n", (t1-t0) / 1000); return 0; } /* returns number of solutions */ int mp_cbrt_cmplx(const double complex *ac, double complex *xc) { double complex y1, y2, p, q, r, s; const double deps = 3.0e-16; if (creal(ac[3]) == 0.0 && cimag(ac[3]) == 0.0) { /* Linear eqn */ if (creal(ac[2]) == 0.0 && cimag(ac[2]) == 0.0) { xc[0] = -ac[0] / ac[1]; return 1; } /* Quadratic eqn */ y1 = -(ac[1] + csqrt(ac[1]*ac[1] - CMPLX(4.0, 0.0)*ac[2]*ac[0])) / (CMPLX(2.0, 0.0)*ac[2]); xc[0] = y1; xc[1] = -ac[1]/ac[2] - y1; return 2; } /* Cubic eqn */ else { s = CMPLX(1.0, 0.0) / (CMPLX(3.0, 0.0)*ac[3]); r = -ac[2]*s; p = ac[1]*s - r*r; q = CMPLX(-1.5, 0.0)*s*(ac[0] + ac[1]*r) + r*r*r; /* Three-fold zero */ if ((cabs(p) < deps) && cabs(q) < deps) { xc[0] = xc[1] = xc[2] = r; return 3; } s = csqrt(q*q + p*p*p) + q; q = (cabs(s) < deps) ? cpow(CMPLX(2.0, 0.0)*q, CMPLX(1.0/3.0, 0.0)) : cpow(s, CMPLX(1.0/3.0, 0.0)); y1 = q - p/q; y2 = CMPLX(-0.5, 0.0)*(y1 + csqrt(-CMPLX(3.0, 0.0)*y1*y1 - CMPLX(12.0, 0.0)*p)); xc[0] = y1 + r; xc[1] = -y1 - y2 + r; xc[2] = y2 + r; return 3; } } \sourceoff Compiliert wurde mit: gcc -Wall -pedantic -Ofast -lm -o mp_cbrt_cmplx mp_cbrt_cmplx.c Prozessor wie gehabt: Intel(R) Core(TM) i3-3240 CPU @ 3.40GHz \quoteon(2021-10-15 21:30 - MontyPythagoras in Beitrag No. 18) Hallo AlphaSigma, ... Und der Beweis, dass der Algorithmus tatsächlich zu nachlässig programmiert wurde, findet sich in Deinem Beitrag #16, im Datensatz 7. Die Lösung wurde dramatisch verfehlt, denn statt dass der Funktionswert null wäre, beträgt er bei zwei Lösungen über $10^{18}$. Wenn ich numerische Lösungen programmiere, dann breche ich erst ab, wenn ich sicher bin, dass ich die gewünschte Genauigkeit erreicht habe. Das würde hier aber bedeuten, dass der Halley-Algorithmus schon einmal mindestens eine Schleife mehr drehen müsste, und dann sieht der Zeitvergleich schon anders aus. ... Ciao, Thomas \quoteoff Ich denke data6 und data7 dienen dazu die Grenzen des Programms zu zeigen. Bei data7 sind die Koeffizienten -1.0e+21 1.00000010000001e+21 -1.00000010000001e+14 1.0e+00, unterscheiden sich also um 21 Größenordnungen und gerechnet wird mit double precision, also eps = 2.22e-16. Wie der Vergleich von data6 mit data6mod oben zeigt, können kleine Abweichungen in den Koeffizienten große Auswirkung auf das Ergebnis haben. Ein Flaschenhals bei der Berechnung in C mit komplexen Zahlen scheint mir die kubische Wurzel cpow(x, 1.0/3.0) zu sein. \sourceon nameDerSprache double complex xc; double x; for (i = 0; i < nrep; i++) xc += cpow(xc,CMPLX(1.0/3.0, 0.0)); for (i = 0; i < nrep; i++) x += pow(x, 1.0/3.0); \sourceoff Der Schritt in der ersten for-Schleife mit cpow benötigt ca. 111 ns ggü. 35 ns in der zweiten mit der rellen pow-Funktion.

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.28, eingetragen 2021-10-20

\quoteon(2021-10-17 15:02 - AlphaSigma in Beitrag No. 27) ... Flaschenhals bei der Berechnung in C mit komplexen Zahlen scheint mir die kubische Wurzel cpow(x, 1.0/3.0) zu sein... cpow benötigt ca. 111 ns ggü. 35 ns in der zweiten mit der rellen pow-Funktion... \quoteoff Wenn 1 einzige (Unter-)Funktion etwa 111 ns benötigt, dann sind die data1, data5 & data6 aber sehr "geschönte Argumente", wenn der komplette Algorithmus mit zig weiteren Funktionen und 3 Ergebnissen teils unter 60 ns liegen... Das habe ich mal mit VC++ und einem i9 4GHz selbst nachgerechnet, jedoch mit float statt double, da durch Fehlerfortpflanzung und Nächerungsfunktionen eh die letzten Stellen falsch sind und durch einen "letzten Newton-Schritt" schnell zu korrigieren sind. Argument wandert in allen Tests von -18...22 mit Schrittweite 3e-8 (keine glatten Data wie in den Beiträgen zuvor); statt 1 Mrd. hier 1073741824 Durchläufe. \sourceon c++ pow(x,0.33333f):8.94 ns statt cpow habe ich (x+y*i)^(1/3)=(y^2 + x^2)^(1/6) (Cos[ArcTan[y/x]/3] + I Sin[ArcTan[y/x]/3]) benutzt: cSqrt3: 30.61 ns Interessant: sin: 6.26 ns cos: 0.24 ns 26 mal schneller?! Sin scheint sich durch die Asymmetrie nicht so gut optimieren zu lassen wie cos... Also Tests mit eigenen sin Nachbildungen, was im float-Bereich und AVX-Befehlen recht gut gelingt: sin256_ps: 1.05 ns Genau: 5.960e-08 sinAVX2 0.45 ns Genau: 1.505e-06 (Bereich -10...10 schon 1.6e-07 genau) Zwar bekommt man die Verbesserung nur, wenn man immer 8 sin parallel in ein AVX-Register schiebt, aber immerhin um etwa 5.2 ns schneller. atan2f 10.91 ns noch nicht optimiert \sourceoff Aber egal wie "explizit" der Algorithmus auch ist: durch die relativ langsamen atan, sin, cos, pow Maschinenbefehle (oder selbst Optimierungen), wird man Halley’s Näherungsfunktion schlecht überholen können. Anders sieht es aus, wenn man spezielle Prozessoren nutzt, die wie FPGAs nur immer die gleiche Anzahl an Takten für 1 Maschinenbefehl benötigen. Allerdings ist dann meist der Befehlssatz weiter eingeschränkt... Korrektur 20.10.2021: !!!!! #### !!!!! ############ Nicht zwischen sin & cos liegt der Geschwindigkeitsfaktor 26, sondern entscheidend ist der Zielspeicher (denn ich werte ja die Ergebnisse auf schlechteste Genauigkeit über kompletten Bereich aus)! Ich hatte testweise sin(x) durch cos(Pi/2-x) ersetzt und war nicht schneller, sondern langsamer! Die selbe cos-Funktion benötigt im selben Programm je nach Speicher1 oder Speicher2 mit der selben for-Schleife mal 0.2 s und mal 6.7 s!! \sourceon c++ Ausrichtung=32;//64 ändert nichts sx2 = (float *) _mm_malloc(uSize_float_double * n, Ausrichtung);//langsam cx1 = (float *) _mm_malloc(uSize_float_double * n, Ausrichtung);//schnell \sourceoff Das klingt nach Adressenausrichtung (https://docs.microsoft.com/de-de/cpp/cpp/align-cpp?view=msvc-160). Zig Compilerschalter & Dummy-Arrays probiert -> noch keine eindeutige Logik gefunden... Noch krasser bei atan2f: mal 10.91 ns mal 0.25 ns (Faktor 43)! Solange das nicht eindeutig ist, ergeben kleine Optimierungen keine logischen Messergebnisse. Es könnte auch an der Zeitmessung oder an der höchsten Geschwindigkeitsoptimierung liegen. Wenn Schleifen mit selben Ergebnissen mehrfach laufen, könnte etwas weggelassen werden...

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.29, eingetragen 2021-10-20

Also an der Zeitmessung liegt es nicht: eine komplett andere Zeitfunktion ist auf die ms exakt gleich. Auch beim Zusehen erkennt man sofort den Faktor 26. Am "Erstzugriff eines frischen Speichers" liegt es auch nicht, denn: - eine extra-Schleife beschreibt alle Arrays bereits zuvor - ohne diese Extra-Schleife würde sich die Zeit nochmals verdoppeln (statt 6,7 s über 11 s egal ob sin oder cos) Auch an der Aufteilung der beiden for-Schleifen liegt es nicht: \sourceon c++ n=1048576, iterations=1024; for (j = 0; j < iterations; j++) for (i = 0; i < n; i++) sx2[i] = cosf(x[i]); \sourceoff Es hätte ja sein können, dass der Compiler erkennt, dass die äußere for Schleife zum Ergebnis nichts beiträgt -> also: n=1048576*2, iterations=1024/2 -> jedoch selbe Zeiten!? Vertauschen der Abarbeitungsreihenfolge: bringt auch nichts! cx1 bleibt das schnellere Array von beiden !?!?!?! Nächster Verdacht: die unterschiedlichen Cache-Typen. Kleine Programme passen in den schnelleren "inneren CPU-Cache" und je größer ein Programm wird, um so langsamer wird es auch. ABER cx1 ist weiter vom anderen Array x[i] entfernt als sx2 jedoch schneller!?!

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.30, eingetragen 2021-10-20

\quoteon(2021-10-20 00:27 - hyperG in Beitrag No. 28) ... Wenn 1 einzige (Unter-)Funktion etwa 111 ns benötigt, dann sind die data1, data5 & data6 aber sehr "geschönte Argumente", wenn der komplette Algorithmus mit zig weiteren Funktionen und 3 Ergebnissen teils unter 60 ns liegen... \quoteoff Hallo hyperG, die kubische Wurzel einer komplexen Zahl wird ja nur für den Fall von drei verschiedenen rellen Nullstellen berechnet. Also nicht bei data1, data5 & data6. data1 bis data7 dienen nur zur Illustration und sind gar nicht für Geschwindigkeitsvergleiche gedacht. data6 und data7 sollen wohl eher die Grenzen des Programms aufzeigen. Siehe die Kommentare zu data6 und data7 in Beitrag No. 1. \quoteon \sourceon c++ ... statt cpow habe ich (x+y*i)^(1/3)=(y^2 + x^2)^(1/6) (Cos[ArcTan[y/x]/3] + I Sin[ArcTan[y/x]/3]) benutzt: ... \sourceoff \quoteoff Das habe ich auch einmal in "meinem" C Programm mit double-Genauigkeit getestet. mp_cbrt_cmplx_2.c \sourceon C if (cabs(s) < deps) { at3 = atan2(cimag(CMPLX(2.0, 0.0)*q), creal(CMPLX(2.0, 0.0)*q)) / 3.0; sqr6 = pow(creal(CMPLX(2.0, 0.0)*q)*creal(CMPLX(2.0, 0.0)*q) + cimag(CMPLX(2.0, 0.0)*q)*cimag(CMPLX(2.0, 0.0)*q), 1.0/6.0); q = CMPLX(sqr6*cos(at3), sqr6*sin(at3)); } else { at3 = atan2(cimag(s), creal(s)) / 3.0; sqr6 = pow(creal(s)*creal(s) + cimag(s)*cimag(s), 1.0/6.0); q = CMPLX(sqr6*cos(at3), sqr6*sin(at3)); } \sourceoff Das bringt hier aber nur ca. 15% Zeitersparnis. \sourceon mp_cbrt_cmplx.c mp_cbrt_cmplx_2.c data2 222 ns 189 ns data3 212 ns 184 ns data4 184 ns 156 ns $ ./mp_cbrt_cmplx < data2 6.0000000000000009e-01 2.0816681711721685e-17 i 2.0816681711721685e-17 3.1225022567582559e-18 i 2.9999999999999982e-01 -5.5511151231257827e-17 i 1.3877787807814457e-17 3.3306690738754672e-18 i 1.0000000000000012e-01 3.1225022567582602e-17 i 1.7347234759768071e-17 3.1225022567582555e-18 i CPU time: 222098 ms $ ./mp_cbrt_cmplx_2 < data2 5.9999999999999998e-01 -2.0816681711721685e-17 i 1.7347234759768071e-17 -3.1225022567582517e-18 i 2.9999999999999993e-01 5.5511151231257827e-17 i 1.3877787807814457e-17 -3.3306690738754672e-18 i 1.0000000000000001e-01 -3.1225022567582546e-17 i 7.8000162747683151e-34 -3.1225022567582544e-18 i CPU time: 189265 ms $ ./mp_cbrt_cmplx < data3 2.0000000000000000e+00 -1.1102230246251565e-16 i -2.4651903288156619e-32 -3.3306690738754696e-16 i 9.9999999999999989e-01 2.2204460492503131e-16 i 2.2204460492503131e-16 -4.4408920985006262e-16 i -9.9999999999999989e-01 -5.5511151231257852e-17 i 6.6613381477509392e-16 -3.3306690738754706e-16 i CPU time: 211852 ms $ ./mp_cbrt_cmplx_2 < data3 2.0000000000000000e+00 -1.1102230246251565e-16 i -2.4651903288156619e-32 -3.3306690738754696e-16 i 9.9999999999999989e-01 2.2204460492503131e-16 i 2.2204460492503131e-16 -4.4408920985006262e-16 i -9.9999999999999989e-01 -5.5511151231257852e-17 i 6.6613381477509392e-16 -3.3306690738754706e-16 i CPU time: 184113 ms $ ./mp_cbrt_cmplx < data4 9.9999999999999964e+00 3.3087224502121107e-24 i -2.2737367544323206e-13 2.1175823681357471e-22 i 2.0000000596046452e+00 -2.2204460492503131e-16 i -2.8421709430404007e-14 2.1175823602471418e-22 i 1.9999999403953561e+00 2.2204460327066996e-16 i -2.8421709430404007e-14 2.1175823365813146e-22 i CPU time: 183983 ms $ ./mp_cbrt_cmplx_2 < data4 9.9999999999999964e+00 3.3087224502121107e-24 i -2.2737367544323206e-13 2.1175823681357471e-22 i 2.0000000596046452e+00 -2.2204460492503131e-16 i -2.8421709430404007e-14 2.1175823602471418e-22 i 1.9999999403953561e+00 2.2204460327066996e-16 i -2.8421709430404007e-14 2.1175823365813146e-22 i CPU time: 155550 ms \sourceoff \quoteon Aber egal wie "explizit" der Algorithmus auch ist: durch die relativ langsamen atan, sin, cos, pow Maschinenbefehle (oder selbst Optimierungen), wird man Halley’s Näherungsfunktion schlecht überholen können. \quoteoff Ja, sieht so aus. Da die Routinen von Deiters et al. nur reell rechnen, ist der Vergleich sowieso fraglich. Wie bereits von MontyPythagoras geschrieben. [Die Antwort wurde nach Beitrag No.28 begonnen.]

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.31, eingetragen 2021-10-20

zu \sourceon c++ double x; for (i = 0; i < nrep; i++) x += pow(x, 1.0/3.0); \sourceoff Da wird ja auch nur "eine Speicherzelle" benutzt. (zwar rekursiv, aber ohne Ausrichtungsprobleme) Wenn ich \sourceon c++ for (i = 0; i < n; i++) sx2[i] = cosf(x[i]); ersetze durch float fh; for (i = 0; i < n; i++) fh = cosf(x[i]); \sourceoff bekomme ich 0 ns! Der optimierte Compiler merkt, wenn Berechnungen zwar angeordnet, aber für eine Ausgabe sinnlos sind. Was ich aber absolut unlogisch finde, ist \sourceon c++ for (i = 0; i < n; i++) sx2[i] = cosf(x[i]); printf(" %.3f s = %.3f s sx2[13]=%f\n",Zeit1,Zeit2,sx2[13]... 6.539 s = 6.539 s sx2[13]=0.939554 ---------------------------------------- for(i=0; i nichts. :-(

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.32, eingetragen 2021-10-20

Ah, ich hab's: Wie in Beitrag No.29 vermutet gibt's noch eine darüberliegende "Zeitschindungsschleife" mit iterations. Der Compiler erkennt diese "nichts zur Berechnung beitragenden Befehle" bei den Arrays unterschiedlich! Erst mit iterations=1 gleichen sich die Zeiten an: \sourceon nameDerSprache sin : 6.2 ns cos : 3.3 ...6.2 ns atanf2:11.5 ns powf : 9.7 ns cpow :29.0 ns \sourceoff zu "nur ca. 15% Zeitersparnis" -> na das ist doch was.

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.33, eingetragen 2021-10-21

Weitere Tests mit AVX2 & FMA Befehlen, wenn man diese "im Block" anwenden kann (gerade bei Wetterdaten oft): \sourceon c++ AVX2 1 Thread Funktion Zeit größte Ungenauigkeit im Wertebereich -18...22 rad sin Tab 0.64 ns 7.749e-07 sin Horn8 0.67 ns 1.511e-06 sin Horn6 0.52 ns 1.511e-06 sin256_ps 1.10 ns 5.960e-08 cos Horn6 0.61 ns 1.628e-06 atan2f FMA 1.36 ns 2.034e-04 \sourceoff Wenn man die Genauigkeit von nur 4 Nachkommastellen toleriert und im Block rechnen kann, können CPUs mit AVX2 & FMA Befehlssatz bei atan2f bis zu 8.5 bei sin bis zu 11.9 mal schneller sein. Noch krasser geht es mit neuem Intel-Compiler & AVX512: Da gibt es den optimierten Befehl __m512 _mm512_sincos_ps (__m512 * mem_addr, __m512 a) \sourceon Basic FOR j := 0 to 15 i := j*32 dst[i+31:i] := SIN(a[i+31:i]) MEM[mem_addr+i+31:mem_addr+i] := COS(a[i+31:i]) Next \sourceoff Also mit 1 Befehl 16 float Inputs je sin & cos zeitgleich berechnen!! Ich habe zwar AVX512, aber nicht diesen neuen C-Compiler mit SVML... In Intel-Compiler Foren findet man aber viele Fragen, da es meist nicht funktioniert... Grüße

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.34, eingetragen 2021-10-23

Das mit _mm512_sincos_ps hat mich nicht losgelassen. Zwar habe ich keinen Intel-C-Compiler, aber es gibt einige Überlegungen im Internet. \sourceon c++ __m512 xin, ysin, ycos, yAtan, yAtan2,y2; xin = _mm512_set_ps(-0.1, 0.042857142857142844, 0.1857142857142857, 0.3285714285714285, 0.4714285714285714, 0.6142857142857143, 0.7571428571428571, 0.9, 1.0428571428571427, 1.1857142857142857, 1.3285714285714285, 1.4714285714285713, 1.614285714285714, 1.7571428571428571, 1.9, 2.0428571428571427); y2 = _mm512_set1_ps(0.4f); sincos512_ps(xin, &ysin, &ycos); yAtan = atan512f_ps(xin); yAtan2=atan2512f_ps(y2, xin);//jede Funktion berechnet 16 float parallel ohne Multitasking for (i = 15; i >=0; i--) printf("x=%.16f sin=%.16f cos=%.16f atan=%.16f, atan2=%.16f\n", xin.m512_f32[i], ysin.m512_f32[i], ycos.m512_f32[i], yAtan.m512_f32[i], yAtan2.m512_f32[i]); ergibt für den Typ float genügend genaue Werte x=-0.1000000014901161 sin=-0.0998334214091301 cos=0.9950041770935059 atan=-0.0996686592698097, atan2=1.8157750368118286 x=0.0428571440279484 sin=0.0428440272808075 cos=0.9990817904472351 atan=0.0428309328854084, atan2=1.4640606641769409 x=0.1857142895460129 sin=0.1846485882997513 cos=0.9828045964241028 atan=0.1836223304271698, atan2=1.1361260414123535 x=0.3285714387893677 sin=0.3226912319660187 cos=0.9465042948722839 atan=0.3174587488174438, atan2=0.8831250667572021 x=0.4714285731315613 sin=0.4541594982147217 cos=0.8909204006195068 atan=0.4405303895473480, atan2=0.7036138176918030 x=0.6142857074737549 sin=0.5763749480247498 cos=0.8171853423118591 atan=0.5508575439453125, atan2=0.5771922469139099 x=0.7571428418159485 sin=0.6868476867675781 cos=0.7268013954162598 atan=0.6480568647384644, atan2=0.4860319495201111 x=0.8999999761581421 sin=0.7833269238471985 cos=0.6216099858283997 atan=0.7328150868415833, atan2=0.4182243943214417 x=1.0428571701049805 sin=0.8638470768928528 cos=0.5037541389465332 atan=0.8063741326332092, atan2=0.3662555515766144 x=1.1857142448425293 sin=0.9267675876617432 cos=0.3756352365016937 atan=0.8701618909835815, atan2=0.3253606557846069 x=1.3285714387893677 sin=0.9708067178726196 cos=0.2398631572723389 atan=0.9255770444869995, atan2=0.2924430072307587 x=1.4714285135269165 sin=0.9950670599937439 cos=0.0992043688893318 atan=0.9738852977752686, atan2=0.2654303908348083 x=1.6142857074737549 sin=0.9990544915199280 cos=-0.0434756726026535 atan=1.0161842107772827, atan2=0.2428953200578690 x=1.7571429014205933 sin=0.9826876521110535 cos=-0.1852699667215347 atan=1.0534031391143799, atan2=0.2238279581069946 x=1.8999999761581421 sin=0.9463000893592834 cos=-0.3232895433902740 atan=1.0863183736801147, atan2=0.2074962258338928 x=2.0428571701049805 sin=0.8906331062316895 cos=-0.4547227025032043 atan=1.1155755519866943, atan2=0.1933578997850418 \sourceoff und analog der letzten Tabelle folgende effektive Geschwindigkeit: \sourceon c++ AVX512 (z.B. i9) Funktion Zeit größte Ungenauigkeit im Wertebereich -18...22 rad sincosAVX512 0.99 ns 5.960e-08 atan2fAVX512 1.10 ns noch verbesserungswürdig \sourceoff Das ist zwar nicht die Geschwindigkeit & Genauigkeit der echten _mm512_sincos_ps, aber immerhin schneller & genauer als je sin & cos optimierte aus Beitrag zuvor... cpow könnte von 29 ns auf (10+1+1=12) ns sinken. Fehlt noch die Optimierung von pow, die mit 9.7 ns recht langsam war... Da der Exponent reell ist, gibt es 3 grobe Möglichkeiten: a) x^(1/6) = exp(log(x)/6) b) Näherung + Iteration c) Intel-c-Compiler mit __m512 _mm512_pow_ps (__m512 a, __m512 b) ...

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.35, eingetragen 2021-10-23

\quoteon(2021-10-14 08:10 - MontyPythagoras in Beitrag No. 9) ... 3. Wie haben sie die cardanischen Formeln umgesetzt? Etwa wie in den Wikipediaartikeln zur kubischen Gleichung (hier) oder den Cardanischen Formeln (hier) beschrieben, mit den ganzen Fallunterscheidungen und $\sin$ und $\cos$ und $\cosh$ und so weiter? Wo man $\Delta=\left(\frac{q}2\right)^2 + \left(\frac{p}3\right)^3 = \frac{27 A^2 D^2+4 B^3D-18 A B C D+4 A C^3-B^2 C^2}{108 A^4}$ berechnet, und die drei Lösungen im komplexen Fall durch Multiplikation mit den dritten Einheitswurzeln berechnet? Geht es denn noch uneffektiver? Dann muss ich mich nicht wundern, wenn die cardanische Lösung langsam ist. ... Daher habe ich meinen Code hier mal reingeschrieben, wo ich eine wirklich effiziente Umsetzung der cardanischen Formel angestrebt habe. Keine Fallunterscheidungen, keine riesigen Formeln wie die oben gezeigte, das wäre alles dramatisch uneffektiv. 9 Zeilen Code mit einer Handvoll algebraischen Anweisungen - und das war's. Ich habe daher arge Zweifel, dass die cardanische Lösung von den Machern des Algorithmus optimiert wurde, und daher würde ich gerne meinen Code umgesetzt sehen. Ich übersetze es auch gerne selbst in C++. Ciao, Thomas \quoteoff Hallo MontyPythagoras, von John Burkardt gibt es eine C-Bibliothek c8lib für komplexe Zahlen mit Routinen zum Lösen von quadratischen, kubischen und quartischen Gleichungen. Er rechnet, soweit es geht, reell und berechnet "die drei Lösungen im komplexen Fall durch Multiplikation mit den dritten Einheitswurzeln". Siehe die Funktion r8poly3_root in c8lib.c. Das scheint aber gar nicht so ineffizient zu sein. Die trig. Funktionen $\arccos$ und $\cos$ sind auf modernen Prozessoren effizient implementiert. \sourceon Zeiten für 1 Iteration: mp_cbrt_cmplx.c mp_cbrt_cmplx_2.c r8poly3_root Halley/Deiters data2 222 ns 189 ns 95 ns 64 ns data3 212 ns 184 ns 104 ns 62 ns data4 184 ns 156 ns 81 ns 29 ns ========================================================================== Zeiten für r8poly3_root für 1000000000 Iterationen: $ ./bh_cbrt_2 < data2 Polynomial coefficients A, B, C, D: A = 1 B = -1 C = 0.27 D = -0.018 Roots: (0.1,0) (0.6,0) (0.3,0) CPU time: 94950 ms $ ./bh_cbrt_2 < data3 Polynomial coefficients A, B, C, D: A = 1 B = -2 C = -1 D = 2 Roots: (-1,0) (2,0) (1,0) CPU time: 104203 ms $ ./bh_cbrt_2 < data4 Polynomial coefficients A, B, C, D: A = 1 B = -14 C = 44 D = -40 Roots: (10,0) (2,0) (2,-0) CPU time: 81137 ms \sourceoff

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.36, eingetragen 2021-10-24

Hallo, auf dieser Internetseite Numerical Analysis & Applications gibt es u.a. Infos und Programme zu Polynominal Root Finders. "Root Finder, finds all zeros (roots) of a polynomial of any degree with either real or complex coefficients using Bairstow's, Newton's, Halley's, Graeffe's, Laguerre's, Jenkins-Traub, Abert-Ehrlich, Durand-Kerner, Ostrowski or Eigenvalue method." Es gibt auch einen interaktiven Rechner zur Bestimmung von Nullstellen Polynomial Root finder. Für einige der Algorithmen gibt es C++-Programme Various Scientific Software Application.

   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 453
  Beitrag No.37, eingetragen 2021-10-24

Hallo, das C++-Programm von Henrik Vestermark Jenkins-Traub C++ für komplexe Polynomkoeffizienten ist für einfache Beispiele mit reellen Koeffizienten ca. eine Größenordnung langsamer als Cardano oder Halley. \sourceon C /* ******************************************************************************* * * * Copyright (c) 2002 * Henrik Vestermark * Denmark * * All Rights Reserved * * This source file is subject to the terms and conditions of the * Henrik Vestermark Software License Agreement which restricts the manner * in which it may be used. * ******************************************************************************* * * * Module name : cpoly.cpp * Module ID Nbr : * Description : cpoly.cpp -- Jenkins-Traub real polynomial root finder. * Translation of TOMS493 from FORTRAN to C. This * implementation of Jenkins-Traub partially adapts * the original code to a C environment by restruction * many of the 'goto' controls to better fit a block * structured form. It also eliminates the global memory * allocation in favor of local, dynamic memory management. * * The calling conventions are slightly modified to return * the number of roots found as the function value. * * INPUT: * opr - double precision vector of real coefficients in order of * decreasing powers. * opi - double precision vector of imaginary coefficients in order of * decreasing powers. * degree - integer degree of polynomial * * OUTPUT: * zeror,zeroi - output double precision vectors of the * real and imaginary parts of the zeros. * to be consistent with rpoly.cpp the zeros is inthe index * [0..max_degree-1] * * RETURN: * returnval: -1 if leading coefficient is zero, otherwise * number of roots found. * -------------------------------------------------------------------------- \sourceoff \sourceon Zeiten für 1 Iteration: cpoly.cpp data1 2203 ns data2 1880 ns data3 2034 ns data4 1661 ns data5 1304 ns data6 1968 ns data7 1220 ns ============================================================== Zeiten für 1000000 Iterationen: $ ./mpoly < data1 Number of roots: 3 -1.0000000000000000 0.0000000000000000 i 0.5000000000000000 0.8660254037844386 i 0.5000000000000000 -0.8660254037844386 i CPU time: 2203 ms $ ./mpoly < data2 Number of roots: 3 0.1000000000000000 0.0000000000000000 i 0.3000000000000002 0.0000000000000000 i 0.5999999999999999 -0.0000000000000000 i CPU time: 1880 ms $ ./mpoly < data3 Number of roots: 3 -0.9999999999999999 0.0000000000000000 i 0.9999999999999997 0.0000000000000000 i 2.0000000000000004 -0.0000000000000000 i CPU time: 2034 ms $ ./mpoly < data4 Number of roots: 3 2.0000000000438583 -0.0000000000198884 i 1.9999999999561346 0.0000000000198790 i 10.0000000000000071 0.0000000000000095 i CPU time: 1661 ms $ ./mpoly < data5 Number of roots: 3 -1.0000000000000018 -0.0000000000000007 i -0.9999999999999989 0.0000000000000002 i -0.9999999999999993 0.0000000000000005 i CPU time: 1304 ms $ ./mpoly < data6 Number of roots: 3 1.0000009999950223 0.0000000000002368 i 1.0000020066058000 -0.0000000072388417 i 0.9999999933991777 0.0000000072386049 i CPU time: 1968 ms $ ./mpoly < data7 Number of roots: 3 1.0000000000000002 0.0000000000000000 i 10000000.0000000000000000 0.0000000009313226 i 100000000000000.0000000000000000 -0.0000000009313226 i CPU time: 1220 ms \sourceoff

   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 2115
  Beitrag No.38, eingetragen 2021-10-25

\quoteon(2021-10-14 19:43 - AlphaSigma in Beitrag No. 14) ... \sourceon method 0 method 1 Cardano Halley (Deiters) data1 59 ns 38 ns data2 120 ns 64 ns data3 121 ns 62 ns data4 133 ns 29 ns data5 31 ns 12 ns data6 87 ns 13 ns data7 133 ns 42 ns \sourceoff Verwendete CPU: \sourceon ... model name : Intel(R) Core(TM) i3-3240 CPU @ 3.40GHz ... \sourceoff ... \quoteoff Mich interessierte mal der Einfluss von fsincos auf Cardano (ohne swap): \sourceon data4 -40 +44 -14 1 Cardano VC2017 | Halley (Deiters) i9 i7 VC2012 fsincos sin(x),cos(x) | ohne AVX AVX2 f:precise f:fast 119.6 ns 85.6 ns | 23.6 ns 20.8 ns 44.5 ns 34.4 ns data3 | 37.5 ns \sourceoff Jetzt ist mir auch klar, warum VC die alten Coprozessorbefehle nicht mehr unterstützt (fsincos musste ich mühsam per ASM -> OBJ einbinden): viel langsamer als einzelne Befehle Sin & cos zusammen! (Coprozessoren waren damals ja auf 80-Bit Genauigkeit ausgelegt, was natürlich mit "einfachen Näherungsfunktionen" nicht erreicht werden konnte; leider unterstützt VC den "long double" nicht mehr für volle 19 Stellen {nur mit ASM-Trick}) Die schnellste Zeit mit AVX2 bezieht sich nur auf den Compilerschalter. Da wurde noch nichts mit 4 oder 8 "pro Block" optimiert, wie es sonst bei AVX der optimierte Fall ist. Solange cpow(z,1/3) etwa so langsam wie der komplette Algorithmus von Halley/Deiters ist und man keine komplexen Ergebnisse benötigt, ist die Variante ohne trigonometrische Funktionen natürlich schneller. Aber nicht 10 mal (die vielen uneffektiven Teile hatten wir ja schon besprochen), sondern höchstens 4 mal. Da man eh durch Fehlerfortpflanzung am Ende eine Newton Iteration zur Genauigkeitskorrektur benötigt, sollte bei cpow eine float Genauigkeit und Näherungs-Algorithmus reichen. So könnte man bis auf Geschwindigkeits-Faktor 3 oder 2 herankommen...

   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
  Beitrag No.39, vom Themenstarter, eingetragen 2021-10-25

Hallo zusammen, Ihr wart ganz schön fleißig in meiner Abwesenheit, da musste ich erst einmal einen Haufen lesen - und verstehen. Was Ihr ausgearbeitet und recherchiert habt, gräbt ja schon sehr tief. Vielen Dank dafür. Ich fasse den aktuellen Stand mal grob zusammen: * Die Aussage im Wikipedia-Artikel gilt vor allem für reelle Koeffizienten und Lösungen. Für komplexe Koeffizienten ist die Aussage zumindest fragwürdig, wenn ich die sehr umfangreichen Informationen richtig gefiltert habe. * In der im Artikel zitierten Cardano-Implementierung wurde noch eine Newton-Schleife hinten angehängt, um die Präzision zu erhöhen und Rundungsfehler zu reduzieren, was die Bewertung verfälscht bzw. die "reine" Cardano-Lösung verlängert. * Im Halley-Deiters-Algorithmus wird (meiner Meinung nach) die Iteration etwas zu früh abgebrochen, wenn die Lösung noch gar nicht sicher gefunden ist. * Es sind insbesondere die trigonometrischen Funktionen und die Potenzfunktion, die in der cardanischen Lösung die Zeit fressen. Beim Potenzieren einer komplexen Zahl mit (1/3) muss der Prozessor letztlich auch trigonometrische Funktionen anwenden. * Es gibt bessere Implementierungen der cardanischen Lösung als die im Artikel zitierte. Etwas undurchsichtig ist für mich, warum eigentlich die trigonometrischen Funktionen und die Potenzfunktion SOO langsam sind. Letzten Endes muss ein Prozessor zum Beispiel zum Wurzelziehen auch einen Algorithmus anwenden, der der Suche nach einer Wurzel eines kubischen Polynoms gleicht. Siehe zum Beispiel das Heron-Verfahren. Wenn das Halleyverfahren, umgesetzt in einer Hochsprache, aus einem allgemeinen (reellen) kubischen Polynom extrem schnell eine Wurzel findet, warum dann nicht ein Prozessor die Wurzel aus dem simplen Polynom $x^3-a=0$? Er könnte intern ja auch das Halleyverfahren anwenden. Sinn macht das für mich irgendwie nicht. Ciao, Thomas

   Profil
-->> Fortsetzung auf der nächsten Seite -->>
Seite 1Gehe zur Seite: 1 | 2 | 3 | 4  

Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2023 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]