Bei der Berechnung von Arcus Sinus machen dessen Singularitäten bei +/-1
gerne Probleme: Die Ableitung von arcsin ist dort nicht endlich, d.h.
Verfahren zur Berechnung von arcsin reagieren u.U. sehr sensibel auf
Rundungsfehler.
Man kann daher folgende Funktion definieren, den "Reduzierten Arcus
Sinus" ra:
ra ist analytisch in 0 mit Konvergenzradius 2.
arcsin ist in einer Umgebung von 1 also i.W. die Quadratwurzel:
Damit lässt sich arcsin in [0,1] einfach berechnen, indem die
Quadratwurzel sowie ra implementiert werden. Quadratwurzel ist Old
School, und ein praktikabler Ansatz für ra ist ein kubischer Lookup:
Um ar auf 15 Nachkomma-Bits genau zu bestimmen genügt ein kubischer
Lookup mit 4 Teilstücken von [0,1] und 16-Bit Arithmetik. Die
Lookup-Tabelle ist in diesem Falle 20 Bytes groß.
Nun die Frage:
In der obigen Darstellung enthält die Quadratwurzel den "üblen" Teil von
arcsin in 1, und ra sowie die Konstanten den glatten, "artigen" Teil.
arcsin ist punktsymmetrisch zum Ursprung, d.h. aus Symmetriegründen geht
eine analoge Darstellung für -1.
Was ich suche ist nun eine Darstellung von arcsin, die sowohl den
Verzweigungspunkt in 1 als auch den in -1 in je eine Quadratwurzel
"auslagert", aber ansonsten nur aus analytischen Teilen besteht.
Dieser analytische Teil ist natürlich punktsymmetrisch zum Ursprung und
muss einen Konvergenzradius echt größer als 1 haben.
Allerdings gelingt es mir nicht, diese Darstellung zu finden; irgendwie
steh ich auf dem Schlauch, wie die beiden Hälften zu einem Ganzen
kombiniert werden können. Eine Division durch
Johann L. schrieb:> Um ar auf 15 Nachkomma-Bits genau zu bestimmen genügt ein kubischer> Lookup mit 4 Teilstücken von [0,1] und 16-Bit Arithmetik. Die> Lookup-Tabelle ist in diesem Falle 20 Bytes groß.
ra ist so langweilig, daß man sogar 4 gleichlange Teilstücke nehmen
könnte.
Dann musst du nicht 5 Wertepaare speichern, sondern nur die 5 ra-Werte -
mithin 10 Byte.
Außerdem genügt vielleicht sogar mit einer erträglichen Anzahl
Stützstellen eine lineare Interpolation; die wäre ggf. schneller?
Johann L. schrieb:> Nun die Frage:
Da verstehe ich die Intention vielleicht nicht richtig:
Willst du für negative x vermeiden, mit dem Betrag den von dir mit ra
beschriebenen Weg zu gehen und das Ergebnis wieder zu negieren?
PS: Dein -1<x<3 in der zweiten Gleichung erschließt sich mir nicht...
Ansonsten bin ich durchaus interessiert, auch wenn ich nichts beitragen
kann - was zum Schluß rauskommt, könnte ich in meine fixpoint-Geschichte
für C++ einbauen ...
Klaus Wachtler schrieb:> Johann L. schrieb:>> Um ar auf 15 Nachkomma-Bits genau zu bestimmen genügt ein kubischer>> Lookup mit 4 Teilstücken von [0,1] und 16-Bit Arithmetik. Die>> Lookup-Tabelle ist in diesem Falle 20 Bytes groß.>> ra ist so langweilig, daß man sogar 4 gleichlange Teilstücke nehmen> könnte.
Nehm ich ja auch. Aber mit einem linearen Lookup reicht das nicht für 15
Nachkomma-Bits. Beispiel: Unterteile [0,1] äquidistant in 4 Intervalle
und berechne mit linearem Lookup ra(7/8).
Mit kubischen Lookup sind die Stützstellen die gleichen wie bei einem
linearen, äquidistanten Lookup, die Interpolante ust jedoch keine Gerade
sondern eine kubische Parabel. Das erfordert eine Tabelle doppelter
Größe.
> Außerdem genügt vielleicht sogar mit einer erträglichen Anzahl> Stützstellen eine lineare Interpolation; die wäre ggf. schneller?
Ja klar, linear ist immer schneller als (pseudo-)kubisch.
Bei arcsin nach obiger Formel geht allerdings schon einige Zeit auf die
Berechnung der 16-Bit Wurzel (als unsigned fract aus stdfix.h) was mit
250 Ticks + 80 Bytes (MUL) bzw. 350 Ticks + 60 Bytes (kein MUL) zu Buche
schlägt. Ausserdem schluckt die Multiplikation von ra mit sqrt ca. 130
Ticks (mit MUL).
Der Code für den Lookup 3. Ordnung hat zudem den Vorteil, daß der auch
für andere Funktionen wie sin verwendet werden kann.
Taylor ist in dem Zusammenhang ein Totalausfall (langsam, Überläufe,
Rundung, ...) und CORDIC ist auch nicht schneller, hat aber zusätzlich
den Nachteil daß er stark verrauscht ist und ebenfalls in
Zwischenergebnissen überlaufen kann. D.h. man muss effektiv mit 24 Bits
rechnen um die Fehler innerhalb von ISO/IEC TC18037 zu halten.
Kubischer Lookup ist monoton per Konstruktion und kann den Code mit
anderen Funktionen teilen. Das einzige, was ausgetauscht werden muss,
ist die Tabelle. Und die ist kleiner als bei CORDIC und kleiner als bei
linearem Lookup.
Angedachter Einsatzzweck der Funktionen ist AVR-Libc, d.h. das Augenmerk
liegt in 1. Linie auf Codegröße und erst in 2. Linie auf Speed. Riesige
Tabellen sind damit ein No-Go.
Wie groß eine sin-Tabelle für 16 Bits Genauigkeit sein muss weiß ich
net, vermutlich um die 20 Einträge, siehe D3 in
Beitrag "Re: Suche speichersparende Möglichkeit für exp10()".
> Johann L. schrieb:>> Nun die Frage:>> Da verstehe ich die Intention vielleicht nicht richtig:> Willst du für negative x vermeiden, mit dem Betrag den von dir mit ra> beschriebenen Weg zu gehen und das Ergebnis wieder zu negieren?
Nein. Für negative x nimmt man natürlich
und beschränt sich auf [0,1] und macht die Berechnung komplett als
unsigned fract, das hat nämlich 1 Digit mehr als signed fract, das 1 Bit
fürs Vorzeichen braucht.
> PS: Dein -1<x<3 in der zweiten Gleichung erschließt sich mir nicht...
Die genannte Darstellung für arcsin gilt für alle x mit |1-x| < 2 und
zwar auch im Komplexen! ra ist ja analytisch in 0, d.h. für alle
komplexen Zahlen mit |x| < 2 konvergiert die o.g. Potenzreihe. Von der
Quadratwurzel nimmt man den Hauptzweig.
Die Frage bzgl. arcsin aus dem OP ist aus theoretischem Interesse und
unabhängig von irgendwelchen Berechnungen oder Algorithmen. Zudem denke
ich nicht, daß eine solche Darstellung einen Vorteil für die Berechnung
von arcsin bringt — zumindest nicht für AVRs.
Der Arcsin, als Gegenkathete zu Radius, ist fuer die winkel gegen 0 & 90
Grad schlecht definiert. Dh die Ableitung des Sinus nach dem winkel ist
nahe Null. Da hilft auch eine andere Approximation nichts.
Sieben von Neun schrieb:> Der Arcsin, als Gegenkathete zu Radius, ist fuer die winkel gegen 0 & 90> Grad schlecht definiert. Dh die Ableitung des Sinus nach dem winkel ist> nahe Null. Da hilft auch eine andere Approximation nichts.
Doch, tut es. Nicht für Rundungsfehler in der Eingabe, aber für
Rundungsfehler, die immer auch innerhalb der Algorithmen selber
auftreten.
Beispiel:
Hier reagiert der Algorithmus sehr sensibel auf kleinste Rundungsfeher
und Auslöschung bei der Wurzelberechnung, wie sie zum Beispiel für
|x| ≈ 1 auftreten. Diese Fehler pflanzen sich in die atan2-Berechnung
fort.
Fehlerfortpflanzung hat man natürlich auch mit ra und Wurzel, aber das
Verfahren ist wesentlich besser konditioniert.
Weiterer Vorteil ist, daß ra + sqrt mit dem Verfahren per Konstruktion
monoton ist, d.h. für größeres x wird das Ergebnis größer (weil arcsin
monoton wächst).
Bei anderen Verfahren ist das i.d.R. nicht der Fall, es sei denn man
schraubt die interne Genauigkeit massiv nach oben. Ansonsten sind die
Werte stark verrauscht und fehlerbehaftet bis zur Unbrauchbarkeit.
Vergleiche z.B. folgende Algorithmen für sin und cos, einmal mit
linearem Lookup und einmal mit CORDIC. Und dabei sind sin und cos noch
nicht mal übellaunig wie arcsin!
http://www.mikrocontroller.net/articles/AVR_Arithmetik/Sinus_und_Cosinus_(CORDIC)#Genauigkeithttp://www.mikrocontroller.net/articles/AVR_Arithmetik/Sinus_und_Cosinus_(Lineare_Interpolation)#Genauigkeit
Eine eigene Frage?
Dieser Thread ist über 14 Jahre alt.
Stell sie in einem eigenen Thread – dort lesen mehr Leute mit, und sie wird öfter und schneller beantwortet.
