mikrocontroller.net

Forum: Mikrocontroller und Digitale Elektronik Näherung für ArcusTangens?


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
Autor: Johann L. (gjlayde) Benutzerseite
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Hi, kennt jemand eine praktikable Näherungsformel für atan (Arcus 
Tangens) im Intervall [0,1]?

Anzunähern ist in double-Präzision, also mit (mindestens) 53 Bits 
Genauigkeit.  Intern wird mit 56 Bits gerechnet.

Bekanntlich hat die Taylor-Reihe um 0 schleichte Konvergenz weil der 
Konvergenzradius nur 1 ist, d.h. in der Nähe von 1 hat man (praktisch) 
keine Konvergenz mehr.  Daher untersuche ich gerade Padé-Approximation, 
also eine Entwicklung als rationale Funktion, die ansonsten die gleichen 
Nebenbedingungen im Entwicklungspunkt hat wie eine Taylor-Reihe: dort 
stimmen die ersten n Ableitungen der Approximante mit den ersten n 
Ableitungen der Funktion überein.

Entwicklungspunkt 0 hat wegen der Symmetrie von atan den Vorteil, dass 
bei gleicher Anzahl Terme der Grad der Approximation doppelt so groß 
ist.  Die Padé Approximation um 0 konvergiert zwar bis jenseits von 1, 
aber man immer noch zu viele Terme.

Daher hab ich das Intervall [0,1] aufgespalten in 2 Teilintervalle: Eine 
Approximation vom Grad [13/14] um 0 (15 Terme) approximiert in [0, 
0.55], und im Intervall [0.55, 1] wird mit Grad [7/7] um 191/256 
approximiert (ebenfalls 15 Terme).

Konkret für die 2. Approximante für [0.55, 1] um 191/256:
atan(191/256) +
(1109278398945277954552108998515656949760x
+3190014013046549743618019835405171425280x^2
+4409271226637051705982201724529121689600x^3
+3567894272956348772810533213265356390400x^4
+1761790565199680283994647188193675313152x^5
+498498185308650458726160978418683871232x^6
+62925645331517674324275193263427682304x^7)
/
(1726764746478277909691505488610409165095
+5793379763241549317328788011779036491520x
+9481302008830900338634583261911194992640x^2
+9412956020430072773739129723727695052800x^3
+6049687104603154491828695341359484108800x^4
+2495103530011474195120698950231953244160x^5
+608983183658054509908458581941561589760x^6
+67669375374053976356851630199244062720x^7)
Für einen ATmega4808 ergeben sich dabei in etwa folgende Kosten:
* 400 Ticks pro Addition
* 600 Ticks pro Multiplikation
* 6000 Ticks pro Division
insgesamt also 2*div + 14*add + 14*mul = ca. 26000 Ticks.  Zu speichern 
sind 31 Terme à 10 Bytes pro Term.

Gibt es da effizientere Alternativen?  Im Netz hab ich nicht wirklich 
viel dazu gefunden.

Cordic hab ich noch nicht angetestet; der liefert pro Schritt 1 Bit; bei 
56 Bits intern braucht's also 56 Schritte. Und man braucht 1 Konstante 
zur Korrektur pro Schritt, also 560 Bytes an Lookup für die Konstanten. 
Außerdem verliert Cordic pro Schritt mindestens 1/2 LSB, insgesamt also 
56/2 LSB ~ 5 Bit.  Daher ist zu erwarten, dass mit Cordic der Fehler 
nicht kleiner als 2^-51 zu bekommen ist.

Padé liefert im schlechtesten Falle (x = 0.55) hingegen einen Fehler von 
2^-55.2 (theoretisch), womit die Genauigkeit für 53 Bit double-Mantisse 
zu erreichen sein sollte.

Ideen?

Autor: minifloat (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Sinustabelle 0...90° und der Observer aus der NXP AppNote hier:
https://www.google.com/search?q=resolver+observer+an3943

Siehe auch Anhang.
Da musst du halt schauen, wieviele Iterationen der braucht, man kann den 
Winkel aber je Quadrant mit einer Schätzung schon vorbereiten, dass die 
Zielwinkelabweichung mit <±45° oder <±22,5° startet. Das Ding lässt sich 
ohne Division aufbauen :)

Autor: A. S. (achs)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Was ist denn die Referenz? Als wieviel Takte braucht atan dort?

Autor: minifloat (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Add1 die Sinustabelle lässt sich Unterpolieren, damit kommst du evtl. 
mit 2k Sinustabelle aus

Add2 beim sin und cos, die du rein fütterst, macht dir die float sowieso 
eine Normierung auf maximal große Zahlen. Kommst du mit 64bit integer 
aus? Dann reicht zur Berechnung eine 128bit Integerarithmetik.

Add3 muss es float64 sein?

Autor: Franko S. (frank_s866)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Gibt es da effizientere Alternativen?  Im Netz hab ich nicht wirklich
> viel dazu gefunden.

Darüber habe ich erst gestern was gelesen. Da stand auch nur dass man 
sich auf 0..1 beschränken solle wenn es geht, besser wäre sein Problem 
nochmal anzuschauen und ob man es ohne arctan lösen kann, 
quadrantenweise Fallunterscheidung, Steigung,... das geht für viele 
Probleme aber eben nicht für alle.

Zum Trost: arcsin, arccos sind noch schlimmer.

Autor: minifloat (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
A. S. schrieb:
> Als wieviel Takte braucht atan dort?

Pro Zyklus Observer
6x Mul
4x Add

Pro Zyklus Sin und Cos mit linearer Interpolation y = 
y1+((x-x1)*((y2-y1)*(1/dx)))

4x Table-Lookup
6x Add
4x Mul

Summe pro Zyklus also
4x Table-Lookup
10x Mul
10x Add

Ergo       10kTicks

Mit Vorbesetzen des Winkels einmalig 4 Fallunterscheidungen und 0-2 
Multiplikationen

Hmmm, sieht gar nicht mal so gut dagegen aus...

Autor: Dirk B. (dirkb2)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
wie ist der arctan denn mit CORDIC ?

Autor: minifloat (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Zu speichern sind 31 Terme à 10 Bytes pro Term.

Der Observer braucht max. 4 Zwischenwerte auf einmal.

Autor: Egon D. (egon_d)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:

> Cordic hab ich noch nicht angetestet; der liefert pro
> Schritt 1 Bit; bei 56 Bits intern braucht's also
> 56 Schritte.

Korrigiere mich, wenn ich mich irre, aber CORDIC ist
doch im Prinzip ein Iterationsverfahren auf Basis der
Intervallhalbierung, oder?

Was wäre mit einem Sender-Jerewan-CORDIC: Im Prinzip
ja, aber erstens verwenden wir eine bessere Start-
näherung, zweitens versuchen wir, mehr als ein Bit
je Schritt zu gewinnen, und drittens verwenden wir
ein Iterationsverfahren von höherer Konvergenzordnung.

Nur mal spekuliert: Startnäherung mittels Taylorreihe
auf z.B. 16bit; das müsste noch recht fix gehen.
Iterationsverfahren mit quadratischer Konvergenz
braucht nochmal zwei Schritte, um in die Nähe von
64 geltenden Bit zu kommen.


> Außerdem verliert Cordic pro Schritt mindestens 1/2 LSB,
> insgesamt also 56/2 LSB ~ 5 Bit.  Daher ist zu erwarten,
> dass mit Cordic der Fehler nicht kleiner als 2^-51 zu
> bekommen ist.

Das Problem entfällt bei einem echten Iterationsverfahren.

Autor: Johann L. (gjlayde) Benutzerseite
Datum:
Angehängte Dateien:

Bewertung
2 lesenswert
nicht lesenswert
minifloat schrieb:
> Add3 muss es float64 sein?

Franko S. schrieb:
> besser wäre sein Problem nochmal anzuschauen und ob man es ohne
> arctan lösen kann,

In dem Fall nicht, Fernziel wäre Implementierung als Teil der 
Standard-Lib im Rahmen von avr-gcc, d.h. als Teil der libgcc und / oder 
avr-libc.  Von daher hab ich auch kein Problem mit salted Code.


A. S. schrieb:
> Was ist denn die Referenz? Als wieviel Takte braucht atan dort?

Als Referenz hab ich momentan nur avr_f64.c/h von

Beitrag "Re: 64 Bit float Emulator in C, IEEE754 compatibel"

Zu der Implementierung kann ich bisher sagen, dass sie Out-of-the-Box 
problemlos funktioniert.  Bei PROGMEM fehlen const und für C++ braucht's 
ein extern "C" Wrapper, das war's.  Die Implementierung ist vanilla C, 
also prinzipiell einfach portabel.  Leider ist die Implementierung 
stellenweise etwas kryptisch; zumindest for mich.

avr_f64.c::f_arctan() belegt folgende Resourcen:
* Ticks: 40000
* Stack: 216 Bytes
* Code: 6 KiB

Für ein paar Funktionen, die ich bereits habe, hab ich mit avr_f64 
verglichen um sicher zu sein, dass ich nicht kompletten Unsinn treibe. 
Meine Implementierung hat nicht Portierbarkeit als Ziel, sonder vor 
allem Effizienz in den 3 wesentlichen Bereichen
* Ausführungszeit
* Code- bzw. Flash-Verbrauch
* Stack- bzw. RAM-Verbrauch
und zwar auch im Hinblich auch die anderen double-Funktionen, d.h. 
Wiederverwendbarkeit anderer Routinen wie Polynom-Auswertung, die vom 
Anwender ebenfalls (implizit) verwendet werden, spielt auch eine Rolle.

siehe dazu die Anhänge.  Übersetzt wird i.W. mit:
    -Os -ffunction-sections -fdata-sections -mcall-prologues -mstrict-X 
-mrelax -Wl,--gc-sections
und für Messung der Codegröße zusätzlich
    -nostartfiles -Wl,--defsym,main=0

f7_t ist die expandierte 10-Byte Darstellung (1 Byte Flags, 7 Byte 
Mantisse, 2 Byte Exponent).  f7 sind Werte für C/C++ IEEE 
double-Wrapper; hier kommen Resourcen für f7_t <-> double hinzu.

Als Resourcen für Padé atan erwarte ich 100-150 Bytes Stack, 3-4 KiB 
Code (bei ansonsten leerem Programm) und 20000-30000 Ticks je nachdem, 
ob 1 oder 2 Divisionen gebraucht werden.  Was die Codegröße betrifft, so 
werden viele Routinen verwendet, die auch von anderen Funktionen 
gebraucht werden: add, sub, mul, div, Vergleiche und Polynom-Auswertung, 
so dass der effektive Codezuwachs in realer Anwendung deutlich kleiner 
ausfallen wird.

Wie asin und acos implementiert werden sollten, weiß ich auch noch 
nicht.  Via atan hätte den Vorteil von Code-Sharing, allerdings kämen 
mindestens 1 Division und 1 Quadratwurzel hinzu, also nochmal 20000 
Ticks...

Franko S. schrieb:
> Zum Trost: arcsin, arccos sind noch schlimmer.

Die stehen dann als nächstes auf der Liste.

Autor: whitespace (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Konkret kann ich nur einen Literaturvorschlag beisteuern: Computer 
Approximations, John F. Hart, Wiley, 1968

Das Buch hatte ich in einer Bibliothek zu solchen Themen einmal in den 
Händen und durchgeblättert. Darin gibt es auch Kapitel über die Näherung 
von trigonometrischen Funktionen und den Inversen. Verschiedene 
Approximationsmethoden werden auch beschrieben. Ist schon etwas älter, 
ob die Methoden noch "state of the art" sind oder ob es mittlerweile 
bessere gibt, weiß ich nicht. Ich habe die darin beschriebenen Methoden 
auch nie selber verwendet und das Buch nur aus Interesse 
durchgeblättert.

Autor: A. B. (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Wenn eine gute (genau und schnell) Quadratwurzelfunktion schon da ist, 
kann man mittels

 arctan(x)=2*arctan(x/(1+sqrt(1+x^2)))

den problematischen Bereich von [0.5, 1.0] so etwa auf [0.25, 0.42] 
reduzieren. Ggf. auch gleich ohne Fallunterscheidung [0, 1.0] auf [0, 
0.42].

Ansonsten würde ich mal Entwicklung nach Legendre- oder 
Chebyscheff-Polynomen versuchen.

Autor: Joggel E. (jetztnicht)
Datum:

Bewertung
-7 lesenswert
nicht lesenswert
Trollalarm !!!

wer benotigt einen Arctan mit 53 bit Praezision ... Leute, die das 
benoetigen haben ueblicherweise eine Floatingpointeinheit.

Autor: Jörg W. (dl8dtl) (Moderator) Benutzerseite
Datum:

Bewertung
4 lesenswert
nicht lesenswert
Joggel E. schrieb:
> Trollalarm !!!

[ ] Du weißt, über wen du gerade sprichst.

Autor: Joggel E. (jetztnicht)
Datum:

Bewertung
-3 lesenswert
nicht lesenswert
Nein. Es ist ja schoen, wenn man das machen moechte. Und dafuer eine 
Library aufblasen moechte. Minimale Codegroesse und Ausfuehrungszeit ist 
auch nett.

Nur soll die auf den ATMega festgelegt werden. Dann kann ich in einen 
grossen ATMega die Library linken plus 2 Operationen und dann ist voll.

Ungeachtet dem, muss man eben alle Moeglichkeiten abklopfen. Taylor, 
stueckweise Taylor, cordic, invers iterativ usw. Und nicht gleich alles 
was ineffizient erscheint verwerfen. Allenfalls kann ein Algorithmus 
aufgrund der Rahmenbedingungen auch beschleunigt werden.

Das gibt's auch noch wissenschaftliche Publikationen. Das Problem wurde 
schon behandelt, sonst waere es nicht in Hardware in einer Floating 
point einheit  implementiert. Der Intel 8087 konnte das schon vor 35 
Jahren.

Als Anwender muss man sich sowieso ueberlegen, ob man diese Genauigkeit 
benotigt... eine Floatlibrary wird ja immer als Ganzes geladen.

Autor: Walter T. (nicolas)
Datum:

Bewertung
3 lesenswert
nicht lesenswert
An dieser Stelle: Jetzt nicht.

In der Allgemeinheit stimmen die Ausführungen zwar - hätte ich nicht den 
Namen des Fragestellers gekannt, wäre mir die angestrebte Genauigkeit 
auch zweifelhaft erschienen.

Aber jetzt und hier können wir davon ausgehen, daß der TO genau weiß, 
was er bezweckt.

Autor: Jörg W. (dl8dtl) (Moderator) Benutzerseite
Datum:

Bewertung
2 lesenswert
nicht lesenswert
Joggel E. schrieb:
> eine Floatlibrary wird ja immer als Ganzes geladen.

Gut, du kennst also nicht nur Johann nicht, sondern du hast auch nicht 
den blassesten Schimmer, wie ein Linker mit Libraries umgeht.

Das disqualifiziert dich in dieser Hinsicht komplett.

Autor: Axel S. (a-za-z0-9)
Datum:

Bewertung
2 lesenswert
nicht lesenswert
Joggel E. schrieb:

[belanglos]

Von Dieter Nuhr gibt es ein Zitat, das paßt auf dich wie die Faust aufs 
Auge. Man möchte fast meinen, er würde dich kennen.

https://www.google.com/search?&q=Nuhr%20Ahnung

Autor: Dieter R. (drei)
Datum:

Bewertung
-2 lesenswert
nicht lesenswert
Joggel E. schrieb:
> Trollalarm !!!
>
> wer benotigt einen Arctan mit 53 bit Praezision ... Leute, die das
> benoetigen haben ueblicherweise eine Floatingpointeinheit.

Abgesehen davon, dass der Trollalarm eher auf Ignoranz beruht:

mich würde auch sehr ehrlich interessieren, wofür man das brauchen kann, 
also für welche konkrete Applikation.

Kann da jemand zu meiner Wissensmehrung beitragen? Danke!

Autor: Jörg W. (dl8dtl) (Moderator) Benutzerseite
Datum:

Bewertung
3 lesenswert
nicht lesenswert
Dieter R. schrieb:
> mich würde auch sehr ehrlich interessieren, wofür man das brauchen kann,
> also für welche konkrete Applikation.

Wenn man (weil man wie Johann halt Mit-"Anbieter" der Toolchain ist) 
eine komplette Standardbibliothek aufbauen möchte, dann muss eben auch 
ein Arcustangens mit dabei sein. Dann ist es recht und billig, sich um 
die Effizienz der Implementierung Gedanken zu machen.

Autor: Jörg W. (dl8dtl) (Moderator) Benutzerseite
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@Johann

Kennst du das hier?

Beitrag "64 Bit float Emulator in C, IEEE754 compatibel"

Keine Ahnung, wie genau oder schnell es ist.

Autor: Dieter R. (drei)
Datum:

Bewertung
-3 lesenswert
nicht lesenswert
Jörg W. schrieb:
> Dieter R. schrieb:
>> mich würde auch sehr ehrlich interessieren, wofür man das brauchen kann,
>> also für welche konkrete Applikation.
>
> Wenn man (weil man wie Johann halt Mit-"Anbieter" der Toolchain ist)
> eine komplette Standardbibliothek aufbauen möchte

Eine Standardbibliothek und die Arbeit daran ist aber nur nützlich, wenn 
es auch Applikationen dafür gibt. Danach hatte (nicht nur) ich gefragt.

Autor: Jörg W. (dl8dtl) (Moderator) Benutzerseite
Datum:

Bewertung
4 lesenswert
nicht lesenswert
Dieter R. schrieb:
> Eine Standardbibliothek und die Arbeit daran ist aber nur nützlich, wenn
> es auch Applikationen dafür gibt.

… oder wenn es eben ein Hobby ist.

Applikationen finden sich immer. :) Ich bin zumindest schon in 
Situationen gelaufen, wo die Rundungsfehler der 32-bit-FP-Arithmetik 
schon deutlich wurden, und ich mir die Option von 64-bit-FP gewünscht 
hätte. Wenn Johanns Implementierung ähnlich effizient ist wie die 
existierende 32-bit-FP (logischerweise zuzüglich Mehrbedarf für die 
höhere Genauigkeit), dann wird sie durchaus ihre Anwender finden. Die 
existierende 32-bit-FP-Arithmetik in der AVR-GCC-Toolchain hat 
jedenfalls aus meiner Sicht sehr eindrücklich gezeigt, dass das alte 
Dogma "Gleitkomma macht man auf einem 8-bit-Prozessor nicht" (*) getrost 
ad acta gelegt werden kann.

(*) Das stimmte ja schon mit Turbo-Pascal auf CP/M nicht mehr. 
Allerdings hat man dort als Kompromiss eine 48-Bit-Implementierung 
gewählt, was übrigens sogar C-Standard-konform als "double" durchgehen 
würde. Nur 32 bit tut's nicht.

: Bearbeitet durch Moderator
Autor: Alexander S. (alesi)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Hi, kennt jemand eine praktikable Näherungsformel für atan (Arcus
> Tangens) im Intervall [0,1]?

Hallo,

wenn man auf https://ch.mathworks.com/ nach "Calculate Fixed-Point 
Arctangent" sucht gelang man zu dieser Seite
https://ch.mathworks.com/help/fixedpoint/examples/calculate-fixed-point-arctangent.html?s_tid=srchtitle
Dort gibt es Infos zu
Calculating atan2(y,x) Using the CORDIC Algorithm
und
Calculating atan2(y,x) Using Chebyshev Polynomial Approximation
und
Calculating atan2(y,x) Using Lookup Tables

Johann L. schrieb:
> Bekanntlich hat die Taylor-Reihe um 0 schleichte Konvergenz weil der
> Konvergenzradius nur 1 ist,

Neben der Taylor-Reihe gibt es auch noch die Euler-Reihe für den arctan
https://www.cambridge.org/core/journals/mathematical-gazette/article/8967-an-elementary-derivation-of-eulers-series-for-the-arctangent-function/3D53E3BC23BEAF049A9F1013E89E9F46

Dann habe ich noch diesen Artikel gefunden:
“Efficient approximations for the arctangent function”, Rajan, S. Sichun 
Wang Inkol, R. Joyal, A., May 2006
http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf

Vielleicht hilft das etwas.

Autor: Alexander S. (alesi)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo nochmal,

in dieser Datei nach "r8_atan" suchen: R8_ATAN evaluates the arc-tangent 
of an R8 argument.

https://people.sc.fsu.edu/~jburkardt/c_src/fn/fn.c

Autor: Johann L. (gjlayde) Benutzerseite
Datum:

Bewertung
2 lesenswert
nicht lesenswert
Jörg W. schrieb:
> @Johann
>
> Kennst du das hier?
>
> Beitrag "64 Bit float Emulator in C, IEEE754 compatibel"
>
> Keine Ahnung, wie genau oder schnell es ist.

Ja, kenn ich.  Dient wie gesagt als Basis für den Vergleich von hier:

Beitrag "Re: Näherung für ArcusTangens?"

Dieter R. schrieb:
> mich würde auch sehr ehrlich interessieren, wofür man das brauchen kann,
> also für welche konkrete Applikation.

Eine Anwendung ist Auswertung von GPS-Daten, wo Genauigkeit von float 
nicht ausreicht.  Üblicherweise sind das Bereichnungen aus der 
sphärischen Geometrie:  Zuerst Winkelfunktionen, Berechnugen auf deren 
Werten und dann eine Arcus-Funktion so wie hier:

Beitrag ""echter" double mit AVR-GCC"

Da wird zwar asin verwendet, zeigt aber dass float nicht ausreicht und 
dass gelegentlich, wenn auch selten, Bedarf für double besteht.

Außerdem ist an avr-gcc prinzipiell anzumäkeln, dass double nicht 
standardkonform ist.  Wäre ein rein theoretisches Problem wenn niemand 
double brauchen würde; aber auch der Autor der o.g. avr_f64.c/h hat 
seine double-Implementierung gemacht weil er sie brauchte, um ein 
konkretes Problem zu lösen.  AFAIR gab es auch bei avrfreaks 
entsprechende Anfragen.

Autor: Dieter R. (drei)
Datum:

Bewertung
-3 lesenswert
nicht lesenswert
Johann L. schrieb:

> Eine Anwendung ist Auswertung von GPS-Daten, wo Genauigkeit von float
> nicht ausreicht.  Üblicherweise sind das Berechnungen aus der
> sphärischen Geometrie

Ah, danke. Endlich mal eine kompetente Auskunft dazu!

Autor: W.S. (Gast)
Datum:

Bewertung
-1 lesenswert
nicht lesenswert
Johann L. schrieb:
> Cordic hab ich noch nicht angetestet; der liefert pro Schritt 1 Bit; bei
> 56 Bits intern braucht's also 56 Schritte. Und man braucht 1 Konstante
> zur Korrektur pro Schritt, also 560 Bytes an Lookup für die Konstanten.

Dann lies erstmal die Literatur zum CORDIC.

Kurz gesagt, wenn du keine Schritte ausläßt, sondern jeden Schritt 
ausführst (mit wechselnder Richtung, so daß der Restwinkel gegen Null 
läuft), dann ist dein Rest nach endlich vielen Schritten in seiner 
Amplitude determiniert (da ja alle Schritte ausgeführt werden und jeder 
einen Faktor 1/cos(d) einbringt).

Dann kannst du den Rest zum Berechnen des Endwertes heranziehen, da du 
ihn mit dem konstanten Produkt aus allen 1/cos(dn) wieder normalisieren 
kannst.

Der Rest enthält so etwa die doppelte Stellenzahl wie das, was du in den 
Drehschritten errechnet hast.

In Wahrheit ist es ein bissel weniger, aber man spart sich trotzdem 2/3 
der Schritte ein. Also brauchst du für 56 Bit nur 19 Cordic-Schritte, 
dann eine Multiplikation des Restes mit o.g. 1/cos.. Korrekturfaktor und 
Addition des Restes.

W.S.

Autor: Egon D. (egon_d)
Datum:

Bewertung
1 lesenswert
nicht lesenswert
W.S. schrieb:

> Dann kannst du den Rest zum Berechnen des Endwertes
> heranziehen, da du ihn mit dem konstanten Produkt aus
> allen 1/cos(dn) wieder normalisieren kannst.
>
> Der Rest enthält so etwa die doppelte Stellenzahl wie
> das, was du in den Drehschritten errechnet hast.

Sehr hübsch.
Hast Du zufällig auch eine ungefähre Quellenangabe zur
Hand?

Autor: Helmut S. (helmuts)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
whitespace schrieb:
> Konkret kann ich nur einen Literaturvorschlag beisteuern: Computer
> Approximations, John F. Hart, Wiley, 1968
>
> Das Buch hatte ich in einer Bibliothek zu solchen Themen einmal in den
> Händen und durchgeblättert. Darin gibt es auch Kapitel über die Näherung
> von trigonometrischen Funktionen und den Inversen. Verschiedene
> Approximationsmethoden werden auch beschrieben. Ist schon etwas älter,
> ob die Methoden noch "state of the art" sind oder ob es mittlerweile
> bessere gibt, weiß ich nicht. Ich habe die darin beschriebenen Methoden
> auch nie selber verwendet und das Buch nur aus Interesse
> durchgeblättert.

Ich habe mir das Buch interessehalber heute mal angeschaut. Das sieht 
gut aus.

arctan(x) = x*P(x^2)/Q(x^2)

Es werden also nur x, x^3, x^5, .... verwendet.

Es gibt auch ein Polynom x*P(x^2) aber nur bis 10 Dezimalstellen.

Wenn die Zahl unter Precision wirklich Dezimalstellen sind, was ich auch 
annehme, dann reichen 6 Koeffizienten im Zähler und 6 Koeffizienten im 
Nenner um die 53bit Genauigkeit zu erreichen. Da wird man aber nicht 
umhin kommen mit 64bit Mantisse zu rechnen. Das Buch liegt z. B. in der 
UB Stuttgart im Magazin.

: Bearbeitet durch User
Autor: Percy N. (vox_bovi)
Datum:

Bewertung
-2 lesenswert
nicht lesenswert
Jörg W. schrieb:
>
> (*) Das stimmte ja schon mit Turbo-Pascal auf CP/M nicht mehr.
> Allerdings hat man dort als Kompromiss eine 48-Bit-Implementierung
> gewählt, was übrigens sogar C-Standard-konform als "double" durchgehen
> würde. Nur 32 bit tut's nicht.

Schon vorher gab es unter cp/m von MS BASIC80 als Compiler und Mbasic 
als Interpreter - zumindest letzterer eher lahm, und die 
"Assembler-Schnittstelle" beschränkte sich (bei Mbasic) auf PEEK und 
POKE, aber IO-Befehle gehörten zur Sprache ebenso wie 8-byte-reals - nur 
leider keine transzendenten Funktionen über 4-byte-reals hinaus, iirc.

Autor: Jörg W. (dl8dtl) (Moderator) Benutzerseite
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Percy N. schrieb:
> Schon vorher gab es unter cp/m von MS BASIC80 als Compiler und Mbasic
> als Interpreter

Mit denen hatte ich nicht viel zu tun. Turbo-Pascal war so viel schöner 
und bequemer und eben – wie der Namen nicht nur versprach – verdammt 
schnell.

Autor: Alexander S. (alesi)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

ganz unten auf dieser Seite ist ein code atan_137 für ca. 13.7 
Dezimalstellen:   http://www.ganssle.com/approx.htm

Autor: Helmut S. (helmuts)
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Hallo,
>
> ganz unten auf dieser Seite ist ein code atan_137 für ca. 13.7
> Dezimalstellen:   http://www.ganssle.com/approx.htm

Dieses Programm verwendet exakt eines der Polynome aus dem bereits 
zitierten Buch "Computer Approximations, John F. Hart, Wiley, 1968".
Die angegebene Genauigleit von 13,7 Stellen dieser Approximation gilt 
aber nur bis zum Winkel pi/12. Gefordert wird aber der Winkelbereich 0 
bis pi/2 -  entspricht arctan(0 bis 1).

Autor: Alexander S. (alesi)
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Helmut S. schrieb:
> Die angegebene Genauigleit von 13,7 Stellen dieser Approximation gilt
> aber nur bis zum Winkel pi/12. Gefordert wird aber der Winkelbereich 0
> bis pi/2 -  entspricht arctan(0 bis 1).

Dafür wird auf den "range reduction code" weiter oben auf der Seite 
verwiesen: "atan_137 computes the arctangent to about 13.7 decimal 
digits of accuracy using a simple rational polynomial. It's input range 
is 0 to Π /12; use the previous range reduction code."

Autor: Helmut S. (helmuts)
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Helmut S. schrieb:
>> Die angegebene Genauigleit von 13,7 Stellen dieser Approximation gilt
>> aber nur bis zum Winkel pi/12. Gefordert wird aber der Winkelbereich 0
>> bis pi/2 -  entspricht arctan(0 bis 1).
>
> Dafür wird auf den "range reduction code" weiter oben auf der Seite
> verwiesen: "atan_137 computes the arctangent to about 13.7 decimal
> digits of accuracy using a simple rational polynomial. It's input range
> is 0 to Π /12; use the previous range reduction code."

OK. Das hatte ich übersehen.
Statt einem "range reduction code" könnte man auch gleich mit 5 
Koeffizienten im Zähler und 5 Koeffizienten im Nenner arbeiten. Damit 
bekäme man 14,5 Stellen.
Der Fragesteller hätte gerne mindestens 16 Stellen (53bit). Dazu 
benötigt man 6 Koefizienten im Zähler und 6 Koeffizienten im Nenner.

Nachtrag zu meiner vorherigen Antwort - dort sollte pi/4 statt pi/2 
stehen.

Autor: Johann L. (gjlayde) Benutzerseite
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> "atan_137 computes the arctangent to about 13.7 decimal
> digits of accuracy using a simple rational polynomial.

Wobei man da aufpassen muss, denn soweit ich sehe sind für atan_137 
absolute Fehler angegeben.  Bei Floats zähl aber der relative Fehler, 
d.h. in der Nähe von 0 muss der absolute Fehler dann nochmal kleiner 
sein.

Jack Ganssle:
> All of the approximations here are polynomials, or ratios of
> polynomials.  All use very cryptic coefficients derived from
> Chebyshev series and Bessel functions.

Wie werden hier Besselfunktionen eingesetzt?

Bzw. wie bekommt man beste (im Sinne der max-Norm) rationale 
Approximationen? Für Polynome ist es Chebyshev, aber für rationale 
Funktionen kenn ich nur das Analogon zu Taylor (Padé), was ist das 
rationale Analogon zu Chebyshev?

Autor: Alexander S. (alesi)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

hier die Ausgabe von
printf("%4.2f %18.16f %18.16f %11.4e\n", x, atan_137(x), atan(x), atan_137(x)-atan(x));
0.00 0.0000000000000000 0.0000000000000000  0.0000e+00
0.05 0.0499583957219434 0.0499583957219428  6.4532e-16
0.10 0.0996686524911621 0.0996686524911620  8.3267e-17
0.15 0.1488899476094955 0.1488899476094973 -1.8041e-15
0.20 0.1973955598498836 0.1973955598498808  2.7756e-15
0.25 0.2449786631268659 0.2449786631268641  1.7208e-15
0.30 0.2914567944778712 0.2914567944778671  4.1078e-15
0.35 0.3366748193867234 0.3366748193867272 -3.7748e-15
0.40 0.3805063771123672 0.3805063771123648  2.3315e-15
0.45 0.4228539261329406 0.4228539261329406  0.0000e+00
0.50 0.4636476090008049 0.4636476090008061 -1.1657e-15
0.55 0.5028432109278610 0.5028432109278608  2.2204e-16
0.60 0.5404195002705838 0.5404195002705842 -3.3307e-16
0.65 0.5763752205911844 0.5763752205911837  6.6613e-16
0.70 0.6107259643892097 0.6107259643892087  9.9920e-16
0.75 0.6435011087932825 0.6435011087932845 -1.9984e-15
0.80 0.6747409422235513 0.6747409422235527 -1.4433e-15
0.85 0.7044940642422213 0.7044940642422178  3.4417e-15
0.90 0.7328151017865063 0.7328151017865068 -4.4409e-16
0.95 0.7597627548757679 0.7597627548757709 -2.9976e-15
1.00 0.7853981633974536 0.7853981633974484  5.2180e-15
und zur Kontrolle den mit bc berechneten atan Werten:
0
.0499583957219427
.0996686524911620
.1488899476094972
.1973955598498807
.2449786631268641
.2914567944778670
.3366748193867271
.3805063771123648
.4228539261329407
.4636476090008061
.5028432109278608
.5404195002705841
.5763752205911836
.6107259643892086
.6435011087932843
.6747409422235526
.7044940642422177
.7328151017865065
.7597627548757708
.7853981633974483
Auf der Seite der Ganssle Group mit den auf Hart, "Computer 
Approximations" basierenden Routinen http://www.ganssle.com/approx.htm 
gibt es auch den source code der Routinen 
http://www.ganssle.com/approx/sincos.cpp
Die Ausgabe habe ich angehängt. Die erste Spalte ist der Winkel in Grad. 
Das wird intern in rad gewandelt. Dabei habe ich nur die fehlende 
Ausgabe der atan Werte korrigiert. g++ gibt noch einige Warnungen aus.

: Bearbeitet durch User
Autor: Alexander S. (alesi)
Datum:

Bewertung
0 lesenswert
nicht lesenswert

Autor: Johann L. (gjlayde) Benutzerseite
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Hab noch dieses Paper von Castellanos und Rosenthal gefunden: "Rational 
Chebyshev Approximations of Analytic Functions"

https://www.fq.math.ca/Scanned/31-3/castellanos.pdf

Danke erstmal für all die Hinweise und Links, ist vorerst mal einiges an 
Hirnfutter...

Autor: Rainer V. (a_zip)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
...ein selten spannender und ergiebiger Faden! Danke.
Gruß Rainer

Autor: Hamburger (Gast)
Datum:

Bewertung
-1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Auf der Seite der Ganssle Group

Ist das in Zeiten von viel RAM und schneller Interpolation noch 
zeitgemäß? Hier hammse schon normale Sinusfunktionen in Tabellen und 
interpolieren, weil es schneller ist ..

Autor: Rainer V. (a_zip)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hamburger schrieb:
> Ist das in Zeiten von viel RAM und schneller Interpolation noch
> zeitgemäß? Hier hammse schon normale Sinusfunktionen in Tabellen und
> interpolieren, weil es schneller ist ..

Ich glaube, du wirfst da was durcheinander....
Gruß Rainer

Autor: Detlef _. (detlef_a)
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Liebe Rechenfreunde, hallo Johann,

ich brauchte ein wenig Spass und habe den atan Algorithmus für den 
Argumentbereich [0 , 1] in Cordic und reiner 64 Bit integer Rechnung 
implementiert.

Beitrag "Atan function implementation"

Die gemessene Genauigkeit sollte besser als 51 Bit sein, mehr kann ich 
nicht messen weil die double Implementierung von atan nicht besser ist. 
Ausgelegt ist der Algo auf 61 Bit damit man mit dem 64 Bit integer Type 
rechnen kann. Ich habe jetzt nicht ausufernd testen können, vllt. kommt 
da noch ein Korken hoch aber so komplex ist die Sache am Ende auch 
nicht, der Weg war etwas steinig gewesen wie man am auskommentierten 
Code erahnen kann.

Jederzeit gerne Kommentare und Anmerkungen.

War grosser Spaß jewesen!
math rulez!
Cheers
Detlef

Autor: Rainer V. (a_zip)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Man oh Man, bin forth-freund und habe neben den üblichen 
Implementierungen, sowas noch nicht gesehen. Aber immerhin.
Gruß Rainer

Autor: Martin Cibulski (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Johann,

ich habe den ArcTan für meine Teleskopsteuerung benötigt, programmiert 
habe ich damals in Assembler, weil es für den AtMega128 in der C-Library 
nur Fließkommazahlen in einfacher Genauigkeit gab.

Basis für die Algorithmen war eine C-Library ('Cephes Mathematical 
Library') von Steve Moshier.
http://www.netlib.org/cephes/index.html und
http://moshier.net/#Cephes

Hier meine Implementation in Asm, die seit Jahren zuverlässig 
funktioniert:
http://martin-cibulski.de/atm/mount_controller_4/mc_4_asm_main/float_atan.asm

Wenn ich mich richtig erinnere, läuft die Routine in ca. einer 
Millisekunde durch; ich habe damals gegen einen PC mit Zufallszahlen 
getestet.
Gruß,
Martin

Autor: W.S. (Gast)
Datum:
Angehängte Dateien:

Bewertung
-1 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Jederzeit gerne Kommentare und Anmerkungen.

Na denn wolle mer mal.

Ich hatte vor ein paar Jahren mir den Cordic mal näher angesehen, um mir 
über die Nützlichkeit des restlichen Drehwinkels nach getaner 
CORDIC-Arbeit klar zu werden.

Ergebnis: siehe Anlagen.
Ursprünglich mit Delphi 6 (jaja, immer noch), hab's grad eben mal nach 
XE3 umgesetzt, deswegen die Umfänglichkeit der EXE.

Also: du verschwendest in deinem Code Rechenzeit und Speicherplatz. Wie 
ich schon weiter oben geschrieben habe, braucht man nur so etwa 1/3 der 
Durchläufe. Also für 30 Bit 10 Runden CORDIC, die restlichen Stellen 
stecken im Restwinkel.

Es ist aber wirklich WICHTIG, im CORDIC nicht bedingt zu drehen, 
sondern immer und die Richtung geht immer in Richtung Abszisse. Dann 
wird jede DrehStreckung ausgeführt und jede fügt für die Zeigerlänge 
einen Faktor 1/cos(drehwinkel) hinzu. Deshalb ist die finale Zeigerlänge 
und eben auch der Restwinkel mit einem simplen Faktor wieder normierbar 
und der Restwinkel kann danach auf das Ergebnis einfach addiert werden.

Ach ja: die Tafel mit all den Konstanten hab ich mir hier verkniffen, 
aber das siehst du ja in der Quelle selber.

W.S.

Autor: Detlef _. (detlef_a)
Datum:
Angehängte Dateien:

Bewertung
2 lesenswert
nicht lesenswert
>>>>>>>
Also: du verschwendest in deinem Code Rechenzeit und Speicherplatz. Wie
ich schon weiter oben geschrieben habe, braucht man nur so etwa 1/3 der
Durchläufe.
<<<<

Nein, das stimmt nicht, die Behauptung wird auch durch nachdrückliches 
Vertreten nicht wahr.

Angehängtes Bild zeigt für Zufallszahlen in blau die Anzahl der 
Drehungen bis der Imaginärteil 0 wird wenn man den Imaginärteil immer 
positiv läßt, in rot wenn man auch ins negative dreht. Beide Zahlen sind 
ungefähr gleich, nix Drittel.

Das liegt daran, dass Zufallszahlen die gleiche mittlere Anzahl von 
'Einsen' als Binärzahl haben, ob negativ oder positiv.

Ausserdem kann ich billig drehen, das kostet in meinem Code zwei 
Additionen und zwei shifts.

math rulez!

Cheers
Detlef


clear
pluss =zeros(1,100);
minuss=zeros(1,100);
for(k=1:100000)
  ph=1;
  zz=1+i*rand(1,1);
  zz1=zz;
  nd=0;
  while(abs(angle(zz1))>1E-16)
      if(angle(zz1*(1-i*ph))>0)
        zz1= zz1*(1-i*ph);
      end;
      nd=nd+1;
      ph=ph/2;
  end;
  pluss(nd)=pluss(nd)+1;

  zz1=zz;
  nd=0;
  ph=1;
  while(abs(angle(zz1))>1E-16)
      ph;
      if(angle(zz1)>0)
        zz1= zz1*(1-i*ph);
      else
        zz1= zz1*(1+i*ph);
      end;
      nd=nd+1;
      ph=ph/2;
  end;
  minuss(nd)=minuss(nd)+1;
end;
plot(1:100,pluss,'b.-',1:100,minuss,'r.-');

Autor: W.S. (Gast)
Datum:

Bewertung
-1 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Nein, das stimmt nicht, die Behauptung wird auch durch nachdrückliches
> Vertreten nicht wahr.

Guck einfach in die Datei, um es zu sehen. Mit nur 8 Drehungen plus Rest 
kommt man auf 5 gültige Nachkommastellen des Winkels (im Gradmaß 0..359 
Grad) - und wieviele Drehungen brauchst du für diese Genauigkeit, he?

Ich habe eher den Eindruck, daß du über was ganz anderes redest als ich.

W.S.

Autor: Johann L. (gjlayde) Benutzerseite
Datum:

Bewertung
1 lesenswert
nicht lesenswert
Aaaalso...

Ein wichtiger Tipp war der Link von Alexander zu Ganssle und der dort 
verwendeten Reduktion des Arguments:

Nachdem das Argument auf [0, 1] reduziert ist, spaltet man das Intervall 
auf in 2 Intervalle [0, a] und [a, 1], und zwar so, dass per 
Addttionstheorem

atan(x) - atan(y) = atan((x - y) / (1 + xy))

das Intervall [a, 1] auf [-a, a] abgebildet wird.  Das ergibt ein 
Gleichungssystem in 2 Variablen, aus dem alle Schikanen ausfallen und 
das liefert:  y = 1/sqrt(3),  atan(y) = pi/6,  a = (1 - y) / (1 + y) ~ 
0.268.

Für dieses reduzierte Intervall genügt dann

(14549535x+29609580x^3+19801782x^5+4813380x^7+307835x^9)
/
(14549535+34459425x^2+28378350x^4+9459450x^6+1091475x^8+19845x^10)

um auf 58 Bits Genauigkeit zu kommen (auf dem PC bestimmt).  Außerdem 
hab ich die Division jetzt in Assembler, so dass Divisionen nicht mehr 
ganz so weh tun was Laufzeit angeht (vorher hatte ich aus Bequemlichkeit 
Newton-Raphson, was ca 3× langsamer war).

Laufzeit ist ca. 18000 Ticks bei 2.5 KiB Code (inclusive aller 
Konstanten und Abhängigkeiten wie Code für Division etc.)  bzw. < 450 
Bytes Code (mit Konstanten aber ohne sonstige Abhängigkeiten).

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.