Arcus Sinus und seine Verzweigungspunkte

[Johann L. (gjlayde)]

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:

$$\text{ra}(x) = \frac{\frac \pi2 - \arcsin (1-x)}\sqrt{2x} = 1 + \frac1{12}x + \frac3{160}x^2 + \cdots$$

ra ist analytisch in 0 mit Konvergenzradius 2.

arcsin ist in einer Umgebung von 1 also i.W. die Quadratwurzel:

$$\arcsin(x) = \frac \pi2 - \sqrt{2-2x}\cdot\text{ra}(1-x) \quad\text{ mit } -1 < x < 3$$

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

$$\sqrt{1-x^2}$$

tut's übrigens nicht.

[Klaus W. (mfgkw)]

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 W. (mfgkw)]

Für Mitleser der Verlauf von ra

[Johann L. (gjlayde)]

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

$$\arcsin(x) = -\arcsin(-x)$$

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.

[Purzel H. (hacky)]

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.

[Klaus W. (mfgkw)]

Sieben von Neun schrieb:
> ...schlecht definiert...

das heißt ja nicht, daß man den asin nicht mehr oder weniger effizient 
berechnen kann.

[Johann L. (gjlayde)]

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:

$$\arcsin(x) = \text{atan2}(x,\sqrt{1-x^2})$$

 
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)#Genauigkeit
http://www.mikrocontroller.net/articles/AVR_Arithmetik/Sinus_und_Cosinus_(Lineare_Interpolation)#Genauigkeit