Arcus Tangens
Dass mathematische Funktionen so ihre Tücken haben können, hat wohl jeder schon erfahren, der einmal in Megamax-C (Version 1.0) versucht tat. den Arkustangens zu berechnen.
Bekanntlich geschieht das mit dem Funktionsaufruf ‘atan(x)’. Das Üble an 'atan(x)’ ist, daß es - je nach dem Wert für x - abwechselnd korrekte, völlig falsche oder - das ist nun wiederum ein entschiedener Vorzug bei der Fehlersuche - gar keine Resultate liefert. Für Ihre eigenen Versuche: alles in Sicherheit bringen und dann atan(x) für
x =5.79/sqrt (7.0)
berechnen. Wenn sich der Rechner nochmal zurückmeldet, dann haben Sie vermutlich eine verbesserte Version erwischt. Ich hatte auch nach einer halben Stunde noch kein Ergebnis! Abbrechen ließ sich das Programm auch nicht - seitdem habe ich eine resetfeste Ram-Disk.
Über die Fehler des Megamax-Entwicklungssystems möchte ich mich an dieser Stelle nicht weiter auslassen. Inzwischen ist ja auch die Version 2.0 erschienen, die so manche Besserung verspricht. Sie sind hier nur der Anlaß, Sie als Leser dazu anzuregen, sich einmal etwas näher mit numerischer Mathematik zu beschäftigen. Denn auch die hat ihre Tücken, denen man nur zu '»eicht zum Opfer fällt. Wenn Sie Lust haben, daran teilzunehmen, so sind Sie herzlich eingeladen.
Sie wissen sicher alle noch aus der Schule, wie der Tangens eines Winkels definiert ist. Am rechtwinkligen Dreieck gilt:
Gegenkathete
tan(x) =...................
Ankathete
Der Winkel x wird normalerweise - entgegen unseren Schulgewohnheiten - nicht in Grad, sondern im Bogenmaß (rad) ausgedrückt. Die Umrechnung ist sehr einfach:
rad =grad * 2 pi/360; ji=3.14159265358.,
grad= rad * 360/(2 * n)
An diese Darstellungsweise werden Sie sich schnell gewöhnen. Selbstverständlich würde man vor der Umrechnung in beiden Gleichungen mit 2 kürzen.
Relativ häufig steht man vor dem Problem, daß der Tangens eines Winkels bekannt ist. Gesucht wird der zugehörige Winkel. Die Funktion, mit der wir diesen Winkel berechnen, heißt Arkustangens oder in C ‘atan()’. Diese Umkehrfunktion des Tangens liefert uns allerdings nur die sog. Hauptwerte des gesuchten Winkels, denn wie Sie sich sicherlich erinnern, gilt:
tan(x)=tan(x+7i)=tan(x+2n:)=...=tan(x+rm
Das heißt, es macht keinen Unterschied, ob wir den Tangens des Winkels x berechnen oder vorher auf x ganze Vielfache von n aufaddieren. Entsprechend kann uns atan(x) nur Winkel im Bereich ±Jt liefern. Das wird allerdings erst dann eine Bedeutung haben, wenn wir diese Funktion anwenden, nicht jedoch bei der Programmierung der Funktion. Ferner gilt:
tan(x) = tan ( atan (tan(x)) )
und, allerdings nur für die Hauptwerte:
atan(x) = atan (tan (atan(x)))
Das heißt: die ‘Umkehrung der Umkehrung’ liefert wieder den Ausgangswert. Atan() kann man nun durch eine sog. Potenzreihe berechnen. Die Formeln hierfür stammen aus dem ‘Taschenbuch der Mathematik’ 1:
Ich gehe einmal davon aus, daß Sie keine Erfahrung im Umgang mit Potenzreihen haben und werde deshalb die einzelnen Teile etwas näher erklären. X ist in beiden Formeln natürlich der Tangens des gesuchten Winkels und wird dem Programm vorgegeben. Die Kreiszahl n kennt auch jeder. Bliebe in Gleichung (1) noch der Ausdruck:
x2n+1
... +...+ (-l)n—-------- ; | x |≤1
2n+1
'+...+’ bedeutet, daß beliebig (unendlich) viele gleichartig gebaute Summanden auftreten, bevor man zum ‘letzten’ (nach +...+) kommt, der dann endlich das Prinzip angibt, nach dem alle vorherigen gebaut sind. Die Variable n steht dabei für die Nummer des Summanden: n=1,2,3... Der Ausdruck (-1)n liefert -1 oder +1, je nachdem ob n ungerade oder gerade ist und soll nur für die wechselnden Vorzeichen sorgen. Für den 17-ten Summanden ergibt sich dann:
n = 17: -1 * x35 / 35; |x|≤1
Die Aufgabe unseres Programms besteht also darin, all diese Quotienten nacheinander zu berechnen und aufzusummieren. Allerdings müssen wir noch das Vorzeichen von x beachten. Da x jedoch nur in ungeraden Potenzen auftritt, genügt es, sich das Vorzeichen zu merken, die ganze Berechnung für positives x durchzuführen und dann das Ergebnis mit dem Vorzeichen zu multiplizieren. Das ist alles. Bleibt nurnoch zu klären, wieviele solcher Quotienten zu berechnen sind, denn unendlich viele schafft auch der schnellste Rechner nicht. Dazu wieder zurück zu unseren Formeln. Beide Ausdrücke haben zwei sehr wesentliche Gemeinsamkeiten:
- Die Summe der Quotienten konvergiert mit wachsendem n gegen einen Grenzwert. Dieser Grenzwert ist der gesuchte Winkel im Bogenmaß. Diese Eigenschaft habe ich bisher stillschweigend vorausgesetzt.
- Auch die einzelnen Summanden konvergieren mit wachsendem n und zwar gegen Null. Man sagt auch: die SuM manden der Reihe bilden eine Nullfolge.
Von diesen Eigenschaften gehen wir im folgenden aus, ohne uns um einen Beweis dafür zu kümmern. Er wurde bereits von anderen für uns erbracht. Konvergieren heißt in diesem Fall übrigens: je mehr Teilbrüche wir berechnen, d.h. je größer n wird, umso näher kommt die Summe dieser Brüche dem Grenzwert, ohne ihn je wirklich exakt erreichen zu können. Dabei muß selbstverständlich auch die Rechengenauigkeit eine erhebliche Rolle spielen. Man macht die Entscheidung über den Abbruch der Berechnung übrigens nicht von n abhängig, sondern benutzt die zweite Eigenschaft, die Konvergenz der Summanden gegen Null: ein Summand, der hinreichend nahe bei Null liegt, leistet für die Endsumme keinen wesentlichen Beitrag. Da die Beträge aller folgenden Summanden noch kleiner sind (Stichwort: Nullfolge), brechen wir die ganze Berechnung dann ab, wenn ein kritischer Wert unterschritten wird. Das hat den Vorteil, daß wir auf jeden Fall alle Winkel mit der gleichen Anzahl gültiger Stellen berechnen. Wie wir noch sehen werden, wäre das bei einer festen Vorgabe von n nicht der Fall. Der maximale Fehler, den man macht, wenn man die Berechnung für ein bestimmtes n einfach abbricht, kann über das (n+1)-te Glied der Reihe geschätzt werden. Als kritischen Wert, den die Quotienten unterschreiten müssen, legen wir 0.0000001 fest. Dann erhalten wir den gesuchten Winkel auf sechs Stellen genau nach dem Komma.
So einfach ist das also
Zunächst zu Gleichung (1). Sie gilt nur für |x| ≤ 1. Die Nenner berechnen wir, indem wir eine gleichnamige Variable bei jedem Schleifendurchlauf um 2 erhöhen. Für die Potenzen von x machen wir von den Regeln für die Potenzrechnung und Multiplikation Gebrauch, nach denen gilt:
x * (-x2) = -x3
(-x3) * (-x2) = x5 ...
oder
(-x) * (-x2) = x3
x3 * (-x2) = -x5 ...
D.h. durch Multiplikation mit -x2 werden nicht nur die verschiedenen Potenzen von x erzeugt, sondern auch gleichzeitig der Vorzeichenwechsel berücksichtigt. Für die Aufsummierung dient die Variable arcus mit dem Anfangswert 0. Die Berechnung wird abgebrochen, sobald ein Quotient den kritischen Wert unterschreitet. Weil die einzelnen Quotienten für diesen Vergleich benötigt werden, summieren wir sie nicht direkt auf, sondern speichern sie zunächst in der Variable quotient. Die Variable, die der Zähler aufnehmen soll, habe ich mit Rücksicht auf Gleichung (2) faktor genannt.
Die zweite Gleichung scheint auf den ersten Blick wesentlich anders gebaut zu sein, als die erste. Die wesentlichen Unterschiede:
- die Potenzen von x treten nur im Nenner auf
- sämtliche Vorzeichen sind vertauscht
- es gibt den zusätzlichen Summanden ±_π, wobei gilt: +_π für x ≥ 0
Trotzdem kann die Berechnung mit der gleichen Schleife erfolgen. Faktor erhält für fl Gleichung (2) allerdings nicht den Wert -x2, sondern -1/x2. Um die Vorzeichen gegenüber Gleichung (1) zu vertauschen, beginnen wir die Folge der Zähler mit -1/x statt mit x und die Summenvariable arcus wird nicht mit 0, sondern mit _π vorbesetzt. Wir rechnen ja nur für den Betrag von x und berücksichtigen das Vorzeichen erst am Schluß. Alles andere bleibt gleich. Das Struktogramm für eine vorläufige Version sieht dann folgendermaßen aus:
Ganz so einfach ist es doch nicht
Wer es nicht erwarten konnte und daraus bereits ein lauffähiges Programm erstellt hat. sollte es erst einmal auf Diskette sichern und dann versuchen, damit den Arkustangens von 1 oder-1 zu berechnen. ... Hat's geklappt ? Klassischer Reinfall, nehme ich an. Was passiert ist ? Ganz einfach: das Gleiche, wie bei meiner Megamax-Version. Ihr Rechner hat gerechnet -und zwar endlos. Um das zu verstehen nochmals ein Blick auf die einzelnen Quotienten. Nehmen wir der Einfachheit halber die Gleichung (1). Wenn hier x = 1, so reduziert sich die ganze Potenzreihe auf den Ausdruck:
1 1 1
atan(1)=1 - - + - - - + (-1)2n+1—;
3 5 7
In unserem Programm hatten wir gefordert, daß der Betrag des letzten Glieds kleiner als 0.0000001 sein soll. Da es um den Betrag geht, können wir das Vorzeichen vergessen und schreiben:
1
--- ≤ 0.0000001
2n+1
-> 2n+1 ≥ 10 000 000.
-> 2n ≥ 9 999 999.
-> n ≥ 4 999 999.5
Ich hoffe, Sie haben die letzten Schritte nachvollziehen können. Nur als Tip: mit Ungleichungen kann man im Prinzip genauso rechnen, wie mit normalen Gleichungen. Nur muß man beachten, daß bei Kehrwertbildung oder Vorzeichenumkehr sich auch die Richtung des ‘≤’-Zeichens umkehrt.
Schlußfolgerung: Bevor unser Programm bei x = 1 das Abbruchkriterium erreichen kann, müssen mindestens 5 Millionen (!) Quotienten berechnet und aufsummiert werden ! Erschwerend kommt noch hinzu, daß bei sehr großen bzw. kleinen Zahlen die Rechengenauigkeit mitspielt, so daß es durchaus sein kann, daß das Abbruchkriterium - aufgrund von Rundungsfehlern -nie erreicht wird. In diesem Fall hätten wir es sogar mit einer 'echten' Endlosschleife zu tun. Probieren Sie das Programm doch noch einmal für Werte von
x= 1.5 . 1.1 . 1.01 , 1.001,1.0001 ..
bzw.
x = 0.5,0.9,0.99,0.999,0.9999 ..
aus. Sie werden feststellen, daß die Berechnung umso länger dauert, je näher x dem kritischen Wert von 1 kommt. Zur Begründung wollen wir uns der Einfachheit halber auf die positiven x-Werte beschränken. Betrachten wir Gleichung (2) für Werte |x|>1. Die Potenzen von x stehen hier im Nenner. Wenn x = 1, so leistet nur der jeweilige Nenner 1, 3, 5, 7 usw. einen Beitrag dazu, den Bruch zu verkleinern (ihn dem Abbruchkriterium näher zu bringen). Ist x > 1, so tragen auch die Potenzen von x dazu bei. Allerdings: Je näher x bei 1 liegt, umso größer muß n gewählt werden, damit der Beitrag von x-n wesentlich wird. So ist:
0.0000001- 1.1-169 -> n = 169
0.0000001 - 1.01-1621 -> n=1621
0.0000001 - 1.001-16127 -> n = 16127
Entsprechendes gilt für Gleichung (1). Wenn Sie soweit alles verstanden haben, wundem Sie sich vermutlich über gar nichts mehr. Wahrscheinlich ahnen Sie so langsam auch, was bei der Megamax-Funktion atan() schiefläuft. Nun kann man zwar ganz und gar nicht davon sprechen, daß der von uns gewählte Algorithmus vollständig versagt, denn auch für den kritischen Wert von x=1 ist die Konvergenz der Potenzreihe - theoretisch - gesichert (den Beweis erbrachte übrigens der Philosoph und Mathematiker Leibniz, und der Grenzwert ist 1/4π). Aber sie ist für diesen Wert am schlechtesten. Nach einem simplen Programmierfehler werden Sie also lange und vergeblich suchen !
Für unser Programm heißt das, daß es keinesfalls genügt, wenn wir nur eine ‘kleine Abfrage' einbauen, mit der der kritische Wertx= 1 behandelt wird. Wir müssen vielmehr für einen größeren Bereich von x die Konvergenz der Potenzreihe verbessern. Das Stichwort dazu heißt Winkelsubstitution.
So geht es auf jeden Fall schneller
Ich habe im vorigen Abschnitt etwas weiter ausgeholt, um zu zeigen, daß das was nun folgt, keineswegs nur einer fragwürdigen ‘Eleganz’ dienen soll, sondern unumgänglich ist! In diesem Abschnitt werde ich Sie nun etwas mit Formeln foltern. Aber ich hoffe, es wird sich in Grenzen halten.
Der Grundgedanke ist eigentlich ganz einfach. Wir berechnen für 0.5<|x|<1.5 den gesuchten Winkel nicht direkt, sondern ersetzen (substituieren) den Tangens dieses Winkels durch den Tangens eines anderen Winkels, von dem wir zumindest soviel wissen, daß er zu dem gesuchten Winkel y in einer besonders einfachen Beziehung steht. Dabei nutzen wir die Tatsache aus. daß sich jede reelle Zahl - also auch der gesuchte Winkel - als Summe zweier anderer reeller Zahlen darstellen läßt, wovon wir eine stets frei wählen können. Der gesuchte Winkel sei y. Tan (y) = x ist vorgegeben. Gesucht wird y = atan(x). Es gilt:
(3) y = a+b
-> x = tan(y) = tan (a+b)
atan(x) = atan(tan(y)) = atan(tan(a+b))
atan(x) = y = atan(tan(a+b))
Ich hoffe, die mehrfachen ‘=’-Zeichen irritieren Sie nicht zu sehr. Solche Kettengleichungen besagen nichts anderes, als daß alle gleichgesetzten Ausdrücke - na, eben gleich sind und ich mir daraus ‘Teilgleichungen’ entnehmen kann. Wenn u = v = w gilt, dann gilt auch u = w,u = v,v = w. Davon werde ich noch mehrfach Gebrauch machen.
Wie gesagt: einen der Werte a oder b können wir willkürlich festlegen. Viel gewonnen ist damit zunächst nicht. Glücklicherweise gibt es aber Umformungen für den Tangens einer Summe von Winkeln. Insbesondere gilt:
tan(a) + tan(b)
(4) tan(a+b) = x = -----------
1- tan(a) * tan (b)
und wegen tan(l/4π) = 1 sind wir so frei, wählen für a =1/4π und erhalten:
tan(_π+b) = x = tan(y) = tan (a+b)
1 + tan (b)
(5) tan (_π+b) = x = tan(y) = ———
1 - tan(b)
Wir wollen tan(b) berechnen. Der Wert von tan(y) ist vorgegeben. Nach entsprechender Umformung von (5) erhält man:
x - 1
(6) tan(b) = ----- ; x <> -1
x + 1
In kurzen und dürren Worten (statt Formeln) heißt das:
- berechne ein tan(b) so, daß Gleichung (6) erfüllt ist.
- berechne für diesen Wert den Arkus tangens durch Reihenentwicklung (vgl. oben).
- berechne aus dem gefundenen Winkel den gesuchten Winkel y durch Addition von 1/4π.
Die Einschränkung x <> -1 ist notwendig, da wir andernfalls durch Null dividieren könnten. Da wir aber nur für positive x rechnen und das Vorzeichen erst wieder am Ende berücksichtigen, stört uns das nicht weiter. Der gesuchte Winkel ergibt sich dann für positive x (vgl, oben):
x-1
(7) y = atan(x) = atan(—) + _π;
x+1
Hier also die endgültige Fassung des Programms kann man nebenstehend sehen.
Wie Sie inzwischen vielleicht festgestellt haben, benötigt das Programm für beliebiges x nie mehr als rund 20 Schleifendurchläufe - gegenüber maximal ~ 5 Millionen in der ersten Fassung. Das ist doch eine ganz beachtliche Temposteigerung, auch ohne jede zusätzliche Hardware, sprich: floating point Coprozessor. Es aßt sich übrigens noch verbessern: durch veränderte Wahl von a (dann aber Gleichung (4) benutzen) und einen anderen Bereich für die Winkelsubstitution. Ohne es selbst weiter nachgerechnet oder ausprobiert zu haben, vermute ich, daß eine ~ 2 bis 5-fache Geschwindigkeitssteigerung durch-aus drin ist. Aber das soll Ihnen überlassen bleiben. Der Nachweis, daß eine Lösung ‘optimal’ ist, ist übrigens nicht mehr so ganz trivial.
Zur Ehrenrettung der Megamax-Programmierer sei noch vermerkt: Diesen ‘Trick’ haben sie wohl auch gekannt, nur ist ihnen dabei ein Fehler unterlaufen. Daher auch der zu Beginn erwähnte krumme Wert von x. Mir wird man wahrscheinlich Vorhalten, daß meine Programmversion nicht unbedingt ein Muster an Übersichtlichkeit und Strukturiertheit ist. Ganz unrecht hätte man damit vielleicht nicht. So wäre es durchaus möglich, die Winkelsubstitution rekursiv zu gestalten und den Rest auf mehrere if-Blöcke zu verteilen. Gewiß, gewiß. Da merkt man halt noch den alten Fortran-Stil. Aber in der Numerik ist es was Effizienz und Genauigkeit betrifft, häufig sehr viel günstiger, einen Ausdruck erst einmal in seine Bestandteile zu zerlegen und dabei auf Modularität und ‘Lesbarkeit’ wenig Rücksicht zu nehmen. Es gibt ja schließlich auch Kommentare. Strukturierung kann keineswegs das Verständnis für den zugrundeliegenden Algorithmus ersetzen. Genau darauf, dieses Verständnis zu vermitteln, kam es mir an, und nicht so sehr darauf, möglichst viele Programme in den Beitrag zu packen.
Zum Schluß vielleicht aber doch noch eines, sozusagen als Programmrätsel des Monats, späte Rache des ehemaligen Fortran-Programmierers an der versammelten C-Gilde oder Zugeständnis an den neu(deutsch)en Stil: Was berechnet die folgende Funktion mit 4 Buchstaben ?
Listing einfügen
Sehr richtig: die Quadratwurzel. Dieser Algorithmus ist übrigens nicht nur außerordentlich effizient, sondern auch ausgesprochen gutmütig: er korrigiert sich nämlich selbst ! Gleichgültig, was Sie für IRGEND_W AS einsetzen - nur größer Null muß es sein - Sie bekommen immer das richtige Ergebnis2.
/*....................allgemeine Defines....................*/
#define abs(x) ((x>0.0)?(x):(-x)) /* Betrag bilden */
#define PI 3.14159265358979323846 /* kennt wohl jeder */
#define EPS 0.0000001 /* Abbruchkriterium */
/*...........................................................*/
/* Berechnung des Arkus-Tangens */
/*...........................................................*/
float atan(x)
float x;
/*.........................Anfangswerte......................*/
{
float x2; /* - x^2 */
float faktor; /* Anfangswert f.Zähler */
float nenner=1.0; /* Anfangswert f.Nenner */
float arcus=0.0; /* Anfangswert f.Summe */
float vorzeichen=1.0; /* */
float quotient=1.0; /* */
if (x<0.0) vorzeichen=-l.0; /* Vorzeichen x abspalten */
x=abs(x); /* und Betrag von x bilden*/
x2=-x*x; /* -x^2 fuer Potenzierung */
faktor=x; /* 'Zaehler', Anfangswert */
/*.......................Winkelsubstitution ?............... */
if ((0.5 <= x) && (x <= 1.5)) /* x im kritischen Bereich? */
{ /* ja */
arcus+=PI*0.25; /* 'a' (vgl. Text) */
x=(x-l.0)/(x+1.0); /* tan(b) */
x2=-x*x; /* nicht vergessen: x ueberall durch */
faktor=x; /* tan(b) ersetzen */
}
/*......................Welche Gleichung gilt ?.................*/
/* */
/* diese Entscheidung erst nach der Winkelsubstitution ! ! */
if (x > 1.0)
{ /* Gleichung (2) */
faktor= -1.0/x; /* Kehrwerte bilden */
x2=1.0/x2; /* auch fuer Potenzierung */
arcus+=PI*0.5; /* Ergebnis + 1/2 PI */
}
/*............................Schleife......................*/
do
{
quotient=faktor/nenner; /* neuer Quotient */
arcus+=quotient; /* aufsummieren */
faktor*=x2; /* x potenzieren */
nenner+=2.0; /* Nenner erhöhen */
} while (abs(quotient)>EPS); /* Abbruchkriterium ? */
return (vorzeichen*arcus); /* Ergebnis nicht vergessen */
} /* Ende 'atan()' */
Manfred Klaus Herrmann