| |
| Seite 2 [1 2 3 4] | 4 Seiten |
| Autor |
 Kubische Gleichung: numerische Näherung schneller als Cardanische Formeln? |
|
AlphaSigma
Senior 
Dabei seit: 23.11.2012
Mitteilungen: 453
 |  Beitrag No.40, eingetragen 2021-10-25
|
 Profil
|
hyperG
Senior 
Dabei seit: 03.02.2017
Mitteilungen: 2117
 | Beitrag No.41, eingetragen 2021-10-26
|
\quoteon(2021-10-25 21:41 - MontyPythagoras in Beitrag No. 39)
...
warum eigentlich die trigonometrischen Funktionen und die Potenzfunktion SOO langsam sind.
...
\quoteoff
Das kann man sich mit folgender Tabelle schnell klar machen:
https://www.ardent-tool.com/CPU/docs/IIT/3c87dx.pdf
\sourceon Maschinenbefehle & nötige Takte Beispiel Coprozessor
Befehl Clock Count
FADD 12
FMUL 17
FMUL dMEM 35
FDIV dMEM 64
FSQRT 48
FCOS 306
FSIN 306
FSINCOS 343
FPTAN 297
FPATAN 337
\sourceoff
Natürlich variiert es von CPU/FPU etwas, aber die groben Verhältnisse ähneln sich.
Bei den alten Coprozessoren war es (wie schon zuvor berichtet) besonders langsam, da sie auch die 80 Bit (long double) unterstützen mussten.
Jede Nachkommastelle an Genauigkeit kostet Takte!
Deshalb hatte ich ja im Beitrag 33 die Genauigkeit mit angegeben.
Dann muss man auch unterscheiden, ob der Inputwert bereits im Rechenregister vorliegt, oder aus dem Speicher (MEM) geholt werden muss!
Bei FMUL habe ich mal beide angegeben: aus dem RAM mehr als doppelt so langsam!
Das einfach aussehende Pow(x,reell) wird meist mit
x^y=2^(log2(x)*y)
\
x^y=2^(log_2(x)*y)
also FYL2X, F2XM1 and FSCALE berechnet:
\sourceon c++
double Pow(double x, double y)
{
__asm
{
fld y
fld x
fabs
fyl2x
lea eax, controlword
fstcw word ptr [eax]
or byte ptr [eax+1], 0Ch
fldcw word ptr [eax]
fld st(0)
fld st(0)
frndint
fsubp st(1), st(0)
f2xm1
fld1
fadd
fscale
fstp x
lea eax, controlword
fldcw word ptr [eax]
}
return x;
}
\sourceoff
... und da kommen 'ne Menge Takte zusammen, weil ja auch die Register alle noch gefüllt werden müssen...
Die neuen Compiler kürzen alle etwas ab und nehmen Ungenauigkeiten in Kauf. Compilerschalter f:precise (möglichst genau) und f:fast (schneller aber etwas ungenauer) kommen noch hinzu.
Eine letzte Newton-Iteration ist also immer nötig, um volle double Genauigkeit zu erreichen.
Man kann bei cpow(z,1/3) (komplexe 3. Wurzel) immer weiter Genauigkeit herunterschrauben, um Geschwindigkeit herauszuholen. 8 Stellen wäre ideal, um durch die letzte Newton-Iteration auf 16 richtige Nachkommastellen zu kommen.
Jede Optimierung bei cpow(z,1/3) {z.B. mit AVX-Befehlen}, muss man natürlich gerechterweise auch bei Halley/Deiters anwenden. Und so werden die Algorithmen mit trogonom. Funktionen immer hinterherhinken, solange man keine besseren Clock Count Werte zu kaufen bekommt.
Ich habe auch schon überlegt, ob ich den Intel-Compiler teste... (Beitrag 33 angedeutet: 16 float Inputs, die parallel sin & cos berechnen).
Oder __m512 _mm512_pow_ps (__m512 a, __m512 b)
16 float-pow Berechnungen...
Das sind aber auch keine "echte Maschinenbefehle" {denn dann hätte ich ihn längst mit MASM selbst testen können}, sondern auch nur hoch optimierter Intel-Code...
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.42, eingetragen 2021-10-26
|
Vergleich Geschwindigkeit zu Genauigkeit bei x^1/3 =cbrt(x) unter C++:
Ideen aus http://metamerist.com/cbrt/cbrt.htm
\sourceon c++
32-bit float tests -- VC 2005 CPU 2008| VC2012 i7 | VC2017 i9 x86 VC2017 i9 AVX2
Algorithmus dt BitGenauigkeit dt Bit dt Bit dt Bit
cbrt_5f 8.8 ns 5 2.2 ns 5 1.4 ns 5 1.4 ns 5
powf 144.5 ns 23 22.7 ns 23 10.1 ns 23 10.1 ns 23
halley x 1 31.8 ns 15 3.8 ns 15 2.6 ns 15 2.7 ns 15
halley x 2 59.0 ns 23 6.7 ns 22 5.0 ns 22 5.2 ns 22
newton x 1 23.4 ns 10 3.2 ns 10 2.2 ns 10 1.9 ns 10
newton x 2 48.9 ns 20 5.7 ns 20 3.8 ns 20 3.6 ns 20
newton x 3 72.0 ns 23 8.6 ns 23 5.9 ns 23 5.5 ns 23
newton x 4 89.6 ns 23 12.5 ns 23 8.3 ns 23 7.5 ns 23
64-bit double tests ----- 64-bit double tests 64-bit double 64-bit double tests
cbrt_5d 23.5 ns 5 6.0 ns 5 5.0 ns 5 5.0 ns 5
pow 151.3 ns 52 38.6 ns 52 25.9 ns 52 26.0 ns 52
halley x 1 56.5 ns 16 16.5 ns 16 12.8 ns 16 12.8 ns 16
halley x 2 95.3 ns 47 27.4 ns 47 20.6 ns 47 20.6 ns 47
halley x 3 122.1 ns 51 37.9 ns 51 28.4 ns 51 28.3 ns 51
newton x 1 48.8 ns 10 14.4 ns 10 11.0 ns 10 11.0 ns 10
newton x 2 78.6 ns 20 23.0 ns 20 17.2 ns 20 15.4 ns 20
newton x 3 104.5 ns 40 32.0 ns 40 23.1 ns 40 21.3 ns 40
newton x 4 134.7 ns 51 40.8 ns 51 29.1 ns 51 25.8 ns 52
\sourceoff
Die 10.1 ns sind im Vergleich zu 9.7 ns (Beitrag 32 ) DESHALB minimal größer,
weil die Argumente erst zur Laufzeit (innerhalb der Schleife) berechnet werden (und nicht schon in einem Array zuvor abgelegt wurden).
AVX2 ist hier auch nur wieder der Compilerschalter und noch nicht die "echte Block-Berechnung", die erst so richtig 8 parallele Berechnungen optimiert.
Grüße Gerd
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
 | Beitrag No.43, vom Themenstarter, eingetragen 2021-10-26
|
Hallo zusammen,
@AlphaSigma:
\quoteon(2021-10-25 21:55 - AlphaSigma in Beitrag No. 40)
da Du gerne Python verwendest, ist evtl. diese Seite On computing roots of quartic and cubic equations in Python für Dich interessant.
\quoteoff
Ich habe mir die Seite angeschaut und auch den dortigen Code für die cardanische Lösung heruntergeladen und mit meiner verglichen. Meinen eigenen Code habe ich noch weiter optimiert und noch einmal einiges herausgeholt. Das ist mein aktueller Stand:
\sourceon Python
\numberson
sqrt3_j=complex(0,3**0.5)
rcpr3=1/3
def cbrt(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=1/(3*a) # (Von hier an echte kubische Gleichung)
r=-b*s
q=r*r
p=c*s-q
q=r*q-1.5*s*(d+c*r)
if p==0:
r+=(2*q)**rcpr3
return r,r,r
s=((q*q+p*p*p)**0.5+q)**rcpr3
q=p/s
s,q=s-q,-0.5*(s-q-(s+q)*sqrt3_j)
return r+s,r+q,r-s-q
\sourceoff
Mein Code ist um den Faktor 1,6 schneller als der Code von der genannten Seite. Und das, obwohl ich auch komplexe Koeffizienten verarbeiten kann, während der dortige Code nur reelle Koeffizienten verarbeitet und noch nicht einmal die bösen Trivialfälle rausfiltert. Perfektes Beispiel für ineffizienten Code... 😁
@HyperG:
Ich finde Dein Wissen über die Abläufe in den Prozessoren sehr beeindruckend und tiefgreifend. Meine Aussage ("warum eigentlich die trigonometrischen Funktionen und die Potenzfunktion SOO langsam sind"), hast Du aber ein wenig falsch verstanden, oder ich habe sie unklar formuliert.
Wenn man die cardanische Lösung umsetzt, muss man an einem Punkt eine dritte Wurzel aus einem Wert ziehen, also die Gleichung $x^3=a$ lösen. Der Halley-Algorithmus in der Hochsprache tut nicht nur das, sondern er berechnet gleich die Nullstelle eines (komplizierteren) kubischen Polynoms $ax^3+bx^2+cx+d=0$, und nicht nur eine dritte Wurzel, und das auch noch schneller, als der Prozessor in der cardanischen Lösung eine dritte Wurzel zieht. Das macht doch keinen Sinn - auf den ersten Blick.
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...
Trotzdem habe ich noch keinen Vergleich gesehen (oder ich habe ihn übersehen), wo die cardanische Lösung in C++ mal ohne Newton und Sortierung am Ende und das Halleyverfahren mit einem Abbruchkriterium von $10^{-15}$ statt $10^{-6}$ gegeneinander antreten. 🙂
Ciao,
Thomas
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.44, eingetragen 2021-10-26
|
Hallo MontyPythagoras,
in den "Numerical Recipes" Section 5.6 "Quadratic and Cubic Equations" sind auch Formeln für das Lösen von Gleichungen 3. oder 4. Grades mit komplexen Koeffizienten angegeben. Die habe ich mehr oder weniger Eins zu Eins in C implementiert. Soweit ich das nach wenigen Tests mit beliebigen komplexen Koeffizienten sagen kann, ist das schneller als die direkte Übersetzung des Python Codes aus Beitrag No. 3 in C. Die Zeiten in ms sind für 1 Mio. Iterationen. Die Zeiten für 1 Iteration ergeben sich einfach durch ersetzen von ms in ns. In beiden Programmen verwende ich cpow(z, 1.0/3.0) bei komplexen Koeffizienten.
\sourceon
$ ./nr_cubic_eqn
1000000
-2.0 3.0
5.0 7.0
-4.0 1.0
1.0 8.0
-1.5578708012718054e-01 -3.6769974519211840e-01 i 4.4408920985006262e-16 0.0000000000000000e+00 i
-1.8961719049962567e-01 -1.1011823952492525e+00 i -8.8817841970012523e-16 4.4408920985006262e-16 i
2.8386580908834469e-01 9.6118983274906322e-01 i 3.1086244689504383e-15 -8.8817841970012523e-16 i
CPU time: 247 ms
$ ./mp_cbrt_cmplx
1000000
-2.0 3.0
5.0 7.0
-4.0 1.0
1.0 8.0
-1.5578708012718065e-01 -3.6769974519211840e-01 i 0.0000000000000000e+00 -8.8817841970012523e-16 i
2.8386580908834469e-01 9.6118983274906311e-01 i 0.0000000000000000e+00 8.8817841970012523e-16 i
-1.8961719049962561e-01 -1.1011823952492525e+00 i -8.8817841970012523e-16 -4.4408920985006262e-16 i
CPU time: 357 ms
$ ./nr_cubic_eqn
1000000
0.7 6.1
-3.3 -5.5
13.6 6.2
-0.01 2.22
-9.6234441549459704e-02 6.8690348224395059e-01 i -6.8833827526759706e-15 -5.3290705182007514e-15 i
4.4145318691477686e-01 -4.0792569009939683e-01 i -9.7699626167013776e-15 -7.9936057773011271e-15 i
-3.1103602691555174e+00 5.8596039264310766e+00 i -6.6613381477509392e-14 -1.6875389974302379e-14 i
CPU time: 284 ms
$ ./mp_cbrt_cmplx
1000000
0.7 6.1
-3.3 -5.5
13.6 6.2
-0.01 2.22
-9.6234441549459815e-02 6.8690348224395015e-01 i -8.8817841970012523e-16 -8.8817841970012523e-16 i
4.4145318691477708e-01 -4.0792569009939639e-01 i -1.7763568394002505e-15 -3.5527136788005009e-15 i
-3.1103602691555166e+00 5.8596039264310766e+00 i -8.8817841970012523e-16 -5.3290705182007514e-14 i
CPU time: 319 ms
\sourceoff
\quoteon(2021-10-26 20:49 - MontyPythagoras in Beitrag No. 43)
Hallo zusammen,
...
Mein Code ist um den Faktor 1,6 schneller als der Code von der genannten Seite. Und das, obwohl ich auch komplexe Koeffizienten verarbeiten kann, während der dortige Code nur reelle Koeffizienten verarbeitet und noch nicht einmal die bösen Trivialfälle rausfiltert. Perfektes Beispiel für ineffizienten Code... 😁
Ciao,
Thomas
\quoteoff
Um wieviel ist der neue Python Code denn schneller als der in Beitrag No. 3? Lohnt es sich den auch noch nach C zu transformieren?
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.45, eingetragen 2021-10-27
|
\quoteon(2021-10-26 20:49 - MontyPythagoras in Beitrag No. 43)
...
\sourceon Python
\numberson
...
s=((q*q+p*p*p)**0.5+q)**rcpr3
q=p/s
s,q=s-q,-0.5*(s-q-(s+q)*sqrt3_j)
return r+s,r+q,r-s-q
\sourceoff
\quoteoff
Nur weil die 4 Zeilen in Python so schön einfach und explizit aussehen, sind sie es nicht!
- ...**0.5 -> muss komplex und somit per trigonometrischen Funktionen und pow berechnet werden. Natürlich wird Python das gut optimiert haben, aber unter 50 ns wird man "krumme Zahlen" eines 2D-Arrays nicht auf 15 Stellen genau berechnen können
- ...**rcpr3 -> analog das Gleiche wie ...**0.5
- Fehlerfortpflanzung für große krumme Werte fehlt mir hier total! Ich vermute es werden nicht mal 12 richtige Stellen sein, denn es müssen 2 atan, 2 sin, 2 cos , 2 pow(z,1/6) berechnet werden!
Die Zeile 19: s,q=s-q,-0.5*(s-q-(s+q)*sqrt3_j)
verstehe ich nicht: Variable s (hängt in der Luft??) , q=s-q (einzige logische Zuweisung), Formel (hängt in der Luft?)
\quoteon(2021-10-26 20:49 - MontyPythagoras in Beitrag No. 43)
Mein Code ist um den Faktor 1,6 schneller als der Code von der genannten
Seite.
\quoteoff
"Code von der genannten Seite" ist zu ungenau. Wir haben sooo viele Codes vorgestellt. Den "alten uneffektiven Cardano"-Code der Original-Seite, den wir alle bereits als uneffektiv angesehen haben -> der ist doch komplett raus, da
data4 133 ns
data7 133 ns
Am besten mit ns Angabe, sonst sehe ich da nicht mehr durch.
\quoteon(2021-10-26 20:49 - MontyPythagoras in Beitrag No. 43)
...
Trotzdem habe ich noch keinen Vergleich gesehen (oder ich habe ihn übersehen), wo die cardanische Lösung in C++ mal ohne Newton und Sortierung am Ende und das Halleyverfahren mit einem Abbruchkriterium von $10^{-15}$ statt $10^{-6}$ gegeneinander antreten. 🙂
...
\quoteoff
Dass ich die Sortierung am Ende immer weggelassen habe, hatte ich bereits genannt. Macht aber nicht mal 1 ns aus & kann vernachlässigt werden.
Die eine Newton-Iteration macht auch nicht mal 1 ns aus und ist vermutlich bei ALLEN Algorithmen immer nötig, um bei großen krummen Zahlenwerten (die Beispiele gehen immer von kleinen gerundeten ganzen Zahlen aus, was "echte double Berechnungen" aber nicht würdig ist) immer die 15 Stellen zu schaffen!
Und wie gesagt brauchen Cardano und Halley/Deiters nur 8 Stellen genau sein, da die letzte Iteration die Genauigkeit schnell verdoppelt.
Eigentlich ist der Optimierungsvorgang fließend:
- vom exakten umständlichen, langsamen Weg über die vielen trig. & Pow (jedoch mit enormer Fehlerfortpflanzung!)
- über zig Optimierungen für Näherungsformeln (auch Python gehört dazu)
- bis hin zur schnellsten Halley/Deiters ohne trig. & pow
Sollte Halley/Deiters jedoch nur 6 Stellen schaffen, bringt die letzte Iteration nur 12 richtige Stellen -> was Schummeln darstellt. denn 8 müssen es schon sein! (vermutlich weitere Iteration nötig; geht aber schnell)
Da hier die "komplexe Lösung" immer so hoch gelobt wird: ich dachte, dass Halley/Deiters nie dafür gedacht war, um komplex zu rechnen.
Man kann doch den Vorteil "bei nicht komplexen Lösungen kommt er ohne komplexe Zwischenergebnisse aus" nicht als Nachteil ansehen...
Eher anders herum: Cardanos Algor. benötigt:
- entweder zig Fallunterscheidungen (alte Algo vom echten Cardano )
- oder beim neuen Cardano muss IMMER komplex (also 2D-Array) gerechnet werden
Ich werde mal "schöne krumme Parameter" heraussuchen, damit wir hier nicht immer "geschönte Geschwindigkeiten" für kleine glatte Werte bekommen!
Gerade die Näherungsfunktionen können damit besser erkannt werden:
im Wertebereich -10...10 schaffen alle locker 8...14 Stellen. Danach wird es interessant! Selbst der genauste ASM-Befehl FSIN(x) bringt bei 10^12
nur noch 8 richtige Stellen, obwohl 19 versprochen wurden!
Ich weiß z.B. nicht, welchen Befehl der Linux gcc Kompiler für sincos benutzt: a) den langsamen genauen FSINCOS des Coprozessors,
oder b) den schnellen ungenauen "der neuen AVX-Welt".
(zeitlich machen das über 40% aus! -> bei mir war der Unterschied nur das 4fache, statt Faktor 10 -> deshalb vermute ich ja, dass der alte Cardano mit FSINCOS absichtlich schlecht gemacht wurde)
|
Profil
|
zippy
Senior
Dabei seit: 24.10.2018
Mitteilungen: 5019
 | Beitrag No.46, eingetragen 2021-10-27
|
\quoteon(2021-10-27 00:08 - hyperG in Beitrag No. 45)
Ich weiß z.B. nicht, welchen Befehl der Linux gcc Kompiler für sincos benutzt: a) den langsamen genauen FSINCOS des Coprozessors,
oder b) den schnellen ungenauen "der neuen AVX-Welt".
\quoteoff
Gibt es in "der neuen AVX-Welt" einen Befehl, der irgendetwas trigonometrisches berechnet?
Die glibc hat bis Version 2.27 für x86 (nicht für x86_64) FSINCOS benutzt, mit Version 2.28 ist das rausgeflogen.
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.47, vom Themenstarter, eingetragen 2021-10-27
|
Hallo zusammen,
@HyperG
\quoteon(2021-10-27 00:08 - hyperG in Beitrag No. 45)
Nur weil die 4 Zeilen in Python so schön einfach und explizit aussehen, sind sie es nicht!
\quoteoff
Also manchmal glaube ich schon, Du hältst mich für ein bisschen dämlich. Glaubst Du, das weiß ich nicht? Meine erste Einschätzung aus dem Threadstart beruhte im wesentlichen auf der Annahme, dass der Algorithmus von Halley-Deiters komplexe Koeffizienten kann und der Vergleich Sinn macht. Macht er aber nicht. Das wissen wir schon seit etlichen Beiträgen, und das habe ich selbst im Detail geschildert. Wir könnten nur
1. komplexe Cardano-Lösung mit komplexem Halley-Algorithmus
2. reelle Cardano-Lösung mit reellem Halley-Deiters-Algorithmus
vergleichen, denn sonst vergleichen wir Äpfel mit Birnen. Natürlich ist es kein Nachteil eines für reelle Koeffizienten gedachten Algorithmus, wenn er nur reell arbeitet, ist doch klar. Und wie ich schon in Beitrag #23 dargelegt habe, besteht eine komplexe Multiplikation in Wirklichkeit aus 4 Multiplikationen und 2 Additionen, eine Division sogar aus 6 Multiplikationen, 3 Additionen und 2 Divisionen. Und natürlich muss der Prozessor beim Potenzieren erst einmal Radius und Argument (Winkel) berechnen, für letzteres trigonometrische Funktionen anwenden (wie sonst), und das anschließend wieder zurück rechnen. Weiß ich doch alles.
\quoteon(2021-10-27 00:08 - hyperG in Beitrag No. 45)
Und wie gesagt brauchen Cardano und Halley/Deiters nur 8 Stellen genau sein, da die letzte Iteration die Genauigkeit schnell verdoppelt.
\quoteoff
Das ist eine weit verbreitete, gravierende Fehleinschätzung. Quadratische Konvergenz (Newton) oder kubische Konvergenz (Halley) sind keine Garantie! Ich werde es mal mit einem Beispiel erläutern. Nachfolgend die Näherung für die Gleichung $x^2-1=0$, Startwert $x_0=2$, Lösung ist $x=1$, dargestellt sind Newton und Halley:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Konvx2-1.pngSehr zügige Konvergenz, wie erwartet verdoppelt sich pro Schritt die Anzahl der korrekten Stellen bei Newton und verdreifacht sich etwa bei Halley. Das entspricht jeweils der erwarteten Konvergenzgeschwindigkeit.
Und hier nun die Lösung der Gleichung $x^2-2x+1=0$, Startwert und Lösung so wie oben:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Konvx2-2x_1.pngUn nu? Nix is mit Verdoppelung oder Verdreifachung der korrekten Stellen pro Iterationsschritt. Das Kriterium $10^{-8}$ erreichen Newton nach 28 und Halley nach 17 Schritten, 16 Stellen Genauigkeit aber erst nach jeweils gut doppelt so viel Iterationsschritten.
Ich könnte sogar ganz bösartig sein und die Lösung der Gleichung $x^2-2x+1.00000000001$ ausrechnen lassen. Wenn Du da zu früh abbrichst, gibt dir der Algorithmus eine Lösung, die gar nicht existiert!
Ich hoffe, damit ist klar geworden, warum ich sage, dass man nicht bei einem Iterationsschritt von $x_n$ zu $x_{n+1}$ von relativ $10^{-6}$ abbrechen darf, wenn man sicherstellen will oder muss, dass die Lösung auf 15, 16 Stellen genau gefunden wurde. Du müsstest innerhalb des Algorithmus aufwendige Checks anstellen, ob nicht die Ableitung und ggf. zweite Ableitung eventuell nahe null sind, um solche Fälle wie oben auszuschließen, und solche Checks fressen auch wieder Zeit. Und wenn die Macher des Halley-Deiters-Algorithmus bei einer relativen Veränderung um $10^{-6}$ abbrechen, dann brechen sie mindestens eine Schleife zu früh ab! Darauf versuche ich schon seit etlichen Beiträgen hinzuweisen.
\quoteon(2021-10-27 00:08 - hyperG in Beitrag No. 45)
"Code von der genannten Seite" ist zu ungenau. Wir haben sooo viele Codes vorgestellt.
\quoteoff
Das bezog sich auf die Seite, die AlphaSigma in seinem Beitrag #40 verlinkt hatte, also auf Python-Code, mit dem ich meinen Algorithmus verglichen hatte. Eine komplexe Lösung lobe ich nur, wenn sie auch für komplexe Koeffizienten gedacht ist. Sonst würde sie keinen Sinn machen.
\quoteon(2021-10-27 00:08 - hyperG in Beitrag No. 45)
Die Zeile 19: s,q=s-q,-0.5*(s-q-(s+q)*sqrt3_j)
verstehe ich nicht: Variable s (hängt in der Luft??) , q=s-q (einzige logische Zuweisung), Formel (hängt in der Luft?)
\quoteoff
Das ist Python-typisch, da hängt nichts in der Luft. Man kann in Python in einer Zeile zwei Werten zwei neue Werte zuordnen: Hier wird der Variablen $s$ der Wert $s-q$ und der Variablen $q$ der Wert nach dem Komma zugeordnet. Wollte ich das nacheinander in zwei Zeilen machen, müsste ich eine Hilfsvariable hinzunehmen, weil sonst die erste Anweisung die zweite Anweisung verändern würde.
@AlphaSigma:
\quoteon(2021-10-26 21:58 - AlphaSigma in Beitrag No. 44)
Um wieviel ist der neue Python Code denn schneller als der in Beitrag No. 3? Lohnt es sich den auch noch nach C zu transformieren?
\quoteoff
Ja, sollte es. Ich bin um den Faktor 1,85 schneller geworden. Wesentlichste Änderung ist, dass ich nach der Berechnung der ersten Lösung die zweite nicht durch Lösung der verbleibenden quadratischen Gleichung, die ich durch Polynomdivision erhalte (wie es ironischerweise auch der Halley-Deiters-Algorithmus tut), berechne, sondern über eine der komplexen Einheitswurzeln. Aber auch weiter oben habe ich ein wenig optimiert und die Anzahl der arithmetischen Operationen reduziert.
Ciao,
Thomas
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.48, eingetragen 2021-10-27
|
\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.
Der Vollständigkeit halber nochmal die Info zur CPU:
\sourceon
model name : Intel(R) Core(TM) i3-3240 CPU @ 3.40GHz
stepping : 9
microcode : 0x21
cpu MHz : 1767.566
cache size : 3072 KB
\sourceoff
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.49, eingetragen 2021-10-27
|
\quoteon(2021-10-27 21:54 - AlphaSigma in Beitrag No. 48)
...
Voraussichtlich mache ich noch ähnliche Versuche für complex.
...
\quoteoff
\sourceon C
/* cpow(z, 1.0/3) */
const double rd3 = 1.0/3.0;
sumc += cpow(rc[i], rd3);
/* im Vergleich zu */
/* atan2, sqrt, cbrt */
at3 = atan2(cimag(rc[i]), creal(rc[i])) / 3.0;
sqr6 = sqrt(cbrt(creal(rc[i])*creal(rc[i]) + cimag(rc[i])*cimag(rc[i])));
sumc += CMPLX(sqr6*cos(at3), sqr6*sin(at3));
\sourceoff
Bei komplexer Arithmetik komme ich mit Standardmitteln nur von ca. 140 ns auf ca. 100 ns.
\sourceon
$ gcc -Wall -pedantic -lm -o rand_cmplx rand_cmplx.c
$ ./rand_cmplx
Method Total CPU time Single CPU time
cpow(z, 1.0/3) 14619 ms 146 ns 857460903.237076 229812969.864184 i
atan2, sqrt, cbrt 10551 ms 105 ns 857460903.237076 229812969.864184 i
$ gcc -Wall -pedantic -O2 -lm -o rand_cmplx rand_cmplx.c
$ ./rand_cmplx
Method Total CPU time Single CPU time
cpow(z, 1.0/3) 14121 ms 141 ns 857460903.237076 229812969.864184 i
atan2, sqrt, cbrt 9180 ms 91 ns 857460903.237076 229812969.864184 i
$ gcc -Wall -pedantic -Ofast -lm -o rand_cmplx rand_cmplx.c
$ ./rand_cmplx
Method Total CPU time Single CPU time
cpow(z, 1.0/3) 14147 ms 141 ns 857460916.794279 229812973.483622 i
atan2, sqrt, cbrt 10130 ms 101 ns 857460916.794279 229812973.483622 i
\sourceoff
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.50, eingetragen 2021-10-28
|
\quoteon(2021-10-26 20:49 - MontyPythagoras in Beitrag No. 43)
...
Mein Code ist um den Faktor 1,6 schneller als der Code von der genannten Seite. Und das, obwohl ich auch komplexe Koeffizienten verarbeiten kann, während der dortige Code nur reelle Koeffizienten verarbeitet und noch nicht einmal die bösen Trivialfälle rausfiltert. Perfektes Beispiel für ineffizienten Code... 😁
...
\quoteoff
\quoteon(2021-10-27 17:51 - MontyPythagoras in Beitrag No. 47)
...
@AlphaSigma:
\quoteon(2021-10-26 21:58 - AlphaSigma in Beitrag No. 44)
Um wieviel ist der neue Python Code denn schneller als der in Beitrag No. 3? Lohnt es sich den auch noch nach C zu transformieren?
\quoteoff
Ja, sollte es. Ich bin um den Faktor 1,85 schneller geworden. Wesentlichste Änderung ist, dass ich nach der Berechnung der ersten Lösung die zweite nicht durch Lösung der verbleibenden quadratischen Gleichung, die ich durch Polynomdivision erhalte (wie es ironischerweise auch der Halley-Deiters-Algorithmus tut), berechne, sondern über eine der komplexen Einheitswurzeln. Aber auch weiter oben habe ich ein wenig optimiert und die Anzahl der arithmetischen Operationen reduziert.
Ciao,
Thomas
\quoteoff
Ok, das bedeutet dann aber, dass der Code aus Beitrag No. 3 in etwa so schnell ist, wie der unter On computing roots of quartic and cubic equations in Python .
Ich habe den Python Code aus Beitrag No. 43 nun auch in C programmiert. Dabei habe ich allerdings für die komplexe Kubikwurzel cpow(z, 1/3) verwendet, was langsamer ist, als die separate Berechnung von atan, sqrt und cbrt.
Kompiliert wurde mit "gcc -Ofast"
Getestet habe ich bisher nur einige wenige Fälle mit beliebigen komplexen Koeffizienten.
mp_cbrt_cmplx.c Implementierung von Python Code aus Beitrag No. 3 in C.
mp_cbrt_cmplx_3.c Implementierung von Python Code aus Beitrag No. 43 in C.
nr_cubic_eqn_2.c Umsetzung der Formeln in "Numerical Recipes" mit atan2, sqrt, cbrt anstatt cpower(z, 1/3) in C.
\sourceon
Programm Total CPU time Single CPU time
mp_cbrt_cmplx.c 244098 ms 244 ns
mp_cbrt_cmplx_3.c 209701 ms 210 ns
nr_cubic_eqn_2.c 174244 ms 174 ns
------------------------------------------------------------------
mp_cbrt_cmplx.c 408291 ms 408 ns
mp_cbrt_cmplx_3.c 377753 ms 378 ns
nr_cubic_eqn_2.c 179496 ms 179 ns
===================================================================
$ ./mp_cbrt_cmplx
1000000000
0.33 -0.11
11.1 5.55
-13.1 4.7
-7.7 0.0022
6.0427687361243565e-01 3.5778249108357352e-01 i -4.6074255521943996e-15 2.2759572004815709e-15 i
-2.0404401459785504e-02 1.9179580570848453e-02 i -2.7339241981394480e-15 5.1347814888913490e-15 i
-2.2853454315871393e+00 2.3294140360392157e-01 i 1.7659484985443896e-15 1.2767564783189300e-15 i
CPU time: 244098 ms
$ ./mp_cbrt_cmplx_3
1000000000
0.33 -0.11
11.1 5.55
-13.1 4.7
-7.7 0.0022
6.0427687361243565e-01 3.5778249108357352e-01 i -4.6074255521943996e-15 2.2759572004815709e-15 i
-2.0404401459785837e-02 1.9179580570848315e-02 i -5.8703042427055152e-15 1.9290125052862095e-15 i
-2.2853454315871393e+00 2.3294140360392168e-01 i 1.0451708942760263e-15 -5.8147930914742574e-15 i
CPU time: 209701 ms
$ ./nr_cubic_eqn_2
1000000000
0.33 -0.11
11.1 5.55
-13.1 4.7
-7.7 0.0022
-2.2853454315871389e+00 2.3294140360392152e-01 i -2.2649417064091182e-14 2.2065682614424986e-15 i
-2.0404401459785393e-02 1.9179580570848093e-02 i 2.9143354396410359e-16 1.5126788710517758e-15 i
6.0427687361243543e-01 3.5778249108357363e-01 i -1.0547118733938987e-15 3.5388358909926865e-15 i
CPU time: 174244 ms
-------------------------------------------------------
$ ./mp_cbrt_cmplx
1000000000
-101.34 0.055
0.33 -22.22
-57.57 0.0
5.9 -100.100
7.7819316926153903e-01 7.2807736920504373e-01 i 2.8421709430404007e-14 1.4495349365262200e-14 i
-1.2584096142931656e-02 -9.0492639189306567e-01 i 4.2632564145606011e-14 -3.3688329903469594e-14 i
-7.3182796117428628e-01 7.4998280330269251e-01 i 5.6843418860808015e-14 -4.2348069495545815e-14 i
CPU time: 408291 ms
$ ./mp_cbrt_cmplx_3
1000000000
-101.34 0.055
0.33 -22.22
-57.57 0.0
5.9 -100.100
7.7819316926153903e-01 7.2807736920504373e-01 i 2.8421709430404007e-14 1.4495349365262200e-14 i
-7.3182796117428628e-01 7.4998280330269229e-01 i 2.1316282072803006e-14 2.1600776722863202e-14 i
-1.2584096142931656e-02 -9.0492639189306556e-01 i -1.4210854715202004e-14 -3.3910374508394625e-14 i
CPU time: 377753 ms
$ ./nr_cubic_eqn_2
1000000000
-101.34 0.055
0.33 -22.22
-57.57 0.0
5.9 -100.100
7.7819316926154458e-01 7.2807736920504762e-01 i 1.8616219676914625e-12 3.7687214460291329e-13 i
-1.2584096142930434e-02 -9.0492639189307122e-01 i 1.8616219676914625e-12 4.0063091732989164e-13 i
-7.3182796117429294e-01 7.4998280330269407e-01 i 1.8687273950490635e-12 3.6976671724531229e-13 i
CPU time: 179496 ms
\sourceoff
Zur Kontrolle geben die beiden rechten Zahlen die komplexen Funktionswerte an den berechneten Nullstellen an. Wie man in den Beiträgen No. 48 u. 49 sieht, liefert -Ofast etwas andere, vmtl. ungenauere, Werte als -O2 oder keine Optimierung. Wenn ich Zeit habe, werde ich evtl. noch einen Newton-Schritt anhängen. Etwas verwunderlich sind die deutlich unterschiedlichen Zeiten für die beiden Sätze von Koeffizienten bei den ersten beiden Codes. Cardano ist ja keine iterative Methode.
Jetzt haben wir einige ganz gut optimierte Implementierungen der Cardano-Methode für komplexe Gleichungen 3. Grades in C. Es fehlt noch eine C-Implementierung der Methode nach Halley für den komplexen Fall.
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.51, eingetragen 2021-10-28
|
@zippy:
Du fragst zu recht, denn wie schon gesagt, gibt es zwar unter
https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#=undefined&techs=SVML
viele Befehle, ABER keine echten Maschinenbefehle, sondern in Wirklichkeit sind das spezielle Befehle, die bei den neusten Intel-C-Compilern mit "SVML" & AVX fehlende Funktionen nachbilden.
Mit schnellen ungenauen "der neuen AVX-Welt" hätte ich besser allg. schreiben sollen: mit Compilern, die per Näherungsfunktion & schnellen AVX-Befehlen die fehlenden ASM-Befehle ersetzen. (das würde dann auch für andere Compiler gelten, dann aber ohne die vielen SVML.)
@MontyPythagoras:
Also ich betrachte zunächst nur die reellen Ergebnisse, da ich ja den Vergleich zu Deiters & Macias-Salinas als Haupt-Thema sehe.
zur Konvergenz:
Ja, es ist bekannt, dass eine Iteration schlecht gegen ein Grenzwert konvergiert, dessen (Teil-)Term genau an dieser Stelle eine Polstelle {abgeleitete Funktion bildet an dieser Stelle ein 1/0} hat.
Deshalb hatte ich ja immer wieder auf große krumme Parameter & Ergebnisse gepocht.
@All:
Um den Vergleich von Deiters & Macias-Salinas & Cardano gerechter zu machen, habe ich 2 Fakten optimiert:
1. Genauigkeit bei krummen Parametern: Data9: 233328643357133 -14029858000 -628199 0.091
Hier schaffte Deiters nur 11 richtige Dezimalstellen bei 1 von 3 Ergebnissen! Ursache: Im Code fehlte für 1 Ergebnis die letzte Iteration!!
\sourceon c++
if (p < 0.0) {
x[i] = aux_newton1(o, p - d);
x[j] = aux_newton1(o, p + d);// hier stand x[j] = p + d; !!!!!!
\sourceoff
2. Der alte Cardano-Code, der zwar statt atan ein acos verwendet, hat bereits ein x^3, was mit dem darauffolgenden x^(1/6) zum schnellen sqrt(x) optimiert werden kann. Dann fehlte ja noch der Beweis, dass durch Fehlerfortzpflanzung nicht alle 15 Dezimalstellen richtig sind. -> also mal mit & mal ohne die letzte Iteration zur Genauigkeitsoptimierung.
Dass man swap (Vertauschung ) weglassen muss, hatten wir ja schon geklärt.
\sourceon c++ mit Genauigkeitsangabe der Dezimalstellen
Cardano 15 richtige mit Newton |Cardano 14 richtige ohne Newton | Deiters & Macias-Salinas 15 richtige | Deiters 11 richtige
data9 59.39 ns | 44.95 ns | 40.53 ns | 31.66 ns
\sourceoff
Wie man sieht, gleichen sich die Zeiten immer mehr an.
Durch weitere Optimierungen (jede trigonom. Funktion durch AVX-Näherung)
könnte man die Zeiten immer weiter annähern...
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.52, eingetragen 2021-10-29
|
\quoteon(2021-10-28 21:31 - AlphaSigma in Beitrag No. 50)
...
Ich habe den Python Code aus Beitrag No. 43 nun auch in C programmiert. Dabei habe ich allerdings für die komplexe Kubikwurzel cpow(z, 1/3) verwendet, was langsamer ist, als die separate Berechnung von atan, sqrt und cbrt.
Kompiliert wurde mit "gcc -Ofast"
Getestet habe ich bisher nur einige wenige Fälle mit beliebigen komplexen Koeffizienten.
...
\quoteoff
Um die Zeiten etwas belastbarer zu machen, habe ich jetzt 100000 Sätze zufälliger komplexer Polynomkoeffizienten erzeugt und für jeden Satz die Nullstellen-Routine 100-Mal aufgerufen.
mp_cbrt_cmplx.c Implementierung von Python Code aus Beitrag No. 3 in C.
mp_cbrt_cmplx_3.c Implementierung von Python Code aus Beitrag No. 43 in C.
nr_cubic_eqn_2.c Umsetzung der Formeln in "Numerical Recipes" mit atan2, sqrt, cbrt anstatt cpower(z, 1/3) in C.
\sourceon
Programm Single iteration CPU time
mp_cbrt_cmplx.c 336 ns
mp_cbrt_cmplx_3.c 282 ns
nr_cubic_eqn_2.c 206 ns
\sourceoff
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.53, vom Themenstarter, eingetragen 2021-10-30
|
Hallo zusammen,
es wird mittlerweile schwierig, Euch beiden gleichzeitig zu folgen. 🙂
@HyperG:
\quoteon(2021-10-28 23:44 - hyperG in Beitrag No. 51)
Also ich betrachte zunächst nur die reellen Ergebnisse, da ich ja den Vergleich zu Deiters & Macias-Salinas als Haupt-Thema sehe.
zur Konvergenz:
Ja, es ist bekannt, dass eine Iteration schlecht gegen ein Grenzwert konvergiert, dessen (Teil-)Term genau an dieser Stelle eine Polstelle {abgeleitete Funktion bildet an dieser Stelle ein 1/0} hat.
\quoteoff
Zum ersten Punkt: das Hauptthema war für mich eigentlich ein genereller Schnelligkeitsvergleich, unter der irrtümlichen Annahme, dass beide >>komplexe<< Koeffizienten verarbeiten. In dem Wikipedia-Artikel ist keine Rede davon, dass die Aussage, der Halley/Deiters-Algorithmus sei schneller als die cardanischen Formeln, nur für reelle Koeffizienten gelte. Diese Unterscheidung mussten wir zwangsläufig machen, erst nachdem sich herausstellte, dass der Halley/Deiters-Algorithmus (ich kürze den mal mit "HD" ab) nur mit reellen Koeffizienten arbeitet - was fraglos in der physikalischen Welt häufiger der Fall ist als mit komplexen.
Zum zweiten Punkt der Konvergenz: wenn das bekannt ist, und ich glaube Dir durchaus, dass Du das wusstest, warum musste ich mir dann den Mund fusselig labern, dass ein Abbruchkriterium von $10^{-6}$ oder $10^{-8}$ nicht reicht? So, wie der HD programmiert wurde, ist er bei $10^{-6}$ lange noch nicht fertig.
\quoteon(2021-10-28 23:44 - hyperG in Beitrag No. 51)
\sourceon c++ mit Genauigkeitsangabe der Dezimalstellen
Cardano 15 richtige mit Newton |Cardano 14 richtige ohne Newton | Deiters & Macias-Salinas 15 richtige | Deiters 11 richtige
data9 59.39 ns | 44.95 ns | 40.53 ns | 31.66 ns
\sourceoff
\quoteoffIch betrachte hier eigentlich den Vergleich "Cardano 14 richtige ohne Newton" mit "Deiters & Macias-Salinas 15 richtige" als am relevantesten. Es wird hier doch der HD mit Cardano verglichen. Ein Nachschieben einer Newton-Iteration bei Cardano zwecks Erhöhung der Präzision verfälscht das Ergebnis, während man einen Näherungsalgorithmus schon so programmieren muss, dass die Genauigkeit zuverlässig vergleichbar ist und der Algo eben nicht zu früh abbricht. Und dann sind wir hier mittlerweile auf eine Differenz runter von nur noch gut 10%.
Vielen Dank für Deine Mühe.
@AlphaSigma:
Was genau ist "nr_cubic_eqn_2"?
\quoteon(2021-10-26 21:58 - AlphaSigma in Beitrag No. 44)
in den "Numerical Recipes" Section 5.6 "Quadratic and Cubic Equations" sind auch Formeln für das Lösen von Gleichungen 3. oder 4. Grades mit komplexen Koeffizienten angegeben. Die habe ich mehr oder weniger Eins zu Eins in C implementiert.
\quoteoffWas ist mit den "Numericsl Recipes" gemeint? Bezieht sich das auf die in Deinem Beitrag #36 verlinkte Seite? Dort finde ich das Stichwort oder ein Kapitel "5.6" nicht.
\quoteon(2021-10-28 21:31 - AlphaSigma in Beitrag No. 50)
Ok, das bedeutet dann aber, dass der Code aus Beitrag No. 3 in etwa so schnell ist, wie der unter On computing roots of quartic and cubic equations in Python .
\quoteoff
Jein. Zumindest in Python war auch schon der Code aus #3 rund 15% schneller als der von der besagten Seite, und zwar selbst unter Verwendung von komplexen Zahlen, während der verlinkte Code mit reellen Koeffizienten arbeitet.
Ich denke, mit #43 habe ich zumindest in Python das Maximum aus der cardanischen Lösung mit komplexen Koeffizienten herausgekratzt. Es sollte auch in C etwa das Optimum sein, da ich viel Gehirnschmalz hineingesteckt habe, die Anzahl der arithmetischen Operationen so weit wie möglich zu reduzieren. Es hängt an den internen trigonometrischen und Potenzfunktionen.
Interessant ist, dass der Unterschied zwischen #3 und #43 in C++ deutlich geringer ausfällt als in Python, was wohl auf die hohe Effizienz der C++-Compiler zurückzuführen ist.
Das bringt mich zu einem neuen Gedanken. Du fandest zwar die von Dir verlinkte Datei "apple_tr_kt32_cuberoot.pdf" heute nicht mehr relevant, aber warum eigentlich? Der Autor geht einen ganz pragmatischen, aber auch durchaus vernünftigen Weg: er bestimmt die Anzahl der arithmetischen Operationen. Natürlich darf man dabei die zeitfressenden trig. und Potenzfunkton nicht außer acht lassen. Wenn man die Einzelzeiten jeder Operation bestimmen kann (was nicht so einfach ist und was HyperG ansatzweise in Beitrag #41 dargestellt hatte), sollte man in der Lage sein, die Geschwindigkeit abzuschätzen. Cardano sollte in C++ immer die ungefähr gleiche Zeit brauchen, weil die Anzahl der Rechenoperationen immer gleich ist (von Trivialfällen abgesehen). Bei HD könnte man nur die Anzahl der Rechenoperationen für die Initialisierung und anschließend pro Schleife zählen. Meint Ihr, ob das ein lohnenswerter Ansatz ist?
Ciao,
Thomas
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.54, eingetragen 2021-10-30
|
\quoteon(2021-10-30 09:30 - MontyPythagoras in Beitrag No. 53)
...
@AlphaSigma:
Was genau ist "nr_cubic_eqn_2"?
\quoteon(2021-10-26 21:58 - AlphaSigma in Beitrag No. 44)
in den "Numerical Recipes" Section 5.6 "Quadratic and Cubic Equations" sind auch Formeln für das Lösen von Gleichungen 3. oder 4. Grades mit komplexen Koeffizienten angegeben. Die habe ich mehr oder weniger Eins zu Eins in C implementiert.
\quoteoffWas ist mit den "Numericsl Recipes" gemeint? Bezieht sich das auf die in Deinem Beitrag #36 verlinkte Seite? Dort finde ich das Stichwort oder ein Kapitel "5.6" nicht.
Ciao,
Thomas
\quoteoff
Ich dachte die Numerical Recipes wären bekannt. Siehe auch Wikipedia Numerical Recipes. Konkret beziehe ich mich auf das Buch Numerical Recipes in C, Second Edition (1992).
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.55, eingetragen 2021-10-30
|
\quoteon(2021-10-30 09:30 - MontyPythagoras in Beitrag No. 53)
...
unter der irrtümlichen Annahme, dass beide >>komplexe<< Koeffizienten verarbeiten.
\quoteoff
Deshalb hatte ich immer von den "alten".. gesprochen:
Die "alten" cardanischen Formeln ...erstmals 1545 von dem Mathematiker Gerolamo Cardano in seinem Buch Ars magna...
Die ganze Optimierung, Reduzierung von Fällen und komplexen Zahlen kam erst später.
\quoteon(2021-10-30 09:30 - MontyPythagoras in Beitrag No. 53)
Und dann sind wir hier mittlerweile auf eine Differenz runter von nur noch gut 10%.
Vielen Dank für Deine Mühe.
\quoteoff
Mich haben die interessanten SVML-Funktionen des Intel-Compilers nicht losgelassen!
So besorgte ich mir eine DLL (nur 18 MB; nicht die kompletten 1 GB!)
und baute sie in mein VC mit ein. Da der Befehl acos vorkommt, programmierte ich ihn mit "normalen acos" nach, um diesen 512-Bit-Befehl (8 double) objektiv zu vergleichen:
\sourceon c++
__m512d _mm512_acos_pd_Alt(__m512d aIn512)
{
__m512d ret; for (int i = 0; i < 8; i++) ret.m512d_f64[i] = acos(aIn512.m512d_f64[i]);
return ret;
}
\sourceoff
Es gab in der DLL nun mehrere versteckte Funktionen, von denen 2 sehr interessant waren:
\sourceon 8er double-Blöcke (natürlich über 300000000, um die ns genau zu bekommen)
Befehl | Zeit pro 1 Befehl | Genauigkeit in Dezimalstellen
normaler acos | 7.153 ns | 15
mm512_acos_pd_ep_z0| 1.897 ns !! | 9..11 richtige
mm512_acos_pd_z0 | 2.185 ns ! | 15
\sourceoff
Bei gleicher Genauigkeit bei selber CPU & selben Compiler alles ohne Multitasking also um 4.968 ns schneller!
(dass man nur 3.3 statt 7.5 mal schneller wurde, zeigt:
- dass es kein echter Maschinenbefehl ist
- dass Intel seine eigene CPU gut optimieren kann)
Wenn man "Blockberechnung" (und das macht man bei vielen Wetterdaten so), zulässt, wird allein schon mit Optimierung eines trigon. Befehls der "alte Cardano" schneller als die Näherungsiteration.
(natürlich hat die "HD-Iteration" ja auch noch ihre AVX-Reserven...)
Und dann gibt es ja noch 2 weitere trigon. Befehle...
Bei sincos "zusammen in 1 Befehl" scheitert jedoch (noch) der Rückgabewert der Intel-DLL mit der Microsoft-Welt...
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.56, eingetragen 2021-10-30
|
Hallo,
nach der Geschwindigkeit habe ich mir jetzt die Genauigkeit angeschaut. Für die Methode nach "Numerical Recipes" habe ich zufällige Polynomkoeffizienten mit Real- und Imaginärteil zwischen -1000 und +1000 erzeugt und dann den Betrag der Funktionswerte an den berechneten Nullstellen ausgegeben. Die x-Achsen der Histogramme sind log10-skaliert. Ein Newton-Schritt bringt eine deutliche Verbesserung. Weitere Newton-Schritte bringen relativ wenig.
https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/35344_logerr_no_newton.png
https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/35344_logerr_1_newton.png
https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/35344_logerr_2_newton.png
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.57, eingetragen 2021-10-31
|
\quoteon(2021-10-14 08:10 - MontyPythagoras in Beitrag No. 9)
...
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
\quoteon(2021-10-16 23:46 - MontyPythagoras in Beitrag No. 23)
...
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.
...
\quoteoff
\quoteon(2021-10-28 23:44 - hyperG in Beitrag No. 51)
@MontyPythagoras:
Also ich betrachte zunächst nur die reellen Ergebnisse, da ich ja den Vergleich zu Deiters & Macias-Salinas als Haupt-Thema sehe.
...
\quoteoff
\quoteon(2021-10-30 09:30 - MontyPythagoras in Beitrag No. 53)
Hallo zusammen,
es wird mittlerweile schwierig, Euch beiden gleichzeitig zu folgen. 🙂
...
\quoteoff
Das liegt wohl daran, dass ich vorrangig an dem Fall mit komplexen Koeffizienten interessiert bin und keine CPU-spezifischen Optimierungen berücksichtige und HyperG, zumindestens erst einmal, am rellen Fall und dabei auch spezifische Optimierungen für Intel-Prozessooren untersucht.
Ein C-Programm mit der Methode nach Cardano für komplexe Koeffizienten habe ich mittlerweile. Z.B. die Umsetzung der NR Formeln mit anschließendem Newton-Schritt. Für einen Vergleich mit Halley auf dem gleichen PC und mit gleichem C-Compiler und Einstellungen fehlt noch eine C-Implementierung von Halley für komplexe Koeffizienten. Laut Wikipedia Halley-Verfahren ist "das Halley-Verfahren eine Methode der numerischen Mathematik zur Bestimmung von Nullstellen f(x) = 0 reeller Funktionen $f\colon \mathbb {R} \to \mathbb {R}$." Es gibt aber z.B. einen Artikel von Lily Yau and Adi Ben-Israel The Newton and Halley Methods for Complex Roots, in dem die wesentlichen Formeln für den komplexen Fall beschrieben werden. Daraus ein C-Programm zu machen, ist aber ein gewisser Aufwand.
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.58, eingetragen 2021-10-31
|
Hallo,
auf der bereits in Beitrag No. 36 verlinkten Seite Polynominal Root Finders habe ich diesen Artikel Practical_Implementation_of_Polynomial_Root_Finders.pdf entdeckt.
Darin sind C++-Programme für mehrere Verfahren und relle oder komplexe Polynome (beliebigen Grades) enthalten. Am Programm für Halley musste ich einige kleine Anpassungen machen, um es unter Linux zu kompilieren und einen kleinen Fehler zu beheben. Die Halley-Methode ist für Polynome mit beliebigem Grad programmiert, d.h. nicht auf kubische Polynome spezialisiert bzw. optimiert. Weiterhin werden in der Hauptroutine mehrere Unterprogramme aufgerufen. Damit will ich sagen, dass evtl. bzgl. der Geschwindigkeit noch Optimierungspotenzial besteht. Ich habe es bisher nur für wenige Koeffizienten getestet. Ich erhalte die identischen Nullstellen, wie unter dem Online Polynomial Root Finder. Die Rechenzeiten hängen stark von der gewählten Compiler-Optimierung ab. Gesamtzeiten für 1000000 Iterationen.
\sourceon
Optimierung Total time Single Time
keine 3580 ms 3580 ns
-O2 1342 ms 1342 ns
-Ofast 650 ms 650 ns
========================================================
$ g++ -Wall -pedantic halley_cmplx.c -o halley_cmplx
$ ./halley_cmplx
(0.707, -33.123)
(-77.77, 0.005)
(1.11, 6.432)
(0.007, -19.91)
(-1.2694948853467996,-1.1666575448955212)
(1.5628141270307716,1.5484697287350955)
(0.0297148592118700,-0.4376766358087487)
Total CPU time: 3580 ms
Single CPU time: 3580 ns
$ g++ -Wall -O2 -pedantic halley_cmplx.c -o halley_cmplx
$ ./halley_cmplx
(0.707, -33.123)
(-77.77, 0.005)
(1.11, 6.432)
(0.007, -19.91)
(-1.2694948853467996,-1.1666575448955212)
(1.5628141270307716,1.5484697287350955)
(0.0297148592118700,-0.4376766358087487)
Total CPU time: 1342 ms
Single CPU time: 1342 ns
$ g++ -Wall -Ofast -pedantic halley_cmplx.c -o halley_cmplx
$ ./halley_cmplx
(0.707, -33.123)
(-77.77, 0.005)
(1.11, 6.432)
(0.007, -19.91)
(-1.2694948853467998,-1.1666575448955214)
(1.5628141270307714,1.5484697287350955)
(0.0297148592118700,-0.4376766358087488)
Total CPU time: 650 ms
Single CPU time: 650 ns
\sourceoff
Das ist selbst mit -Ofast erst einmal (etwas) langsamer als Cardano, aber wie geschrieben: das Programm ist nicht auf kubische Polynome spezialisiert und vermutlich nicht voll auf Geschwindigkeit optimiert.
Prozessor wie in Beitrag No. 52 und 48.
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.59, eingetragen 2021-11-01
|
Hallo,
analog zu den Beiträgen No. 52 und 56 für Cardano habe ich nun auch die Zeiten und Fehler für Halley nach dem C++-Programm von Henrik Vestermark Polynominal Root Finders für zufällig erzeugte Polynomkoeffiziente bestimmt.
\sourceon
Optimierung Total time Single Time
keine 4472 ms 4472 ns
-O2 1762 ms 1762 ns
-Ofast 725 ms 725 ns
\sourceoff
Real- und Imaginärteil der zufälligen Polynomkoeffizienten liegen wieder zwischen -1000 und +1000. Ich habe wieder den Betrag der Funktionswerte an den berechneten Nullstellen ausgegeben. Die x-Achsen der Histogramme sind log10-skaliert.
https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/35344_h_logerr_no_newton.png
https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/35344_h_logerr_1_newton.png
https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/35344_h_logerr_2_newton.png
Um direkt den Fehler der Nullstellen zu bestimmen, müsste man die Routinen zusätzlich mit höherer Genauigkeit implementieren oder nur Polynome mit "glatten" Nullstellen wählen.
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.60, vom Themenstarter, eingetragen 2021-11-01
|
Hallo AlphaSigma,
\quoteon(2021-10-31 11:29 - AlphaSigma in Beitrag No. 57)
Laut Wikipedia Halley-Verfahren ist "das Halley-Verfahren eine Methode der numerischen Mathematik zur Bestimmung von Nullstellen f(x) = 0 reeller Funktionen $f\colon \mathbb {R} \to \mathbb {R}$."
\quoteoff
Warum sollte das so sein. Der danach von Dir zitierte Artikel macht viel Lärm um nichts. Schau Dir die Herleitung zum hyperbolischen Verfahren an, und zwar hier. Warum sollte da eine Einschränkung auf $\mathbb {R}$ gelten? Wir suchen die Nullstelle $x$ von $ax^3+bx^2+cx+d=0$. Wir haben ein ungefähres $x_0$, es gilt $x=x_0+\varepsilon$ (alles $\in\mathbb C$). Dann gilt unzweifelhaft:
$$a(x_0+\varepsilon)^3+b(x_0+\varepsilon)^2+c(x_0+\varepsilon)+d=0$$Multipliziert man das aus und vernachlässigt $\varepsilon^3$:
$$ax_0^3+bx_0^2+cx_0+d+(3ax_0^2+2bx_0+c)\varepsilon+(6ax_0+2b)\varepsilon^2=0$$Alles ganz simpel algebraisch. Vernachlässigt man auch noch $\varepsilon^2$, erhält man das Newtonverfahren, löst man dagegen die quadratische Gleichung, erhält man das (parabolische) Halleyverfahren, und mit dem Trick
$$\left((3ax_0^2+2bx_0+c)+(6ax_0+2b)\varepsilon\right)\left((3ax_0^2+2bx_0+c)-(6ax_0+2b)\varepsilon\right)\approx\left((3ax_0^2+2bx_0+c)\right)^2$$wegen $\varepsilon^2\approx0$ auch das hyperbolische. Alles "very basic", rein algebraisch. Die Einschränkungen, die zum Tragen kommen können ($\varepsilon$ nicht vernachlässigbar gegen eine Ableitung, weil die zufällig null ist), kann man sich an fünf Fingern abzählen, ohne in die Tiefen der komplexen Analysis abtauchen zu müssen.
Schwierig wird es beim Deiters-Algorithmus allenfalls dort, wo mithilfe der Laguerre-Samuelson-Ungleichung Schranken der Lösungen gefunden werden müssen. Komplexe Zahlen lassen sich so schlecht per "<" und ">" vergleichen... 😉
Was die Newton-Iteration nach der cardanischen Lösung angeht: ich glaube, da betreibst Du Selbsttäuschung. Was Du dort siehst, ist im Grunde die Verteilung Deiner Koeffizienten, nicht die der Lösungen. Wenn Dein Zufallsgenerator zufällig "große" Koeffizienten ausspuckt, ist auch der Funktionswert größer, der aufgrund von Rundungsfehlern übrig bleibt. Den Funktionswert für die Abschätzung heranzuziehen ist nicht zielführend.
Du musst meines Erachtens die genaue Lösung kennen, und die Näherungsformel oder Cardano dann damit vergleichen. Ich würde also eher folgendes vorschlagen, um die Genauigkeit zu untersuchen:
1. Lege per Zufallsgenerator 3 komplexe x-Werte fest.
2. Streue gelegentlich Spezialfälle ein, wie zum Beispiel doppelte oder dreifache Nullstellen (z.B. 100:10:1 oder 10000:100:1).
3. Lege per Zufallsgenerator noch ein $a$ fest, und berechne dann daraus erst die Koeffizienten gemäß $ax^3+bx^2+cx+d=a(x-x_1)(x-x_2)(x-x_3)$. Mir ist klar, dass bei der Berechnung der Koeffizienten nun auch schon Rundungsfehler passieren.
4. Vergleiche die per Algorithmus gefundene Lösung mit der vorher konstruierten Lösung, wobei ich nur den komplexen Abstand $\left|x_{1,algo}-x_1\right|$ im Verhältnis zu $\left|x_1\right|$ betrachten würde.
Nur so bekommen wir einen halbwegs sinnvollen Vergleich. Du musst bedenken, dass der Rechner zwar Rundungsfehler macht, aber an der gleichen Stelle immer die gleichen. Nur weil der Rechner glaubt, dass der Funktionswert jetzt besser oder näher an null ist, muss das noch lange nicht stimmen.
Ich habe selbst auch mit einer Newton-Schleife nach der Erstberechnung experimentiert. Manchmal bekomme ich eine Verbesserung, aber manchmal auch eine Verschlechterung. Folgendes Zahlenbeispiel:
Rechne mal mit folgenden Koeffizienten:
\sourceon Python
A,B,C,D=1+1j,-1.8-4.2j,0.12+4.92j,0.616-1.656j
\sourceoff
Die einzige, exakte Lösung ist $1+0\mathord,4i$. Wenn ich den Funktionswert berechne, erhalte ich
\sourceon
f(1+0.4j)= 2.220446049250313e-16
\sourceoff
Aber z.B. bei den Werten
\sourceon
x=1.0000000000032483+0.4j
x=0.9999976741835636+0.40000402839089727j
x=1.000000000000924+0.39999999999961733j
\sourceoff
ist laut Python der Funktionswert exakt 0,0. Was soll da Newton helfen? Wenn zum Beispiel der genaue Funktionswert eigentlich negativ wäre, er aber aufgrund von Rundungsfehlern positiv ist, stimmt ja schon das Vorzeichen des Deltas, das zu $x$ hinzuaddiert wird, nicht, und es kommt zu einer Verschlechterung statt einer Verbesserung.
Ciao,
Thomas
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.61, vom Themenstarter, eingetragen 2021-11-02
|
... 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
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.62, vom Themenstarter, eingetragen 2021-11-02
|
Hallo zusammen,
mit der beschriebenen Methode habe ich die Genauigkeitsverteilung der Cardano-Methode ohne Newton-Schleife berechnet. Da ich eine logarithmische Darstellung benutze, habe ich, wenn die Lösung exakt getroffen wurde, den Wert auf -18 gesetzt. "n=20.000" bedeutet, dass ich 20.000 kubische Gleichungen gelöst habe, was also 60.000 Lösungswerte ergab, die ich hier klassifiziert habe.
Nachfolgend das Bild bei drei unterschiedlichen Lösungen:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Cardano3.png
Viel besser geht es nicht. Bei zwei Lösungen, also einer doppelten Nullstelle, ergibt sich folgendes Bild:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Cardano2.png
Noch schlechter wird es bei dreifachen Nullstellen:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Cardano1.png
Im Normalfall ist die Genauigkeit bei einer FLOAT-Zahl ungefähr $10^{-16}$. Da früher oder später die dritte Wurzel gezogen wird, führt das bei einer Dreifach-Nullstelle dazu, wenn eigentlich die dritte Lösung aus null gezogen werden sollte, die dritte Wurzel aus einem Wert in der Nähe von $10^{-16}$ gezogen wird, was zu einer relativen Abweichung von rund $10^{-5\mathord,3}$ führt. Und das in der überwiegenden Anzahl der Fälle.
Der Fall "dreifache Nullstelle" ist allerdings auch der Fall, bei dem Halley und Newton versagen werden, weil sich dort auch die Konvergenz dramatisch verlangsamen wird.
Jedenfalls macht es in meinen Augen keinen Sinn, bei Cardano eine Newton-Schleife hinten dranzubasteln, wenn ich mir das erste Bild ansehe.
Ciao,
Thomas
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.63, eingetragen 2021-11-02
|
Für Interessierte:
unter hier im Forum
habe ich einen tollen schnellen Befehl gefunden:
\sourceon c++
__m256 _mm256_csqrt_ps (__m256 a)
\sourceoff
der die komplexe Wurzel nicht nur "mit einem Schlag",
sondern auch noch in 4-facher Ausführung berechnet!
Grüße
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.64, vom Themenstarter, eingetragen 2021-11-04
|
Hallo zusammen,
ich habe nun mit der in #61 beschriebenen Methode einen Vergleich der Genauigkeit zwischen der Cardano-Lösung und dem Halley-Algorithmus (jeweils komplexe Koeffizienten und Lösungen) angestellt. Datenbasis sind jeweils 100.000 kubische Gleichungen.
Zunächst für drei verschiedene Nullstellen:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Halley_Cardano_three_null.pngTrotz eines Abbruchkriteriums von nur $10^{-6}$ erreicht der Halley-Algorithmus sogar eine geringfügig bessere Genauigkeit als die Cardano-Lösung. Allerdings rechtfertigt ein Unterschied von rund 0,15 Nachkommastellen nicht, dass man der Cardano-Lösung noch eine Newton-Scheife nachschiebt, zumal eine Verbesserung der Genauigkeit ohnehin nicht zu erwarten ist. Mindestens 14 Stellen sind praktisch immer korrekt.
Anders sieht es aus bei doppelten Nullstellen:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Halley_Cardano_dplnull.pngDie "Standard"-Cardano-Lösung schafft im Normalfall eine Genauigkeit von rund 15 bis 16 Stellen bei der einfachen Nullstelle und rund 8 Stellen bei der Doppel-Nullstelle. Halley kommt nur auf rund 13 bei der Einfach- und 6 bis 7 Stellen bei der Doppel-Nullstelle.
Ich habe darüber hinaus noch einen "optimierten" Cardano-Algorithmus entwickelt (in diesem Bild orange dargestellt), der immer auf 14 bis 16 Stellen Genauigkeit kommt. Der Trick ist simpel: wenn ich in der Cardano-Lösung feststelle, dass zwei Nullstellen sehr dicht beieinander liegen (das ist nur eine sehr einfache if-Abfrage), unterstelle ich, dass in Wirklichkeit eine "absichtliche" Doppelnullstelle vorliegt und treffe eine solche Annahme, indem ich einen Wert, der sehr nahe null ist, gleich null setze. Der Funktionswert ändert sich dadurch tatsächlich nicht, weil die Steigung der Funktion an dieser Stelle nun einmal null ist. Im Halley-Algorithmus ist eine solche heuristische Entscheidung nicht ohne größeren Aufwand möglich.
Kommen wir zur dreifachen Nullstelle:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Halley_Cardano_trplnull.pngHier liegen interessanterweise die einfache Cardano-Lösung und Halley dicht beieinander, die optimierte Cardano-Lösung schlägt die anderen beiden jedoch um Längen. Sie erreicht praktisch immer 16 exakte Stellen, die anderen beiden liegen bei nur etwas mehr als 5.
Bei allen obigen Untersuchungen habe ich den Halley-Algorithmus bei einem Kriterium von $10^{-6}$ abgebrochen, was ich für verfrüht hielt, siehe meine früheren Posts. Nachfolgend der Vergleich zwischen $10^{-6}$ und $10^{-9}$ bei drei unterschiedlichen Nullstellen:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Halley_criterion.pngPraktisch deckungsgleich. Das war hier zu erwarten, da bei drei unterschiedlichen Nullstellen die Ableitungen in der Nähe der Nullstelle nicht null sind und dadurch saubere kubische Konvergenz vorliegt.
Aber eine interessante Erkenntnis ist, dass ich hier keinen Vergleich unterschiedlicher Abbruchkriterien für Doppel- und Dreifach-Nullstellen zeigen kann. Bei einem Abbruchkriterium von $10^{-9}$ konvergiert Halley nämlich in beiden Fällen gar nicht. Das ist ohne weiteres damit erklärbar, dass die Genauigkeit in diesen Fällen ja, wie oben gezeigt, nur rund 8 bzw. 5 Stellen beträgt. $10^{-9}$ relative Distanz von einem Iterationswert zum nächsten wird daher in aller Regel nie erreicht bzw. unterschritten. Selbst beim Kriterium von $10^{-6}$ ist es quasi ein Wunder, dass der Algorithmus überhaupt konvergiert. Allerdings muss ich dazusagen, dass es mir hier nicht um Geschwindigkeit ging, sondern nur die Genauigkeit. Daher habe ich hier "faul" programmiert und immer bis zum Abbruchkriterium durchgezogen, während in dem Algorithmus von Deiters eine explizite Abfrage eingebaut wurde ob eine Dreifach-Nullstelle vorliegt (nicht jedoch für eine Doppelnullstelle).
Wer es in Python nachvollziehen möchte, ich habe den Cardano-Algorithmus noch ein wenig weiter optimiert.
Hier ist der "einfache" (und schnellere):
\sourceon Python
\numberson
rcpr3=1/3
eps1=complex(-0.5,0.75**0.5)
eps2=complex(-0.5,- 0.75**0.5)
def cbrt(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
p,q=c*s-q,r*q-1.5*s*(d+c*r)
s=(q*q+p*p*p)**0.5
s=q-s if q.real<0 else q+s
if s==0: return r,r,r
s**=rcpr3
q=s*eps1
t=s*eps2
return r+s-p/s,r+q-p/q,r+t-p/t
\sourceoff
Nachfolgend der optimierte. Man vergleiche die einfachen if-Abfragen (hier Zeilen 17 und 19), die den Unterschied zum obigen Algorithmus darstellen:
\sourceon Python
\numberson
rcpr3=1/3
eps1=complex(-0.5,0.75**0.5)
eps2=complex(-0.5,- 0.75**0.5)
def cbrt_opt(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
if abs(s)<2e-5*abs(q): s=0
s=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
return r+s-p/s,r+q-p/q,r+t-p/t
\sourceoff
Die if-Abfragen brauchen ein wenig Zeit, weil dafür halt Absolutwerte von komplexen Zahlen berechnet werden müssen. Dieser Algorithmus ist daher rund 10% langsamer als der obige.
Und der Vollständigkeit halber, der Halley:
\sourceon Python
\numberson
def cbrt_halley(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
x1=b/(-3*a)
dx=x1
while abs(dx)>1e-6*abs(x1):
f0=((a*x1+b)*x1+c)*x1+d
f1=(3*a*x1+2*b)*x1+c
f2=6*a*x1+2*b
p=(0.5*f0*f2-f1*f1)
if p==0: break
dx=f0*f1/p
x1+=dx
p=a*x1+b
q=p*x1+c
x2=(p+(p*p-4*a*q)**0.5)/(-2*a)
return x1,x2,-b/a-x1-x2
\sourceoff
Unter dem Strich muss man sich fragen, ob eine Optimierung hinsichtlich doppelter oder dreifacher Nullstellen lohnt. Ich würde sagen
- ja, wenn aus welchen guten Gründen auch immer mit dem Auftreten von mehrfachen Nullstellen zu rechnen ist (z.B. bei trigonometrischen Berechnungen)
- nein, wenn die Koeffizienten absolut zufällig verteilt sind. Wenn ich Doppel- und Dreifachnullstellen in die oben praktizierte zufällige Verteilung der Nullstellen einbeziehe, treten sie trotzdem nicht wirklich in Erscheinung, weil die Wahrscheinlichkeit, dass zwei oder sogar drei zufällige Nullstellen sehr nah beieinander liegen, extrem unwahrscheinlich ist.
Ciao,
Thomas
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.65, vom Themenstarter, eingetragen 2021-11-04
|
Hallo zusammen,
ich habe nun in Python auch einen Geschwindigkeitsvergleich der Algorithmen für komplexe Koeffizienten gemacht. Hier ist Cardano (einfach) um den Faktor 2,4 bis 3,5 schneller als Halley.
Allerdings ist der Startwert der Halley-Näherung hier nicht optimiert wie in der reellen Variante von Deiters.
Ich habe mich noch nicht damit beschäftigt, ob die Laguerre-Samuelson-Ungleichung in der Variante von Laguerre sinngemäß auch in $\mathbb C$ gilt, also anstatt in $\mathbb R$
$$-\frac{a_1}{n}-\sqrt{n-1}\,b\leq x_j \leq -\frac{a_1}{n}+\sqrt{n-1}\,b$$mit
$$b=\frac{1}{n}\sqrt{(n-1)a_1^2 - 2na_2}$$könnte es in $\mathbb C$ lauten:
$$\left|x_j+\frac{a_1}{n}\right|\leq \frac{1}{n}\sqrt{\left|(n-1)^2a_1^2 - 2n(n-1)a_2\right|}$$Aber selbst wenn das so sein sollte, dann ist ein großer Kreis um den Wendepunkt in der komplexen Zahlenebene jetzt auch nicht gerade eine gravierende Einschränkung des Suchbereichs. Eine schlechte Wahl des Startpunktes könnte immerhin ein lustiges Mandelbrotbildchen ergeben, wenn man es denn visualisiert. Ansonsten äußert sich eine schlechte Wahl des Startpunktes nur in einer nicht enden wollenden Dauerschleife. Wie man in meinem Halley-Algorithmus oben zumindest sehen kann, ist der Wendepunkt als Startpunkt zumindest immer erfolgreich - jedenfalls ist er unter 100.000 Versuchen nicht hängengeblieben. Im Durchschnitt braucht obiger Algorithmus knapp 4,8 Schleifen:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Halley_Schleifen.pngDas scheint ziemlich viel zu sein, und so ist es kein Wunder, dass er langsamer ist als Cardano. Im Reellen bei Halley-Deiters scheint er oft nach nur ein oder zwei Schleifen am Ziel zu sein.
Ciao,
Thomas
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.66, eingetragen 2021-11-04
|
\quoteon(2021-10-31 22:30 - AlphaSigma in Beitrag No. 58)
... diesen Artikel Practical_Implementation_of_Polynomial_Root_Finders.pdf entdeckt.
Darin sind C++-Programme für mehrere Verfahren und relle oder komplexe Polynome (beliebigen Grades) enthalten. Am Programm für Halley musste ich einige kleine Anpassungen machen, um es unter Linux zu kompilieren und einen kleinen Fehler zu beheben....
\quoteoff
Hallo @AlphaSigma,
da sind ja viele grobe Fehler drin.
Meine ersten schnellen Versuche mit
a) der fehlenden Funktion {gefunden auf anderen Internetseiten}
\sourceon c++
int complexdeflation( const register int n, complex a[], const complex z )
{
register int i;
double r, u;
r = -2.0 * z.real();
u = z.real() * z.real() + z.imag() * z.imag();
a[ 1 ] -= r * a[ 0 ];
for( i = 2; i < n - 1; i++ )
a[ i ] = a[ i ] - r * a[ i - 1 ] - u * a[ i - 2 ];
return n - 2;
}
\sourceoff
b) der inkompatiblen Parameter {denn fw = horner(n, a, wz, &fwz); passt zu keiner der 3 Funktionen in der PDF!}
\sourceon c++
//Original: double horner(const int n, const double a[], const double r, double *fz)
double horner(const int n, const complex a[], const complex r, complex *fz)
{
double fval;
fval = a[0].real();
for (int i = 1; i <= n; i++)
fval = fval * r.real() + a[i].real();
return fval;
}
\sourceoff
c) falscher Rückgabeparameter:
\sourceon c++
double horner( const int n, const complex a[], const complex z)
hatte ich passend zur Return-Variablen
complex fval; // geändert nach:
complex horner( const int n, const complex a[], const complex z)
\sourceoff
sind beim 1. Versuch natürlich falsch, denn Ergebnisse stimmen alle nicht:
\sourceon Ergebnisse
In : -0.437676635808749-19.910000000000000 i * x^3
In : 1.110000000000000+6.432000000000000 i * x^2
In : -77.769999999999996+0.005000000000000 i * x^1
In : 0.707000000000000-33.122999999999998 i * x^0
Erg: 0.000000000000000+0.000000000000000 i
Erg: 0.263604382472102-0.055864451969175 i
Erg: 0.000000000000000+0.000000000000000 i
\sourceoff
... und das richtige Ergebnis des einen komplexen Parameters ist
zum Input des ersten reellen Parameters gewandert!?
Bevor ich jetzt langwierige Fehlersuche mache, kannst Du vielleicht schnell aushelfen, da Du ja auch schon Anpassungen vornehmen musstest.
@MontyPythagoras:
Du schreibst bei allen Deinen Python Ergebnissen nie absolute Zeiten hin...
Konkret Beitrag 64: wie schnell ist es denn Deine Berechnungszeit bezüglich 1 komplexen Kubischen Berechnung (also in ns)
mit den Parametern aus Beitrag No. 58:
\sourceon nameDerSprache
In : 0.00700000000000-19.910000000000000 i * x^3
In : 1.110000000000000+6.432000000000000 i * x^2
In : -77.7700000000000+0.005000000000000 i * x^1
In : 0.707000000000000-33.12300000000000 i * x^0
\sourceoff
??
Grüße
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.67, vom Themenstarter, eingetragen 2021-11-04
|
Hallo HyperG,
\quoteon(2021-11-04 22:03 - hyperG in Beitrag No. 66)
Du schreibst bei allen Deinen Python Ergebnissen nie absolute Zeiten hin...
Konkret Beitrag 64: wie schnell ist es denn Deine Berechnungszeit bezüglich 1 komplexen Kubischen Berechnung (also in ns)
mit den Parametern aus Beitrag No. 58:
\sourceon nameDerSprache
In : 0.00700000000000-19.910000000000000 i * x^3
In : 1.110000000000000+6.432000000000000 i * x^2
In : -77.7700000000000+0.005000000000000 i * x^1
In : 0.707000000000000-33.12300000000000 i * x^0
\sourceoff
\quoteoff
Warum auch? Python ist bekanntermaßen langsam, Welten langsamer als C++. Mich interessieren erst einmal nicht die absoluten Zeiten. Aber wenn es Dich so sehr interessiert:
\sourceon
runfile('F:/Privates/Programmieren/Python/Cubicroot.ipy', wdir='F:/Privates/Programmieren/Python')
((0.007-19.91j), (1.11+6.432j), (-77.77+0.005j), (0.707-33.123j))
((0.029714859211870226-0.43767663580874894j), (1.5628141270307714+1.5484697287350953j), (-1.2694948853467996-1.1666575448955214j))
1.705 µs ± 60.622 ns per loop (mean ± std. dev. of 10 runs, 500000 loops each)
\sourceoff
Die erste Zeile sind die Koeffizienten, die zweite die Lösungen. 1.705ns auf einem i7-Prozessor 8. Generation in Spyder.
Ich habe mir derweil mal den von AlphaSigma in Beitrag #14 eingestellten C++ Code der von Deiters und Co. programmierten Cardano-Lösung angesehen. Meine Fresse, ist der schlecht programmiert! Den Aufruf der "POW"-Funktion, um eine sechste Wurzel zu berechnen, habe ich schon komplett wegrationalisiert. Eine dritte Wurzel ist auch überflüssig. Wiederholte Multiplikationen kann man sinnvoll zusammenfassen. Auf Newton kann man verzichten. Wenn "d==0" ist, wird trotzdem zweimal die dritte Wurzel gezogen. Ich werde morgen mal den Deiters-Code in optimierter Form hier reinstellen für die reelle Lösung. Heute Abend wird das nix mehr. 🙃
Ciao,
Thomas
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.68, eingetragen 2021-11-04
|
\quoteon(2021-11-04 22:30 - MontyPythagoras in Beitrag No. 67)
...
1.705 µs ± 60.622 ns ...
\quoteoff
Hi, na das ist doch mal was messbares! ("Welten langsamer als C++" ist für mich zu schwammig)
Ohne jegliche Optimierung habe ich das mal mit dem alten i7 und VC2012
1:1 (ohne Sonderfallbehandlung) übersetzt:
\sourceon c++
void cbrt_opt(const complex coeff[], complex res[])
{
complex rcpr3=complex(0.3333333333333333333,0.0);
complex s=rcpr3/coeff[0]; //echte kubische Gleichung)
complex r=-coeff[1]*s;
complex q=r*r;
complex t=q*r;
complex p=coeff[2]*s-q;
q=t-s*(coeff[3]+coeff[2]*r)*1.5;
s=sqrt((q*q+p*p*p));
if(q.real()<0.0) s=q-s;
else s+=q;
s=pow(s,0.33333333333333333333);//s**=rcpr3
q=s*complex(-0.5,0.8660254037844386467637);
t=s*complex(-0.5,-0.8660254037844386467637);
res[0]=r+s-p/s;res[1]=r+q-p/q;res[2]=r+t-p/t;
}
...
Erg: 0.029714859211870-0.437676635808748 i
Erg: 1.562814127030772+1.548469728735095 i
Erg: -1.269494885346799-1.166657544895522 i
in 360 ns +/- 10 ns {also etwa nur 4,7 mal }
\sourceoff
Der Vorteil: nur 1 if -> was sich dann gut zu AVX-Blöcken optimieren lässt.
Übrigens war ich von der __m256d _mm256_csqrt_pd (__m256d a)
nicht begeistert:
- nur 51 % schneller -> das hat Intel schlecht umgesetzt!
- keine AVX512 dabei
... wo doch die anderen Befehle 3...5 mal schneller waren.
Mal sehen, wann ich mal Zeit zum Optimieren habe, denn da steckt noch viel
Potential drin...
Grüße
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.69, vom Themenstarter, eingetragen 2021-11-05
|
... falls noch jemand Zweifel hatte, ob eine Newton-Schleife nach der Cardano-Lösung erhöhte Genauigkeit bringt:
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Cardano_Newton.pngDie durchschnittliche Genauigkeit erhöht sich von 15,69 auf 15,76 Nachkommastellen.
Ciao,
Thomas
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.70, eingetragen 2021-11-05
|
\quoteon(2021-11-04 22:03 - hyperG in Beitrag No. 66)
...
Hallo @AlphaSigma,
da sind ja viele grobe Fehler drin.
Meine ersten schnellen Versuche mit
a) der fehlenden Funktion {gefunden auf anderen Internetseiten}
...
b) der inkompatiblen Parameter {denn fw = horner(n, a, wz, &fwz); passt zu keiner der 3 Funktionen in der PDF!}
...
c) falscher Rückgabeparameter:
...
sind beim 1. Versuch natürlich falsch, denn Ergebnisse stimmen alle nicht:
... und das richtige Ergebnis des einen komplexen Parameters ist
zum Input des ersten reellen Parameters gewandert!?
Bevor ich jetzt langwierige Fehlersuche mache, kannst Du vielleicht schnell aushelfen, da Du ja auch schon Anpassungen vornehmen musstest.
Grüße
\quoteoff
Hallo hyperG,
den einen falschen Aufruf der Funktion horner in Halley(...) habe ich analog den anderen Aufrufen von horner so ersetzt:
\sourceon C
complex horner( const int n, const complex a[], const complex z);
// fw = horner(n, a, wz, &fwz);
fwz = horner(n, a, wz); fw = abs( fwz );
\sourceoff
Die im Artikel nicht existierende Funktion complexdeflation habe ich durch die Funktion compositedeflation für complex Polynomial and complex root ersetzt:
\sourceon C
int compositedeflation(const int n, complex a[], complex z );
// n = complexdeflation(n, a, z);
n = compositedeflation(n, a, z);
\sourceoff
Weiterhin habe ich
\sourceon C
delete [] a, a1, a2;
\sourceoff
durch
\sourceon C
delete [] a;
delete [] a1;
delete [] a2;
\sourceoff
ersetzt.
Speziell für die GNU-CC (unter Linux) habe ich die #define Variablen im Aufruf der pow Funktion
\sourceon C
_DBL_RADIX, -DBL_MANT_DIG
\sourceoff
durch
\sourceon C
__FLT_RADIX__, -__DBL_MANT_DIG__
\sourceoff
ersetzt.
Beim Aufrufen der Funktion Halley ist natürlich dieser Kommentar zu beachten:
\sourceon C
// Notice that a[0] is an, a[1] is an-1 and a[n]=a0
\sourceoff
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.71, eingetragen 2021-11-05
|
Danke AlphaSigma,
das hat mir einige Such-Zeit erspart.
Die entschiedene Stelle, warum mein erstes Ergebnis nicht kam, lag an:
\sourceon c++
res[0] = z; //res[n] = z;
n=compositedeflation(n, a, z);//n = complexdeflation(n, a, z);
\sourceoff
Und da n "von hinten" beginnt, kam nie bei res[0] was an...
Oder die quadratic(n, a, res);
am Ende bahandelt Index n anders...
\sourceon Ergebnisse & Zeiten
cbrt_opt:
In : 0.007000000000000-19.910000000000000 i * x^3
In : 1.110000000000000+ 6.432000000000000 i * x^2
In : -77.769999999999996+0.00500000000000 i * x^1
In : 0.707000000000000-33.122999999999998 i * x^0
Erg: 0.029714859211870- 0.437676635808748 i
Erg: 1.562814127030772+ 1.548469728735095 i
Erg: -1.269494885346799-1.166657544895522 i
Zeit: 387.50 ns (i7 3.5 GHz VC2012 Win7; exe wurde etwas größer -> daher 10 ns langsamer)
Halley:
In : 0.007000000000000 -19.910000000000000 i * x^3
In : 1.110000000000000 6.432000000000000 i * x^2
In : -77.769999999999996 0.005000000000000 i * x^1
In : 0.707000000000000 -33.122999999999998 i * x^0
Erg: 0.029714859211870 -0.437676635808749 i
Erg: -1.269494885346800 -1.166657544895521 i
Erg: 1.562814127030771 1.548469728735095 i
Zeit: 3218.75 ns (extrem langsam)
ENTER...
\sourceoff
|
Profil
|
AlphaSigma
Senior
Dabei seit: 23.11.2012
Mitteilungen: 453
| Beitrag No.72, eingetragen 2021-11-06
|
\quoteon(2021-11-05 23:54 - hyperG in Beitrag No. 71)
\sourceon Ergebnisse & Zeiten
...
Halley:
In : 0.007000000000000 -19.910000000000000 i * x^3
In : 1.110000000000000 6.432000000000000 i * x^2
In : -77.769999999999996 0.005000000000000 i * x^1
In : 0.707000000000000 -33.122999999999998 i * x^0
Erg: 0.029714859211870 -0.437676635808749 i
Erg: -1.269494885346800 -1.166657544895521 i
Erg: 1.562814127030771 1.548469728735095 i
Zeit: 3218.75 ns (extrem langsam)
ENTER...
\sourceoff
\quoteoff
Hallo hyperG,
ist das mit Optimierung? Beim g++ macht die Optimierung einen großen Unterschied (siehe Beitrag No. 58):
\sourceon
Optimierung Total time Single Time
keine 3580 ms 3580 ns
-O2 1342 ms 1342 ns
-Ofast 650 ms 650 ns
\sourceoff
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.73, eingetragen 2021-11-06
|
\quoteon(2021-11-06 00:09 - AlphaSigma in Beitrag No. 72)
...
ist das mit Optimierung? Beim g++ macht die Optimierung einen großen Unterschied (siehe Beitrag No. 58):
...
\quoteoff
Ja, alles auf maximale Geschwindigkeit & Gleitkomma fp:fast.
Zum Vergleich der i9 mit VC2017
\sourceon nameDerSprache
Algorithmus i9 VC2017 SSE2 i9 VC2017 AVX2
cbrt_opt: 228.75 ns 226.25 ns nur 25 Zeilen
Halley: 2140.00 ns 2187.50 ns etwa 300 Zeilen
\sourceoff
Einen alten g++ (mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev2\mingw64\bin\mingw32-make.exe von 2017, und Du vermutlich echtes Linux?) könnte ich ja mal probieren, aber da er kein AVX kann,
interessiert es mich weniger. (es kann auch sein, dass ich noch in den 300 Zeilen was optimieren könnte, aber fast Faktor 10 mit 25 Zeilen ist doch unschlagbar)
Der cbrt_opt lässt sich noch gut mit AVX blockweise optimieren...
(habe schon 4 Vergleiche allein für die komplexen Wurzeln laufen...)
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.74, vom Themenstarter, eingetragen 2021-11-06
|
Hallo zusammen,
ich hatte, wie in #67 bereits geschrieben, die Umsetzung von Deiters bzw. der Uni Köln als "schlecht programmiert" kritisiert. Wer kritisiert, muss auch besser machen. Daher habe ich die Cardano-Lösung für reelle Koeffizienten und Lösungen in C++ umgesetzt:
\sourceon C++
\numberson
void cbrt_real(const double coeff[], double res[])
{
double rcpr3 = 0.3333333333333333333;
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;
if ((p == 0.0) and (q == 0.0)) {
res[0] = r; res[1] = r; res[2] = r;
return;
}
s = q * q + p * p * p;
if (s > 0.0) {
(q < 0.0) ? s = - pow(sqrt(s) - q, rcpr3) : pow(sqrt(s) + q, rcpr3);
res[0] = r + s - p / s;
}
else {
double t = sqrt(-p);
sincos(rcpr3 * asin(q / (p * t)), &s, &p);
//Alternative zu voriger Zeile: 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
Der Code ist nicht debugged, Trivialfälle habe ich rausgelassen. Vielleicht möchte / kann es jemand ausprobieren. In nachfolgender Tabelle sieht man die Unterschiede:
(die Zahlen geben die Anzahl der arithmetischen Operationen bzw. Funktionsaufrufe an)
Deiters/Uni Köln
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Deiters_kub.png
Monty
https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_Monty_kub.png
Bei der kubischen Gleichung mit drei Lösungen habe ich gegenüber Deiters/Uni Köln 6 Multiplikationen, 3 Divisionen, 1 sqrt-Aufruf und 1 pow-Aufruf gespart, bei der mit nur einer Lösung 1 Multiplikation, 1 Division, 1 sqrt und 1 pow).
Ich benutze außerdem asin, Deiters acos. Das sollte keinen Unterschied machen. sincos wird jeweils einmal aufgerufen.
Und einen solchen Code: "0.5 * p * 2.0 * c" in Deiters' Cardano-Lösung muss man wohl nicht weiter kommentieren.
Die Newton-Schleife aus der Deiters-Programmierung habe ich wegen Sinnlosigkeit komplett gestrichen, dadurch sind noch weitere 15 Additionen, 21 Multiplikationen und 3 Divisionen entfallen.
Hier der Code in Python:
\sourceon Python
\numberson
import math as m
sqrt3 = m.sqrt(3)
rcpr3 = 1/3
def cbrt_real(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
if s < 0:
return ()
else:
s = m.sqrt(s)
return c + s, c - s
s = rcpr3 / a # (Von hier an echte kubische Gleichung)
r = -b * s
q = r * r
p, q = c * s - q, r * q - 1.5 * s * (d + c * r)
if p == q == 0: return r, r, r
s = q * q + p * p * p
if s > 0:
s = -(m.sqrt(s) - q) ** rcpr3 if q < 0 else (m.sqrt(s) + q) ** rcpr3
return (r + s - p / s, )
else:
t = m.sqrt(-p)
phi = rcpr3 * m.asin(q / (p * t))
s = m.sin(phi)
q = t * (sqrt3 * m.cos(phi) - s)
s *= 2 * t
return r + s, r + q, r - s - q
\sourceoff
Nachfolgend der Zeitunterschied zu meiner "komplexen" Variante:
\sourceon
Koeffizienten:
(1, 3, 6, 1)
Lösungen komplexer Algorithmus
(-0.18226832611317645, (-1.4088658369434117+1.8712332344895732j), (-1.4088658369434117-1.8712332344895732j))
Lösungen reeller Algorithmus
-0.18226832611317645
Zeit komplex
1.4071 µs ± 37.506 ns per loop (mean ± std. dev. of 10 runs, 500000 loops each)
Zeit reell
904.77 ns ± 6.3087 ns per loop (mean ± std. dev. of 10 runs, 500000 loops each)
\sourceoff
Ciao,
Thomas
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.75, eingetragen 2021-11-06
|
\quoteon(2021-11-06 13:58 - MontyPythagoras in Beitrag No. 74)
...
\sourceon
Koeffizienten:
(1, 3, 6, 1)
Lösungen komplexer Algorithmus
(-0.18226832611317645, (-1.4088658369434117+1.8712332344895732j), (-1.4088658369434117-1.8712332344895732j))
Lösungen reeller Algorithmus
-0.18226832611317645
Zeit komplex
1.4071 µs ± 37.506 ns per loop (mean ± std. dev. of 10 runs, 500000 loops each)
Zeit reell
904.77 ns ± 6.3087 ns per loop (mean ± std. dev. of 10 runs, 500000 loops each)
\sourceoff
...
\quoteoff
Hi, diesen Parametersatz darf man natürlich NICHT vergleichen, da hier der Spezialfall "nur 1 reelle Lösung", also nur 1 statt 3 Ergebnisse berechnet werden!
Meine Referenz im Beitrag 51 (Cardano 14 richtige ohne Newton mit 44.95 ns) war der Parametersatz:
\sourceon c++
A + B*x + C*x^2 + D*x^3=0; Anzahl (1...1000000000):800000
Wiederholungen=800000
Parameter:233328643357133 -14029858000 -628199 0.091
233328643357133.0000000000+-14029858000.0000000000*x+-628199.0000000000*x^2+0.0910000000*x^3=0
-33319.878552630354534 Err:1.25874127359790e+02
11111.494876193013624 Err:1.05657008805782e+02
6925494.097962150350213 Err:-9.44501767163936e+03
CPU time: 36.0 ms also 44.95 ns
MontyPythagoras Beitrag No.74
233328643357133.0000000000+-14029858000.0000000000*x+-628199.0000000000*x^2+0.0910000000*x^3=0
11111.494876183569431 Err:3.69674058567250e+02
6925494.097962150350213 Err:-9.44501767163936e+03
-33319.878552620299160 Err:4.08808726274875e+02
CPU time: 32.0 ms also 40.00 ns
\sourceoff
(Err ist hier mal nicht die Differenz zur Ideallösung, sondern bei eqn_cubic war es das Ergebnis nach dem Einsetzen in's Polynom, welches ja 0 sein soll)
Also 5 ns schneller (aber nicht immer volle 15.5 richtige Dezimalstellen).
Hinweis zu 3 Anpassungen:
- statt AND -> &&
- Parameter kompatibel zu anderen: const Polynomial& coeff, Vector& res
- Coeff Reihenfolge an vorhandenes Programm angepasst.
Im Komplexen (also 6 Ergebnisse) bin ich gerade bei 200 ns (i9 8er Block Optimierung noch ohne echte AVX; Beitrag 71 i7 war noch 387.50 ns)
. Später mehr...
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.76, vom Themenstarter, eingetragen 2021-11-06
|
Hallo HyperG,
\quoteon(2021-11-06 20:03 - hyperG in Beitrag No. 75)
Also 5 ns schneller (aber nicht immer volle 15.5 richtige Dezimalstellen).
\quoteoff
Das heißt, wir haben bei den reellen Koeffizienten und Lösungen den Halley-Deiters-Algorithmus eingeholt oder sogar knapp überholt. Sehe ich das richtig?
Was den Vergleich in #74 angeht: Vergleichen kann ich alles mögliche. Wenn mich nur die reellen Lösungen interessieren, z.B. in der Trigonometrie, warum sollte ich die konjugiert komplexen Lösungen berechnen? Aber prinzipiell hast Du natürlich recht, bei einer Gleichung mit drei Lösungen fällt der Unterschied geringer aus.
Noch eine wichtige Anmerkung zu meinem Code in #74:
Ich habe dort in Code-Zeile 22 eine Alternative zur Zeile 21 angegeben. Das hat folgenden Hintergrund:
Theoretisch sind wir hier bei $s\leq0$, was $q^2+p^3\leq0$ bedeutet. Daraus folgt eigentlich, dass $|q|\leq|p\sqrt{-p}|$ und somit $\left|\frac q{p\sqrt{-p}}\right|\leq1$ ist. Damit sollte es in Zeile 21 kein Problem sein, den $\arcsin$ zu berechnen. Allerdings habe ich Versuche gemacht mit Zufallszahlen für $q$, habe dann $p=-\sqrt[3]{q^2}$ berechnet, und dann $s$ und den Arkussinus berechnet. Man glaubt es kaum, aber Rundungsfehler sind teuflisch. Es kommt in 0,063% der Fälle vor, dass $s<0$ ist, und das Argument für den Arkussinus trotzdem größer als 1! Ich hatte mich gewundert, warum bei Halley-Deiters eine Beschneidung auf das Intervall $[-1; 1]$ stattfindet. (Vielleicht waren sie aber auch nur übervorsichtig). Der Code in Zeile 22 umgeht das Problem elegant. Die zeitliche Auswirkung kann ich jedoch nicht beurteilen. Gegenüber Zeile 21 wird in Zeile 22 eine Multiplikation und eine Division eingespart, aber ein weiterer Aufruf der sqrt-Funktion hinzugefügt. Auch weiß ich nicht, wie sich der "atan2" im Verhältnis zu "asin" verhält. In Python jedenfalls ist Zeile 22 langsamer als 21 - aber sicherer.
Welche Zeile hast Du verwendet in Deinem Vergleich?
Ciao,
Thomas
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 2117
| Beitrag No.77, eingetragen 2021-11-07
|
\quoteon(2021-11-06 21:24 - MontyPythagoras in Beitrag No. 76)
..bei den reellen Koeffizienten und Lösungen den Halley-Deiters-Algorithmus eingeholt oder sogar knapp überholt. Sehe ich das richtig?
...
Welche Zeile hast Du verwendet in Deinem Vergleich?
...
\quoteoff
Ich hatte die asin-Version genommen, weil mich die zusätzliche sqrt Funktion abgeschreckt hatte.
... geht gleich weiter -> muss kurz weg
... wieder da ->Fortsetzung:
Ja, wir haben durch Optimierung der "expliziten Funktionen" den Halley-Deiters-Algorithmus mindestens eingeholt.
Und ja, natürlich gibt es Sonderfälle, wo die Iterationsalgorithmen noch schlechter werden (wo die Ableitung/Anstieg sehr klein ist).
Gerechterweise muss man aber auch sagen, das wir fast nur auf der einen Seite optimiert haben. Da ich sehr hardware-nahe programmiere, kenne ich auch Tricks, wie man selbst einen Prozessor mal gut & mal schlecht dastehen lassen kann (ein alter AMD hatte eine einen Maschinenbefehl mit internem 128 Bit Integer-Zwischenergebnis, der jede neue Intel-CPU alt aussehen lassen kann...)
Wie ich bereits sagte, müsste man eigentlich mehr Wert auf Genauigkeit legen. Das gilt nicht nur für den 64-Bit-Datentyp double, wo oft schon 14 richtige Dezimalstellen als "OK" toleriert werden, sondern bei noch mehr Nachkommastellen (aber in die Richtung wolltest Du ja nicht).
Bei meiner Cardano-Version (45 ns) hatte ich damals das "fmax(-1.0, fmin(1.0, h))" übernommen, was bei uns hier eigentlich nicht nötig ist. Ohne diese Überprüfung kommt man auch nochmals um 5 ns runter auf 40 ns.
Hier nun einige interessante Überraschungen:
\sourceon c++ mit atan2 + sqrt
233328643357133.0000000000+-14029858000.0000000000*x+-628199.0000000000*x^2+0.0910000000*x^3=0
11111.494876193348318 Err:9.62975333681458e+01
6925494.097962151281536 Err:-5.37523331293836e+03
-33319.878552630543709 Err:1.20545579608691e+02
CPU time: 36.0 ms also 36.00 ns
\sourceoff
Das hätte ich nicht gedacht, dass atan2 + sqrt nochmals 4 ns schneller als die asin Funktion ist!
Manchmal liegt es auch am Compiler, dass die " / (p * t)" uneffektiv umgesetzt wurde.
Noch ein Beispiel, mit dem ich nicht gerechnet hatte:
die Optimierung der komplexen expliziten Version.
i7 VC2012 -> i9 VC2017: 387 ns -> 228.75 ns (Compiler SSE2 oder AVX kaum Unterschiede)
Dann wollte ich auf AVX umstrukturieren und bereitete 8er Blöcke vor:
jeder Befehl also 8 mal
for (i = 0; i < 8; i++) s[i] = rcpr3 / pKomplexPar4mal8[i].par[0];
Testweise habe ich das mal gemessen und kam auf 200 ns (also 28 ns schneller)!
Vermutlich konnte nun der Compiler die AVX2 besser anwenden & schon so einen Erfolg.
Dann Optimierungsversuch der komplexen Wurzel (complex):
for (i = 0; i < 8; i++) s[i] = sqrt(q[i] * q[i] + p[i] * p[i] *p[i]);
durch
die 51% schnellere _mm256_csqrt_pd komplexe Wurzel:
\sourceon c++
for (i = 0; i < 4; i++)
{
i20 = i * 2; i21 = i20 + 1;
ReIm256.m256d_f64[0]= (merk=q[i20] * q[i20] + p[i20] * p[i20] * p[i20]).real();
ReIm256.m256d_f64[1]= merk.imag();
ReIm256.m256d_f64[2] =(merk=q[i21] * q[i21] + p[i21] * p[i21] * p[i21]).real();
ReIm256.m256d_f64[3]= merk.imag();
Erg256=_mm256_csqrt_pd(ReIm256);//kann 2 komplexe Wurzeln
s[i20] = complex(Erg256.m256d_f64[0], Erg256.m256d_f64[1]);
s[i21] = complex(Erg256.m256d_f64[2], Erg256.m256d_f64[3]);
}
\sourceoff
Aber 243 ns (247 ns), also 43 ns langsamer?!
Das zeigt, dass die "Umschiftung" der vielen complex Variablen
nach AVX-Variablen (und zurück!) mehr Zeit verschlingt, als man denkt!
Es bringt also erst Zeitgewinn, wenn die komplette cbrt_opt auf AVX umgestellt wurde.
Da es die _csqrt_pd nur in der 256bit Version gibt und nur 51% (statt 200% der sin...) bringt, werde ich vermutlich bei der atan2 Umgehung landen.
Grüße Gerd
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.78, vom Themenstarter, eingetragen 2021-11-07
|
Hallo HyperG,
das ist in der Tat mal ein spannendes Ergebnis. Scheint, als hängen wir Halley langsam ab. Daher habe ich auch gleich noch eine weitere Optimierung:
\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 = - pow(sqrt(s) - q, rcpr3) : pow(sqrt(s) + q, rcpr3);
res[0] = r + s - p / s;
} else if (p == 0.0) {
res[0] = r; res[1] = r; res[2] = r;
} 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
Im C++ Code in Beitrag #74, Zeilen 10 bis 13, hatte ich frühzeitig die Abfrage drin, ob eventuell gleichzeitig p und q gleich null sind, also eine dreifache Nullstelle vorliegt. Die Abfrage muss man früher oder später machen, damit man nicht durch null teilt, aber jedes Mal sehr früh fragen, obwohl der Fall äußerst selten ist, war Zeitverschwendung. Und dann fiel mir auf, dass man diese Abfrage elegant in die nachfolgende if-Schleife integrieren und sich dadurch die Abfrage q==0 sparen kann. Also teste bitte noch einmal diese Variante.
Theoretisch ist die Zeile 12 sogar auch überflüssig, man könnte sich unabhängig vom Vorzeichen von q für eine der beiden Varianten entscheiden, beide würden funktionieren. Diese Zeile reduziert nur den Rundungsfehler, wenn die Diskriminante nahe null ist.
Um nach Ausreizen der Optimierungen einen Geschwindigkeitsvergleich mit Halley zu machen, müssten wir uns mal eine sinnvolle und vor allem realitätsnahe Test-Routine ausdenken. Wie mir scheint, haben Deiters und Co. 7 Datensätze erzeugt, und daran den Test durchgeführt. Der Faktor 10, den Halley angeblich teilweise schneller sein sollte als Cardano, war vermutlich nur genau dann der Fall, wenn der Wendepunkt, der als erstes bestimmt wird, zufällig auch genau eine Dreifach-Nullstelle ist. Aber was bringt ein Algorithmus, der nur in extremst seltenen Fällen so schnell ist, einfach weil er dann ja fast nichts rechnen muss, und in allen "normalen" Fällen langsamer ist?
Ich würde mir vorstellen, dass man per Zufallsgenerator ein paar Hunderttausend realitätsnahe Koeffizientensätze erzeugt, die entweder drei oder eine reelle Lösung enthalten: weit streuend, dicht beieinander usw.. Zusätzlich noch ein paar Hundert "Spezialfälle" wie Doppel- und Dreifachnullstellen. Und dann schaut man, welcher Algorithmus das Paket am schnellsten abarbeitet.
Ciao,
Thomas
|
Profil
|
MontyPythagoras
Senior
Dabei seit: 13.05.2014
Mitteilungen: 3417
Wohnort: Werne
| Beitrag No.79, vom Themenstarter, eingetragen 2021-11-08
|
... weiter geht das muntere Wegrationalisieren.
\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 = - pow(sqrt(s) - q, rcpr3) : pow(sqrt(s) + q, rcpr3);
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
Dank des sehr toleranten "atan2" kann man aus dem vorigen Beitrag die Zeilen 14 und 15 ersatzlos streichen, was ich hier gemacht habe. Laut Dokumentation gibt der atan2 in dem Fall, dass beide Argumente null sind, einfach null zurück. Bei Verwendung der asin-Funktion müsste man Division durch null vermeiden und dementsprechend zwingend p==0.0 abfragen, aber das müssen wir hier nicht. Auch bei p==0 stimmt diese Berechnung auf jeden Fall. Es dauert zwar ein bisschen länger, bis man die Lösung einer dreifachen Nullstelle gefunden hat, aber das würde ich opfern, da wie gesagt dieser Fall extrem selten sein wird und man eher Wert auf die schnelle Berechnung der allgemeinen, "zufälligen" Koeffizienten legen sollte.
Dieser Code ist also vorläufig das Optimum. Man könnte noch Genauigkeitsvergleiche anstellen, wieviel man verliert, wenn man Zeile 12 einfach kompromisslos durch
\sourceon C++
s = pow(sqrt(s) + q, rcpr3);
\sourceoff
ersetzt und sich eine weitere Abfrage spart. Die Ergebnisse würden dann größere Rundungsfehler aufweisen, wenn p klein und q negativ wäre.
Ciao,
Thomas
Ciao,
Thomas
|
Profil
|
|
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2023 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die
Nutzungsbedingungen,
die Distanzierung,
die Datenschutzerklärung und das Impressum.
[Seitenanfang]
|
Home / Seite ohne Frame
Wie unterstütze ich den Matheplaneten mit Geld?
