Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Kubische Gleichung: numerische Näherung schneller als Cardanische Formeln?
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: 2827
Wohnort: Werne
  Themenstart: 2021-10-13 17:46

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: 298
  Beitrag No.1, eingetragen 2021-10-13 20:40

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: 2827
Wohnort: Werne
  Beitrag No.2, vom Themenstarter, eingetragen 2021-10-13 21:03

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: 2827
Wohnort: Werne
  Beitrag No.3, vom Themenstarter, eingetragen 2021-10-13 21:23

...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: 298
  Beitrag No.4, eingetragen 2021-10-13 22:48

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: 2056
Wohnort: Thüringen
  Beitrag No.5, eingetragen 2021-10-13 23:24

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: 298
  Beitrag No.6, eingetragen 2021-10-13 23:37

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: 211
  Beitrag No.7, eingetragen 2021-10-14 06:49

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: 2056
Wohnort: Thüringen
  Beitrag No.8, eingetragen 2021-10-14 07:26

@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: 2827
Wohnort: Werne
  Beitrag No.9, vom Themenstarter, eingetragen 2021-10-14 08:10

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: 2056
Wohnort: Thüringen
  Beitrag No.10, eingetragen 2021-10-14 08:27

"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: 2827
Wohnort: Werne
  Beitrag No.11, vom Themenstarter, eingetragen 2021-10-14 10:04

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: 2827
Wohnort: Werne
  Beitrag No.12, vom Themenstarter, eingetragen 2021-10-14 10:10

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: 211
  Beitrag No.13, eingetragen 2021-10-14 18:17

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: 298
  Beitrag No.14, eingetragen 2021-10-14 19:43

\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: 2827
Wohnort: Werne
  Beitrag No.15, vom Themenstarter, eingetragen 2021-10-14 20:39

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: 298
  Beitrag No.16, eingetragen 2021-10-14 22:04

\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: 298
  Beitrag No.17, eingetragen 2021-10-15 18:07

\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: 2827
Wohnort: Werne
  Beitrag No.18, vom Themenstarter, eingetragen 2021-10-15 21:30

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: 298
  Beitrag No.19, eingetragen 2021-10-15 21:56

\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: 2827
Wohnort: Werne
  Beitrag No.20, vom Themenstarter, eingetragen 2021-10-15 22:10

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: 298
  Beitrag No.21, eingetragen 2021-10-16 22:51

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: 298
  Beitrag No.22, eingetragen 2021-10-16 23:04

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: 2827
Wohnort: Werne
  Beitrag No.23, vom Themenstarter, eingetragen 2021-10-16 23:46

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: 1512
  Beitrag No.24, eingetragen 2021-10-17 00:52

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: 2827
Wohnort: Werne
  Beitrag No.25, vom Themenstarter, eingetragen 2021-10-17 07:39

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
Folgende Antworten hat der Fragesteller vermutlich noch nicht gesehen.
Er/sie war noch nicht wieder auf dem Matheplaneten
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1512
  Beitrag No.26, eingetragen 2021-10-17 09:03

\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: 298
  Beitrag No.27, eingetragen 2021-10-17 15:02

\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: 1512
  Beitrag No.28, eingetragen 2021-10-20 00:27

\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...


   Profil

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-2021 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]