Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

atan durch cos und sin ???

1,581 views
Skip to first unread message

Harald Ploch

unread,
Jul 16, 1999, 3:00:00 AM7/16/99
to
Hi,

tan = sin / cos

Aber wie läßt sich der atan ausdrücken??

mfg mike

Horst Kraemer

unread,
Jul 16, 1999, 3:00:00 AM7/16/99
to
On Fri, 16 Jul 1999 11:31:08 +0200, "Harald Ploch"
<harald...@xpoint.at> wrote:

> Hi,
>
> tan = sin / cos
>
> Aber wie läßt sich der atan ausdrücken??

Im wesentlichen nur durch atan...

Jedenfalls nicht durch eine endliche "Formel mit" sin,cos,exp +,-,*,/

Bei Uebergang ins Komplexe kannst Du ihn ueber log ausdruecken, aber
zur konkreten Berechnung des komplexen log brauchst Du wieder atan ;-)

MfG
Horst


Hans Steih

unread,
Jul 16, 1999, 3:00:00 AM7/16/99
to
Im Artikel <7mmu3f$g5n$1...@fleetstreet.Austria.EU.net>, "Harald Ploch"
<harald...@xpoint.at> schreibt:

>
>Hi,
>
>tan = sin / cos
>
>Aber wie läßt sich der atan ausdrücken??
>

>mfg mike
>
>

Hallo !

...als inverse trigonometrische Funktion durch andere inverse trigonometrische
Funktionen:

arctan(x) = arcsin(x/sqrt(1+x^2)) = pi/2 - arccos(x/sqrt(1+x^2))

MfG
Hans

______________________________
Hans Steih || HSt...@aol.com
D-47533 Kleve, Germany
"Ich hoffe, es wird niemanden befremden, dass ich den Homer und Virgil zu
Asymptoten gemacht habe" (Lichtenberg, Vom Nutzen der Mathematik)


Ralf Muschall

unread,
Jul 17, 1999, 3:00:00 AM7/17/99
to
Hans Steih schrieb:

> ...als inverse trigonometrische Funktion durch andere inverse trigonometrische
> Funktionen:

> arctan(x) = arcsin(x/sqrt(1+x^2)) = pi/2 - arccos(x/sqrt(1+x^2))

Geht im reellen (nicht global nachgeprüft, ob die Vorzeichen
immer stimmen).

Über den komplexen Zahlen ist die offizielle Definition

arctan x= -i * artanh (i x)
mit
artanh x = [log(1+x)-log(1-x)]/2
log x = log |x| + i*phase(x)

phase x = Winkel von der positiven reellen Achse aus gemessen,
Wertebereich: -pi<phase(x)<=pi.

Systeme, die eine negative Null unterstützen (IEEE), liefern
ohase x= -pi, wenn der Realteil von x negativ und der Imaginärteil
minus Null ist.

Andere, "äquivalent" aussehende Formeln müssen auf richtige
Lage der Unstetigkeiten überprüft werden.

Ralf

Harald Ploch

unread,
Jul 19, 1999, 3:00:00 AM7/19/99
to
Hi,

Danke für Antworten, leider nützen sie mir nicht sehr viel. Mein Problem
besteht darin das atan eine zweitaufwändige Funktion ist und ich sie daher
ersetzen muss. Cos und Sin habe ich vorberechnet und in eine Tabelle
gestellt (jeden halben Grad ein Eintrag). Leider ist das mit acos und asin
nicht so leicht möglich.

Wüsste einer eine Näherung zu atan, die sich im reellen mit elemtarer
Mathematik lösen läßt.

Danke,

MFG Mike

Horst Kraemer

unread,
Jul 19, 1999, 3:00:00 AM7/19/99
to
On Mon, 19 Jul 1999 16:15:54 +0200, "Harald Ploch"
<harald...@xpoint.at> wrote:

> Danke für Antworten, leider nützen sie mir nicht sehr viel. Mein Problem
> besteht darin das atan eine zweitaufwändige Funktion ist und ich sie daher
> ersetzen muss. Cos und Sin habe ich vorberechnet und in eine Tabelle
> gestellt (jeden halben Grad ein Eintrag). Leider ist das mit acos und asin
> nicht so leicht möglich.

Warum nicht ?

Wie so etwas zu loesen ist, laesst sich nicht allgemeingueltig
beantworten.

Entscheidend ist oft, _wozu_ der arctan eigentlich benutzt werden
soll. Vielleicht braucht man ihn gar nicht. Wenn er aus cos|sin
abgeleitet werden soll, die ihrerseits aus einer Tabelle x->cos(x)
oder x->sin(x) bestimmt werden, hat man oft kein Problem, da man
einfach eine Tabelle c->arccos(c) bzw. s->arcsin(s) erstellen kann,
die in c,s gleichabstaendig gequantelt ist und daher einen
Direktzugriff erlaubt.

> Wüsste einer eine Näherung zu atan, die sich im reellen mit elemtarer
> Mathematik lösen läßt.

Vielleicht breitest Du Dein konkretes Problem etwas aus.
Moeglicherweise hast Du Dich auf einen bestimmten Weg festgelegt, weil
Dir kein besserer eingfallen ist und ein anderer Weg fuehrt leichter
zum Ziel (Erfahrungstatsache aus langjaehriger Praxis ;-).

MfG
Horst


Harald Ploch

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to
Hi Horst,

Tja mein problem besteht darin das ich die beiden Richtungsvektoren einer
Liene habe und daraus den Winkel errechnen muss. Ganz einfach - eigentlich.

winkel = atan (vy/vx)
EinheitsvektorX = cos(winkel) --> Cos und Sin sind in einer Tabelle
vorberechnet
Einheitsvektory = sin(winkel)

Leider ist wie gesagt atan zu langsam um meine Bedürfnisse zu erfüllen.

Vielleich fällt dir ja eine leichtere Lösung ein, wäre nicht schlecht.

PS.: Das Problem mit den "Hügeln und Vektoren" habe ich gelöst. Eigentlich
ganz simpel:

Die neuen Richtungsvektoren ergeben sich aus der Vektoraddition.
Die neue Geschwindigkeit ergibt sich aus folgender Schlußrechnung:

Alte_Vektorlänge : Neuer_Vektorlänge = AlteSpeed² : NeuerSpeed²

Danke für deine Ansätze in dieser Richtung.

mfg,
Mike

Wolfram Hinderer

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to
In article <7n1bhg$fqn$1...@fleetstreet.austria.eu.net>,
"Harald Ploch" <harald...@xpoint.at> writes:

> Tja mein problem besteht darin das ich die beiden Richtungsvektoren einer
> Liene habe und daraus den Winkel errechnen muss. Ganz einfach - eigentlich.
>
> winkel = atan (vy/vx)
> EinheitsvektorX = cos(winkel) --> Cos und Sin sind in einer Tabelle
> vorberechnet
> Einheitsvektory = sin(winkel)

Falls Du den Vektor (vx,vy) normieren willst, versuche mal (vx/d,vy/d)
mit d=sqrt(vx^2+vy^2). Sollte ziemlich schnell sein, wird aber nur
funktionieren, wenn Du den Winkel nicht noch fuer was anderes brauchst.

(Bist Du Dir sicher, dass das die richtige Stelle ist, um das Programm
zu optimieren?)

Wolfram

Harald Ploch

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to

Wolfram Hinderer schrieb in Nachricht
<7n1eub$t0m$1...@news.rz.uni-karlsruhe.de>...

Hi,

Tut mir leid dich das fragen zu müssen, aber: Hast du eigentlich die
verhergehenden Messages gelesen?

Nun um es klaren auszudrücken, ich benötige die Vektoren nicht mehr (die
habe ich bereits anderwertig bekommen). Aber ich benötige den Winkel
unbedingt, und zwar mit hilfe der beiden Parameter (vx,vy). Daran führt kein
Weg vorbei.

Was deine Frage der Optimierung betrifft. Das Programm ist ein Spiel und
soll auf einem P75 mit 16 MB RAM und Win95 mit DirectX 6.1 laufen (und zwar
flüssig). Da ich den Atan an bestimmten Stellen bis zu 20 oder 30 mal pro
Sekunde aufrufe erscheint es mir doch vernünftig hier zu optimieren.
Weil wir gerade von Optimierung sprechen, du weißt nicht zufällig eine
möglichkeit fabs und abs durch schnellere Funktionen zu ersetzen.

Bye,
mike

Wolfram Hinderer

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to
In article <7n1g1g$hjr$1...@fleetstreet.austria.eu.net>,
"Harald Ploch" <harald...@xpoint.at> writes:


> Tut mir leid dich das fragen zu müssen, aber: Hast du eigentlich die
> verhergehenden Messages gelesen?

Schon, aber die Darstellung deines Problems liess mich daran zweifeln,
ob ich verstanden hatte, was dem Problem eigentlich zugrunde liegt.


> Was deine Frage der Optimierung betrifft. Das Programm ist ein Spiel und
> soll auf einem P75 mit 16 MB RAM und Win95 mit DirectX 6.1 laufen (und zwar
> flüssig). Da ich den Atan an bestimmten Stellen bis zu 20 oder 30 mal pro
> Sekunde aufrufe erscheint es mir doch vernünftig hier zu optimieren.
> Weil wir gerade von Optimierung sprechen, du weißt nicht zufällig eine
> möglichkeit fabs und abs durch schnellere Funktionen zu ersetzen.

Ich benutze hier einen K6/266, wo pro Sekunde 1 Mio. atan drin sind.
Ich glaube, auch auf Deinem Rechner sind weniger als 100 Aufrufe
keine grosse Angelegenheit. Zu fabs und abs kann ich leider nichts
sagen, glaube aber nicht, dass dabei viel Zeit verlorengeht
(es sei denn natürlich, Du hast das ausprobiert).

Die bereits von Horst vorgeschlagene Tabelle waere doch eine Idee zur
Beschleunigung. Ich vermute aber, dass die Zeit an anderer Stelle,
verlorengeht typischerweise passiert das bei irgendwas was mit dem
Zeichnen zu tun hat. Eine kurze Erläuterung, was alles in der
Hauptschleife passiert, wuerde vielleicht was nuetzen, aber da sind
wir hier falsch.

Tip: Alles, was mit Zeichnen zu tun hat, aus dem Programm werfen,
schauen, wie schnell es läuft. Dann genau umgekehrt. So merkt man, wo
es klemmt.
(Mit Profiler geht das natuerlich viel schöner.)

Wolfram

Harald Ploch

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to
>Die bereits von Horst vorgeschlagene Tabelle waere doch eine Idee zur
>Beschleunigung.

Für Cos und Sin arbeite ich schon mit Tabellen, bei ATan ist mir nicht ganz
klar wie das funktinieren soll. Wie soll ich eine solche Tabelle indizieren.
Bei Sinus und Cosinus funktioniert das über Grade (lassen sich ganz einfach
in Radianten und zurück rechnen).

>Ich vermute aber, dass die Zeit an anderer Stelle,
>verlorengeht typischerweise passiert das bei irgendwas was mit dem
>Zeichnen zu tun hat.

Das Zeichnen ist nicht das Problem, bei DirectX habe ich im Debug - Modus 40
Frames / Sekunde.
Ich denke das reicht völlig.

Die einzigen Optimierungsmöglichkeiten liegen bei der Mathematik und einigen
Arrays die ich verwende. Da es sich dabei allerdings schon um fixe Array
handelt, läßt sich hier wohl nicht mehr viel machen.

Ein zusätzliches Problem ist das ein Vergleich sehr schwer fällt, da ich auf
einem PII/300 arbeite und das fertige Spiel auf einem P75 laufen soll.

mfg,
mike

Horst Kraemer

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to
On Tue, 20 Jul 1999 10:22:03 +0200, "Harald Ploch"
<harald...@xpoint.at> wrote:

> Hi Horst,



> Tja mein problem besteht darin das ich die beiden Richtungsvektoren einer
> Liene habe und daraus den Winkel errechnen muss. Ganz einfach - eigentlich.
>
> winkel = atan (vy/vx)
> EinheitsvektorX = cos(winkel) --> Cos und Sin sind in einer Tabelle
> vorberechnet
> Einheitsvektory = sin(winkel)

Irgendwie glaube ich, dass wir Sprachmissverstaendnisse haben. Sofern
wir von der Ebene reden:

cos(winkel) ist _kein_ _Einheitsvektor_.

Das Paar (cos(winkel),sin(winkel)) ist _der_ _Einheitsvektor_ in
Richtung "Winkel". Ein _Einheitsvektor_ ist ein Vektor, dessen Betrag
1 ist, und dies ist bei cos(winkel) oder sin(winkel) allein i.A. nicht
der Fall.

Ich bin wieder mal frech: Wozu brauchst Du eigentlich den Winkel ?

Wenn Du ihn dazu brauchst, um die neue Richtung eines sich geradlinig
bewegenen Objekts o,.ae. zu bestimmen, das an einer Wand, dem
Spielfeldrand oder an einem Kreis o.ae. reflektiert wird - dazu
braucht man keinen Winkel. Das macht man allein mit Multiplikationen,
Additionen und Divisionen auf (pardon) Einheitsvektoren und zwar
voellig exakt...
(Das bestgehuetete mathematische Geheimnis der Flipper-Programmierer
;-)

[Sorry, falls ich schief liege, aber dies ist _der_ Standardfehler,
den angehende Spieleprogrammierer immer wieder machen: Sie glauben
fuer Reflektionen a la "Einfallspinsel=Ausfallspinsel" atan oder atan2
bemuehen zu muessen. Vielleicht bist Du ja einen Ausnahme...]

MfG
Horst


Thomas Richard

unread,
Jul 20, 1999, 3:00:00 AM7/20/99
to
Harald Ploch schrieb:

> Für Cos und Sin arbeite ich schon mit Tabellen, bei ATan ist mir nicht ganz
> klar wie das funktinieren soll. Wie soll ich eine solche Tabelle indizieren.

Heißt das, es läßt sich nicht vorhersagen, in welchem Intervall das atan-
Argument liegt?
Ansonsten sehe ich hier kein Problem, atan ist doch streng monoton.

> Bei Sinus und Cosinus funktioniert das über Grade (lassen sich ganz einfach
> in Radianten und zurück rechnen).

Jaja, aber hat das mit der Tabellen-Erstellung zu tun?

> Ein zusätzliches Problem ist das ein Vergleich sehr schwer fällt, da ich auf
> einem PII/300 arbeite und das fertige Spiel auf einem P75 laufen soll.

Takt reduzieren?
Ihr Spiele-Bubis (sorry) schraubt doch immer daran rum. ;-)

MfG
Richie
--
Thomas....@post.rwth-aachen.de
Thomas_...@ac2.maus.de (max. 16kB pro Mail)

Ernst-Udo Wallenborn

unread,
Jul 21, 1999, 3:00:00 AM7/21/99
to
"Harald Ploch" <harald...@xpoint.at> writes:

> >Die bereits von Horst vorgeschlagene Tabelle waere doch eine Idee zur
> >Beschleunigung.
>

> Für Cos und Sin arbeite ich schon mit Tabellen, bei ATan ist mir nicht ganz
> klar wie das funktinieren soll. Wie soll ich eine solche Tabelle indizieren.

> Bei Sinus und Cosinus funktioniert das über Grade (lassen sich ganz einfach
> in Radianten und zurück rechnen).
>

Felix Pütsch hat das hier vorgeschlagen:

|x|<1: atan x=x-x^3/3+x^5/5-...
x>1 : atan x=pi/2-1/x+1/(3x^2)-1/(5x^5)+...
x<-1 : atan x=-pi/2-1/x+-...


Also die Idee ist: x=vy/vx. Dann machst du eine Tabelle für x
zwischen 0 und 1. Werte für x außerhalb des Intervalls evaluierst
du mit den beiden Reihenentwicklungen. Profis würden den atan außerhalb
von [-1,1] stückweise durch Polynome annähern, was vielleicht ein
wenig Rechenzeit spart (in den guten Implementationen von erf(x)
wird das zum Beispiel so gemacht).

--
Ernst-Udo Wallenborn
Laboratorium fuer Physikalische Chemie
ETH Zuerich

Dirk Thierbach

unread,
Jul 21, 1999, 3:00:00 AM7/21/99
to
Horst Kraemer <horst....@snafu.de> wrote:
> On Tue, 20 Jul 1999 10:22:03 +0200, "Harald Ploch"
> <harald...@xpoint.at> wrote:
>> Tja mein problem besteht darin das ich die beiden Richtungsvektoren einer
>> Liene habe und daraus den Winkel errechnen muss. Ganz einfach - eigentlich.
>>
>> winkel = atan (vy/vx)
>> EinheitsvektorX = cos(winkel) --> Cos und Sin sind in einer Tabelle
>> vorberechnet
>> Einheitsvektory = sin(winkel)

> Irgendwie glaube ich, dass wir Sprachmissverstaendnisse haben. Sofern
> wir von der Ebene reden:

> cos(winkel) ist _kein_ _Einheitsvektor_.

Das wurde mit ihm in dem Thread "Schiefe Ebene / Vektorrechnung (wichtig!!!)"
schon mal diskutiert. Siehe dort fuer seine Definition von
Einheits- und Richtungsvektor.

[...]

> Ich bin wieder mal frech: Wozu brauchst Du eigentlich den Winkel ?

Ich gehe mal davon aus, dass es sich immer noch um das Golfspiel
handelt. Ich hatte ihm in <7mf8is$j2e$1...@sun27.hrz.tu-darmstadt.de>
(selber Thread wie oben) schon mal vorgeschlagen, auf die Winkel doch
zu verzichten, aber offenbar will er das nicht. Wenn er mit den
Winkeln gluecklicher wird, ist mir das dann auch recht - nur darf er
nicht erwarten, dass sein Programm dann schneller oder fehlerfreier
wird. Sinuesse, Cosinuesse und atans sind immer langsamer als eine
einfache Multiplikation oder Division. Aber das ist nicht mein
Problem :-)

> [Sorry, falls ich schief liege, aber dies ist _der_ Standardfehler,
> den angehende Spieleprogrammierer immer wieder machen: Sie glauben
> fuer Reflektionen a la "Einfallspinsel=Ausfallspinsel" atan oder atan2
> bemuehen zu muessen. Vielleicht bist Du ja einen Ausnahme...]

Glaube ich nicht.

- Dirk

Horst Kraemer

unread,
Jul 21, 1999, 3:00:00 AM7/21/99
to
On 21 Jul 1999 10:09:58 +0200, Ernst-Udo Wallenborn
<walle...@phys.chem.ethz.ch> wrote:


> Felix Pütsch hat das hier vorgeschlagen:

> |x|<1: atan x=x-x^3/3+x^5/5-...
> x>1 : atan x=pi/2-1/x+1/(3x^2)-1/(5x^5)+...
> x<-1 : atan x=-pi/2-1/x+-...


Das ist die bekannte Reihenentwicklung, die fuer |x|~1 herzlich
schlecht konvergiert.




> Also die Idee ist: x=vy/vx. Dann machst du eine Tabelle für x
> zwischen 0 und 1. Werte für x außerhalb des Intervalls evaluierst
> du mit den beiden Reihenentwicklungen.

Naivere Gemueter wuerden fuer x>1 einfach pi/2-arctan(1/x) aus der
Tabelle ablesen...

> Profis würden den atan außerhalb
> von [-1,1] stückweise durch Polynome annähern, was vielleicht ein
> wenig Rechenzeit spart (in den guten Implementationen von erf(x)
> wird das zum Beispiel so gemacht).

Ich glaube, dass Realwelt-Implementierungen eher den atan fuer 0<=x<=1
stueckeweise oder wie auch immer durch Polynome annaehern und fuer x>1
schlicht pi/2-arctan(1/x) berechnen.


In der Praxis kann man auch folgende Formel benutzen

x-d
arctan(x) - arctan(d) = arctan -----
1+x*d

Berechnet man z.B. vor

d := sqrt(2)-1 ad:=arctan(d)


So erhaelt man fuer d<=x<=1 die Formel

x-d
arctan(x) = arctan ----- + ad
1+x*d

x-d
Hier gilt "zufaellig" dass 0<= ----- <=d
1+x*d


wenn d<=x<=1

so dass d eine gemeinsame obere Schranke fuer x in der obigen
Reihenentwicklung ist.

Bei diesem einfachen Splitting ist der relative Fehler des Resultats
bei der nach x^5/5 abgebrochenen Reihenentwicklung bereits fuer alle
0<=x<=1 unter 1/1000. Man kann dieses Splitting beliebig weit treiben
und jedes Intervall wiederum im Verhaeltnis sqrt(2)-1 splitten. Durch
jedes binaere Splitting verringert sich die obere Schranke fuer den
relativen Maximalfehler um den Faktor d^6 = 0.005, wenn ich richtig
gepeilt habe.


MfG
Horst


Roland Franzius

unread,
Jul 21, 1999, 3:00:00 AM7/21/99
to
Hallo
Ich hab mal vor 10 Jahren eine 80bit Arithmetik auf 68k CPU
geschrieben. Für die Approximationspolynome habe ich
>>Hart, Computer Approximations<< benutzt, auch eine gute Quelle für unbekannte Funktionalformeln. So weit ich erinnere: Der atan wird unterteilt in [0,1] und [1,\infty]. Die Spiegelung beruht auf dem Additionstheorem atan x +- atan y = atan (x+-y)/(1-+ x y) mit y=1/x.
Für eine 18 stellige Genauigkeit muß der Wertebereich [0,Pi/2]
in 12 gleiche Intervalle geteilt werden, d.h das Intervall
[0,1] wird durch
tan (n Pi/24) ungleichmäßig unterteilt. Der Rest ist durch
eine rationale Funktion minimal 5. Grades zu bekommen.
actan x - arctan (tan n pi/24 ) = actan ((x - tan(n pi/24))/(1
+ x tan (n pi/24))
Mit Hash-Table ist das neben Gamma in etwa die langsamste
Funktion im All.


Gruss roland


Horst Kraemer schrieb:
>
>[snip]

> x-d
> arctan(x) - arctan(d) = arctan -----
> 1+x*d
> Berechnet man z.B. vor
>
> d := sqrt(2)-1 ad:=arctan(d)
>
> So erhaelt man fuer d<=x<=1 die Formel
>
> x-d
> arctan(x) = arctan ----- + ad
> 1+x*d
>
> x-d
> Hier gilt "zufaellig" dass 0<= ----- <=d
> 1+x*d
>
> wenn d<=x<=1
> so dass d eine gemeinsame obere Schranke fuer x in der obigen
> Reihenentwicklung ist.
>

= Bei diesem einfachen Splitting ist der relative Fehler des
Resultats
= bei der nach x^5/5 abgebrochenen Reihenentwicklung bereits
fuer alle
= 0<=x<=1 unter 1/1000. Man kann dieses Splitting beliebig
weit treiben
= und jedes Intervall wiederum im Verhaeltnis sqrt(2)-1
splitten. Durch
= jedes binaere Splitting verringert sich die obere Schranke
fuer den
= relativen Maximalfehler um den Faktor d^6 = 0.005, wenn ich
richtig
= gepeilt habe.

--
Roland Franzius
Theor. Physik FB Physik, Univ. Osnabrueck
+++ exactly <<n>> lines of this message have value <<FALSE>>
+++

Ernst-Udo Wallenborn

unread,
Jul 22, 1999, 3:00:00 AM7/22/99
to

horst....@snafu.de (Horst Kraemer) writes:
> Naivere Gemueter wuerden fuer x>1 einfach pi/2-arctan(1/x) aus der
> Tabelle ablesen...

duh. Klarer Fall von Wald vor lauter Bäumen nicht gesehen.


> In der Praxis kann man auch folgende Formel benutzen
>

> x-d
> arctan(x) - arctan(d) = arctan -----
> 1+x*d
>
>
> Berechnet man z.B. vor
>
> d := sqrt(2)-1 ad:=arctan(d)
>
>
> So erhaelt man fuer d<=x<=1 die Formel
>
> x-d
> arctan(x) = arctan ----- + ad
> 1+x*d
>
> x-d
> Hier gilt "zufaellig" dass 0<= ----- <=d
> 1+x*d
>
>
> wenn d<=x<=1
>
> so dass d eine gemeinsame obere Schranke fuer x in der obigen
> Reihenentwicklung ist.
>

> Bei diesem einfachen Splitting ist der relative Fehler des Resultats

> bei der nach x^5/5 abgebrochenen Reihenentwicklung bereits fuer alle

> 0<=x<=1 unter 1/1000. Man kann dieses Splitting beliebig weit treiben

> und jedes Intervall wiederum im Verhaeltnis sqrt(2)-1 splitten. Durch

> jedes binaere Splitting verringert sich die obere Schranke fuer den

> relativen Maximalfehler um den Faktor d^6 = 0.005, wenn ich richtig

> gepeilt habe.


Wow. Das kannte ich nicht. Laß mal sehen: Nach dem ersten Schritt
hat das Splitting eine Fallunterscheidung mit sieben Fällen
(x<-1/d, x<-1, x<-d, x<d, x<1, x<1/d, x>1/d) und für jeden
Splittingschritt kommt etwa ein Faktor zwei bei der Zahl
der Fallunterscheidungen dazu. Demgegenüber kann man eine simple
Interpolationstabelle ohne großen Mehraufwand (fast) beliebig
groß machen. Benchmarken wir das mal:

atan(x) für (x=-10;x>=10;x+=0.1), 201 Punkte, 10^6 Durchläufe,
compiliert mit cc -O3 auf einer R10000 SGI O^2:

atan2 294 sec
interpolation 200 sec |max err| 8.1e-6 mean 4.1e-6
split 1 134 sec |max err| 2.6e-4 mean 3.7e-5

Die Interpolationstabelle hatte hier grade 100 Einträge,
wurde mit atan2 gefüllt und linear interpoliert.

Der Kern der Splitting-Implementation sieht dabei so aus:

#define PI2 1.570796326794896619231321692
#define D1 0.414213562373095048801688724
#define D1I 2.414213562373095048801688724
#define AD1 0.392699081698724154807830423

double atan0 (double x) {
return x * ( 1 - x*x * ( 1/3.0 - x*x/5.0 ) );
}

double atan (double x) {
double r;
int i;
if (x<-D1I) {
r = -PI2 - atan0(1/x);
} else if (x<-1) {
r = -PI2 + atan0((1+x*D1)/(D1-x)) + AD1;
} else if (x<-D1) {
r = -atan0((x+D1)/(D1*x-1)) - AD1;
} else if (x<D1) {
r = atan0(x);
} else if (x<1.0) {
r = atan0((x-D1)/(1+D1*x)) + AD1;
} else if (x<D1I) {
r = PI2 - atan0((1-x*D1)/(x+D1)) - AD1;
} else {
r = PI2 - atan0(1/x);
}
return r;
}

modulo bugs. Das ist mehr als doppelt so schnell wie atan2,
und bei minimalem overhead liegt der maximale relative Fehler
unter 10^-3.

0 new messages