Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Kubische Gleichung: numerische Näherung schneller als Cardanische Formeln?
Thema eröffnet 2021-10-13 17:46 von MontyPythagoras
Seite 3   [1 2 3]   3 Seiten
Autor
Kein bestimmter Bereich Kubische Gleichung: numerische Näherung schneller als Cardanische Formeln?
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.80, vom Themenstarter, eingetragen 2021-11-08 14:24

Hallo zusammen, nachdem ich ein wenig in der C++-Dokumentation gekramt habe, habe ich festgestellt, dass mit C++11 die dritte Wurzel als eigene Funktion eingeführt wurde: siehe hier. Warum quälen wir uns hier eigentlich mit der pow-Funktion rum? Laut verlinkter Seite ist die Genauigkeit sogar höher als mit der pow-Funktion, und das Vorzeichen berechnet sie auch richtig! Fantastisch. 😁 Und ich glaube kaum, dass man die dritte Wurzel einführt als eigene Funktion und dann so schlecht umsetzt, dass die pow-Funktion schneller ist. Da können wir also weiter optimieren: \sourceon C++ \numberson void cbrt_real(const double coeff[], double res[]) { const double rcpr3 = -0.3333333333333333333; const double sqrt3 = 1.7320508075688772935; double s = rcpr3 / coeff[0]; double r = coeff[1] * s; double q = r * r; double p = coeff[2] * s + q; q = q * r + s * (coeff[3] + coeff[2] * r) * 1.5; s = q * q - p * p * p; if (s > 0.0) { (q < 0.0) ? s = cbrt(q - sqrt(s)) : s = cbrt(q + sqrt(s)); res[0] = r + s + p / s; } else { double t = sqrt(p); sincos(rcpr3 * atan2(q, sqrt(-s)), &s, &p); q = t * (sqrt3 * p - s); s *= 2 * t; res[0] = r + s; res[1] = r + q; res[2] = r - s - q; } } \sourceoff EDIT: Ich habe noch ein paar Vorzeichen umgedreht, was im Laufe des Codes ein paar erzwungene Vorzeichenwechsel erspart. Ich hoffe, dabei ist nichts schief gegangen. 😉 Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.81, vom Themenstarter, eingetragen 2021-11-08 16:34

Einen habe ich noch: in Zeile 18 multipliziere ich "t" mit 2. Das wird vermutlich als Fließkomma-Multiplikation kompiliert. Man kann das exakt Doppelte eines Fließkommawertes aber auch durch direkte Bitmanipulation der Fließkommazahl umsetzen. Dafür stellen die C++-Entwickler in ihrer grenzenlosen Weisheit die Funktion "ldexp" zur Verfügung. In Zeile 18 könnte man daher alternativ verwenden: \sourceon C++ \numberson s *= ldexp(t,1); \sourceoff Allerdings scheint das nicht immer schneller zu sein als die direkte Multiplikation mit 2, je nach Compiler oder Prozessor. Kommt auf einen Versuch an. Damit bin ich aber auch erstmal am Ende mit Optimiererei. @HyperG: Jetzt hast Du erst einmal ein bisschen zu tun... 😁 Ciao, Thomas


   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 322
  Beitrag No.82, eingetragen 2021-11-08 17:11

\quoteon(2021-11-08 14:24 - MontyPythagoras in Beitrag No. 80) Hallo zusammen, nachdem ich ein wenig in der C++-Dokumentation gekramt habe, habe ich festgestellt, dass mit C++11 die dritte Wurzel als eigene Funktion eingeführt wurde: siehe hier. Warum quälen wir uns hier eigentlich mit der pow-Funktion rum? Laut verlinkter Seite ist die Genauigkeit sogar höher als mit der pow-Funktion, und das Vorzeichen berechnet sie auch richtig! Fantastisch. 😁 Und ich glaube kaum, dass man die dritte Wurzel einführt als eigene Funktion und dann so schlecht umsetzt, dass die pow-Funktion schneller ist. ... Ciao, Thomas \quoteoff Hallo MontyPythagoras, darauf hatte ich bereits in Beitrag No. 48 hingewiesen. Ein Pendant für complex scheint es aber nicht zu geben. \quoteon(2021-10-27 21:54 - AlphaSigma in Beitrag No. 48) \quoteon(2021-10-26 20:49 - MontyPythagoras in Beitrag No. 43) ... Das eigentliche Problem an der Stelle scheint mir zu sein, dass der Prozessor keine Routine hat, um eine dritte Wurzel zu ziehen. Es gibt zwar den Prozessorbefehl FSQRT (also zweite Wurzel ziehen), aber keinen für die dritte Wurzel. Man muss sich daher bei der Programmierung damit behelfen, dass man mit $\tfrac13$ potenziert, was vermutlich über den Umweg $x^{\frac13}=2^{{\frac13}\cdot\log_2x}$ bewerkstelligt wird und daher eine Ewigkeit länger dauert. Das ist meines Erachtens der Hauptgrund, warum der Halleyalgorithmus schneller ist. Vielleicht kann man letzteren überholen, wenn man eine "Dritte-Wurzel"-Funktion in Assembler programmiert. Der FSQRT-Befehl beruht meines Erachtens auf dem Heron-Verfahren, welches man auch auf beliebige Wurzeln verallgemeinern könnte... Ciao, Thomas \quoteoff Hallo, seit C99 (glaube ich) gibt es für double eine cbrt Funktion, aber nicht für complex. Ich habe versucht die Zeiten für pow(x, 1.0/3) und cbrt auf meinem PC zu bestimmen. Ich schreibe extra "versucht", da das manchmal gar nicht so einfach ist. Im Internet habe ich noch diese apple_tr_kt32_cuberoot.pdf historisch interssante, aber heute wohl nicht mehr relevante, Routine gefunden und in den Vergleich mit aufgenommen. \sourceon $ gcc -Wall -pedantic -lm -o rand rand.c $ ./rand Method Total CPU time Single CPU time float cbrt(x) 3889 ms 38 ns 268435456.000000 float pow(x, 1.0/3) 6535 ms 65 ns 268435456.000000 float ap_cbrt(x) 4327 ms 43 ns 268435456.000000 double cbrt(x) 2856 ms 28 ns 749732938.328323 double pow(x, 1.0/3) 6241 ms 62 ns 749732938.328323 $ gcc -Wall -pedantic -O2 -lm -o rand rand.c $ ./rand Method Total CPU time Single CPU time float cbrt(x) 2811 ms 28 ns 268435456.000000 float pow(x, 1.0/3) 6224 ms 62 ns 268435456.000000 float ap_cbrt(x) 2644 ms 26 ns 268435456.000000 double cbrt(x) 2704 ms 27 ns 749732938.328323 double pow(x, 1.0/3) 6039 ms 60 ns 749732938.328323 $ gcc -Wall -pedantic -Ofast -lm -o rand rand.c $ ./rand Method Total CPU time Single CPU time float cbrt(x) 2984 ms 29 ns 268435456.000000 float pow(x, 1.0/3) 5661 ms 56 ns 268435456.000000 float ap_cbrt(x) 2604 ms 26 ns 268435456.000000 double cbrt(x) 91 ms 0? ns 749732938.318400 double pow(x, 1.0/3) 91 ms 0? ns 749732938.318400 \sourceoff float ist nicht schneller als double. cbrt ist deutlich schneller als pow(x, 1.0/3). -O2 ist für float schneller als ohne Optimierung. Den Zeiten mit -Ofast für double traue ich noch nicht ganz. Ich habe die Argumente für die kubische Wurzel mit rand() erzeugt und in einem Array abgelegt. Damit der Compiler nicht alles wegoptimiert, habe ich alle Ergebnisse aufsumiert und am Ende ausgeben lassen. Die Zufallszahlen variieren von 0 bis 1000.0. Es sind insgesamt 100000*1000 = 100000000 Iterationen. Bei ap_cbrt() ist das Flag ITERATE auf 1 gesetzt. Voraussichtlich mache ich noch ähnliche Versuche für complex. ... \sourceoff \quoteoff


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.83, vom Themenstarter, eingetragen 2021-11-08 18:18

Hallo AlphaSigma, sorry, das ist zumindest bei mir untergegangen. Dabei hatte ich ja sogar die verlinkte PDF-Datei wahrgenommen und kommentiert. Wenn Du magst, könntest Du ja schon einmal den Code aus #80 (ggf. mit der Alternative aus #81) testen im Vergleich zu Halley von Deiters/Uni Köln als auch zu deren Cardano-Umsetzung. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.84, eingetragen 2021-11-08 20:02

\quoteon(2021-11-07 21:21 - MontyPythagoras in Beitrag No. 78) Hallo HyperG, ... die Abfrage q==0 sparen kann. Also teste bitte noch einmal diese Variante. ... Ciao, Thomas \quoteoff Dieses if ((p == 0.0) && (q == 0.0)) wollte ich schon damals nicht, weil man gerade bei "double" durch Rundungsfehler diesen Sonderfall fast nur mit Vielfachen von 2^x erreichen kann. Die Messung zeigt auch keinen messbaren Unterschied, da ein IF kaum Zeit benötigt: beide 36.0 ns zu "s *= ldexp(t,1);" + Jetzt hast Du erst einmal ein bisschen zu tun... s *= ldexp(t, 1); 56.00 ns -> 20 ns langsamer!!! Das ist also nur graue Theorie, die Praxis zeigt was anderes, da *2 bereits vom Compiler gut optimiert wird (z.B. t+t -> ADD ist super schnell)! Nachtrag: Außerdem bringt er den bereits auf AVX2 optimierten Code durcheinander. zu cbrt: Die hatte ich auch bereits in meinem Cardano, welches ich ja auf 40 ns angegeben hatte. Es spielt aber keine Rolle bei meinen Messungen, da es nur den Abkürzungs- Fall "1 reelle Lösung" betrifft, den ich absichtlich ausgeschlossen hatte. (wir wollen doch keine Benchmarks mit Fällen voller Abkürzungen) Noch was kurioses: Bei meinen Optimierungen mit AVX512 hatte ich gerade die Wandlung vom complex Array zum AVX512-Array als Test vor der eigentlichen Zeitmessung eingebaut... ...und schon stieg die Rechenzeit der alten (auf knapp über 200 ns optimierten komplexen Berechnung der 3 komplexen Ergebnisse) auf das DOPPELTE an!!!!! Zunächst dachte ich an Einflüsse wie: - größere Programme passen nicht mehr in den kleinen CPU-First-Level-Cache... - falsche Speicherausrichtung (Alignment) Aber als die AVX-Befehle an das Ende verschoben wurden -> war alles wieder schnell!!! Es scheint also hier dran zu liegen: Intel-Chip drosselt seine Frequenz , wenn er AVX-Code ausführt: https://www.zdnet.de/88381431/scharfe-kritik-an-intel-avx-512/ Das muss ich mal mit AVX256 testen... Grüße


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.85, vom Themenstarter, eingetragen 2021-11-08 20:17

Hallo HyperG, also zwischen #80 und #74 mit atan2 besteht fast kein Unterschied? Krass. Du kannst statt s *= ldexp(t,1) natürlich auch gerne s *= (t+t) versuchen. Dass der ldexp nicht zwangsläufig schneller ist als die Multiplikation mit 2 hatte ich schon gelesen. Dass er drastisch langsamer sein soll, ist dann doch überraschend, wenn man bedenkt, dass bis zum Endergebnis eigentlich doch das gleiche passiert. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.86, eingetragen 2021-11-08 20:56

\quoteon(2021-11-08 20:17 - MontyPythagoras in Beitrag No. 85) Hallo HyperG, ... drastisch langsamer sein soll, ist dann doch überraschend, wenn man bedenkt, dass bis zum Endergebnis eigentlich doch das gleiche passiert. ... \quoteoff Außerdem bringt er den bereits auf AVX2 optimierten Code durcheinander. Mit s *= t + t; -> teilweise noch 1 ns schneller: 35...36 ns


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.87, eingetragen 2021-11-08 23:11

Verrückt, welche AVX-Befehle den x86-Folgecode um welchen Faktor reduzieren: \sourceon nameDerSprache AVX-Befehl Faktor _mm256_unpacklo_pd 1 _mm256_permutex2var_pd 1.6 _mm512_permutex2var_pd 2 \sourceoff Das bedeutet, wenn man einmal mit AVX angefangen hat, sollte man bei AVX bleiben!? Oder: erst x86 Code und dann AVX-Code laufen lassen. Das macht eine Optimierung extrem kompliziert... 🤔


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.88, vom Themenstarter, eingetragen 2021-11-09 20:24

Hallo HyperG, mach doch jetzt bitte noch einmal einen Vergleich zwischen Halley und Cardano. Die ursprüngliche Frage war ja, welcher Algorithmus der schnellere ist. Die Optimierung der Prozessorbefehle hin, absolute Zeiten her - welcher ist schneller? Ich würde gerne ein paar realistische Koeffizientensätze verglichen sehen, einmal für kubische Gleichungen mit drei verschiedenen reellen Lösungen und einmal mit einfachen Lösungen. Auch für allgemeine, komplexe Koeffizienten würde mich ein finaler Vergleich interessieren. Mein letzter Stand, wenn ich nichts überlesen habe, ist, dass Cardano dabei erst einmal der schnellste Algorithmus ist (in Python auf jeden Fall sehr deutlich). Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.89, eingetragen 2021-11-10 22:07

Also hier noch eine extra EXE zum Vergleich der: a) langsamer gemachten "Deiters & Macias", weil beim Original von den 3 Ergebnissen 1 nicht mit der letzten Iteration berechnet wurde und zu wenig Genauigkeit hatte. b) Beitrag 80 in c++ (also ohne den extrem seltenen If-Fall; mit atan2) \sourceon Data9: 233328643357133 -14029858000 -628199 0.091 Algorithmus | i7-3770K 3.5GHz SSE2 | i9 4GHz SSE2 | AVX2 genauer Deiters & Macias | | | 15 Stellen | 72.50 ns | 48.75 ns | 40.00 ns -------------------------+----------------------+--------------+------- Cardano reell 14.5 Stell.| 53.75 ns | 37.50 ns | 35.00 ns \sourceoff Beide EXE habe ich unter https://www.lamprechts.de/exe/eqn_cubic2017.ZIP (gegen 22:29 Uhr korrigierte Version hochgeladen) hochgeladen. Wer keine neue AVX2 CPU hat, sollte nur die SSE2 exe testen! (eine Überprüfung der CPU-Kompatibilität habe ich aus Zeitgründen mal weggelassen; bei meinem i7 stürzt sie unter Win7 einfach ab, weil die neuen Maschinenbefehle unbekannt sind) \sourceon i7 SSE2 (Eingabe der Parameter mit Leerzeichen) A + B*x + C*x^2 + D*x^3=0; Anzahl (1...1000000000):800000 Wiederholungen=800000 Parameter (z.B. Data4: -40 44 -14 1 Data7:-1.0e21 1.00000010000001e21 -1.0000001 0000001e14 1.0e00 ) Parameter (z.B. Data8: -2 37 - 35 8; Data9: 233328643357133 -14029858000 -628199 0.091):233328643357133 -14029858000 -628199 0.091 Deiters & Macias :233328643357133.0000000000+-14029858000.0000000000*x+-628199.0000000000*x^2+0.0910000000*x^3=0 -33319.878552634829248 Err:-3.12500000000000e-02 11111.494876196791665 Err:3.12500000000000e-02 6925494.097962152212858 Err:-2.98128125000000e+03 CPU time: 58.0 ms also 72.50 ns cbrt_real_ohneIfP0:233328643357133.0000000000+-14029858000.0000000000*x+-628199.0000000000*x^2+0.0910000000*x^3=0 11111.494876191020012 Err:1.61375000000000e+02 6925494.097962151281536 Err:-2.99446875000000e+03 -33319.878552628681064 Err:1.73062500000000e+02 CPU time: 43.0 ms also 53.75 ns \sourceoff Die fehlenden msvcp140.dll, vcruntime140.dll, vcomp140.dll kann man sich alle bei Bedarf kostenlos aus dem Internet herunterladen. (exe meldet fehlende DLL) Grüße Gerd


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.90, vom Themenstarter, eingetragen 2021-11-11 21:02

Hallo Gerd, ich habe mir Dein Programm heruntergeladen und ein wenig damit herumprobiert. Ich denke, man kann sagen, dass im Falle von drei reellen Lösungen die Cardano-Lösung schneller ist, typischerweise um den Faktor 1.3, während im Falle von nur einer reellen Lösung Halley-Deiters manchmal schneller zu sein scheint - aber nicht immer. Halley-Deiters ist vor allem bei Trivialfällen deutlich schneller. Das ist er deshalb, weil er gezielt danach sucht bzw. damit beginnt, aber ich betrachte das als bedeutungslos. Der eigentliche Halleysche Algorithmus kommt ja nur zum Tragen, wenn eben kein Trivialfall vorliegt. Und zumindest bei meinen bisherigen Codes ist bei den komplexen Koeffizienten und Lösungen Cardano erheblich schneller als Halley. Ich betrachte das als schönen Erfolg. Gerade die Behauptung des Wikipediaartikels (Zitat: "bei sorgfältiger Implementierung") ist ziemlich absurd, denn die Cardano-Lösung wurde alles andere als sorgfältig implementiert. Eigentlich ist es nun Zeit, einen entsprechenden Artikel zu schreiben, und den Wikipedia-Artikel umzuschreiben... Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.91, eingetragen 2021-11-12 19:25

\quoteon(2021-11-11 21:02 - MontyPythagoras in Beitrag No. 90) ... heruntergeladen und ein wenig damit herumprobiert... um den Faktor 1.3,... \quoteoff Wenn Du hier keine absoluten Zahlen nennen willst, kannst Du sie mir auch per "private Nachricht" zukommen lassen. Mich interessiert das schon (natürlich nur die mit 3 reellen Ergebnissen), da nur so auch folgende Vergleiche möglich werden: - Python zu c++ - Deine CPU (welche genau) zu meinen beiden - lief die AVX2.exe, & wenn JA, um wieviel schneller ? Ich hatte dann noch versucht, meine acos-Version unter 40 ns zu bekommen, aber einige Varianten hatten sogar (vermutlich wegen Speicherverschiebung oder Vergrößerung des nötigen CPU-Caches) negativen Einfluss auf die schnellste 35 ns Version (dann 38 ns)... Grüße Gerd


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.92, vom Themenstarter, eingetragen 2021-11-14 11:56

Hallo HyperG, warum sollte ich damit ein Problem haben. Also meinetwegen gern: SSE: https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_i7_SSE.png AVX: https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_i7_AVX.png Rechner: Surface Pro 6, i7/16GB Prozessordaten hier Energieoptionen: "Beste Leistung", aber im Akkubetrieb (Frequenz schwankend um 3,9GHz) Aber nochmal: mich interessieren die absoluten Zahlen nicht. Es ging mir um die generelle Aussage des Wikipediaartikels, dass der Halley-Deiters-Algorithmus schneller sei als die direkte Berechnung mit Cardano. Natürlich sollte man das nicht unbedingt mit einer Interpretersprache wie Python messen, sondern mit einer schnellen Compilersprache wie C++ oder Fortran. Aber Vergleiche von Absolutwerten sagen mehr über den Prozessor als den Algorithmus aus, und sind daher für mich nicht interessant. Ich habe übrigens noch einige Genauigkeitsvergleiche gemacht und einen verbesserten Algorithmus auf Basis des multidimensionalen Newtonverfahrens gemacht. Anstatt also jede der drei Nullstellen durch die lokale Anwendung des einfachen Newtonverfahrens zu verbessern, was nachgewiesenermaßen wenig bringt, kann man durch die simultane Optimierung aller drei Nullstellen noch deutlich Genauigkeit rausholen und den Halley sogar noch übertreffen. Bei zufällig verteilten Lösungen sind die Genauigkeiten ("Stellenanzahl") wie folgt: 15,59 Cardano ohne Optimierung 15,68 Cardano mit Newton lokal 15,81 Halley-Deiters 15,93 Cardano mit multidimensionalen Newton Dazu später mehr. Ciao, Thomas


   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 322
  Beitrag No.93, eingetragen 2021-11-14 13:54

\quoteon(2021-11-14 11:56 - MontyPythagoras in Beitrag No. 92) ... Ich habe übrigens noch einige Genauigkeitsvergleiche gemacht und einen verbesserten Algorithmus auf Basis des multidimensionalen Newtonverfahrens gemacht. Anstatt also jede der drei Nullstellen durch die lokale Anwendung des einfachen Newtonverfahrens zu verbessern, was nachgewiesenermaßen wenig bringt, kann man durch die simultane Optimierung aller drei Nullstellen noch deutlich Genauigkeit rausholen und den Halley sogar noch übertreffen. Bei zufällig verteilten Lösungen sind die Genauigkeiten ("Stellenanzahl") wie folgt: 15,59 Cardano ohne Optimierung 15,68 Cardano mit Newton lokal 15,81 Halley-Deiters 15,93 Cardano mit multidimensionalen Newton Dazu später mehr. Ciao, Thomas \quoteoff \quoteon(2021-11-02 11:38 - MontyPythagoras in Beitrag No. 61) ... um das von mir vorgeschlagene Verfahren zur Untersuchung der Genauigkeit noch einmal aufzugreifen, würde ich folgendes vorschlagen: Wenn man die Lösungen vorher festlegt und erst danach die Koeffizienten berechnet, macht man bei der Berechnung der Koeffizienten schon Rundungsfehler, siehe Punkt 3 in meinem vorigen Beitrag. Das kann man aber elegant umgehen! 1. Wir nehmen Lösungen, deren Real- und Imaginärteil jeweils im Bereich von -10 bis 10 liegen, aber in Schritten von 0,1. Das lässt sich ja relativ leicht bewerkstelligen: Man nimmt einfach ganzzahlige komplexe Zahlen von -100 bis +100 und teilt durch 10. 2. Nun berechnet man die Koeffizienten der kubischen Gleichung. Dabei entstehen Rundungsfehler, aber man weiß, dass die "genauen Koeffizienten" jeweils Vielfache von 0,1, bei B von 0,01, bei C von 0,001 und bei D von 0,0001 sind. Man rundet die Koeffizienten jeweils auf diese Schrittweite. 3. Man berechnet per Halley oder Cardano die Lösungen. Indem man diese dann wieder auf das nächste Zehntel rundet, kennt man die ursprünglich vorgesehene, genaue Lösung. Sprich: der Abstand des Real- und des Imaginärteils zum nächsten Zehntel ist die tatsächliche Ungenauigkeit. Ciao, Thomas \quoteoff Hallo MontyPythagoras, ich nehme an, dass die Genauigkeiten in Beitrag No. 92 für den Wertebereich der Koeffizienten aus Beitrag No. 61 gelten, d.h. Real- und Imaginärteil jeweils im Bereich von -10 bis 10. Wenn man mit begrenzter Double-Genaugkeit (ca. 1.0e-16) rechnet, hängt der Fehler natürlich vom Wertebereich der Koeffizienten ab. Für den reellen Fall mit Koeffizienten zwischen -1000 und +1000 erhalte ich deutlich größere Fehler. Damit man den Überblick behält, sollten wir möglichst immer auch den Wertebereich der Koeffizienten, ob es um den rellen oder komplexen Fall geht und die Anzahl der Eingangsdaten mit angeben. Allgemein stellt sich die Frage, welcher Wertebereich der Koeffizienten sinnvoll ist. Das hängt natürlich vom Anwendungsfall ab, aber -10 bis +10 bei Verwendung von double halte ich für relativ klein.


   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 322
  Beitrag No.94, eingetragen 2021-11-14 14:13

Hallo, Routinen zur Berechnung der Nullstellen eines Polynoms dritten Grades wurden eigentlich genug vorgestellt. Trotzdem hier noch der Hinweis auf die GSL - GNU Scientific Library, siehe auch GNU Scientific Library. Die GSL enthält die Routine solve_cubic.c für den reellen Fall, und die Routine zsolve_cubic.c für den komplexen Fall. Der Ansatz ist ähnlich, wie in den Numerical Recipes (siehe die Beiträge No. 52 und No. 54).


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.95, vom Themenstarter, eingetragen 2021-11-14 18:19

Hallo zusammen, \quoteon(2021-11-14 14:13 - AlphaSigma in Beitrag No. 94) Routinen zur Berechnung der Nullstellen eines Polynoms dritten Grades wurden eigentlich genug vorgestellt. \quoteoff Ja, und praktisch alle waren Scheiße. Beispiele gefällig? Deiters und Konsorten (Ausschnitt): \sourceon C++ \numberson 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; \sourceoff Davon abgesehen, dass man die sechste Wurzel überhaupt nicht braucht, warum multipliziert man hier in Zeile 1 die pow-Funktion mit 2 (ausdrücklich als Fließkommazahl), wenn man danach bei jeder einzelnen Verwendung des Wertes $p$ wieder mit 0,5 multipliziert? Wenngleich hier die Multiplikation mit Potenzen von 2 vermutlich keine Rundungsfehler hinzufügt, sind das vier unnütze Multiplikationen. Bei einem Algorithmus aus einer Bibliothek interessiert mich nicht mehr die Lesbarkeit, sondern vor allem die Anzahl der arithmetischen Operationen, da diese unmittelbaren Einfluss haben auf a) Geschwindigkeit und b) Präzision (wg. Rundungsfehlern). Noch ein Beispiel, aus dem soeben von Dir zitierten Code solve_cubic.c. Zähl da nur mal durch, wie oft er "a/3" berechnet. Warum zum Teufel speichert er den Wert nicht einmal ab, z.B. unter der Variablen "a_3", und fertig? Und alle Codes, die ich gesehen habe, berechnen die drei Lösungen nach Bilderbuch wie im Wikiartikel, in diesem Beispiel über die Winkel 0°, 120° und -120° (das machen Deiters und Co. tatsächlich besser). Die trigonometrischen Funktionen sind zeitraubend, genauso wie die Potenzfunktion. Warum macht sich niemand zunutze, dass die Summe der drei Lösungen $-\frac ba$ ergibt, man also durch eine simple, hundert mal schnellere Addition/Subtraktion die dritte Lösung berechnen kann, wenn man die ersten beiden Lösungen schon hat? Und dass man mit einem einzigen Aufruf von atan2 und sincos die 3 Lösungen bzw. im Falle nur einer reellen Lösung durch einen einzigen Aufruf von cbrt berechnen kann, scheint auch keiner auf die Kette zu kriegen. Das ist alles schlecht, durch die Bank! Nebenbei: der Code zsolve_cubic.c ist nicht für komplexe Koeffizienten, sondern nur (konjugiert) komplexe Lösungen von kubischen Gleichungen mit reellen Koeffizienten. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.96, vom Themenstarter, eingetragen 2021-11-14 21:09

Hallo zusammen, ich habe noch weiter darüber nachgedacht, wie man die Präzision der Algorithmen weiter verbessern kann. Hier geht es wieder um die allgemeine Lösung für komplexe Koeffizienten. Außer bei Doppel- und Dreifach-Nullstellen war der Halley-Algorithmus bisher am genauesten (siehe mein Beitrag #64). Deiters et al. haben bei der Implementierung des Cardano-Algorithmus noch jeweils eine Newtonschleife angehängt, um die Präzision zu verbessern, und zwar nach dem bekannten Prinzip $$x_1=x_0-\frac{f(x_0)}{f'(x_0)}$$wobei $$f(x)=ax^3+bx^2+cx+d$$ist, $x_0$ die vorher bestimmte, "ungenaue" Nullstelle und $x_1$ die "verbesserte". Wie ich schon beschrieben hatte, bringt das Newton-Verfahren auf diese Art und Weise nicht viel, weil in der Nähe einer "echten" Nullstelle die Funktionswerte im Rundungsfehlerrauschen untergehen. Tatsächlich bringt das Verfahren nur eine Verbesserung der Genauigkeit gegenüber "Cardano einfach" von grob 0,1 Stellen. Aber man kann durch das mehrdimensionale Newton-Verfahren die Genauigkeit doch noch recht deutlich erhöhen. Seien $x_i$ mit $i=1..3$ die vorher mit den Cardano-Verfahren gefundenen, "ungenauen" Lösungen. Es soll nun jeweils ein $\varepsilon_i$ bestimmt werden, so dass $x_i+\varepsilon_i$ eine genauere Lösung darstellt. Mit dem einfachen Newton-Verfahren wäre daher $$\varepsilon_i=-\frac{f(x_i)}{f'(x_i)}$$Diese Variante nutzt den Funktionswert als Optimierungsziel. Beim mehrdimensionalen Verfahren ziehe ich stattdessen das möglichst genaue Treffen der Koeffizienten heran, denn es gilt ja: $$ax^3+bx^2+cx+d=a\prod_i\left(x-(x_i+\varepsilon_i)\right)$$Multipliziert man das aus und linearisiert, indem man Produkte $\varepsilon_i\varepsilon_j$ vernachlässigt, erhält man folgendes Gleichungssystem: $$\left(\begin{matrix}1 & 1 & 1 \\x_2+x_3 & x_1+x_3 & x_1+x_2\\x_2x_3 & x_1x_3 & x_1x_2 \end{matrix}\right)\vec\varepsilon=\left(\begin{matrix}-\frac ba-(x_1+x_2+x_3) \\\frac ca-(x_1x_2+x_1x_3+x_2x_3)\\-\frac da-x_1x_2x_3\end{matrix}\right)=\left(\begin{matrix}d_1\\d_2\\d_3\end{matrix}\right)$$Im Idealfall, also wenn man die genaue Lösung gefunden hat, sollten $d_1,d_2,d_3=0$ sein. Wenn man das lineare Gleichungssystem löst, erhält man nach einiger Rechnerei: $$\varepsilon_i=\frac{d_1x_i^2-d_2x_i+d_3}{(x_i-x_j)(x_i-x_k)}$$$j$ und $k$ sind jeweils die anderen beiden Indizes, die nicht $i$ sind. Vergleicht man das mit dem lokalen Newton, dann kann man dies auch darstellen als $$\varepsilon_i=-\frac{f(x_i)}{f'(x_i)+2d_1x_i-d_2}$$Man beachte somit den Unterschied im Nenner. Allerdings ist diese Gleichung rundungsfehlertechnisch schlechter als die vorige Gleichung. Umgesetzt habe ich das wie folgt (Code versteckt, nur um Platz zu sparen): \showon \sourceon Python \numberson rcpr3 = 1 / 3 eps1 = complex(-0.5, 0.75 ** 0.5) eps2 = complex(-0.5, - 0.75 ** 0.5) def cbrt_newt_3(a, b, c, d): if a == 0: # (Trivialfall quadratische Gleichung) if b == 0: return (-d / c, ) # (Trivialfall lineare Gleichung) c /= (-2 * b) # (Von hier an echte quadratische Gleichung) s = (c * c - d / b) ** 0.5 return c + s, c - s s = -rcpr3 / a # (Von hier an echte kubische Gleichung) r = b * s q = r * r t = q * r p, q = c * s + q, t + 1.5 * s * (d + c * r) s = (q * q - p * p * p) ** 0.5 s = q if abs(s) < 2e-5 * abs(q) else q - s if q.real < 0 else q + s if abs(s) <= 1e-14 * abs(t): return r, r, r s **= rcpr3 q = s * eps1 t = s * eps2 x0, x1, x2 = s + p / s + r, q + p / q + r, t + p / t + r p = b + a * (x0 + x1 + x2) q = c - a * (x0 * (x1 + x2) + x1 * x2) r = d + a * x0 * x1 * x2 x01 = x0 - x1 x12 = x1 - x2 x20 = x2 - x0 s = a * x01 * x20 if s != 0: x0 += ((p * x0 + q) * x0 + r) / s s = a * x01 * x12 if s != 0: x1 += ((p * x1 + q) * x1 + r) / s s = a * x20 * x12 if s != 0: x2 += ((p * x2 + q) * x2 + r) / s return x0, x1, x2 \sourceoff \showoff Das hat folgende Wirkung: https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Precision_cubic_algorithm.png Gegenüber der reinen Cardano-Lösung hat dieser Algorithmus die Präzision um rund 0,5 Stellen verbessert, während die einfache, lokale Newtonvariante nur rund 0,1 Stellen schafft. In dieser Auswertung wurden die Lösungen zufällig zusammengesetzt im Bereich von -100 bis +100, jeweils für reellen Anteil und imaginären Anteil, und 0,1 Schrittweite, was $2001\times2001=4004001$ mögliche Werte für die konstruierten Lösungen $x_i$ bedeutet. Man beachte rechts am Rand, wie viel häufiger der optimierte Algorithmus die genaue Lösung trifft als die anderen. Die Zahlen in der Legende geben die Durchschnittswerte an. \quoteon(2021-11-14 13:54 - AlphaSigma in Beitrag No. 93 Wenn man mit begrenzter Double-Genaugkeit (ca. 1.0e-16) rechnet, hängt der Fehler natürlich vom Wertebereich der Koeffizienten ab. Für den reellen Fall mit Koeffizienten zwischen -1000 und +1000 erhalte ich deutlich größere Fehler. Damit man den Überblick behält, sollten wir möglichst immer auch den Wertebereich der Koeffizienten, ob es um den reellen oder komplexen Fall geht und die Anzahl der Eingangsdaten mit angeben. Allgemein stellt sich die Frage, welcher Wertebereich der Koeffizienten sinnvoll ist. Das hängt natürlich vom Anwendungsfall ab, aber -10 bis +10 bei Verwendung von double halte ich für relativ klein. \quoteoffJa, kann sein. Daher habe ich den Wertebereich hier vergrößert. Hast Du bei Deinem Feld von -1000 bis +1000 eine gewisse Schrittweite verwendet? Ich würde es nicht größer machen wollen, als ich habe. Bedenke bitte, dass ich hier die Lösung für komplexe Koeffizienten berechne. Die Koeffizienten, die ich aus den vorgegebenen Lösungen konstruiere, werden berechnet, indem ich den komplexen Lösungen ganze Zahlen für Real- und Imaginärteil zuweise, aber dadurch werden die Zahlen automatisch in double umgewandelt, denn Python kennt keine ganzzahligen komplexen Zahlen. Bei der Berechnung der Koeffizienten muss ich daher sichergehen, dass die double Multiplikationen und Additionen keine Rundungsfehler erzeugen. Am Ende lese ich den Real- und Imaginärteil noch einmal aus und runde sie jeweils auf eine Ganzzahl, und teile dann durch 10, so dass der Koeffizient korrekt sein sollte. Ich könnte mir natürlich die Mühe machen, und "zu Fuß" nur über Ganzzahlwerte die Koeffizienten berechnen, aber das war mir bislang zu aufwendig. Besonders bösartige Konstruktionen habe ich bislang außen vor gelassen, da ich eher auf "real world" Anwendungen hin optimieren möchte. (Klar kann man sich fragen, wann denn in der "real world" jemals komplexe Koeffizienten auftreten... 😂) Ciao, Thomas


   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 322
  Beitrag No.97, eingetragen 2021-11-14 21:34

\quoteon(2021-11-14 21:09 - MontyPythagoras in Beitrag No. 96) ... Hast Du bei Deinem Feld von -1000 bis +1000 eine gewisse Schrittweite verwendet? ... Ciao, Thomas \quoteoff Ich habe die Funktion drand48 verwendet. Die erzeugt gleichmäßig verteilte double Werte im Intervall [0; 1[ mit 48-Bit Integer-Arithmetik. Die Werte habe ich dann mit (max - min)*x + min auf [-1000: 1000[ transformiert. Ich erzeuge die Nullstellen zufällig und berechne daraus die Polynomkoeffizienten mit erhöhter Genauigkeit (30 Dezimalstellen). Das funktioniert momentan bei mir nur für den Fall mit rellen Koeffizienten und reellen Nullstellen. Den Fehler berechne ich in double Genauigkeit als Absolutbetrag der Differenz zwischen den erzeugten und den berechneten Nullstellen.


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.98, vom Themenstarter, eingetragen 2021-11-14 22:10

Hallo zusammen, ich möchte das gerne weiterhin in Python umsetzen, wo es diese Funktion so 1:1 nicht gibt. Ich kann aber in Python das Package mpmath einbinden, welches grundsätzlich fast unbegrenzte Präzision ermöglicht. Andere Variante: wir erzeugen eine ASCII- oder wie auch immer geartete Binärdatei, die man einlesen kann und wo die Koeffizienten und "genauen" Lösungen von einer großen Anzahl an kubischen Gleichungen im (complex) double Format abgelegt sind. So würden wir alle gemeinsam mit dem gleichen Datensatz arbeiten. Die Datei könntest Du ja per C++ erzeugen, und ich dann auch in Python einlesen. Ciao, Thomas


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 2948
  Beitrag No.99, eingetragen 2021-11-14 22:33

\quoteon(2021-11-14 18:19 - MontyPythagoras in Beitrag No. 95) \quoteon(2021-11-14 14:13 - AlphaSigma in Beitrag No. 94) Routinen zur Berechnung der Nullstellen eines Polynoms dritten Grades wurden eigentlich genug vorgestellt. \quoteoff Ja, und praktisch alle waren Scheiße. \quoteoff Wenn ein Beitrag mit so einer "alle sind blöd außer mir"-Aussage anfängt, sollte man vielleicht lieber die Finger von einer Antwort lassen, aber auf einen Punkt möchte ich doch eingehen. \quoteon(2021-11-14 18:19 - MontyPythagoras in Beitrag No. 95) Zähl da nur mal durch, wie oft er "a/3" berechnet. Warum zum Teufel speichert er den Wert nicht einmal ab, z.B. unter der Variablen "a_3", und fertig? \quoteoff Es gibt keinen Grund, sich um derartige Optimierungen zu kümmern, weil das der Compiler macht. Die Code-Fragmente, die du ansprichst, haben ja etwa die folgende Form: \sourceon C void fun (double a, double *x0, double *x1, double *x2) { *x0 = 1 - a / 3; *x1 = 2 - a / 3; *x2 = 3 - a / 3; } \sourceoff Daraus macht der Compiler: \sourceon Disassembly of section .text 0000000000000000 : 0: f2 0f 59 05 00 00 00 mulsd 0x0(%rip),%xmm0 # 8 7: 00 8: f2 0f 10 0d 00 00 00 movsd 0x0(%rip),%xmm1 # 10 f: 00 10: f2 0f 5c c8 subsd %xmm0,%xmm1 14: f2 0f 11 0f movsd %xmm1,(%rdi) 18: f2 0f 10 0d 00 00 00 movsd 0x0(%rip),%xmm1 # 20 1f: 00 20: f2 0f 5c c8 subsd %xmm0,%xmm1 24: f2 0f 11 0e movsd %xmm1,(%rsi) 28: f2 0f 10 0d 00 00 00 movsd 0x0(%rip),%xmm1 # 30 2f: 00 30: f2 0f 5c c8 subsd %xmm0,%xmm1 34: f2 0f 11 0a movsd %xmm1,(%rdx) 38: c3 retq \sourceoff Wie du siehst, wird der Ausdruck $a/3$ nicht für jede Zeile neu, sondern nur einmal am Anfang berechnet. Die von dir vorgeschlagene Einführung einer Hilfsvariablen a_3 = a / 3 führt zu exakt dem gleichen Ergebnis und bringt somit nichts. [Die Antwort wurde nach Beitrag No.97 begonnen.]


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.100, vom Themenstarter, eingetragen 2021-11-15 10:11

Hallo zusammen, \quoteon(2021-11-14 22:33 - zippy in Beitrag No. 99) Wenn ein Beitrag mit so einer "alle sind blöd außer mir"-Aussage anfängt, sollte man vielleicht lieber die Finger von einer Antwort lassen, aber auf einen Punkt möchte ich doch eingehen. \quoteoff Danke, Mama. Ich habe nicht gesagt, dass ich der weltbeste Programmierer bin, sondern dass zumindest die verlinkten Codes schlecht waren. Ich wüsste auch nicht, wo ich irgendjemanden als "blöd" bezeichnet hätte und bitte um ein entsprechendes Zitat. Ich würde sogar Geld drauf wetten, dass sowohl Du als auch AlphaSigma es besser programmieren könnten, wenn Ihr es selbst machen würdet. Ich bin ein Feind davon, vorgekaute Lösungen kritiklos zu übernehmen, nur weil sie in irgendeine Bibliothek übernommen wurden, die das Wort "scientific" trägt. Ist es deshalb automatisch gut? Natürlich bügelt ein guter Compiler manche Unzulänglichkeiten der Programmierer weg. Trotzdem darf man diese Unzulänglichkeiten ja wohl kritisieren, denn es gibt so etwas wie ein "Best Practice"-Prinzip. Eine Interpretersprache wie Python zum Beispiel wird hier x-mal hintereinander "a/3" berechnen. Das mit dem "a/3" ist eine Unzulänglichkeit, die der Compiler problemlos ausbügeln kann. Wie schön. Und was ist mit den anderen Beispielen? Dann erleuchte mich doch mal, ob Dein Compiler auch in dem Deiters Code die 4 unsinnigen Multiplikationen ausbügelt, oder die schon prinzipiell unnötige Berechnung der sechsten Wurzel, oder den Kosinus von drei unterschiedlichen Winkeln, um die drei Lösungen zu berechnen, wenn es auch zwei oder sogar eine einzige Berechnung tut. Offenbar ist Dir entgangen, dass der von mir programmierte Code aus Beitrag 80 tatsächlich erheblich schneller ist als zumindest der Code von Deiters, da er bei kubischen Gleichungen mit drei Lösungen den Halley-Deiters hinter sich lässt. Also erzähl hier bitte nicht, dass es der Compiler schon richten wird, wo das offenkundig nicht der Fall ist. Oberstes Prinzip sollte sein (zumindest war es das mal), dass man erst einmal eine effiziente Lösung umsetzt, und den Rest macht der Compiler. Denn egal wie effizient der Compiler arbeitet: eine sechste Wurzel, die man gar nicht braucht, sollte man ihn erst gar nicht berechnen lassen. Ciao, Thomas


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 2948
  Beitrag No.101, eingetragen 2021-11-15 10:15

Der Ton deiner Antwort bestätigt leider meine Bedenken. Ich hätte nicht antworten sollen. \quoteon(2021-11-14 22:33 - zippy in Beitrag No. 99) Wenn ein Beitrag mit so einer "alle sind blöd außer mir"-Aussage anfängt, sollte man vielleicht lieber die Finger von einer Antwort lassen \quoteoff


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.102, vom Themenstarter, eingetragen 2021-11-15 12:27

Hallo zippy, ich entschuldige mich für das "Danke, Mama". Das war unangebracht. An dem Rest halte ich der Sache nach fest. Es geht hier darum, einen möglichst schnellen Algorithmus zu finden, denn das war der Aufhänger des ganzen Threads. Wir haben den Cardano durchoptimiert, so dass er jetzt manchmal schneller, manchmal nur knapp langsamer ist als Halley-Deiters. Wir könnten uns nun den Halley-Deiters vornehmen und auch durchoptimieren, am Ende ist er vielleicht wieder durchgängig schneller als Cardano. Das wäre lustig. Aber es ist unbestreitbar, dass die verlinkten Codes nicht auf Geschwindigkeit optimiert waren. Das muss nicht an mangelndem Können der Programmierer gelegen haben, vielleicht lag der Schwerpunkt auf Lesbarkeit oder Didaktik. In diesem Beispiel: \sourceon C++ \numberson 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; \sourceoff wird der Compiler "merken", dass dreimal "0.5*p" berechnet wird, wird es nur in Zeile 2 berechnen und den Wert "0.5*p" für Zeile 3 und 4 abspeichern. Aber es würde mich sehr wundern, wenn er die Sinnlosigkeit der ersten Multiplikation mit 2.0 und aller anderen danach erkennen würde. Statt vier Multiplikationen, die hier eine Interpretersprache macht, macht der C++ Compiler also nur zwei - vermute ich. Aber notwendig sind null. Von den anderen Themen ganz zu schweigen. Ich gebe zu, dass es mich ärgert, wenn jemand sagt "Es ist in einer Bibliothek auf Github, es trägt den Namen 'scientific', also ist es gut". Das kann es nicht sein, das ist innovationsfeindlich. Wir sind ja hier nicht in der Grundschule, wo der Schüler gelobt werden muss, auch wenn er eine schlechte Arbeit abgeliefert hat. Es ist eine objektiv messbare Tatsache, dass der Code aus #80 schneller ist als der Cardano-Algorithmus bei Deiters. Vielleicht gibt es irgendwo einen Algorithmus, Halley oder sonst wie, der noch schneller ist. Ich möchte nur die Aussage in dem Wiki-Artikel prüfen, und am Ende ist sie wahr oder falsch. Und wenn sie falsch ist, dann sollten wir den Artikel korrigieren. Offenkundig interessiert Dich das Thema ja. Daher würde ich mich freuen, wenn Du Deine unbestreitbaren Fachkenntnisse mit einbringen würdest, anstatt mit schwachen Argumenten schlechte Programmiertechniken zu verteidigen. Wen das Thema nicht genug interessiert oder es für verschenkte Zeit hält, weil es ja schon so viel Code gibt, der muss ja nicht mitmachen. Ich denke aber, wir haben schon einige interessante Fakten zu Tage gefördert. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.103, vom Themenstarter, eingetragen 2021-11-15 17:09

Hallo AlphaSigma, \quoteon(2021-11-14 21:34 - AlphaSigma in Beitrag No. 97) Ich habe die Funktion drand48 verwendet. Die erzeugt gleichmäßig verteilte double Werte im Intervall [0; 1[ mit 48-Bit Integer-Arithmetik. Die Werte habe ich dann mit (max - min)*x + min auf [-1000: 1000[ transformiert. Ich erzeuge die Nullstellen zufällig und berechne daraus die Polynomkoeffizienten mit erhöhter Genauigkeit (30 Dezimalstellen). Das funktioniert momentan bei mir nur für den Fall mit rellen Koeffizienten und reellen Nullstellen. Den Fehler berechne ich in double Genauigkeit als Absolutbetrag der Differenz zwischen den erzeugten und den berechneten Nullstellen. \quoteoff Für Polynome mit drei deutlich verschiedenen Lösungen wird das funktionieren, aber bei Polynomen mit einer doppelten oder dreifachen Nullstelle erzeugst Du eine vermeidbare Ungenauigkeit. Du hast Polynomkoeffizienten, die 30 Dezimalstellen haben, schneidest sie dann aber auf double precision, also auf 16 Dezimalstellen, ab (wir reden hier ja eigentlich von binärer Darstellung, und auch in double precision sind die Dezimalstellen ab Stelle 17 nicht wirklich null, aber Du weißt schon, was ich meine). Die Lösungen, die Du vorher erzeugt hattest, sind jetzt aber nicht mehr die Lösungen >>dieses<< Polynoms mit diesen double-Koeffizienten, sondern nach wie vor eines anderen Polynoms mit 30-stelligen Koeffizienten. Kleine Abweichungen bei den Koeffizienten, und sei es nur der Rundungsfehler, erzeugen eine neue, exakte Lösung, die deutlich mehr als nur den typischen Rundungsradius von der ursprünglichen Lösung abweicht, wenn die Ableitung der Funktion nahe null ist. Man kann einem mit double precision arbeitenden Algorithmus nicht anlasten, dass er die Dezimalstellen 17 bis 30 nicht ahnt. Der Weg muss daher meiner Meinung nach sein, dass man die double-Koeffizienten als die tatsächlich "genauen" Koeffizienten betrachtet, dass also die Stellen 17 bis 30 tatsächlich null sind. Für diese beschnittenen double-Koeffizienten berechnen wir mit einem "Gott"-Algorithmus mit unbegrenzter Präzision die exakten Lösungen, und runden diese auf den nächsten double precision Wert. Das ist dann der bestmögliche Wert, den ein double precision Algorithmus finden könnte. Realistisch sollte man also mit mindestens dreifacher Stellenanzahl eines double-Wertes die Lösungen berechnen, da im Falle einer Dreifach-Nullstelle die Präzision gedrittelt wird, siehe mein Beitrag #62. Da ich mit dem Python Package "mpmath" eine beliebige Genauigkeit festlegen kann, könnte ich ja problemlos die kubische Gleichung mit 100 Stellen Genauigkeit lösen. Natürlich wäre das langsam, aber das ist egal, das macht man ja nur einmal. Wahrscheinlich gibt es solche Librairies wie "mpmath" für C++ auch, aber da bin ich zu lange raus. Mathematica und andere CAS könnten das mit Sicherheit auch. Ciao, Thomas


   Profil
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 322
  Beitrag No.104, eingetragen 2021-11-15 18:47

\quoteon(2021-11-15 17:09 - MontyPythagoras in Beitrag No. 103) ... Wahrscheinlich gibt es solche Librairies wie "mpmath" für C++ auch, aber da bin ich zu lange raus. Mathematica und andere CAS könnten das mit Sicherheit auch. Ciao, Thomas \quoteoff Hallo MontyPythagoras, leicht off-topic: Bereits dc auf einer PDP-11 und sein Nachfolger bc können mit beliebiger Genauigkeit rechnen, d.h. die Genauigkeit ist nur durch den Speicher limitiert. dc gab es schon vor C auf Unix. Die Syntax von dc ist wohl nur für Leute, die vorher mit Lochkarten programmiert haben (oder FORTH Programmierer) komfortabel. Hier ein Beispiel zur Berechnung von e auf 2570 Stellen: Rosetta Code - dc. Das ging schon ca. 1970. Heute gibt es dafür u.a. die GNU Multiple Precision Arithmetic Library und darauf basierend die GNU MPC für komplexe Zahlen mit beliebiger Genauigkeit. Auch mpmath für Python basiert zum Teil auf GMP/MPIR. Zitat: "mpmath internally uses Python's builtin long integers by default, but automatically switches to GMP/MPIR for much faster high-precision arithmetic if gmpy is installed or if mpmath is imported from within Sage." von mpmath.org. Ich verwende momentan bc. Wenn ich etwas mehr Zeit habe, werde ich GNU MPC testen.


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.105, vom Themenstarter, eingetragen 2021-11-15 20:02

Hallo AlphaSigma, ja, ich war überzeugt, dass es so etwas für C++ gibt. Ich habe das letzte Mal aktiv mit C++ programmiert im Jahre 1996. Von den Lochkarten war ich also gar nicht so weit entfernt. 😂 Mit Programmiersprachen ist es wie mit Fremdsprachen: wenn man lange nicht gesprochen hat, ist der passive Wortschatz noch da, aber es hapert am aktiven. Ich kann den C++ Code nach wie vor problemlos lesen, aber selber schreiben ist arg mühselig. Daher freue ich mich auch, dass Ihr, Du und HyperG, mir da so aktiv helft, falls das nicht rübergekommen sein sollte. Was hältst Du von dem Gedanken in #103? Ich hoffe, ich konnte mich verständlich ausdrücken. Ich denke, wir sollten ein paar Testdateien für unterschiedliche Zwecke schreiben: a) Eine Datei für komplexe Koeffizienten und Lösungen b) Eine Datei für reelle Koeffizienten und drei reelle Lösungen c) Eine Datei für reelle Koeff. und nur eine reelle Lösung Wobei die Frage ist, ob man b) und c) nicht in eine Datei scheibt. Aber dann folgt sofort die Frage, zu welchen Anteilen. Und dann wäre da auch noch die Frage nach den ungünstigen Konstellationen, also Doppel- und Dreifachnullstellen. Ignoriert man sie? Wenn sie zufällig vorkommen, ist gut, und wenn nicht, auch gut? Widmet man ihnen sogar eine eigene Testdatei? Und auch hier: wie gewichtet man sie gegenüber den unproblematischen? Vielleicht sollte man alles getrennt betrachten und es dem jeweiligen Leser oder Anwender überlassen, auf was er Wert legt. Ich muss da noch ein wenig drauf herumdenken. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.106, eingetragen 2021-11-16 18:01

\quoteon(2021-11-15 18:47 - AlphaSigma in Beitrag No. 104) ... die GNU Multiple Precision Arithmetic Library und darauf basierend die GNU MPC für komplexe Zahlen mit beliebiger Genauigkeit... \quoteoff Das dachte ich anfangs auch. Also habe ich mal "nur einen einfachen asin-Wert" berechnet (ohne komplexen Anteil) -> und schon nach etwa 300000 Dezimalstellen war der Rest falsch! Hier ist also IMMER eine Validierungsberechnung nötig. Ich konnte zwar noch etwas herausholen (oder mit anderen Funktionen umgehen), aber alle Funktionen, die komplizierter als die Wurzel sind, kommen nicht über 800000 "echte Stellen". Mathematica kann zwar mehr, aber beim Ziegenfaktor kam ich auch nicht über 1,1 Mrd., da der RAM von 32 GB die Begrenzung darstellte. Richtig schnell und mehr als 1 Mrd. Stellen -> da kenne ich nur YMP vom Pi-Weltmeister. Nachteil: alles komplizierter als Wurzel muss man selbst programmieren. Kennt jemand schnelle asin- oder atan-Algorithmen, die schneller als die bekannten Fakultäts-Reihen konvergieren? (also z.B. wie die quadratisch konvergierende Iterationen von Newton) Grüße P.S.: @MontyPythagoras: Du hattest doch mal eine Jacobi-Iteration vorgestellt - oder? Kann man die Gleichung ArcSin[z] == Pi/2 - InverseJacobiCD[z, 0] dafür umstellen?


   Profil
MontyPythagoras
Senior Letzter Besuch: im letzten Monat
Dabei seit: 13.05.2014
Mitteilungen: 2864
Wohnort: Werne
  Beitrag No.107, vom Themenstarter, eingetragen 2021-11-17 14:26

Hallo zusammen, @HyperG, zu Deinem P.S.: in meinem Artikel über nichtlineare Biegung kommen die Jacobischen elliptischen und Theta-Funktionen vor, aber nur am Rande. Außerdem habe ich mal in einem anderen Zusammenhang darauf hingewiesen, dass es sehr schnell konvergierende Algorithmen für die elliptischen Integrale gibt. Dazu verwendet man die quadratische Transformation der symmetrischen Integrale nach Carlson, siehe hier. Ich würde das Thema aber nicht so gern in diesem Thread behandeln. Das können wir gerne per PM machen, oder mach dafür bitte einen eigenen Thread auf. Wenn mir dazu was einfällt, werde ich mich gern beteiligen. Ich werde drüber nachdenken. Ich habe wie in #105 Punkt a) beschrieben vier Dateien mit komplexen Koeffizienten erzeugt. Vier deshalb, weil ich für den ersten Koeffizienten $a$ und die drei Lösungen $x_i$ drei unterschiedliche Wertebereiche getestet habe, nämlich $\pm10$, $\pm10^3$, $\pm10^5$ und $\pm10^7$, jeweils reeller und imaginärer Anteil. Jede der Dateien ist knappe 11MB groß, weil sie jeweils 100000 Datensätze mit je 7 komplexen Zahlen enthalten. Die ersten 4 Zahlen sind die Koeffizienten $a$ bis $d$, dann folgen 3 Zahlen für $x_i$, jeweils complex double. Die $x_i$ sind wie beschrieben die Lösungen, die den exakten Lösungen in der complex double Präzision am nächsten kommen. Sie enthalten keine Spezialfälle wie Doppel- oder Dreifachnullstellen. Ergebnisse unten. Ich habe inzwischen festgestellt, dass der von mir in #96 beschriebene Algorithmus bekannt ist unter mehreren Namen, nämlich "Weierstraß" und "Durand-Kerner", siehe hier, wenngleich die Herleitung bei Wiki anders ist als meine. Dort ist im Zähler nur der Funktionswert $f(x_i)$ angegeben. Im Hinblick auf die Minimierung der Rundungsfehler ist jedoch die Darstellung $$\varepsilon_i=\frac{d_1x_i^2-d_2x_i+d_3}{(x_i-x_j)(x_i-x_k)}$$wie ich sie hergeleitet hatte, empirisch besser, vermutlich da sie die tatsächlichen Abweichungen $d_i$ von den genauen Koeffizienten berücksichtigt, und nicht den Funktionswert. Die Ergebnisse sind wie folgt: https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_1_Prec_cubic_1e1.pngBeim Lösungsbereich $\pm10$ decken sich die Ergebnisse ziemlich gut mit #96. Der einfache Newton bringt gegenüber Cardano so gut wie nix, Halley ist ein bisschen besser, Durand-Kerner deutlich besser mit einem Durchschnitt >16. https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_1_Prec_cubic_1e3.pngBeim Lösungsbereich $\pm1000$ bringt der einfache Newton einen recht deutlichen Sprung, die Genauigkeit ist besser als mit Halley, aber trotzdem mehr als 0,3 Stellen schlechter als Durand-Kerner, der wieder über 16 liegt. https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_1_Prec_cubic_1e5.pngBeim Lösungsbereich $\pm100000$ fallen Cardano und Halley schon deutlich ab, Durand erreicht wieder >16. https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_1_Prec_cubic_1e7.pngDer Trend setzt sich fort, Cardano+1x Newton stabil bei 15,7 Stellen und Cardano+1x Durand bei >16. Halley und Cardano fallen ab. Vor allem hat Halley hier schon einen recht hohen Anteil mit unter 15 genauen Stellen. Außerdem habe ich mit dem großen Wertebereich untersucht, ob man den Halley noch verbessern kann. Eine Variante war, nach Erreichen des Abbruchkriteriums noch eine Halley-Schleife nachzuschieben ("Halley +1"), sowie ebenfalls nach Abbruch dem Halley noch eine Newton-Schleife anzuhängen ("Halley + 1x Newton"). https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Prec_cubic_1e7_Halley.pngDie zusätzliche Halley-Schleife bringt nichts, die Genauigkeit mit und ohne nachgeschobene Schleife ist nahezu identisch (0,007 Stellen Unterschied). Allerdings bringt die Newton-Schleife viel, es hebt den Halley auf das Genauigkeitsniveau von Cardano mit Newton, was insofern interessant ist, als Newton ja nur quadratische, Halley dagegen kubische Konvergenz aufweist. Auch wenn es hier um komplexe Koeffizienten ging, dürfte es bei den reellen genauso aussehen. Der Vergleich von Deiters et al. hinkt daher wohl ein wenig. Im größeren Wertebereich müsste man dem Halley noch eine Newtonschleife anhängen, wenn man es bei Cardano auch tut. Oder eben bei beiden weglassen. Ich werde in meiner "persönlichen" Python-Bibliothek für kubische Gleichungen daher wohl einen schnellen, "unpräzisen" Algorithmus (Cardano ohne alles - klingt wie 'ne Pizza 🙃) und einen etwas langsameren, aber präziseren Algorithmus (Cardano + 1x Durand-Kerner) mit "Mehrfachnullstellen-Erraten-Heuristik" führen. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.108, eingetragen 2021-11-17 20:08

\quoteon(2021-11-17 14:26 - MontyPythagoras in Beitrag No. 107) ... nach Carlson, siehe hier. ... mach dafür bitte einen eigenen Thread auf. Wenn mir dazu was einfällt, werde ich mich gern beteiligen... \quoteoff Das mit Carson hatten wir hier schon. Ich bekomme aber die 4 Funktionen nicht zusammen nach ArcSin umgestellt: ArcSin[z] == Pi/2 - InverseJacobiCN[z, 0] ArcSin[z] == InverseJacobiSN[z, 0] EllipticK[z] == InverseJacobiSN[1, z] EllipticK[z] == InverseJacobiDN[Sqrt[1 - z], z] Auch "eigenen Thread" gibt's schon: hier asin aber alle nicht so richtig geeignet ...


   Profil
Folgende Antworten hat der Fragesteller vermutlich noch nicht gesehen.
AlphaSigma
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 23.11.2012
Mitteilungen: 322
  Beitrag No.109, eingetragen 2021-11-26 18:30

Hallo, hier noch ein interessanter, etwas älterer, Artikel von William Kahan zur Lösung kubischer Gleichungen (mit reellen Koeffizienten): To Solve a Real Cubic Equation (Cubic.pdf) Auf S. 16 werden Koeffizienten für kritische Fälle angegeben (9. Some Trial Data for Cubic Equation Solvers).


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1553
  Beitrag No.110, eingetragen 2021-11-26 19:47

\quoteon(2021-11-12 19:25 - hyperG in ...) ...schnellste 35 ns Version ... \quoteoff hier im anderen Beitrag des selben Forums konnte ich nun zeigen, dass man die langsame sin + cos Berechnung nochmals durch AVX-Befehlsoptimierung & durch Multithreading um Faktor 43 beschleunigen kann. Zusammen mit dem Python/c++ Faktor 4,7 -> sind es dann wirklich Welten :-) Grüße Gerd


   Profil
Seite 3Gehe zur Seite: 1 | 2 | 3  

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]