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.
von Johann L. (gjlayde) Benutzerseite


Bewertung
6 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?

von minifloat (Gast)


Angehängte Dateien:

Bewertung
1 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 :)

von A. S. (achs)


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

von minifloat (Gast)


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?

von Franko S. (frank_s866)


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.

von minifloat (Gast)


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...

von Dirk B. (dirkb2)


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

von minifloat (Gast)


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.

von Egon D. (egon_d)


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.

von Johann L. (gjlayde) Benutzerseite


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.

von whitespace (Gast)


Bewertung
1 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.

von A. B. (Gast)


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.

von Joggel E. (jetztnicht)


Bewertung
-13 lesenswert
nicht lesenswert
Trollalarm !!!

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

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


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

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

von Joggel E. (jetztnicht)


Bewertung
-7 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.

von Walter T. (nicolas)


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.

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


Bewertung
3 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.

von Axel S. (a-za-z0-9)


Bewertung
1 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

von Dieter R. (drei)


Bewertung
-3 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!

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


Bewertung
4 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.

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


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.

von Dieter R. (drei)


Bewertung
-5 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.

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


Bewertung
5 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
von Alexander S. (alesi)


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.

von Alexander S. (alesi)


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

von Johann L. (gjlayde) Benutzerseite


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.

von Dieter R. (drei)


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!

von W.S. (Gast)


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.

von Egon D. (egon_d)


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?

von Helmut S. (helmuts)


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
von Percy N. (vox_bovi)


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.

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


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.

von Alexander S. (alesi)


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

von Helmut S. (helmuts)


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).

von Alexander S. (alesi)


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."

von Helmut S. (helmuts)


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.

von Johann L. (gjlayde) Benutzerseite


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?

von Alexander S. (alesi)


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
von Alexander S. (alesi)


Bewertung
0 lesenswert
nicht lesenswert

von Johann L. (gjlayde) Benutzerseite


Bewertung
2 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...

von Rainer V. (a_zip)


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

von Hamburger (Gast)


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 ..

von Rainer V. (a_zip)


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

von Detlef _. (detlef_a)


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

von Rainer V. (a_zip)


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

von Martin Cibulski (Gast)


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

von W.S. (Gast)


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.

von Detlef _. (detlef_a)


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.-');

von W.S. (Gast)


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.

von Johann L. (gjlayde) Benutzerseite


Bewertung
2 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).

von Johann L. (gjlayde) Benutzerseite


Bewertung
3 lesenswert
nicht lesenswert
Franko S. schrieb:
> Zum Trost: arcsin, arccos sind noch schlimmer.

Hab ich inzwischen auch, war garnicht sooo schlimm :-)

Man beginnt mit der Beobachtung, dass acos in 1 i.w. die Quadratwurzel 
ist.  Die Singularität kann man also einfach rausfaktorisieren:

acos(1-x) = sqrt(2x) * a(x)

wobei sich a() als rationale McLaurin-Reihe darstellen lässt. 
Konvergenzradius von a() ist 2 und a(0) = 1 (d.h. a(0) ist ungleich 0), 
also schon mal nicht schlecht.

Ist eigentlich auch klar: 1 - x^2/2 ist Asymptote an cos in (0,1), 
folglich ist sqrt(2-2x) Asymptote an acos in (1,0).

Implementiert ergibt das:
Genauigkeit: ca. 52-53 Bits
Code:        ca. 3 KiB (inclusive aller Abhängigkeiten und Konstanten)
Laufzeit:    ca. 30000 Ticks (worst case), 20000 Ticks (Median)

Quadratwurzel ist momentan Newton-Raphson, die frisst schon 10000 Ticks 
:-/

Als Referenzen hab ich avr_f64.c/h (C):
Genauigkeit:  ca. 54 Bits
Code:         ca. 6 KiB
Ticks:        80000 / 80000

sowie fp64lib (Assembler):
Genauigkeit:  ca. 47-48 Bits
Code:         ca. 4.5 KiB
Ticks:        20000 / 20000

Evtl. ergibt sich eine vorteilhaftere Näherung für a() wenn man nicht um 
0 entwicklelt sondern um ein x0 > 0 irgendwo im Intervall (0, 1/4).

: Bearbeitet durch User
von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Helmut S. schrieb:
> 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.

Johann L. schrieb:
> 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.

Im Buch "Computer Approximations" von Hart, auf das sich auch 
http://www.ganssle.com/approx.htm bezieht, gibt es u.a. 
arctan-Approximationen für 13.70 (2/3), 16.06 (3/3) bis 25.51 (5/5) 
Dezimalstellen
im Bereich [0, tan pi/12]. (N/M) bezeichnet den Grad der Polynome (in 
x^2) in Zähler und Nenner. Mit "range reduction" geht auch [0, tan pi/4] 
oder weiter. Die arctan-Approximationen basieren auf dem relativen 
Fehler.
Die angegebene "precision", z.B. 16.06 bezeichnet aber absolute 
Dezimalstellen. Das 3/3-Polynom für 16.06 Stellen liefert z.B. für x = 
0.05 einen relativen Fehler von
(0.0499583957219428-0.0499583957219427)/0.0499583957219427 = .0000000000000020
Die Werte für x = 0.05 bis 2.0 für atan_13.70 und atan_16.06 und die 
Polynomkoeffizienten (ohne Gewähr) für diese beiden Approximationen 
finden sich im Anhang. Die Fehler in atan_double.txt beziehen sich auf 
die Werte der C libmath.

von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
2 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Im Buch "Computer Approximations" von Hart, auf das sich auch
> http://www.ganssle.com/approx.htm bezieht, gibt es u.a.
> arctan-Approximationen für 13.70 (2/3), 16.06 (3/3) bis 25.51 (5/5)
> Dezimalstellen im Bereich [0, tan pi/12].

Ich hab das
p00 = 12.82297680061262709335
p01 = 16.2428078151342213662
p02 =  4.923443373635526899
p03 =   .205608723964456163
q00 = 12.82297680061262821474
q01 = 20.517133415336848884
q02 =  9.197892485655692716
q03 =  1.0

mal eingebaut, sieht ganz gut aus.  Die Genauigkeit[1] mit der neuen 
Funktion "swift" ist etwas schlechter als meine aktuelle "old" (Jeweils 
die untere Teilgrafik).  Umgerechnet in die interne Darstellung hab ich 
mit Python/Decimal, ich hoff mal ich hab da keinen Fehler gemacht.
#if defined (FOR_NUMERATOR)
// 12.82297680061262709335 + 16.2428078151342213662x + 4.923443373635526899x^2 + 0.205608723964456163x^3
F7_CONST_DEF (X, 0,  0xcd,0x2a,0xe9,0xb8,0xbf,0xf7,0x96, 3)
F7_CONST_DEF (X, 0,  0x81,0xf1,0x45,0x39,0x49,0xb7,0xb2, 4)
F7_CONST_DEF (X, 0,  0x9d,0x8c,0xd9,0x1e,0x2f,0x1e,0xb1, 2)
F7_CONST_DEF (X, 0,  0xd2,0x8b,0x17,0xe4,0xcc,0x6b,0xfd, -3)
#elif defined (FOR_DENOMINATOR)
// 12.82297680061262821474 + 20.517133415336848884x + 9.197892485655692716x^2 + x^3
F7_CONST_DEF (X, 0,  0xcd,0x2a,0xe9,0xb8,0xbf,0xf7,0x9b, 3)
F7_CONST_DEF (X, 0,  0xa4,0x23,0x16,0xd8,0x14,0x53,0x0c, 4)
F7_CONST_DEF (X, 0,  0x93,0x2a,0x91,0x4f,0xa0,0x3e,0xdf, 3)
F7_CONST_DEF (X, 0,  0x80,0x00,0x00,0x00,0x00,0x00,0x00, 0)
#else
F7_CONST_DEF (1_sqrt3,  0,  0x93,0xCD,0x3A,0x2C,0x81,0x98,0xE2,  -1)
F7_CONST_DEF (pi_6,  0,  0x86,0x0a,0x91,0xc1,0x6b,0x9b,0x2c,  -1)
#endif

Argument-Reduktion ist wie gehabt auf 2-sqrt3 ~ 0.26795.

Hast du auch die [5/5] Näherung?

[1] Wird monemtan bestimmt aus tan(atan(x)) - x, kommt also noch die 
Rundungsfehler von tan dazu.

von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Hast du auch die [5/5] Näherung?

Ich habe das Buch vor mir liegen. Ich habe jetzt auch noch die 
Koeffizienten  für atan_20.78 (4/4) und atan_25.51 (5/5) abgetippt. In C 
reicht auch "long double" auf x86 nicht mehr. Für Kontrollen hänge ich 
auch noch die mit bc und "scale=26" berechneten Werte an. Ich habe die 
abgetippten Koeffizienten einmal zur Kontrolle gelesen, aber weiterhin 
"keine Gewähr".

von Schakeline (Gast)


Bewertung
-1 lesenswert
nicht lesenswert
Hallo Jörg W.,

Jörg W. schrieb:
> Joggel E. schrieb:
>> Trollalarm !!!
>
> [ ] Du weißt, über wen du gerade sprichst.

Beitrag "wer is Gjlayde"

hat mir auch nicht weitergeholfen.

Wer ist er und was macht er?

Ich hätte auf GCC-Programmierer getippt, kann aber seinen Vornamen nicht 
auf der Contributor-Seite finden.

von Jörg W. (dl8dtl) (Moderator) Benutzerseite


Bewertung
2 lesenswert
nicht lesenswert
AVR-GCC-Programmierer

Warum Johann nicht bei deren Contributors auftaucht, musst du wohl die 
GCC-Leute fragen.

: Bearbeitet durch Moderator
von K. L. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Was ich mal gebrauchen könnte, wäre eine Implementierung, die ohne 
Real-Divisionen auskommt. Lässt sich der Nenner gfs so umformen, dass 
nur Binärdivisionen auftauchen?

von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Johann L. schrieb:
>> Hast du auch die [5/5] Näherung?
>
> Ich habe das Buch vor mir liegen. Ich habe jetzt auch noch die
> Koeffizienten  für atan_20.78 (4/4) und atan_25.51 (5/5) abgetippt.

Vielen Dank für die (Tipp)Arbeit!  Ich meinte eigentlich das [4/4], also 
5 Terme in Zähler und Nenner.  Mit dem (4/4) komme ich auf die 
Genauigkeit meiner bisherigen Näherung (11 Terme, d.h. [9/10] für atan 
bzw. [4/5] für atan(sqrt)/sqrt).  Sieht also so aus als hätte sich kein 
Tippteufel eingeschlichen :-)

Bei Gelegenheit werd ich mal nen anderen Test als f_inv(f(x)) - x 
einbauen um Genauigkeit besser beurteilen zu können.

Klaus E. schrieb:
> Was ich mal gebrauchen könnte, wäre eine Implementierung, die ohne
> Real-Divisionen auskommt. Lässt sich der Nenner gfs so umformen, dass
> nur Binärdivisionen auftauchen?

Was meinst du denn mit "Binärdivision"? Im Nenner steht ja ein Polynom.

von K. L. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Was meinst du denn mit "Binärdivision"? Im Nenner steht ja ein Polynom.

Das ist das Problem. Etwas MicroController-unfreundlich.

Die Taylor-Entwicklung erlaubt nach meiner Erinnerung die Entwicklung um 
jeden Punkt. Könnte man das nicht irgendwie umformen?

von Antoni Stolenkov (Gast)


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

Wie die Zahl Pi schon vermuten läßt, gibt es im Gegensatz zum Quadrat 
auch mathematisch keinen perfekten Kreis.

von Antoni Stolenkov (Gast)


Bewertung
-2 lesenswert
nicht lesenswert
Ja gut, vielleicht algebraisch.

von Johann L. (gjlayde) Benutzerseite


Bewertung
4 lesenswert
nicht lesenswert
Klaus E. schrieb:
> Johann L. schrieb:
>> Was meinst du denn mit "Binärdivision"? Im Nenner steht ja ein Polynom.
>
> Das ist das Problem. Etwas MicroController-unfreundlich.

Na, Division ist ja kein Hexenwerk. Für die Reduktion des Arguments 
braucht's eh Division(en).

Und man kann auch dividieren ohne zu dividieren:

https://en.wikipedia.org/wiki/Newton-Raphson_division#Newton–Raphson_division

> Die Taylor-Entwicklung erlaubt nach meiner Erinnerung die Entwicklung um
> jeden Punkt. Könnte man das nicht irgendwie umformen?

Man nimmt deshalb rationale Approximationen weil die i.d.R. besser 
funktionieren als Taylor bzw. weniger Terme brauchen. Du kannst auch 
Taylor nehmen und um 0 entwickeln: x - x^3/3 + x^5/5 - x^7/7 ... 
Konvergenzradius ist 1, d.h. die Summanden werden nur kleiner weil man 
höhere Potenzen vo x nimmt; die Koeffizienten leisten praktisch keinen 
Beitrag.  Hat man das Intervall wie oben beschrieben auf 2-sqrt3 ~ 1/4 
reduziert, dann gewinnt man an der Intervallgrenze ca. 4 Bit pro Term. 
Für double-Genauigkeit (53Bits) brauch man also mindestens ca. 14 Terme: 
bis -x^27/27.

Du kannst auch um andere Punkte entwickeln, wobei du aber Symmetrie 
verlierst, d.h. es ist nicht jeder 2. Term der Entwicklung 0 wie bei der 
Entwicklung um 0 (wegen atan(-x) = -atan(x)).  atan hat Singularitäten 
bei i und -i, der Konvergenzradius bei Entwicklung um 1 ist also sqrt2. 
Falls du damit Werte in [0.5, 1.5] nähern willst, dann gewinnst du pro 
Term ca. 1.5 Bit.  Für double-Genauigkeit braucht's alo schon 36 Terme 
an der Intervallgrenze...  und da ist Rundung und Auslöschung noch nicht 
drinne.

Sinnvoll ist um Unendlich zu entwickeln — das macht man effektiv per 
atan(x) + atan(1/x) = sgn(x)*pi/2  wenn man Argumente |x|>1 auf |x|<1 
reduziert.

: Bearbeitet durch User
von K. L. (Gast)


Bewertung
-3 lesenswert
nicht lesenswert
aha(?)

von Antoni Stolenkov (Gast)


Bewertung
-2 lesenswert
nicht lesenswert
Klaus E. schrieb:
> aha(?)

Die Perfektion der Mathematik existiert nur für den Betrachter außerhalb 
der mathematischen Welt. Innerhalb der mathematischen Welt existiert 
kein Betrachter.

Der Betrachter eines Gegenstandes existiert nur innerhalb dieser Welt. 
Außerhalb dieser Welt existiert kein Betrachter.

Zur Welt der Rechenmaschinen fällt mir jetzt grade nichts Anständiges 
ein.

von Antoni Stolenkov (Gast)


Bewertung
-3 lesenswert
nicht lesenswert
Antoni Stolenkov schrieb:
> Außerhalb dieser Welt existiert kein Betrachter.

Jetzt aber keine falschen Schlüsse ziehen. Das Universum ist scheiss 
groß.

von Sven B. (scummos)


Bewertung
0 lesenswert
nicht lesenswert
Blöde Frage, warum funktioniert denn an mehreren Punkten 
taylorentwickeln und dann den passenden auswählen nicht?

von S. Z. (ceving)


Bewertung
2 lesenswert
nicht lesenswert

: Bearbeitet durch User
von Christoph db1uq K. (christoph_kessler)


Bewertung
2 lesenswert
nicht lesenswert
Die Numerical Recipies wurden noch nicht erwähnt, ich dachte, das ist 
ein Standardwerk:
https://de.wikipedia.org/wiki/Numerical_Recipes
Zum Lesen muss man allerdings eine Software installieren, die ich noch 
nicht kannte. "Empanel" oder "OpenFile"

Damit
http://numerical.recipes/oldverswitcher.html
bin ich noch auf diese ausführliche Formelsammlung gestoßen:
http://www.nrbook.com/abramowitz_and_stegun_html/page_81.htm
ein paar Näherungsreihen für arctan

: Bearbeitet durch User
von Walter T. (nicolas)


Bewertung
1 lesenswert
nicht lesenswert
Christoph db1uq K. schrieb:
> Die Numerical Recipies wurden noch nicht erwähnt, ich dachte, das ist
> ein Standardwerk:

Wobei bei den "Numerical Recipes" die Nutzungsrechte der Rezepte nicht 
unproblematisch sind.

von Johann L. (gjlayde) Benutzerseite


Bewertung
4 lesenswert
nicht lesenswert
Sven B. schrieb:
> Blöde Frage, warum funktioniert denn an mehreren Punkten
> taylorentwickeln und dann den passenden auswählen nicht?

Zunächst mal weil, egal wo man entwickelt, der Konvergenzradius immer 
endlich ist.  Heißt: Um ganz R zu überdecken bräuchte man oo viele 
Taylor-Reihen.  Außerdem verliert man die Symmetrie wenn man nicht um 
(0,0) entwickelt, d.h. in einer Umgebung von 0 hat man doppelt zu viele 
Terme auszuwerten.

Wenn man sich auf der reellen Achse von (0,0) entfernt wird der 
Konvergenzradius zwar größer (der Radius ist sqrt(1+|x0|) wenn x0 der 
Entwicklungspunkt ist), aber man muss eben nen ganzen Zoo von Reihen 
speichern.

Wenn die geforderte Präzision nicht allzu groß ist, und man z.B. mit 
linearen oder kubischen Teilstücken hinkommt, entspricht das im 
Endeffekt einer linearen Approximation bzw. einer mit kubischen Splines, 
wobei Splines und lineare Approximation den Vorteil haben, dass die 
stetig sind und im Idealfall den Fehler im Teilintervall minimieren (was 
Taylor nicht macht).

Hier noch die Reihe um 1 (Konvergenzradius ist wie gesagt sqrt2):

atan(x+1) - atan(1) = 1/2x - 1/4x^2 + 1/12x^3 - 1/40x^5 + 1/48x^6 - 
1/112x^7 + 1/288x^9 - 1/320x^10 + 1/704x^11 - 1/1664x^13 + 1/1792x^14 - 
1/3840x^15 + 1/8704x^17 - 1/9216x^18 + 1/19456x^19 ...

Für x=-1 ergibt das -atan(1) = -pi/4.  Wertet man bis zu der gezeigten 
Ordnung aus, ergibt sich:

-130023561/165541376 ~ -pi/4 - 4.6E-5

das sind 14.5 Bits. Warum auch immer verschwindet jeder 4. Term, d.h. 
man braucht nur ca. 3N/4 Terme für Ordnung N.

: Bearbeitet durch User
von Sven B. (scummos)


Bewertung
0 lesenswert
nicht lesenswert
Danke für die Erklärung!

von Alexander S. (alesi)


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.

Erstmal Danke für den Literaturhinweis. Wie die Beiträge weiter oben 
zeigen, sind die Approximationen in dem Buch schon noch "state of the 
art". Mit heutigen Computer-Algebra-Systemen (CAS) kann man sich die 
Koeffizienten, passend zu den Anforderungen, auch berechnen lassen. Für 
das CAS Maple ist das z.B. hier beschrieben:
Maplesoft Online Help numapprox minimax rational approximation
https://www.maplesoft.com/support/help/Maple/view.aspx?path=numapprox/minimax
Das Abtippen der Koeffizienten, wie ich das gemacht habe, ist eigentlich 
nicht mehr zeitgemäß. In den freien CAS, wie Axiom oder Maxima habe ich 
aber noch nicht nach Funktionen wie "minimax" geschaut.

von Helmut S. (helmuts)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hier mal die Lösung mit 6 Koefitienten im Zähler und im Nenner. Der Plot 
zeigt die Differenz zwischen dem approximierten Wert und der normalen 
atan()-Funktion von Cpp. die normale atan()-Funktion hat vermutlich nur 
"double" und kein "long double".

/*    atan(x) x-range 0...1
      y = x*P(x^2)/Q(x^2)
      Approximations, John F. Hart, Wiley, 1968
*/

Es ist schon verblüffend, dass es bereits 1968 so genaue Approximationen 
gab.

: Bearbeitet durch User
von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Helmut S. schrieb:
> die normale atan()-Funktion hat vermutlich nur
> "double" und kein "long double"

Hallo, daher habe ich Dein Programm mal so modifiziert, dass es die mit 
bc (bc - An arbitrary precision calculator language) auf 22 
Dezimalstellen berechneten atan-Werte in ein array einliest und dann die 
Differenz zur Approximation mit Polynomgraden (in x^2) 6/6, d.h. 7/7 
Koeffizienten für das Intervall [0; tan pi/4] mit precision 17.24 im 
Buch von Hart ausgibt. Für die gewählten x-Werte stimmt die precision 
ganz gut.

von Helmut S. (helmuts)


Bewertung
0 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Helmut S. schrieb:
>> die normale atan()-Funktion hat vermutlich nur
>> "double" und kein "long double"
>
> Hallo, daher habe ich Dein Programm mal so modifiziert, dass es die mit
> bc (bc - An arbitrary precision calculator language) auf 22
> Dezimalstellen berechneten atan-Werte in ein array einliest und dann die
> Differenz zur Approximation mit Polynomgraden (in x^2) 6/6, d.h. 7/7
> Koeffizienten für das Intervall [0; tan pi/4] mit precision 17.24 im
> Buch von Hart ausgibt. Für die gewählten x-Werte stimmt die precision
> ganz gut.

Hallo Alexander,
vielen Dank für die Tabelle mit den genauen Werten.
Der Fehlerplot sieht jetzt richtig gut aus. Besonders gut finde ich, 
dass der absolute Fehler ungefähr proportional mit dem Wert von x 
steigt. Damit ist auch der relative Fehler nicht viel schlechter. Hut ab 
vor John F. Hart.

von Helmut S. (helmuts)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Alexander,
ich habe jetzt deinen Code genommen und die Approximation atan1(x) mit 
der atan(x)-Standard-Funktion verglichen. Die Approximation ist wirklich 
Spitze verglichn mit der eingebauten atan()-Funktion.

PS:
#include <stdio.h>
Diese Zeile hat mich eine halbe Stunde gekostet.
Mit stdio.h kamen nur xxx-317 Werte heraus, egal ob printf() oder 
fprintf(). Irgendwann fiel mir auf, dass ich <iostream> verwendet hatte. 
Nach dieser Änderung hat alles wieder funktioniert - 
main3_long_double.cpp

//#include <stdio.h>
#include <iostream>

Beitrag #5990454 wurde von einem Moderator gelöscht.
von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Helmut S. schrieb:
> Hallo Alexander,
> ich habe jetzt deinen Code genommen und die Approximation atan1(x) mit
> der atan(x)-Standard-Funktion verglichen. Die Approximation ist wirklich
> Spitze verglichn mit der eingebauten atan()-Funktion.

Kannst du auch mal den relativen Fehler plotten? Da hat die rote Kurve 
wohl konstante Amplitude.

Hat eigentlich jemand ein CA System verfügbar, mit dem man solche 
rationalen minimax Approximationen bestimmen kann?  Die Numerik ist ja 
nicht ganz trivial, siehe z.B. https://arxiv.org/pdf/1705.10132

Konkret würde mich arcsin(sqrt(x/2)) / sqrt(x/2) im Intervall [0, 1/2] 
mit relativem Fehler -55 Bits oder kleiner interessieren.

von Helmut S. (helmuts)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Alexander,
der screenshot zeigt den relativen Fehler.
Die Frage ist jetzt ob der konstante Anstieg eventuell mit der endlichen 
Genauigkeit von long double (80bit) zusammenhängt.

: Bearbeitet durch User
von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Kannst du auch mal den relativen Fehler plotten? Da hat die rote Kurve
> wohl konstante Amplitude.

Helmut S. schrieb:
> Hallo Alexander,
> der screenshot zeigt den relativen Fehler.

Helmut war schneller. Trotzdem hier noch meine Variante (siehe 
Diagramm).
        if (i>0) {
              e = (y-z)/z;
              fprintf(file2,"%28.20Le %28.20Le\n",x,e);
        }

Helmut S. schrieb:
> #include <stdio.h>
> Diese Zeile hat mich eine halbe Stunde gekostet.

Entschuldigung. Aber eigentlich ist der Code komplett C und nicht C++.

Johann L. schrieb:
> Hat eigentlich jemand ein CA System verfügbar, mit dem man solche
> rationalen minimax Approximationen bestimmen kann?

Ich habe unter Linux Axiom, Maxima und Octave u.v.m. Habe aber noch 
nicht geschaut, ob es dafür schon die Routinen gibt.

Etwas offtopic: Ich überlege gerade mir einen RaspberyPi4 zu holen. 
Dafür gibt es Mathematica für privaten Gebrauch kostenlos.

Hier noch ein Link auf "Function Approximation" bei Wolfram:
https://reference.wolfram.com/language/FunctionApproximations/tutorial/FunctionApproximations.html

: Bearbeitet durch User
von Alexander S. (alesi)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Konkret würde mich arcsin(sqrt(x/2)) / sqrt(x/2) im Intervall [0, 1/2]
> mit relativem Fehler -55 Bits oder kleiner interessieren.

Hallo Johann, ich habe mir mal erlaubt Deine Frage auf dem Matheplaneten 
www.matheplanet.de im Unterforum "Mathematica" zu stellen. Der Benutzer 
HyperG hat mit Mathematica dieses Polynom erzeugt:
(0.999999999999985 - 0.795192955982834 x + 0.180219557964123 x^2 - 0.0103245962130046 x^3) /
(1 - 0.87852628931959 x + 0.234680082201635 x^2 - 0.0189892607523874 x^3 + 0.000185818094781771 x^4)
Wieder abgetipt, d.h. ohne Gewähr. Die Genauigkeit ist ca. 3E-14.

https://www.matheplanet.de/matheplanet/nuke/html/viewtopic.php?topic=243737&post_id=1775475

von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
2 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Johann L. schrieb:
>> Konkret würde mich arcsin(sqrt(x/2)) / sqrt(x/2) im Intervall [0, 1/2]
>> mit relativem Fehler -55 Bits oder kleiner interessieren.
>
> Hallo Johann, ich habe mir mal erlaubt Deine Frage auf dem Matheplaneten
> www.matheplanet.de im Unterforum "Mathematica" zu stellen. Der Benutzer
> HyperG hat mit Mathematica dieses Polynom erzeugt:
>
> (0.999999999999985 - 0.795192955982834 x + 0.180219557964123 x^2 - 
> 0.0103245962130046 x^3) /
> (1 - 0.87852628931959 x + 0.234680082201635 x^2 - 0.0189892607523874 x^3 
> + 0.000185818094781771 x^4)
> 
> Wieder abgetipt, d.h. ohne Gewähr. Die Genauigkeit ist ca. 3E-14.
>
> 
https://www.matheplanet.de/matheplanet/nuke/html/viewtopic.php?topic=243737&post_id=1775475

Hi, vielen Dank für dein Engagement!

Inzwischen hab ich meine eigene Approximation zusammengeklöppelt weil 
mich das Rumgesuche genervt hat...

Just bestimmt und ausgetestet ist folgende rationale [5/4] MiniMax für 
o.g. Funktion in [0, 1/2]:
p(x) =
+ 0.99999999999999999664796559276984386465
- 1.0349496317294863228413401035536833883x
+ 0.35268117737125871147808800094621749993x^2
- 0.043287377507331303209307905646261090153x^3
+ 0.0012535143820048849783737210746210873110x^4
+ 0.0000084478020529822158900907931511410767987x^5

q(x) =
+ 1
- 1.1182829650628212672605440708760882999x
+ 0.42712142445992813936913160697952988905x^2
- 0.063493381095966407972729080832358862541x^3
+ 0.0028776495966594259144201280184595390407x^4

Relativer Fehler ist besser als 7E-18.

Damit sind asin und acos besser als mit der [7/7] Padé.

Math Rules!

: Bearbeitet durch User
von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
2 lesenswert
nicht lesenswert
...zurück zu atan.

Mit meiner MiniMax hab ich folgende [3,4] für a(x) = atan(sqrt(x)) / 
sqrt(x) bestimmt:
p(x) = 
+ 0.999999999999999998080178351225003952632
+ 1.51597040589722809277543885223133087789x
+ 0.636928974763539784860899545005247736093x^2
+ 0.0638944455799886571709448345524306512048x^3

q(x) = 
+1
+ 1.84930373923056200945453682584178320362x
+ 1.05336355450697082895016644607427716580x^2
+ 0.188012025422931152047197803304030906006x^3
+ 0.00507310235929401206762490897042543192106x^4

Damit ergibt sich folgender Algorithmus zur Berechnung von atan:

Argument-Reduktion

1) auf [0,oo]
2) auf [0,1]
3) auf [-w,w] mit w = 2 - sqrt3 ~ 0.268

wie ausgeführt in Beitrag "Re: Näherung für ArcusTangens?"

In [-w^2,w^2] verwendet man dann g(x) = p(x) / q(x) als Näherung für 
a(x) und bestimmt atan gemäß

atan(x) = x·a(x^2) ~ x·g(x^2); x in [-w,w]

d.h. a() wird genähert in [-w^2,w^2] mit w^2 = 7 - 4sqrt3 ~ 0.072.

Der Plot im Anhang zeigt den relativen Fehler von x·g(x^2) zu atan(x) in 
[-w,w].  Der Graph ist nicht komplett symmetrisch weil rechts noch für 
paar Werte > w geplottet wurde.

Der Fehler ist (fast) genau wie erwartet. (Bei 0 und w sind keine Maxima 
weil a() auf einem minimal größeren Intervall als [-w^2,w^2] genähert 
wurde, d.h. bei 0 und w hat der Fehler von g() kein Maximum.)

p.s: Der linke Plot mit x in [0..0.3] kann gelöscht werden.

: Bearbeitet durch User
von Jens G. (jensg)


Bewertung
-3 lesenswert
nicht lesenswert
Zwischenfrage: Weshalb werden hier Näherungsreihen so wesentlich 
diskutiert. Ich würde CORDIC verwenden um an das Ziel zu kommen. Ein zu 
hoher RAM-Verbrauch? .. Eine zu hohe Rechenzeit?

von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Christoph db1uq K. schrieb:
> Die Numerical Recipies wurden noch nicht erwähnt, ich dachte, das ist
> ein Standardwerk:
> https://de.wikipedia.org/wiki/Numerical_Recipes
> Zum Lesen muss man allerdings eine Software installieren, die ich noch
> nicht kannte. "Empanel" oder "OpenFile"

Hallo, auch wenn die ursprüngliche Frage bereits beantwortet ist, hier 
ein Kommentar zu den Numerical Recipes.
Die 2. Aufl. der Numerical Recipes in C kann man sich als flash-Datei 
hier anschauen:  http://numerical.recipes/oldverswitcher.html
In dem Buch gibt es im Abschnitt 5.13 Rational Chebyshev Approximation 
auf S. 206 die Routine ratlsq, die f(x) mit R(x) = 
(p0+p1*x+...)/(1+q1*x+...)
approximiert.
Für die Funktion f(x) = asin(sqrt(x/2)) / sqrt(x/2) erhalte ich mit Grad 
3 / Grad 4 und Intervall [0:0.5] eine Näherung mit einer Polstelle bei x 
= 0.3.
Bei Einschränkung auf das Intervall [0:0.25] kommt eine ganz brauchbare 
Näherung heraus.
 ratlsq iteration=  1 max error=  4.821e-06
 ratlsq iteration=  2 max error=  2.017e-06
 ratlsq iteration=  3 max error=  7.303e-04
 ratlsq iteration=  4 max error=  3.086e-04
 ratlsq iteration=  5 max error=  1.670e-02
[0.000001 : 0.500000]
 0   1.000038095717036946
 1  202244301.401464492082595825
 2  -735806100.350590467453002930
 3  186583373.451535642147064209
 4  202244323.917425632476806641
 5  -752661211.566265940666198730
 6  245531641.993888765573501587
 7  -7564721.532575584016740322

 ratlsq iteration=  1 max error=  2.618e-08
 ratlsq iteration=  2 max error=  4.140e-10
 ratlsq iteration=  3 max error=  1.186e-07
 ratlsq iteration=  4 max error=  3.762e-08
 ratlsq iteration=  5 max error=  1.565e-05
[0.000001 : 0.250000]
 0   1.000001887400718825
 1  8925421822.914939880371093750
 2  369595154.398652076721191406
 3  -870760035.829437494277954102
 4  8925421821.023189544677734375
 5  -374189577.714669704437255859
 6  -1006942402.749973297119140625
 7  41261209.496030598878860474
Siehe die beiden angehängten Diagramme.

Für den atan in [0:pi/12] habe ich keine brauchbare Näherung erhalten.
sin scheint auch ganz gut zu gehen.

von Helmut S. (helmuts)


Bewertung
0 lesenswert
nicht lesenswert
Jens G. schrieb:
> Zwischenfrage: Weshalb werden hier Näherungsreihen so wesentlich
> diskutiert. Ich würde CORDIC verwenden um an das Ziel zu kommen. Ein zu
> hoher RAM-Verbrauch? .. Eine zu hohe Rechenzeit?

Wenn der Prozessor einen Hardwaremultiplizierer mit großer Bitbreite 
hat, dann ist man mit Polynomen schneller als mit CORDIC.

von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
3 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Helmut S. schrieb:
>> die normale atan()-Funktion hat vermutlich nur
>> "double" und kein "long double"
>
> Hallo, daher habe ich Dein Programm mal so modifiziert, dass es die mit
> bc (bc - An arbitrary precision calculator language) auf 22
> Dezimalstellen berechneten atan-Werte in ein array einliest und dann die
> Differenz zur Approximation mit Polynomgraden (in x^2) 6/6, d.h. 7/7
> Koeffizienten für das Intervall [0; tan pi/4] mit precision 17.24 im
> Buch von Hart ausgibt. Für die gewählten x-Werte stimmt die precision
> ganz gut.

"ganz gut"? Fast ideal würde ich sagen!  Hab's mal nachgerechnet, 
relativer Fehler siehe Anhang.  Der größte Max-Abstand ist 5.809E-18, 
der kleinste 5.413E-18.  Deren Quotient ist 1.073 und erreicht damit 
fast das theoretische Minimum von 1.

Irgendwo muss sich bei dir nen Fehlerteufel eingeschlichen haben; so 
ansteigende Fehlerkurven sind ein Indiz dafür, dass die interne 
Genauigkeit dem Ende zugeht.

Zur Referenz hier nochmal die Polynome (Style=Gnuplot):
p(x) = 1209.7470017580907217240715 + 3031.0745956115083044212807*x + 2761.7198246138834959053784*x**2 + 1114.1290728455183546172942*x**3 + 192.5792014481559613474286*x**4 + 11.3221594116764655236245*x**5 + 0.09762721591717633036983*x**6
q(x) = 1209.7470017580907287514197 + 3434.3235961975351716547069*x + 3664.5449563283749893504796*x**2 + 1821.6003392918464941509225*x**3 + 423.0716464809047804524206*x**4 + 39.9178842486537981501999*x**5 + x**6
und als korrekt gerundete IEEE double (Style=C, little Endian, 11-Bit 
Exponent, 53-Bit Mantisse):
uint64_t p[] = { 0x4092e6fcee076437, 0x40b7ae2631655fb5, 0x40b593708cda0ef9, 0x409168842bac0936, 0x40781288d179b406, 0x4036a4f214127f80, 0x3fb8fe18e3905564 };
uint64_t q[] = { 0x4092e6fcee076437, 0x40bad4a5ae669b0c, 0x40bca1170484103d, 0x409c7666bf57e3f0, 0x407a712576c7c5c1, 0x4053f57d3b26bda5, 0x3ff0000000000000 };


Jens G. schrieb:
> Zwischenfrage: Weshalb werden hier Näherungsreihen so wesentlich
> diskutiert. Ich würde CORDIC verwenden um an das Ziel zu kommen. Ein zu
> hoher RAM-Verbrauch? .. Eine zu hohe Rechenzeit?

Polynome bzw. rationale Funktionen sind eben ziemlich allgemein, sie 
können auf unterschiedliche Anforderungen angepasst werden 
(darzustellende Funktion, Intervall, Genauigkeit), und sie sind einfach 
auszuwerten.

Momentan ist meine auf rationalen Approximationen basierende 
Implementierung aber auch besser als die Konkurrenten, als da wären 
avr-f64 und fp64lib von Beitrag "64 Bit float Emulator in C, IEEE754 compatibel"

Von der Genauigkeit bin ich bei asin, atan in etwa wie avr-f64 (Cordic) 
und besser als fp64lib (ebenfalls Cordic).  Siehe untere Zeile von 
(meine Implementierung hat Labels f7) 
https://www.mikrocontroller.net/attachment/430317/test-f7.png

Und Codegröße ist auch besser, siehe 
https://www.mikrocontroller.net/attachment/430869/size.png

Laufzeit meiner Implementierung (C+Asm) ist etwas schlechter als fp64lib 
(Asm).

: Bearbeitet durch User
von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Irgendwo muss sich bei dir nen Fehlerteufel eingeschlichen haben; so
> ansteigende Fehlerkurven sind ein Indiz dafür, dass die interne
> Genauigkeit dem Ende zugeht.

Hallo Johann,
das Diagramm atan3.png im Beitrag 
Beitrag "Re: Näherung für ArcusTangens?" mit dem relativen 
Fehler für die atan-Näherung auf 17.24 Dezimalstellen zeigt den Fehler 
der Polynomberechnung in C mit long double, d.h. nur 18 bis 19 
Dezimalstellen. Die Referenz für die Fehlerberechnung wurde mit bc auf 
22 Stellen genau berechnet.

Das der Anstieg mit der begrenzten Genauigkeit von long double 
zusammenhängt, hat Helmut bereits vermutet.
Helmut S. schrieb:
> Die Frage ist jetzt ob der konstante Anstieg eventuell mit der endlichen
> Genauigkeit von long double (80bit) zusammenhängt.

Wenn ich das rationale Näherungspolynom (Grad 6 / Grad 6) auch in bc mit 
26 Dezimalstellen Genauigkeit berechne, erhalte ich den gleichen 
relativen Fehler wie Du. Siehe atan_17_24.png

von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Hallo Johann,
> das Diagramm atan3.png im Beitrag
> Beitrag "Re: Näherung für ArcusTangens?" mit dem relativen
> Fehler für die atan-Näherung auf 17.24 Dezimalstellen zeigt den Fehler
> der Polynomberechnung in C mit long double, d.h. nur 18 bis 19
> Dezimalstellen. Die Referenz für die Fehlerberechnung wurde mit bc auf
> 22 Stellen genau berechnet.

Ah, ok.  Kann sogar die Kurve nachvollziehen, dazu genügt es, die 
Koeffizienten der Polynome auf float(Python) zu casten.  Gerechnet wird 
mit höherer Präzision.
>>> import sys
>>> print (sys.float_info.mant_dig)
53

float(Python) entspricht also dem üblichen double(C/C++).

> Das der Anstieg mit der begrenzten Genauigkeit von long double
> zusammenhängt, hat Helmut bereits vermutet.
> Helmut S. schrieb:
>> Die Frage ist jetzt ob der konstante Anstieg eventuell mit der endlichen
>> Genauigkeit von long double (80bit) zusammenhängt.

Sieht eher danach aus, dass da irgendwie (teilweise?) mit double 
gerechnet wurde, ansonsten kann ich mir nicht vorstellen, dass ich die 
Kurve so gut mit 53-Bit Mantisse nachvollziehen kann.

> Wenn ich das rationale Näherungspolynom (Grad 6 / Grad 6) auch in bc mit
> 26 Dezimalstellen Genauigkeit berechne, erhalte ich den gleichen
> relativen Fehler wie Du. Siehe atan_17_24.png

Gibt's da eigentlich nen Hinweis, wie die Näherungen bestimmt wurden? 
Mit rationalem Remez kann ich atan in [0,1] annähern per 5/4 oder 4/5, 
bei 5/5 gibt es kaum noch Konvergenz, und für 6/6 ist das Verfahren so 
instabil, dass es nichtmehr funktioniert.

Eine Möglichkeit könnte Carathéodory-Féjer sein, d.h. Approximation auf 
dem Einheitskreis in C und dann Teil des Randes auf [a,b] in R abbilden. 
Alles sehr interessant, aber meine ursprüngliche Frage ist ja durch die 
tollen Beiträge in diesem Thread zu meiner Zufriedenheit beantwortet. 
Und rational Remez zu implementieren war einfach, die meiste Arbeit war 
QR-Zerlegung.

von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Für den atan in [0:pi/12] habe ich keine brauchbare Näherung erhalten.

Wenn man die Symmetrie von atan berücksichtigt und atan(sqrt(x)) / 
sqrt(x) approximiert, erhält man auch mit ratlsq aus den Numerical 
Recipes "brauchbare" Näherungen für den
atan x = x * (p0 + p1*x^2 + p2*x^4 + ...) / (q0 + q1*x^2 + q2*x^4 + ...)
Z.B. mit Grad 4 / Grad 4 in [0:tan pi/12]
 ratlsq iteration=  1 max error=  3.617e-06
 ratlsq iteration=  2 max error=  6.183e-07
 ratlsq iteration=  3 max error=  1.008e-03
 ratlsq iteration=  4 max error=  1.937e-05
 ratlsq iteration=  5 max error=  1.340e-01
[0.000001 : 0.518000]
 0   1.000013216745146760
 1  132475496.082246452569961548
 2  -559127738.893114686012268066
 3  -455648453.819097280502319336
 4  -33228001.094260059297084808
 5  132475502.845357134938240051
 6  -514969853.801637709140777588
 7  -653785760.688392877578735352
 8  -129374889.419161975383758545
Der absolute Fehler in [0:tan pi/12] ist kleiner 5e-09, der relative 
eine Größenordnung höher. Im Bild ist der absolute Fehler gezeigt.

von Alexander S. (alesi)


Bewertung
1 lesenswert
nicht lesenswert
Johann L. schrieb:
> Gibt's da eigentlich nen Hinweis, wie die Näherungen bestimmt wurden?

John F. Hart gibt in seinem Buch "Computer Approximations" nicht für 
jede der zahlreichen Näherungen Quellen an. Am Anfang vom Abschnitt 6.5 
zu arcsin und arctan schreibt er: "Coefficients for expansions of 
arcsin(x) and arctan(x) in a series of Chebyshev polynomials are given 
to 20D by Clenshaw [1962]; Vionnet [1960] gives 24D values.
Clenshaw [1962] "Chebyshev series for mathematical functions," National 
Physical Laboratory Mathematical tables, 5.
https://www.worldcat.org/title/chebyshev-series-for-mathematical-functions/oclc/60945199
Vionnet [1959] Approximation de Tchebycheff d'ordre des functions sin 
x,x-1 sin x, cos x, et exp x.

von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Alexander S. schrieb:
>> Für den atan in [0:pi/12] habe ich keine brauchbare Näherung erhalten.

Soweit nachvollziehbar; mit Nullstelle im Approximationsintervall 
funktioniert Optimieren auf relativen Fehler i.d.R. nicht.

> Wenn man die Symmetrie von atan berücksichtigt und atan(sqrt(x)) /
> sqrt(x) approximiert, erhält man auch mit ratlsq aus den Numerical
> Recipes "brauchbare" Näherungen für den
> [...]
> Z.B. mit Grad 4 / Grad 4 in [0:tan pi/12]
> [code]
>  ratlsq iteration=  1 max error=  3.617e-06
>  ratlsq iteration=  2 max error=  6.183e-07
>  ratlsq iteration=  3 max error=  1.008e-03
>  ratlsq iteration=  4 max error=  1.937e-05
>  ratlsq iteration=  5 max error=  1.340e-01

Mit brauchbar in Gänsefüßchen :-/  Problem ist wohl, einen Startpunkt 
zu finden, von dem aus die Iteration konvergiert.

Übrigens reicht hier bis tan²pi/12 zu approximieren: g(x) = 
atan(sqrt(x)) / sqrt(x) wird bei Reduktion des atan-Arguments auf [-z,z] 
mit z = tan pi/12 = 2-sqrt3 nur im Intervall [0,z²] mit z² < 0.0718 
ausgewertet.

In diesem Fall hat ein [3/3] Proximum Abstand < 8.75E-17 zu g. (Höhere 
Grade schafft meine Implementierung auch nicht, eben wegen des Problems 
einen geeigneten Startpunkt für die Iteration zu finden.)

von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
3 lesenswert
nicht lesenswert
Zum Abschluss noch meine Näherung für asin und acos:
p(x) =
+ 0.9999999999999999944249107313502758620
- 1.0352340338921976278427312087167692142 x
+ 0.3529020623298151981342259189772057401 x^2
- 0.0433348317064168570561235180136569466 x^3
+ 0.0012557428614630796315205218507940285 x^4
+ 0.0000084705471128435769021718764878041 x^5

q(x) =
+ 1
- 1.1185673672255329236623716486696411533 x
+ 0.4273660095987244885409833401675833351 x^2
- 0.0635558848496317165994214838980137828 x^3
+ 0.0028820878185134035637440105959294542 x^4

a(x) = p(x) / q(x) in [0, 1/2]

x in [0, 1/2]  =>  Berechne asin(x) = x·a(2x²)
x in [1/2, 1]  =>  Berechne acos(x) = sqrt(2-2x)·a(1-x)

alle anderen Fälle von asin / acos lassen sich so auf diese beiden Fälle 
zurückführen, dass nur Vorzeichenwechel und/oder Addition eines 
ganzzahligen Vielfaches von pi/2 benötigt werden.

Diese Formel gefällt mir besser als das, was ich so in der Literatur 
gefunden habe, denn dort war die Wurzel immer im Argument von asin bzw. 
acos und oft sogar im Nenner.

Relative Fehler siehe Grafiken.

: Bearbeitet durch User
von S. R. (svenska)


Bewertung
0 lesenswert
nicht lesenswert
Super Arbeit, danke!
Für solche Lichtblicke liebe ich dieses Forum. :-)

von Alexander S. (alesi)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Kann sogar die Kurve nachvollziehen, dazu genügt es, die
> Koeffizienten der Polynome auf float(Python) zu casten.  Gerechnet wird
> mit höherer Präzision.

Hallo Johann, dazu passt dieser Artikel, insbesondere Fig. 2, gut:
Frithjof Blomquist, "OPTIMIZATIONOF RATIONAL APPROXIMATIONS BY CONTINUED 
FRACTIONS", SerdicaJ. Computing 1 (2007), 433{442
http://sci-gems.math.bas.bg:8080/jspui/bitstream/10525/355/1/sjc033-vol1-num4-2007.pdf

von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
2 lesenswert
nicht lesenswert
Johann L. schrieb:
> Hat eigentlich jemand ein CA System verfügbar, mit dem man solche
> rationalen minimax Approximationen bestimmen kann?  Die Numerik ist ja
> nicht ganz trivial, siehe z.B. https://arxiv.org/pdf/1705.10132

Johann L. schrieb:
> ...zurück zu atan.
>
> Mit meiner MiniMax hab ich folgende [3,4] für a(x) = atan(sqrt(x)) /
> sqrt(x) bestimmt:
> p(x) =
> + 0.999999999999999998080178351225003952632
> + 1.51597040589722809277543885223133087789x
> + 0.636928974763539784860899545005247736093x^2
> + 0.0638944455799886571709448345524306512048x^3
>
> q(x) =
> +1
> + 1.84930373923056200945453682584178320362x
> + 1.05336355450697082895016644607427716580x^2
> + 0.188012025422931152047197803304030906006x^3
> + 0.00507310235929401206762490897042543192106x^4
>
> Damit ergibt sich folgender Algorithmus zur Berechnung von atan:
>
> Argument-Reduktion
>
> 1) auf [0,oo]
> 2) auf [0,1]
> 3) auf [-w,w] mit w = 2 - sqrt3 ~ 0.268
>
> wie ausgeführt in Beitrag "Re: Näherung für ArcusTangens?"
>
> In [-w^2,w^2] verwendet man dann g(x) = p(x) / q(x) als Näherung für
> a(x) und bestimmt atan gemäß
>
> atan(x) = x·a(x^2) ~ x·g(x^2); x in [-w,w]

Johann L. schrieb:
> Mit rationalem Remez kann ich atan in [0,1] annähern per 5/4 oder 4/5,
> bei 5/5 gibt es kaum noch Konvergenz, und für 6/6 ist das Verfahren so
> instabil, dass es nichtmehr funktioniert.

Hallo Johann,

falls es noch von Interesse ist: ich habe den atan jetzt mit Mathematica 
in [0;0.268] mit Grad 3 durch Grad 4 und mit Grad 4 durch Grad 5 
approximiert. Erst habe ich die Koeffizienten für ftmp(x) = 
atan(sqrt(x))/sqrt(x) bestimmt und dann atan durch atapp(x) = 
x*ftmp(x^2) approximiert.
 MiniMaxApproximation[ArcTan[Sqrt[x]]/Sqrt[x], {x, {0.000001, 0.0718}, 3, 4}, WorkingPrecision -> 30]
(0.999999999999999999618764963077 + 
   1.51583745397465997042783323587 x + 
   0.636800164486774455102278527599 x^2 + 
   0.0638720393258589120865977766757 x^3)/(1 + 
   1.84917078730799260207599292367 x + 
   1.05319042692298370820594274876 x^2 + 
   0.187958500338027650749845385678 x^3 + 
   0.00507089366452151377474182193301 x^4)

MiniMaxApproximation[ArcTan[Sqrt[x]]/Sqrt[x], {x, {0.000001, 0.0718}, 4, 5},WorkingPrecision -> 30]
(0.999999999999999999999992792996 + 
   2.01329437516665831944251893197 x + 
   1.32777681009864294722106494504 x^2 + 
   0.316766593024067524774032227514 x^3 + 
   0.0197308622074731662156569516898 x^4)/(1 + 
   2.34662770849999165275513365783 x + 
   1.90998604626530684129767865583 x^2 + 
   0.626960209602979190758465977611 x^3 + 
   0.0708418081635094594343411559030 x^4 + 
   0.00124972244476534140880571464115 x^5)
Die Koeffizienten für die 3/4 Approximation liegen sehr nahe an Deinen 
Koeffizienten.
Mit 3/4 liegt der max. relative Fehler bei 4e-19 und mit 4/5 bei 7e-24 
(siehe Diagramme). D.h. Mathematica scheint einen guten Algorithmus 
implementiert zu haben.

von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Relative Fehler siehe Grafiken.

Spontan wäre ich versucht, diese sehr "sinusförmigen" Fehler 
wegzurechnen, indem man den gefundenen Wert noch nachprozessiert. Sieht 
auf den ersten Blick nach leicht komprimierter (quadratischer) X-Achse 
aus, für die man einen Fehlerterm sin(k(x*x)) bestimmen könnte.

?

Auch der A-TAN-rel error direkt hier drüber schreit fast nach weiterer 
Reduktion.

von Johann L. (gjlayde) Benutzerseite


Bewertung
3 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Hallo Johann,
>
> falls es noch von Interesse ist: ich habe den atan jetzt mit Mathematica
> in [0;0.268] mit Grad 3 durch Grad 4 und mit Grad 4 durch Grad 5
> approximiert. Erst habe ich die Koeffizienten für ftmp(x) =
> atan(sqrt(x))/sqrt(x) bestimmt und dann atan durch atapp(x) =
> x*ftmp(x^2) approximiert.

Für atan verwende ich momentan die [6/6] in [0,1], die du oben gepostet 
hattest.

> Mit 3/4 liegt der max. relative Fehler bei 4e-19 und mit 4/5 bei 7e-24
> (siehe Diagramme).

Die [6/6] hat ausreichende Genauigkeit und den Vorteil, dass 
Argument-Reduktion auf [0,1] ausreicht.

Für Genauigkeit siehe in der 1. Grafik ganz unten rechts in 
Beitrag "Re: 64 Bit float Emulator in C, IEEE754 compatibel"

Zur Auswertung eines Polynoms vom Grade n braucht's n Additionen und n 
Multiplikationen.  Nach Reduktion auf [0,1] ergeben sich also folgende 
Kosten bei aktuell folgenden mittleren Laufzeiten in Ticks (ATmega4808):

Add = 400
Mul = 600
Div = 2400 = 4×Mul

[6/6] kostet: 12×Add + 12×Mul + 1×Div = 14400

Weil Reduktion à la atan(x) - atan(y) = atan((x - y) / (1 + xy)) 3×Add + 
1×Mul + 1×Div kostet, ergibt sich:

[3/4] kostet: 10×Add + 8×Mul + 2×Div = 13600

[6/6] - [3/4] kostet: 2×Add + 4×Mul - 1×Div = 2×Add = 800

Mit [3/4] wäre atan also etwa 6% schneller als mit [6/6].

[6/6] braucht 14 Konstanten à 10 Bytes, [3/4] hingegen nur 11 Stück, 
also 30 Bytes weniger, was 15 AVR-Instruktionen entspricht.  [6/6] 
schneidet hier wegen des einfacheren Codes also besser ab, denn ein 
Funktionsaufruf kostet bereits mindestens 8 Bytes: 3×MOVW + 1×RCALL. 
Hinzu kommen bei [3/4] Vergleich und Fallunterscheidung vor Reduktion 
und Fallunterscheidung danach.

von Alexander S. (alesi)


Angehängte Dateien:

Bewertung
1 lesenswert
nicht lesenswert
Johann L. schrieb:
> Für atan verwende ich momentan die [6/6] in [0,1], die du oben gepostet
> hattest.

Ok.

Johann L. schrieb:
> Zum Abschluss noch meine Näherung für asin und acos:
> p(x) =
> + 0.9999999999999999944249107313502758620
> - 1.0352340338921976278427312087167692142 x
> + 0.3529020623298151981342259189772057401 x^2
> - 0.0433348317064168570561235180136569466 x^3
> + 0.0012557428614630796315205218507940285 x^4
> + 0.0000084705471128435769021718764878041 x^5
>
> q(x) =
> + 1
> - 1.1185673672255329236623716486696411533 x
> + 0.4273660095987244885409833401675833351 x^2
> - 0.0635558848496317165994214838980137828 x^3
> + 0.0028820878185134035637440105959294542 x^4
>
> a(x) = p(x) / q(x) in [0, 1/2]
>
> x in [0, 1/2]  =>  Berechne asin(x) = x·a(2x²)
> x in [1/2, 1]  =>  Berechne acos(x) = sqrt(2-2x)·a(1-x)

Zum Vergleich hier die von Mathematica berechneten Werte:
mmlist = MiniMaxApproximation[ArcSin[Sqrt[x/2]]/Sqrt[x/2], {x,{0.000001, 0.5}, 5, 4},WorkingPrecision -> 30]
(0.999999999999999994423051152461 - 
   1.03523410980612966901810389181 x + 
   0.352902121366787757171917578389 x^2 - 
   0.0433348444131093207428523454952 x^3 + 
   0.00125574345990782812575309803332 x^4 + 
   8.47055324771348630110392083335*10^-6 x^5)/(1 - 
   1.11856744313946496503561317403 x + 
   0.427366074961858057240866145321 x^2 - 
   0.0635559015798658404743877188507 x^3 + 
   0.00288208900921310395594554338889 x^4)

nnlist = MiniMaxApproximation[ArcSin[Sqrt[x/2]]/Sqrt[x/2], {x, {0.000001, 0.5}, 6, 5},WorkingPrecision -> 30]
(0.999999999999999999998235405895 - 
   1.29755117580508799919981566022 x + 
   0.609382277207200613471378303921 x^2 - 
   0.124147959189632584750727145119 x^3 + 
   0.0100825882195035338921196573235 x^4 - 
   0.000200565507799759878742583970468 x^5 - 
   9.77108854328832620498643897024*10^-7 x^6)/(1 - 
   1.38088450913842133342518233504 x + 
   0.705705986302069133137423029030 x^2 - 
   0.162645563977985911123266529565 x^3 + 
   0.0162113551808012428119610086569 x^4 - 
   0.000517006584158213720581410722557 x^5)
und den relativen Fehlern im Anhang. Die ersten Koeffizienten sind 
nahezu identisch zu Deinen.

von Geert H. (geerth)


Bewertung
0 lesenswert
nicht lesenswert
I passed this problem to a befriended mathematician.
His quick answer was:
an effcient accurate formula for arctan not too close to 1 is the arctan 
Bailey Borwein Plouffe formula.

von Helmut S. (helmuts)


Bewertung
0 lesenswert
nicht lesenswert
Geert H. schrieb:
> I passed this problem to a befriended mathematician.
> His quick answer was:
> an effcient accurate formula for arctan not too close to 1 is the arctan
> Bailey Borwein Plouffe formula.

This formula has a slowly converging series. Thus I think it's a bad 
choice for compuation.

von Alexander S. (alesi)


Bewertung
1 lesenswert
nicht lesenswert
Nach der Formel für arctan(1/b) mit b aus N in
https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
konvergiert die Reihe linear. Z.B. arctan(1/2)
 1 0.50000000000000000000 +3.635239e-02
 2 0.45833333333333331483 -5.314276e-03
 3 0.46458333333333334814 +9.357243e-04
 4 0.46346726190476189577 -1.803471e-04
 5 0.46368427579365079083 +3.666679e-05
 6 0.46363988658910532115 -7.722412e-06
 7 0.46364927661314381258 +1.667612e-06
 8 0.46364724210793545334 -3.668929e-07
 9 0.46364769089584906281 +8.189504e-08
10 0.46364759050907894400 -1.849173e-08
11 0.46364761321561026586 +4.214804e-09
12 0.46364760803259769117 -9.682084e-10
13 0.46364760922469056004 +2.238845e-10
14 0.46364760894874312847 -5.206297e-11
15 0.46364760901297230600 +1.216617e-11
16 0.46364760899795093296 -2.855163e-12
17 0.46364760900147866662 +6.725741e-13
18 0.46364760900064716509 -1.589639e-13
19 0.46364760900084378559 +3.768359e-14
20 0.46364760900079715622 -8.957166e-15
21 0.46364760900080825845 +2.134233e-15
22 0.46364760900080559392 -5.096469e-16
23 0.46364760900080626005 +1.219466e-16
24 0.46364760900080609352 -2.923270e-17
25 0.46364760900080614903 +7.019400e-18
26 0.46364760900080609352 -1.688200e-18
27 0.46364760900080609352 +4.065000e-19
28 0.46364760900080609352 -9.810000e-20
29 0.46364760900080609352 +2.360000e-20
30 0.46364760900080609352 -5.800000e-21
31 0.46364760900080609352 +1.300000e-21
   0.46364760900080609352

von Johann L. (gjlayde) Benutzerseite


Bewertung
2 lesenswert
nicht lesenswert
Jürgen S. schrieb:
> Johann L. schrieb:
>> Relative Fehler siehe Grafiken.
>
> Auch der A-TAN-rel error direkt hier drüber schreit fast nach weiterer
> Reduktion.

Ich höre da noch nichtmal ein Flüstern :o)

Mal angenommen wir hätten eine (einfach zu berechnende) Näherung h(x) 
des relativen Fehlers r(x).  Dann gilt:

r = (g-f)/f ≈ h  <=>  g ≈ f·(1+h)

In Worten: Um h dazu zu verwenden, eine bessere Näherung g für f zu 
finden, muss man f kennen.

> Sieht auf den ersten Blick nach leicht komprimierter (quadratischer)
> X-Achse aus, für die man einen Fehlerterm sin(k(x*x)) bestimmen
> könnte.

Meine OFR [*] sagt: Polynom :-)

> Spontan wäre ich versucht, diese sehr "sinusförmigen" Fehler
> wegzurechnen, indem man den gefundenen Wert noch nachprozessiert.

Der Fehler ist ja bereits unter dem erforderlichen Limit.  Eine 
Nachkorrektur wäre also dann sinnvoll, wenn der Fehler größer wäre, etwa 
wenn man eine Approximation kleineren Grades wählt.

In diesem Falle hätte man folgende Gemängelage:

Für den (relativen oder absoluten) Fehler eines Proximums vom Grade N 
gilt der Alternantensatz.

Proximum: Hier ein Polynom vom Grade N, das im Sinne der Maximumsnorm 
möglichst dicht (proximal) an der Zielfunktion liegt.  Für rationale 
Funktionen gilt das Gleiche, nur dass der Grad sich aus dem Grad von 
Nenner- und Zählerpolynom ergibt.

Alternantensatz: Der Absolutwert des Proximum-Fehlers hat genau N+2 
lokale Maxima.  Sind x_k die N+2 x-Werte zu diesen lokalen Maxima mit 
x_k < x_{k+1}, dann gilt für den Fehler: r(x_k) = E·(-1)^k.  Der Fehler 
oszilliert also mit konstanter "Amplitude" um 0.  Es kann der seltene 
Fall auftreten, dass es mehr als N+2 solcher Maxima gibt:  Dann gilt E = 
0, d.h. das Proximum stimmt mit der anzunähernden Funktion überein.  Das 
tritt zum Beispiel auf, wenn man ein Polynom vom Grad M durch Polynome 
vom Grade <= M annähert.

Folgerung:  Hat man ein gutes Näherungspolynom vom Grad N gefunden, dann 
hat der Fehler die Alternenteneigenschaft.  Will man den Fehler durch 
ein Polynom vom Grad N+1 darstellen — oder gleichbedeutend damit, ein 
(besseres) Näherungspolynom vom Grad N+1 finden — dann sollte dieses 
ebenfalls die Alternanteneigenschaft haben.  Polynome mit dieser 
Eigenschaft sind die Tschebyschow-Polynome 1. Art T_n:

https://de.wikipedia.org/wiki/Tschebyschow-Polynom#Tschebyschow-Polynome_erster_Art

Im Intervall [-1,1] haben diese folgende kompakte Darstellung:

T_n(x) = cos(n·arccos(x))

Will man also den Fehler als Polynom darstellen, dann ist das nächste 
T_n eine gute Wahl.  Und nach Korrektur ist der Fehler immer noch 
oszillierend, und zwar ungefähr wie

cos((n+1)·arccos(x))

Wenn die Approximation also nicht genau genug ist, dann geht man einfach 
einen Grad höher usw. und fertig ist der Lack.

Analog für rationale Approximationen, nur dass man hier die Wahl hat, 
Nenner- und/oder Zähler-Grad zu erhöhen.  Die besten Ergebnisse erhält 
man i.d.R für grad(Nenner) ≈ grad(Zähler).

Und wie immer: cum granum salis.  Dass der Fehler fast wie 
cos(n·arccos(x)) aussieht, gilt natürlich nur dann, wenn die 
Zielfunktion leicht zu approximieren ist, d.h. bei Entwicklung nach T_n 
werden die Koeffizienten der T_n rasch kleiner mit n, so dass höhere T_n 
nur einen sehr kleinen Beitrag leisten.  Wenn die Zielfunktion nicht 
gutartig ist, z.B. nicht glatt wie |x| in 0, dann ist der Fehler zwar 
immer noch eine Alternante, sieht aber ganz anders aus als die obigen 
Grafiken vermuten lassen.

[*] Optical Function Recognition

von Johann L. (gjlayde) Benutzerseite


Bewertung
3 lesenswert
nicht lesenswert
Alexander S. schrieb:
> Nach der Formel für arctan(1/b) mit b aus N in
> https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
> konvergiert die Reihe linear. Z.B. arctan(1/2)

Wegen der 4 gewinnt man 2 Bit pro Summand.

Die Crux bei BBP ist aber die Art wie die Formel angewandt wird.  Man 
kann mit ihr atan(1/2) geschätzt auf mehrere Hundert oder noch mehr 
Stellen (Basis 4) berechnen, und zwar nur mit 32-Bit float-Arithmetik.

Ganz analog zur Berechnung in 4000 Stellen von Pi mit ATtiny2313.

Beitrag #6030395 wurde vom Autor gelöscht.
von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Wenn die Approximation also nicht genau genug ist, dann geht man einfach
> einen Grad höher usw. und fertig ist der Lack.

Das ist klar, allerdings impliziert ein höherer Grad bei den hier 
vorgestellten Methoden immer auch mehr Divisionen. Ich suche z.B. nach 
Implementierungen, die mit einfachen Divisionen / ohne Wurzelbildungen 
auskommen. Da sind oszillierende Korrekturen sehr viel einfacher zu 
applizieren. Ich mache das z.B. bei Berechnungen nach CORDIC.

von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Jürgen S. schrieb:
> Johann L. schrieb:
>> Wenn die Approximation also nicht genau genug ist, dann geht man einfach
>> einen Grad höher usw. und fertig ist der Lack.
>
> Das ist klar, allerdings impliziert ein höherer Grad bei den hier
> vorgestellten Methoden immer auch mehr Divisionen.

Nicht mit rationalen Approximationen.  Da hat man immer genau eine 
Division: Polynom durch Polynom.  Der Grad von Nenner- und Zählerpolynom 
ist dabei unerheblich.

Einen Grad höher zu wählen kostet also 1 Multiplikation, 1 Addition und 
das Speichern 1 Konstante mehr.

> Ich suche z.B. nach Implementierungen, die mit einfachen
> Divisionen / ohne Wurzelbildungen auskommen.

Also etwa Polynom (z.B. Taylor oder MiniMax) oder rationale Funktion 
(z.B. Padé oder rationale MiniMax).

Bei arccos und arcsin wirst du um Quadratwurzel schwerlich herumkommen, 
denn die Singularitäten in -1 und 1 sind i.W. Quadratwurzeln, und diese 
bekommst du anders nicht vernünftig approximiert.

> Da sind oszillierende Korrekturen sehr viel einfacher zu
> applizieren. Ich mache das z.B. bei Berechnungen nach CORDIC.

??? Aber wenn du zur Korrektur zum Beispiel ein Polynom vom Grad N 
brauchst, dann kannst die Funktion doch direkt per Polynom vom Grad N 
approximieren und CORDIC in die Tonne kloppen.  Ich seh da den Sinn und 
Zweck von CORDIC nicht.  Außerdem hat CORDIC wohl kaum diesen hübschen 
Fehlerverlauf, sondern meine Erfahrung mit CORDIC ist, dass die 
Ergebnisse ziemlich verrauscht sind.

: Bearbeitet durch User
von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Jürgen S. schrieb:
Außerdem hat CORDIC wohl kaum diesen hübschen
> Fehlerverlauf, sondern meine Erfahrung mit CORDIC ist, dass die
> Ergebnisse ziemlich verrauscht sind.

https://www.mikrocontroller.net/attachment/431573/atan_02.png
Ich verwende CORDIC mit 63-bit. Nach 32 Iterationen kann ich das 
Ergebniss direkt ausrechnen. Hintergrund dazu. Nach der Hälfte der 
Iterationen ändert sich der Radius nicht mehr. Das nennt man 
Nano-Stepping.

von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Jens G. schrieb:
> Johann L. schrieb:
>> Außerdem hat CORDIC wohl kaum diesen hübschen
>> Fehlerverlauf, sondern meine Erfahrung mit CORDIC ist, dass die
>> Ergebnisse ziemlich verrauscht sind.
>
> https://www.mikrocontroller.net/attachment/431573/atan_02.png
> Ich verwende CORDIC mit 63-bit. Nach 32 Iterationen kann ich das
> Ergebniss direkt ausrechnen.

Wenn ich die Heuristik richtig lese, ist das Maximum des absoluten 
Fehlers bei 1E-14, d.h. nicht besser als 47 Bits?

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
Da möchte ich doch nochmal in aller Bescheidenheit auf diese, meine 
Implementation des atan in Cordic hinweisen, die mehr als 51 Bit macht:

Beitrag "Atan function implementation"

Cheers
Detlef

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Zur Info: Ich habe den relativen Fehler in meiner Grafik. Aber gut, 
beinahe jede Funktion meiner "Sammlung" sieht ähnlich aus. Mein 
Zahlenformat kann zuverlässig "nur" maximal 31 Bit darstellen.

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Jens G. schrieb:
> Wenn ich die Heuristik richtig lese, ist das Maximum des absoluten
> Fehlers bei 1E-14, d.h. nicht besser als 47 Bits?

Ich habe diesen Randgebieten, weniger als "1 von 1000 Berechnungen", 
bisher keine Beachtung geschenkt.

Mehr Grafiken auch zu anderen Funktionen findet man hier.
Beitrag "Taschenrechner: ATmega1284p - 15x 7-Segmentdisplay"

: Bearbeitet durch User
von Johann L. (gjlayde) Benutzerseite


Angehängte Dateien:

Bewertung
1 lesenswert
nicht lesenswert
Jens G. schrieb:
> Johann L. schrieb:
>> Jens G. schrieb:
>> Wenn ich die Heuristik richtig lese, ist das Maximum des absoluten
>> Fehlers bei 1E-14, d.h. nicht besser als 47 Bits?
>
> Ich habe diesen Randgebieten, weniger als "1 von 1000 Berechnungen",
> bisher keine Beachtung geschenkt.

D.h. man muss nur für viele genügend Werte testen, damit der 
Maximalfehler untem Teppich landet :o)

> Mehr Grafiken auch zu anderen Funktionen findet man hier.
> Beitrag "Taschenrechner: ATmega1284p - 15x 7-Segmentdisplay"

Die Histogramme sagen mir nicht viel; ob die Verteilung von log(Fehler) 
eine Normalverteilung ist oder Poisson-verteilt oder eine spezielle 
Gamma-Verteilung bringt m.E. keinen Erkenntnisgewinnn.

Die wesentliche Metrik, wenn es um Genauigkeit geht, ist da der maximale 
Fehler — je nach Anwendung relativ (wie bei der hier diskutierten 
Floating-Point Anwendung) oder absolut.  Interessant sind vielleicht 
noch Mittelwert oder Median und evtl. noch Quartile wenn man 
verschiedene Implementationen miteinander vergleicht.

Wenn man eine Berechnug hat wie z.B. atan(a-b) und wissen will, ob das 
Ergebnis einer geforderten Genauigkeit genügt, dann gehen in diese 
Betrachtung die Maxmalfehler ein.  Dass die Genauigkeit "im 
Durchschnitt" oder "meistens" eingehalten wird, reicht da nicht.

Und wenn man die Genauigkeit verbessern will, ist oft zielführend, sich 
nen bösen Buben anzuschauen der einen großen Fehler verursacht, um zu 
sehen wo es knirscht in der Berechnung.

Anbei ein recht aktueller Plot der auch das Fehlerbild verschiedener 
Implementierungen zeigt (Grafiken der 3. Zeile).  Die Box-and-Whiskers 
Plots dieser Zeile sind so erstellt, dass die Whiskers bis zu den 
Extremwerten reichen, d.h. es gibt keine Ausreißer jenseits der Whisker 
(im Gegensatz zur Ausführungszeit, wo Ausreißer (z.B. bei Input 0 oder 
NaN) nicht wirklich interessieren).  Die Boxen umfassen den Bereich vom 
25%- bis zum 75%-Quantil, der dicke Strich ist das 50%-Quantil (Median).

Schönes Beispiel sind die Genauigkeiten der sin-Implementierungen: 
avr_f64 etwa hat einen sehr guten Mittelwert, aber ein paar Extremwerte 
verhageln die Qualität.

Interessanterweise gibt es bereits 4 Implementierungen von atan, die 
CORDIC verwenden, aber keine erreicht die Werte der 
Polynom-Implementierung (f7_t bzw. f7), was dan auch deine Frage von 
oben beantwortet:

Jens G. schrieb:
> Zwischenfrage: Weshalb werden hier Näherungsreihen so wesentlich
> diskutiert. Ich würde CORDIC verwenden um an das Ziel zu kommen. Ein zu
> hoher RAM-Verbrauch? .. Eine zu hohe Rechenzeit?

avr_f64:  Genauigkeit ist seht gut, aber Ausführungszeit ca 50% länger, 
Codegröße ca 100% mehr.

fp64lib:  Genauigkeit ca. 7 Bits schlechter; Codegröße und 
Ausführungszeit etwa gleich.  Uwe arbeitet an besserer Genauigkeit, die 
i.W. durch 64-Bit Fixpunkt-Arithmetik begrenzt ist.  Ausführungszeit und 
Codegröße werden aber steigen, und Uwe zieht in Betracht, komplett auf 
Polynome umzusteigen (dann auch für sin, asin, etc.).

Taschenrechner ATmega1284p: Design-Ziel ist nicht mehr als 31 Bits 
Genauigkeit.

Implementierung von Detlef:  Genauigkeit wird limitiert durch 64-Bit 
Fixed-Point, d.h. Problem ähnlich wie bei fp64lib:  In der Nähe von 
Nullstellen kann die für double erforderliche relative Genauigkeit nicht 
erreicht werden.  Codegröße und Ausführungszeit hab ich mir nicht 
angeschaut.

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Taschenrechner ATmega1284p: Design-Ziel ist nicht mehr als 31 Bits
> Genauigkeit.
Besten Dank für die Zusammenfassung der Unterschiede. Mein verwendetes 
Zahlenformat hat offensichtlich die Eigenschaft einen statistischen 
Fehler nach oben (genauer, etwa +2 Dezimalstellen) oder nach unten 
(ungenau, etwa -4 Dezimalstellen) zu haben.

In meiner Anwendnung ist es ein Vorteil. Mit einem Byte mehr als double 
(8 Byte - "54 bit in der Mantisse") erreiche ich eine etwas höhere 
Grundgenauigkeit von etwa "61 bit in der Mantisse" .. mit etwas mehr 
statistischen Fehlern.

Mit welchem Programm wurden die Statistiken erzeugt?

: Bearbeitet durch User
von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Jens G. schrieb:
> Johann L. schrieb:
>> Taschenrechner ATmega1284p: Design-Ziel ist nicht mehr als 31 Bits
>> Genauigkeit.
> Besten Dank für die Zusammenfassung der Unterschiede. Mein verwendetes
> Zahlenformat hat offensichtlich die Eigenschaft einen statistischen
> Fehler nach oben (genauer, etwa +2 Dezimalstellen) oder nach unten
> (ungenau, etwa -4 Dezimalstellen) zu haben.

Das Zahlenformat? Oder die verwendeten Algorithmen?

> In meiner Anwendnung ist es ein Vorteil.

???

Wieso ist Rauschen ein Vorteil?

> Mit welchem Programm wurden die Statistiken erzeugt?

Das Programm ist in C/C++ und wird simuliert mit avrtest.  Auf dem PC 
bestimmt dann ein Python-Skript die Genauigkeiten der Ergebnisse und 
erzeugt eine Tabelle, die dann von R gelesen und dargestellt wird. 
Leider ist R ziemlich langsam und ein Flaschenhals, so dass es nicht 
praktikabel ist, größere Stickproben zu nehmen.  Liegt aber vermutlich 
daran, dass ich irgendein Anti-Pattern verwende.  Mit den R Interna bin 
ich nicht vertraut.

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Jens G. schrieb:
> Wieso ist Rauschen ein Vorteil?
Weil im Ergebniss auch in die richtige Richtung gerundet wird.
Beispiel: asind( acosd( atand( tand( cosd( sind(9)))))) - 9 = 0

> Das Zahlenformat?
Mantisse: Zähler_31_bit / Nenner_31_bit + expo_7_bit
Die Mantisse hat den Bereich 0,3..3,0

Hier eine Idee - Bereichsreduktion: arctan (x) = arctan(c) + arctan((x - 
c) / (1 + x*c))

Mit c=-1/3, c=0, und c=1/3 lässt sich der Bereich zur Berechnung auf 
0,3..0,7 eingrenzen, und damit sind in meinem Fall auch die Fehler 
größer 1e-15 weg.

: Bearbeitet durch User
von Rainer V. (a_zip)


Bewertung
0 lesenswert
nicht lesenswert
Hi, bin mir nicht sicher, ob ich was Erhellendes sagen kann. Als 
"Forthler" rechne ich "eh" mit Ganzzahlen, egal wie skaliert. Und wenn 
ich einen ARCT berechnen muß, dann mache ich das nach Anforderung. Mir 
ist klar, dass das nicht euer Ansatz ist, aber warum für jede Anwendung 
immer das Optimum bereitstellen??
Klar, nicht jeder möchte ARCT aktuell in seiner Anwendung jeweils neu 
berechnen...aber dennoch...
Gruß Rainer und ich habe mit Respekt und Vergnügen mitgelesen!

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
https://www.mikrocontroller.net/attachment/438093/atan_04.png
Ich habe die Berechnung in CORDIC auf 0,3..0,7 begrenzt. Ergebniss: Eine 
leichte Verbesserung in der Genauigkeit, weniger Streuung im Ergebniss.

: Bearbeitet durch User
von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Jens G. schrieb:
> Johann L. schrieb:
>> Jens G. schrieb:
>> Wieso ist Rauschen ein Vorteil?
> Weil im Ergebniss auch in die richtige Richtung gerundet wird.
> Beispiel: asind( acosd( atand( tand( cosd( sind(9)))))) - 9 = 0

Also asind(...) = 9.

Für asin verwende ich wie üblich den Hauptwert, d.h. asin in [-pi/2, 
pi/2].

Oder rechnest du in Altgrad?

mit x = 9*pi/180 und y = asin(acos(atan(tan(cos(sin(x))))))

bekomme ich |x - y| ~ 1.18e-16 ~ 2^{-53}  bei 56-Bit Mantisse für diese 
9 Operationen (inclusive *, /, -).

Wie gesagt ist mein Ziel einen möglichst kleinen Maximalfehler zu 
erreichen.  Ob eine singuläre Berechnung mathematisch exakt ist, ist 
weder kein Design-Zeil noch Design-Kriterium.  Interessant wäre so eine 
Berechnung dann, wenn sie den etablieren Maximalfehler überschreiten 
würde: dann wäre es ein Testfall, anhand dessen die Algorithmen 
verbessert werden könnten.

Hier zeigt sich ein weiterer Vorteil von MiniMax:  Die erhaltenen 
Polynome bzw. rationalen Funktionen haben bekannten Abstand (in der 
Maximumsnorm) von den zu approximierenden Funktionen.  D.h. egal für 
welchen Wert approximiert wird, der Maximalfehler ist nie größer als 
durch den MiniMax-Algorithmus bestimmt.  Die einzige verbleibende Quelle 
für Genauigkeitsverlust ist die verwendete Grundarithmetik (+, -, *, /) 
und im Falle von asin und acos ggf. noch Quadratwurzel.

>> Das Zahlenformat?
> Mantisse: Zähler_31_bit / Nenner_31_bit + expo_7_bit
> Die Mantisse hat den Bereich 0,3..3,0

Also ein Float, der (effektiv) noch durch eine ganze Zahl geteilt wird?
Was bringt das? Nicht-eindeutige Zahlenformate sind doch lästig bei 
Vergleichen, und man vergeudet Darstellungen:

a / b = n·a / n·b

a / b · 2^n = 2a / b · 2^{n-1} = a / (2b) · 2^{n+1}  etc.

> Hier eine Idee - Bereichsreduktion: arctan (x) = arctan(c) + arctan((x -
> c) / (1 + x*c))
>
> Mit c=-1/3, c=0, und c=1/3 lässt sich der Bereich zur Berechnung auf
> 0,3..0,7 eingrenzen, und damit sind in meinem Fall auch die Fehler
> größer 1e-15 weg.

Mit 2-maliger Reduktion kommt man bis auf |x| <= 2-sqrt3 < 0.268

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

Bringt in meinem Falle aber kein Vorteil gegenüber 1-maliger Reduktion 
auf 0 <= x <= 1:

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

: Bearbeitet durch User
von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Jens G. schrieb:
> mit x = 9*pi/180 und y = asin(acos(atan(tan(cos(sin(x))))))
>
> bekomme ich |x - y| ~ 1.18e-16 ~ 2^{-53}  bei 56-Bit Mantisse für diese
> 9 Operationen (inclusive *, /, -).
Das ist ein sehr guter Wert. Wesentlich besser als es die Routinen in 
"C" sind. Excel bekommt bei diesem Test eine schlechte Note.
|x - y| ~ 5.39e-10

: Bearbeitet durch User
von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Hintergrund der Berechnung von .. asind( acosd( atand( tand( cosd( 
sind(9)))))) - 9
.. ist die "Calculator Forensik": 
http://www.rskey.org/~mwsebastian/miscprj/forensics.htm mit Berechnung 
in 360° (360 Degres) Einheitskreis.
Wer mit radian rechnet muss so vorgehen.
http://www.rskey.org/~mwsebastian/miscprj/radians.htm

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Hier gibt es noch einen Algorithmus mit (4,4). Wertebereich 0..0,66. Das 
lässt sich leicht umrechnen.

https://golang.org/src/math/atan.go

Es sind in diesem Dokument noch Verweise zu "C"-Quellen.
// ACCURACY:
//                      Relative error:
// arithmetic   domain    # trials  peak     rms
//    DEC       -10, 10   50000     2.4e-17  8.3e-18
//    IEEE      -10, 10   10^6      1.8e-16  5.0e-17

: Bearbeitet durch User
von Jens G. (jensg)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
CORDIC ist genauer als ihr Ruf - Ich habe mein Algorithmus mal direkt 
getestet. Nach 32 von 63 Iterationen berechne ich direkt das Ergebniss.

Hinweis dazu: Nach 32 Iterationen ändert sich der Radius nicht mehr.

: Bearbeitet durch User
von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Jens G. schrieb:
> Hintergrund der Berechnung von .. asind( acosd( atand( tand( cosd(
> sind(9)))))) - 9
> .. ist die "Calculator Forensik":
> http://www.rskey.org/~mwsebastian/miscprj/forensics.htm

Nach Umwandlung in Altgrad liefert mein Algorithmus

9.000000000000006661338147750939242541790008544921875
= 9 + 15·2^{-51}
= 0x9.000000000001e

Was sagt mir das?

Jens G. schrieb:
> CORDIC ist genauer als ihr Ruf - Ich habe mein Algorithmus mal direkt
> getestet. Nach 32 von 63 Iterationen berechne ich direkt das Ergebnis.

Also auf 58 Bit genau, was man von CORDIC erwarten würde: 63-ld(32) = 
58.

Ist die Genauigkeit relativ, d.h. gilt sie so auch für Float?

: Bearbeitet durch User
von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Johann L. schrieb:
> Jens G. schrieb:
> 9.000000000000006661338147750939242541790008544921875
> = 9 + 15·2^{-51}
> = 0x9.000000000001e
>
> Was sagt mir das?
Hier eine kleine Gegenüberstellung.
x = asind( acosd( atand( tand( cosd( sind(9)))))) - 9

PlanMaker 2018   -3,1247e-08                          
Excel 2013       -1,6732E-10
snc98            -2,4181e-11                          
Genie 92SC        3,1389e-13                          
xxx               6,6613e-15
Win-98 Calc       3,3056e-32
Bei "normaler Berechnung verliert man über Alles etwa 6-7 
Dezimalstellen.

In deiner Rechnung muss irgendwie intern mit höherer Genauigkeit 
gerechnet worden sein. ... weil das Ergebniss dem entspricht, was 
Float64 hergibt (53-bit).

>> CORDIC ist genauer als ihr Ruf - Ich habe mein Algorithmus mal direkt
>> Nach 32 von 63 Iterationen berechne ich direkt das Ergebnis.
>
> Also auf 58 Bit genau, was man von CORDIC erwarten würde: 63-ld(32) =
> 58.
>
> Ist die Genauigkeit relativ, d.h. gilt sie so auch für Float?
Hier die Kommandozeile für x = 0,7500.
spigot -d25 atan(18750/25000)/(5935250132529022095/9223372036854775802)-1 >> atan_cordic_01.txt
5935250132529022095 ist das int64 Ergbeniss direkt aus CORDIC.
spigot ist ein Kommandozeilenprogramm: 
https://www.chiark.greenend.org.uk/~sgtatham/spigot/

von Johann L. (gjlayde) Benutzerseite


Bewertung
2 lesenswert
nicht lesenswert
Die Implementierung hab ich vorerst mal abgeschlossen, für atan wird wie 
gesagt [6/6] MiniMax in [0,1] verwendet.

Zur Integration hab ein paar Makefile-Schnipsel und Wrapper drum gemacht 
und alles in die libgcc eingebaut.  Wer's ausprobieren möchte:

https://sourceforge.net/projects/mobilechessboar/files/avr-gcc%20snapshots%20%28Win32%29/avr-gcc-10.0.0_2019-12-16_-mingw32.zip/download

ist ein ZIP mit avr-gcc.exe für Windos.  Es braucht nichts installiert 
werden, einfach ZIP auspacken und ./bin-Verzeichnis zu PATH hinuifügen, 
oder avr-gcc / avr-g++ per absolutem Pfad aufrufen.  Das Verzeichnis 
kann nach dem Auspacken verschoben oder umbenannt werden.

Es gibt 2 neue Optionen

-mlong-double=32/64 (default: 64)
-mdouble=32/64 (default: 32)

so dass das Layout von (long) double zur Compilezeit gewählt werden 
kann.

Wer sich die Toolchain selbst aus den Quellen generieren will, zum 
Beispiel für Linux:

1) ZIP runterladen und auspacken.

2) Patches aus ./patch auf die Quellen von GCC trunk bzw. avr-libc trunk 
anwenden.

3) Executables x-Flag setzen für:

   ./libgcc/config/avr/libf7/f7renames.sh (gcc)
   ./libgcc/config/avr/libf7/f7wraps.sh (gcc)
   ./devtools/mlib-gen.py (avr-libc)

4) Tools wie üblich aus den Quellen generieren.

Der double-Support ist unvollständig, es fehlen: lround, lrint, atan2, 
Gamma- und Besselfunktionen.  Außerdem gibt es keine 
Konvertierungsfunktionen von / nach String für 64-Bit Floats.

Der Support von 64-Bit Floats in der avr-libc beschränkt sich darauf, 
dass math.h korrekte Prototypen hat, keine double Aliase für 64-Bit 
double definiert werden, es Multilib-Varianten für double = 64 und long 
double = 32 gibt, und xprintf für 64-Bit Floats 8 Bytes überliest (und 
"?" ausgibt).

Ausgabe von 64-Bit Floats wird im avrtest Core-Simulator unterstützt, 
und zwar für double (LOG_DOUBLE), long double (LOG_LDOUBLE), sowie für 
uint64_t, dessen Bits als 64-Bit IEEE 754 double interpretiert werden 
(LOG_D64: 11-Bit Exponent, 52-Bit reduzierte Mantisse, little Endian).

https://sourceforge.net/p/winavr/code/HEAD/tree/trunk/avrtest/

avrtest kann also auch mit avr-Tools verwendet werden, die keinen 
nativen 64-Bit Float support haben.  Zum Beispiel mit Uwes fp64lib oder 
mit Detlefs Implementierung von 
Beitrag "64 Bit float Emulator in C, IEEE754 compatibel"

: Bearbeitet durch User
von M. W. (elektrowagi78) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Schau an: Es gibt doch noch produktive und brauchbare Beiträge hier!

Eine Frage hätte ich:

Jens G. schrieb:
> CORDIC ist genauer als ihr Ruf
Cordic kann beliebig genau sein, wenn man die hinterlegte Tabelle 
entsprechend umfangreich macht. Die Frage wäre aber, wie es mit dem 
Zeitbedarf und dem Aufwand aussieht. Wie würdest du deine Methode 
gegenüber der üblichen Cordic-Leistung / Geschwindigkeit einschätzen?

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Da es ein beim meinem Algorithmus und dem von mir verwendetem 
Zahlenformat ein anderes Entwicklungsziel gibt, möchte ich die Antwort 
auf meinem Thread geben.
Beitrag "Re: Taschenrechner: ATmega1284p - 15x 7-Segmentdisplay"

von Johann L. (gjlayde) Benutzerseite


Bewertung
1 lesenswert
nicht lesenswert
Markus W. schrieb:
> Jens G. schrieb:
>> CORDIC ist genauer als ihr Ruf
> Cordic kann beliebig genau sein, wenn man die hinterlegte Tabelle
> entsprechend umfangreich macht. Die Frage wäre aber, wie es mit dem
> Zeitbedarf und dem Aufwand aussieht. Wie würdest du deine Methode
> gegenüber der üblichen Cordic-Leistung / Geschwindigkeit einschätzen?

Interessant ist natürlich auch der Vergleich mit anderen Methoden, zum 
Beispiel mit rationalen Funktionen.

Meine Implementierung von atan mit Float-Darstellung mit 56-Bit Mantisse 
hat eine Worst-Case Genauigkeit von etwa 55 Bits relativ.

Berechnung dauert ca. 20000 Ticks, das entspricht ca. 1.7 ms bei 12MHz 
(20000 / 12000000).

Jens G. schrieb:
> Zur Zeit verwende ich CORDIC "nur" im Rotationsmodus für jede
> Winkelfunktion: sin, cos, tan, asin, acos, atan. Und das in Radiant
> (Vollkreis = 2 x Pi) und Degrees (Vollkreis = 360°). Um eine andere
> Funktion verwenden zu können wird der Vollkreis intern auf 1
> umgerechnet.
>
> Die Berechnungen dauert etwa 120ms bei einer 8-bit CPU bei 12Mhz Takt.
> Gefühlsmäßig ist das etwas schneller als etwa bei beim TI-30.

Aus der Heuristik ergibt sich eine Worst-Case Genauigkeit von ca. 57 
Bits, also 2 Bit besser, allerdings dauert eine Berechnung ca. 120ms 
also 0.120 * 12000000 Ticks = 1440000 Ticks.  Das wäre über 70 mal 
langsamer als mit rationalen Funktionen für 2 zusätzliche Bits an 
relativer Genauigkeit.

Bist du sicher, dass es da keinen Messfehler gibt?  Dass z.B. Zeiten für 
Ausgabe / String-Konvertierung mit gemessen wurden?

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Bei mir ist alles dabei. Sring-Konvertierung zur Ausgabe an die serielle 
Schnittstelle, ein Betriebssystem auf Basis von Arduino und die interne 
Umrechnung auf mein verwendetes Zahlenformat (Continued Fraction).

Wie man mit professionellen Werkzeugen umgeht - um diese Ticks zu 
messen, ist mir völlig unbekannt. Mein Vorschlag - meinen Quelltext 
holen und selber prüfen.

Die Tastatur wird abgefragt, die LED werden angesteuert, und es gibt ein 
Interrupt 3800 mal in der Sekunde.

: Bearbeitet durch User
von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Damit ich auf die hohe Grundgenauigkeit komme, muss bei der 
Zahleneingabe zu CORDIC mit höherer Genauigkeit gerechnet werden. In 
meiner Arithmetic war eine Multiplikation mit einem Zwischenergebniss 
von 96-bit dabei.
http://www.naughter.com/int96.html

von Johann L. (gjlayde) Benutzerseite


Bewertung
2 lesenswert
nicht lesenswert

: Bearbeitet durch User
von Joggel E. (jetztnicht)


Bewertung
-1 lesenswert
nicht lesenswert
Besten Dank Johann. Super Arbeit.

@Jens.
Die einfachste Methode die Ticks zu messen waere mit einem Oszilloskop 
oder Zaehler. Vorher einen Pin wackeln und nachher wieder. Moderne 
Zaehler koennen gleich ein Histogramm von Zeitdauern anzeigen.

von Rainer V. (a_zip)


Bewertung
0 lesenswert
nicht lesenswert
Jau, tolle Arbeit!!
Gruß Rainer

von Jens G. (jensg)


Bewertung
0 lesenswert
nicht lesenswert
Joggel E. schrieb:
> @Jens.
> Die einfachste Methode die Ticks zu messen waere mit einem Oszilloskop
> oder Zaehler. Vorher einen Pin wackeln und nachher wieder. Moderne
> Zaehler koennen gleich ein Histogramm von Zeitdauern anzeigen.
Die Rechenzeit kann ich bei Arduino bestimmen und auch ausgeben.

Der Messfehler ergibt sich aus der Umwandlung der Rational_64bit auf 
Rational_32bit. Da gibt es erhebliche Streuung in der Rechenzeit -- ist 
jetzt nicht Thema hier bei diesem Thread.

Bitte hier nachschauen: 
Beitrag "Re: Taschenrechner: ATmega1284p - 15x 7-Segmentdisplay"

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.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

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