Gospers Naeherung an n!
Peter Luschny
Hugo Pfoertner machte dieser Tage hier eine Bemerkung, die
bei mir hängen blieb.
Auf Mathworld http://mathworld.wolfram.com/StirlingsApproximation.html
wird eine Näherung von n! hergeleitet, und zwar
n! approx sqrt(2PI) (n/e)^n sqrt(n)
und dann gesagt:
Gosper has noted that a better approximation to n! is given by
n! approx sqrt(2PI) (n/e)^n sqrt(n+1/6)
Meine Frage ist: wie geht Gospers Näherung weiter?
Genauer meine ich damit folgendes:
Stirlings Näherung lässt sich bekanntlich so weiterentwickeln:
n! approx sqrt(2PI) (n/e)^n sqrt(n) (1 + 1/(12n))
Und bei Gosper?
n! approx sqrt(2PI) (n/e)^n sqrt(n+1/6) (1 + ?)
Irgendwelche Vorschläge?
Gruss Peter
Peter Luschny
> Stirlings Näherung lässt sich bekanntlich so weiterentwickeln:
> n! approx sqrt(2PI) (n/e)^n sqrt(n) (1 + 1/(12n))
>
> Und bei Gosper?
> n! approx sqrt(2PI) (n/e)^n sqrt(n+1/6) (1 + ?)
>
> Irgendwelche Vorschläge?
Hier meine erste Lösung, in Maple Parlance:
> stirling0 := proc(n) sqrt(2*Pi)*sqrt(n)*n^n*exp(-n) end:
> stirling1 := proc(n) sqrt(2*Pi)*sqrt(n)*n^n*exp(-n)*(1+(1/(12*n))) end:
> gosper0 := proc(n) sqrt(2*Pi)*sqrt(n+1/6)*n^n*exp(-n) end:
> gosper1 := proc(n) sqrt(2*Pi)*sqrt(n+1/6)*n^n*exp(-n)*(1+(1/(12*n))^2) end:
> approx := proc(n) local prec; prec := 20;
> print("Stirling",n, evalf(stirling0(n)/n!,prec), evalf(stirling1(n)/n!,prec));
> print("Gosper",n, evalf(gosper0(n)/n!,prec), evalf(gosper1(n)/n!,prec)); end:
> seq(approx(10^i),i=1..4);
Formel n, 0-te Naeherung , 1-te Naeherung
-------------------------------------------------------------------
Stirling, 10, 0.99170403955606148635, 0.99996823988569533204
Gosper, 10, 0.99993408971442072976, 1.0000035295817620090
Stirling, 100, 0.99916701656784300019, 0.99999965574831620272
Gosper, 100, 0.99999930910413435918, 1.0000000035480990149
Stirling, 1000, 0.99991667014156998579, 0.99999999653074844997
Gosper, 1000, 0.99999999305910486512, 1.0000000000035492614
Stirling, 10000, 0.99999166670139157014, 0.99999999996528074841
Gosper, 10000, 0.99999999993055910485, 1.0000000000000035493
======================================================================
Noch einmal für das Protokoll, die Gosper-Luschny Näherung an n!.
n / 1 \
n! ~ sqrt(2Pi) sqrt(n + 1/6) (n / e) |1 + ------ |
| 2 |
\ 144 n /
======================================================================
Also klar ist, es geht weiter, und zwar so, dass es dem guten alten
Stirling die Socken auszieht. Jetzt will ich natürlich die ganze Story!
Die nächsten drei Terme bitte, und dann der Beweis.
Gruss Peter
Markus Sigg
> Meine Frage ist: wie geht Gospers Näherung weiter?
Vielleicht hat Gosper selbst etwas dazu geschrieben? Ich weiß nicht, ob
er die zitierte Näherung veröffentlicht hat, und falls ja, wo. Du könntest
es mit dieser Arbeit versuchen:
Gruß,
Markus
Christopher Creutzig
> Gosper has noted that a better approximation to n! is given by
Dazu habe ich leider keine Quelle, aber wenn Dich auch noch andere
Reihen interessieren, emofehle ich wärmstens Lanzcos, A Precision
Approximation of the Gamma Function, Journal of the Society for
Industrial and Applied Mathematics: Series B, Numeric Analysis, Vol. 1
(1964), p. 86-96,
http://links.jstor.org/sici?sici-0887-459X%281964%291%3C86%3AAPAOTG%3E2.0.CO%3B2-F
(ich hoffe, die URL passt, copy&paste ging aus irgendeinem Grund nicht).
Da stehen auch Formeln mit exakten rationalen Koeffizienten drin, die
im Gegensatz zu Stirlings Formel sogar für alle z mit Re(z)>0 konvergieren.
Gruß,
Christopher
--
Manche Aussagen sind nur deshalb nicht falsch, weil sie den dazu
erforderlichen Präzisionsgrad nicht erreichen. (Lothar Schmidt)
Stefan Wehmeier
>
> ======================================================================
>
> Noch einmal für das Protokoll, die Gosper-Luschny Näherung an n!.
>
> n / 1 \
> n! ~ sqrt(2Pi) sqrt(n + 1/6) (n / e) |1 + ------ |
> | 2 |
> \ 144 n /
>
> ======================================================================
>
> Also klar ist, es geht weiter, und zwar so, dass es dem guten alten
> Stirling die Socken auszieht. Jetzt will ich natürlich die ganze Story!
> Die nächsten drei Terme bitte, und dann der Beweis.
>
Hallo Peter,
es wird Dich enttäuschen, aber da Du schon mit dem Computer die numerischen
Werte ausrechnest, könntest Du ihn natürlich auch für Deine Frage
verwenden:
series(n!/(n^n/exp(n)*sqrt(2*PI)*sqrt(n+1/6)), n=infinity)
[in Maple: Pi statt PI]
Modifikation der Aufgabe:
Stattdessen könnte man natürlich auch zusätzliche Terme im Argument der
Wurzel hinzufügen, also
sqrt(n + 1/6 + (1/72)*1/n + ...)
Ich verzichte jetzt mal bewusst darauf, die Lösung einfach abzutippen, um
den anderen den Spass nicht zu verderben. Irgendwie schade, dass die blöden
Rechner alle diese schönen Aufgaben kaputtmachen ...
Gruss,
Stefan
--
Stefan Wehmeier
ste...@math.uni-paderborn.de
Boudewijn Moonen
> Dazu habe ich leider keine Quelle, aber wenn Dich auch noch andere
> Reihen interessieren, emofehle ich wärmstens Lanzcos, A Precision
> Approximation of the Gamma Function, Journal of the Society for
> Industrial and Applied Mathematics: Series B, Numeric Analysis, Vol. 1
> (1964), p. 86-96,
> http://links.jstor.org/sici?sici-0887-459X%281964%291%3C86%3AAPAOTG%3E2.0.CO%3B2-F
> (ich hoffe, die URL passt, copy&paste ging aus irgendeinem Grund nicht).
Der Link hat bei mir nicht funktionoert.
> Da stehen auch Formeln mit exakten rationalen Koeffizienten drin, die
> im Gegensatz zu Stirlings Formel sogar für alle z mit Re(z)>0 konvergieren.
Siehe dazu
http://en.wikipedia.org/wiki/Lanczos_approximation
Dort gibt es auch weiterfuehrende Verweise.
Gruss Boudewijn
Boudewijn Moonen
> Ich verzichte jetzt mal bewusst darauf, die Lösung einfach abzutippen, um
> den anderen den Spass nicht zu verderben.
Klar, die haben ja alle Maple oder Mathematica.
Gruss Boudewijn
Peter Luschny
>> Christopher Creutzig wrote:
>> Dazu habe ich leider keine Quelle,
Ich habe auch schon vergeblich danach gefahndet,
unter anderem in dem legendären Hakmem Report,
ftp://publications.ai.mit.edu/ai-publications/pdf/AIM-239.pdf
bislang ohne Erfolg.
http://mathworld.wolfram.com/StirlingsApproximation.html
gibt keine Referenz. Online-gossip about Gosper? Und der
Rest der Interwelt scheint von Mathworld abzuschreiben.
>> aber wenn Dich auch noch andere
>> Reihen interessieren, emofehle ich wärmstens Lanzcos, A Precision
>> Approximation of the Gamma Function, Journal of the Society for
>> Industrial and Applied Mathematics: Series B, Numeric Analysis, Vol. 1
>> (1964), p. 86-96,
Klar, ein weiterer Klassiker in diesem Feld.
>> [... http://links.jstor.org/sici?sici ...]
> Der Link hat bei mir nicht funktionoert.
Der Link kostet Geld :-( Aber dieser Link ist umsonst und auch noch
hübsch zu lesen, mit Bildchen und so..
http://laplace.physics.ubc.ca/~matt/Doc/ThesesOthers/Phd/pugh.pdf
Gruss Peter
Roland Franzius
Das ist doch dasselbe, aber nur bis in der Näherung 1/n.
sqrt(n+1/6)= sqrt(n) sqrt(1+1/(6n))
= sqrt(n)(1+1/(12n) -1/(288 n^2) + 1 /(3456 n^3)+-..)
Die asymptotische Stirlingreihe ist definitiv besser für Werte über ca
6, denn
Gamma(x)=x^(x-1/2) e^(-x) sqrt(2 pi) (1+1/(12 x)) + 1/(288 z^2) -
139/(51840 z^3) - 571/(2488320 z^4)
Das falsche Vorzeichen im quadratischen Glied ist verheerend.
Die Asymptotische Reihe bei x=0 zu benutzen, wie bei Mathworld
angegeben, beweist nur absolute Unkenntnis der Sachlage. Bei
asymptotische Reihen existiert zu jedem eps>0 ein n(eps) und ein r, so
dass für alle |x|>r die Reihe bei genau n abgebrochen weniger als eps
vom Wert der Funktion entfernt ist. Sowohl anderes n als auch |x|<r ist
das Fehlerglied >0.
Für die Approximation der Gammafunktion im Bereich 1<x<2 muss man für
gegebenes Fehlerepsilon hochgradige Pade-Approximationen nehmen. Die
Funktion ist "noch transzendenter" als Exp- und Winkelfunktionen und
besitzt außer der Multiplikationseigenschaft keine lokalen Symmetrien.
--
Roland Franzius
Roland Franzius
Nachtrag:
Es gibt einen Zeilberger-Gosper-Algorithmus zur hypergeometrischen
Summation zB für Mathematica
http://www.fmf.uni-lj.si/~petkovsek/aisb.m
der intensive Bemühungen zur Behandlung der Fakultät beinhaltet.
--
Roland Franzius
Christopher Creutzig
> Der Link hat bei mir nicht funktionoert.
Copy&Paste:
http://links.jstor.org/sici?sici=0887-459X%281964%291%3C86%3AAPAOTG%3E2.0.CO%3B2-F
Gruß,
Christopher
Peter Luschny
>> Peter Luschny wrote:
>> ======================================================================
>> Noch einmal für das Protokoll, die Gosper-Luschny Näherung an n!.
>>
>> n / 1 \
>> n! ~ sqrt(2Pi) sqrt(n + 1/6) (n / e) |1 + ------ |
>> | 2 |
>> \ 144 n /
>> ======================================================================
> es wird Dich enttäuschen, aber da Du schon mit dem Computer die numerischen
> Werte ausrechnest, könntest Du ihn natürlich auch für Deine Frage verwenden:
> series(n!/(n^n/exp(n)*sqrt(2*PI)*sqrt(n+1/6)), n=infinity)
Yes. Ich bin überhaupt nicht enttäuscht. Mir gefällt es, dass diese Dinger
das können. Was ich hoffe eines Tages noch zu finden, ist eine /gerade/
asymptotische Entwicklung für n!.
> Modifikation der Aufgabe:
> Stattdessen könnte man natürlich auch zusätzliche Terme im Argument der
> Wurzel hinzufügen, also
> sqrt(n + 1/6 + (1/72)*1/n + ...)
Yupp. Nächste gute Idee.
> Ich verzichte jetzt mal bewusst darauf, die Lösung einfach abzutippen, um
> den anderen den Spass nicht zu verderben.
Da habe ich keine Sorge. Ich sehe das als eine 'Benchmark' und will
den sportlichen Ehrgeiz bei unseren Youngsters wecken, noch bessere
Formeln zu finden. -- Nicht Jungs? Aber hallo!
> Irgendwie schade, dass die blöden
> Rechner alle diese schönen Aufgaben kaputtmachen ...
Du siehst das zu akademisch ;-)
Ich referiere mal den Zwischenstand. Wir haben drei Formel-Familien
zur asymptotischen Näherung an n! - Stirling, Gosper und, well, du
verzeihst für den Lauf des Spiels die Verwendung deines Namens - die
Wehmeier Formeln. Höherer Index bedeutet höhere Genauigkeit.
> kern := proc(n) sqrt(2*Pi)*n^n*exp(-n) end:
>
> stirling0 := proc(n) kern(n)*sqrt(n) end:
> stirling1 := proc(n) kern(n)*sqrt(n)*(1+(1/12)/n) end:
> stirling2 := proc(n) kern(n)*sqrt(n)*(1+(1/12)/n+(1/288)/n^2) end:
>
> gosper0 := proc(n) kern(n)*sqrt(n+1/6) end:
> gosper1 := proc(n) kern(n)*sqrt(n+1/6)*(1+(1/144)/n^2) end:
> gosper2 := proc(n) kern(n)*sqrt(n+1/6)*(1+(1/144)/n^2-(23/6480)/n^3) end:
>
> wehmeier0 := proc(n) kern(n)*sqrt(n+1/6) end:
> wehmeier1 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n) end:
> wehmeier2 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n)*(1-(31/12960)/n^3) end:
Welches ist die beste? Halt, hier meldet sich in letzter Minute noch
ein Außenseiter an: Rumpelstilz! Wer ist Rumpelstilz?
> rumpelstilz := proc(x) local a; a := x+x+1;
> exp((log(2*Pi)+log(a/2)*a-a-(1-7/(30*a*a))/(6*a))/2) end:
Das Benchmarkprogramm ermittelt die Anzahl der gültigen dezimalen
Stellen. Hier das Resultat:
Formel n! '0' '1' '2'
-----------------------------------------
Stirling 00010! 2.1 4.5 5.6
Gosper 00010! 4.2 5.5 7.7
Wehmeier 00010! 4.2 5.6 8.5
Rumpelstilz 00010! 8.2
Stirling 00100! 3.1 6.5 8.6
Gosper 00100! 6.2 8.5 11.9
Wehmeier 00100! 6.2 8.6 12.4
Rumpelstilz 00100! 13.1
Stirling 01000! 4.1 8.5 11.6
Gosper 01000! 8.2 11.4 15.9
Wehmeier 01000! 8.2 11.6 16.3
Rumpelstilz 01000! 18.1
Stirling 10000! 5.1 10.5 14.6
Gosper 10000! 10.2 14.4 19.9
Wehmeier 10000! 10.2 14.6 20.3
Rumpelstilz 10000! 23.1
Ein Überraschungssieger! Rumpelstilz hat knapp die Nase vorne!
Gruss Peter
P.S. Hier ist das Benchmarkprogramm:
> decdig := proc(f,n,nf) -log[10](evalf(abs(1.0-f(n)/nf),30)) end:
> approx := proc(n) local prec,nf; nf := n!;
> printf("Stirling %.5d! %6.1f %6.1f %6.1f \n",
> n,decdig(stirling0,n,nf),decdig(stirling1,n,nf),decdig(stirling2,n,nf));
> printf("Gosper %.5d! %6.1f %6.1f %6.1f \n",
> n,decdig(gosper0,n,nf),decdig(gosper1,n,nf),decdig(gosper2,n,nf));
> printf("Wehmeier %.5d! %6.1f %6.1f %6.1f \n",
> n,decdig(wehmeier0,n,nf),decdig(wehmeier1,n,nf),decdig(wehmeier2,n,nf));
> printf("Rumpelstilz %.5d! %20.1f \n",n,decdig(af,n,nf)) end:
HP
allen hier nochmal ein Diskussion zum Nachlesen empfehlen, die ich in
der "numerical analysis" Newsgroup vor ueber 4 Jahren mit David W.
Cantrell ueber dessen Naeherungen der Gammafunktion gefuehrt hab. Siehe
die bei Diskussionsfolge nach
http://groups.google.com/group/sci.math.num-analysis/msg/521fa1a6fb98a300
Peter Luschny
> [...] Diskussion zum Nachlesen empfehlen, die ich in
> der "numerical analysis" Newsgroup vor ueber 4 Jahren mit David W.
> Cantrell ueber dessen Naeherungen der Gammafunktion gefuehrt hab.
Danke für diesen Link, den ich recht anregungsreich fand.
Ich musste etwas schmunzeln darüber, dass du Cantrell sofort auf
den Algorithmus von John L. Spouge angesprochen hast. 5 Jahre zuvor
hatte Cantrell in einem seiner Bücher (Topics in Advanced Scientific
Computation) diesem Algorithmus einen längeren Absatz gewidmet.
Zur Vervollständigung dieses Threads hier, kopiere ich dein Zitat herüber:
John L. Spouge: Computation of the Gamma, Digamma, and Trigamma
Functions, SIAM J. Numer. Anal. Vol. 31, No. 3, pp. 931-944, June 1994
Page 934 contains the following (cited verbally):
<<
The next theorem is a curiosity.
THEOREM 1.3.2. The relative error in the approximation
(13) z! ~ ((z+1/2)/e)^(z+1/2) * (2*pi)^(1/2)
is less than (2*pi*e)^(-1/2)*(ln(2))*[Re(z+1/2)]^(-1).
Since (2*pi*e)^(-1/2) * pi^(-1)*(ln(2)) = 0.053..., equation (13) is not
only simpler than Stirling's approximation, but nearly twice as accurate
for any z>0.
>>
Dies ist aber /nicht/ die Gospersche Formel (ich weiß, das hat auch
niemand behauptet). Wird also in die Sammlung zur Benchmark neu mit
aufgenommen. Aber das "but nearly twice as accurate for any z>0" kann
ich nicht nachvollziehen. Bei der Gosper-Formel trifft dies aber exakt zu!
(Und deswegen hatte mich diese Bemerkung beim ersten Lesen irregeführt.)
Jetzt kommen wir zu Lanczos (Gleichung 21).
Hmm, als erstes sehen wir, dass die kürzeste Form davon gerade die
obige Formel von Spouge ist. Lanczos schrieb vor Spouge, also
taufen wir jetzt die Formel von Spouge in Lanczos0 um.
Und geben noch Lanczos1 und Lanczos2 an.
Wobei ich gleich schon in der Codierung hervorhebe, dass die Grundidee dabei
ist, Gamma(n+1/2) statt Gamma(n) zu approximieren. Cantrell:
I think that (*) clearly indicates that it is more "natural" to define
a function shifted by 1/2, making it midway between z! and Gamma(z),
lanczos0 := proc(n) local N; N:=n+1/2; kern(N) end;
lanczos1 := proc(n) local N; N:=n+1/2; kern(N)*(1-(1/24)*(1/(n+1))) end:
lanczos2 := proc(n) local N; N:=n+1/2; kern(N)*(1-(1/24)*(1/(n+1))
-(23/1152)*(1/((n+1)*(n+2)))) end:
Dabei ist, wie in meinem anderen Posting, kern := proc(n) sqrt(2*Pi)*n^n*exp(-n) end.
In der Benchmark werden wir sehen, dass diese Formeln, mit Ausnahme der 0-ten
Näherung, schlechter abschneiden als die entsprechenden Stirlingschen!
Kommt Cantrell. Will genau das in Lanczos Formeln reparieren. Hier seine
Formeln in den drei einfachsten Formen.
cantrell0 := proc(n) local N; N:=n+1/2; kern(N) end:
cantrell1 := proc(n) local N; N:=n+1/2; kern(N)/(1+1/(24*N-1/2)) end:
cantrell2 := proc(n) local N; N:=n+1/2; kern(N)/(1+1/(24*N-1/2+1/((1440/2021)*N))) end:
So, jetzt sind die Dinge langsam sortiert...
Aber wie sieht die Benchmark aus? Können die Wehmeier Formeln gegen
die Cantrells bestehen? Wir machen eine kleine Pause für die Werbung,
gleich geht's weiter.
<Werbung>
Es gibt Fälle, in denen asymptotische Formeln nicht mehr weiter
helfen. Wenn man es ganz genau haben will. Besuchen Sie uns dann auf
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
Wir sagen Ihnen, wie es geht.
</Werbung>
Also jetzt zum versprochenen Update der Benchmark. Das
Benchmarkprogramm ermittelt die Anzahl der gültigen
dezimalen Stellen. Hier das Resultat:
Formel n! '0' '1' '2'
-----------------------------------------
Lanczos 00010! 2.4 3.8 4.7
Stirling 00010! 2.1 4.5 5.6
Gosper 00010! 4.2 5.5 7.7
Wehmeier 00010! 4.2 5.6 8.5
Cantrell 00010! 2.4 5.7 8.3
Rumpelstilz 00010! 8.2
Lanczos 00100! 3.4 5.7 7.6
Stirling 00100! 3.1 6.5 8.6
Gosper 00100! 6.2 8.5 11.9
Wehmeier 00100! 6.2 8.6 12.4
Cantrell 00100! 3.4 8.6 13.2
Rumpelstilz 00100! 13.1
Lanczos 01000! 4.4 7.7 10.6
Stirling 01000! 4.1 8.5 11.6
Gosper 01000! 8.2 11.4 15.9
Wehmeier 01000! 8.2 11.6 16.3
Cantrell 01000! 4.4 11.6 18.2
Rumpelstilz 01000! 18.1
Lanczos 10000! 5.4 9.7 13.6
Stirling 10000! 5.1 10.5 14.6
Gosper 10000! 10.2 14.4 19.9
Wehmeier 10000! 10.2 14.6 20.3
Cantrell 10000! 5.4 14.6 23.2
Rumpelstilz 10000! 23.1
Erste Analyse:
In der Kategorie '0' ist Gospers Formel einsame Spitze.
In der Kategorie '1' ist 'wehmeier1' und 'cantrell1'
auf gleicher Augenhöhe.
In der Kategorie '2' laufen 'cantrell2' und 'rumpelstilz'
Kopf an Kopf den anderen voraus. Hier müssen weitere
Überlegungen und/oder Kriterien den Sieger ermitteln.
Gruss Peter
Alfred Heiligenbrunner
Christopher Creutzig schrieb:
> ... emofehle ich wärmstens Lanzcos, A Precision
> Approximation of the Gamma Function, Journal of the Society for
> Industrial and Applied Mathematics: Series B, Numeric Analysis, Vol. 1
> (1964), p. 86-96,
> http://links.jstor.org/sici?sici-0887-459X%281964%291%3C86%3AAPAOTG%3E2.0.CO%3B2-F
^
.../sici?sici=...
Dann funktioniert's.
HP
> >HP schrieb:
>
> > [...] Diskussion zum Nachlesen empfehlen, die ich in
> > der "numerical analysis" Newsgroup vor ueber 4 Jahren mit David W.
> > Cantrell ueber dessen Naeherungen der Gammafunktion gefuehrt hab.
>
> Danke für diesen Link, den ich recht anregungsreich fand.
>
> Kommt Cantrell. Will genau das in Lanczos Formeln reparieren. Hier seine
> Formeln in den drei einfachsten Formen.
>
> cantrell0 := proc(n) local N; N:=n+1/2; kern(N) end:
> cantrell1 := proc(n) local N; N:=n+1/2; kern(N)/(1+1/(24*N-1/2)) end:
> cantrell2 := proc(n) local N; N:=n+1/2; kern(N)/(1+1/(24*N-1/2+1/((1440/2021)*N))) end:
>
Hallo Peter & Mathefreunde,
erst mal ganz grossen Dank fuer Deine phantastische Zusammenstellung.
Zu den Naeherungen von David W. Cantrell noch ein kleiner Nachtrag: Ich
hab ihn in einer privaten Mail (ca. vor zwei Jahren) mal gefragt, ob er
nicht die Zaehler und Nenner seiner Kettenbruchentwicklung als Folgen
in Neil Sloanes OEIS einreichen moechte (bevor ich das fuer ihn als
Mitglied des njas'schen Sammeltrupps) uebernehme. Hintergedanke war
natuerlich, David noch zu etwas ausfuehrlicheren Informationen darueber
zu verfuehren, wie er zu den Zahlen in der Kettenbruchentwicklung
gekommen ist. Mein letzter Informationsstand ist naemlich, dass er nie
das zugehoerige Verfahren veroeffentlicht hat.
Soweit ich die Antwort noch in Erinnerung habe, war das aber noch auf
einem hinteren Warteplatz in der Pipeline.
Vielleicht sollten wir mal versuchen, David ueber die aktuelle
Diskussionsfolge zu informieren und ihn dann nochmal auf die Herleitung
anzusprechen. Ich bin mir nicht sicher, aber eventuell kann David sogar
Deutsch lesen, was angesichts dessen, dass er als Saenger in
Bach-Oratorien aufgetreten ist, nicht auszuschliessen ist.
Viel Spass
Hugo Pfoertner
Peter Luschny
>> Peter Luschny wrote:
>>>> HP schrieb:
>>>> [...] Diskussion zum Nachlesen empfehlen, die ich in
>>>> der "numerical analysis" Newsgroup vor ueber 4 Jahren mit David W.
>>>> Cantrell ueber dessen Naeherungen der Gammafunktion gefuehrt hab.
>> Kommt Cantrell. Will genau das in Lanczos Formeln reparieren.
> Mein letzter Informationsstand ist naemlich, dass er nie
> das zugehoerige Verfahren veroeffentlicht hat.
Ich habe auch danach gesucht und nichts gefunden.
> Vielleicht sollten wir mal versuchen, David ueber die aktuelle
> Diskussionsfolge zu informieren und ihn dann nochmal auf die Herleitung
> anzusprechen.
Das ist keine schlechte Idee. Oder wir swappen den Thread in eine
englischsprachige Gruppe. Das würde ich dann gerne dir überlassen.
Ich versuche zu einer *ersten vorläufigen* Bewertung zu kommen.
Dabei setzte ich voraus, dass Cantrell einen Konvergenzbeweis
für sein Verfahren hat.
"Cantrells Formel sind die wichtigsten seit den
Stirlingschen und machen diese gleichzeitig obsolet."
Was meine ich genauer damit? Wir müssen bei den Näherungen ja zwei
Klassen unterscheiden, die asymptotischen und die konvergenten.
Die Stirlingsche Entwicklung ist eine asymptotische Formel.
Die Cantrellsche ist eine konvergente Formel.
Für viele wichtige numerische Größen gibt es ja beide Sorten
von Näherung, und sie können glücklich nebeneinander leben.
Will man eine einfache, 'schnelle' Näherung und begnügt sich
mit einer geringen Genauigkeit, so greift man zu der asymptotischen.
Will man eine genaue Näherung und scheut einen eventuellen Aufwand
bei der Berechnung nicht, greift man zu der konvergenten. (Grob
gesprochen.)
So sehe ich zum Beispiel das Verhältnis zwischen den
Stirlingschen und den Formeln von Lanczos.
Nun gibt es aber auch Extremfälle. Die zerfallen wieder in zwei Klassen.
In der einen Klasse sind die asymptotischen Näherungen, die für
alle Zwecke des Rechnens genau genug sind. Ein Beispiel dafür
ist die Hardy-Ramanujan-Rademacher Entwicklung der Anzahl der
Zahlenpartitionen. Jeder, der p(n) ausrechnen möchte, wird zu
der asymptotischen Formel von Hardy-Ramanujan greifen, weil sie
einfacher ist und genau genug um an den ganzzahligen Wert heran zureichen.
Die Rademacher Entwicklung ist dann nur mehr das ästhetische und
theoretische I-Tüpfelchen auf einer schier unglaublichen Formel.
Im anderen Extremfall sind die konvergenten Formeln so gut, dass
die Verwendung asymptotischer keinen Sinn macht. Dies scheint
mir in Bezug auf n! jetzt durch die Formeln von Cantrell erreicht
worden zu sein. Cantrells Formeln liefern bei vergleichbarem
Aufwand gleich gute oder bessere Näherungen als die Stirlingschen.
Also die Nachricht wäre:
Die Stirlingschen Formeln haben das Zeitliche gesegnet
und sind nur mehr von historischem Interesse.
Ich sage 'wäre', weil dies ausdrücklich eine erste vorläufige
Bewertung ist, es kann leicht sein, dass ich wesentliche
Einschränkungen der Formeln übersehen habe.
Wirklich klären wird man dies erst können, wenn Cantrell sein
Verfahren und seinen Konvergenzbeweis veröffentlicht hat.
Ihn da freundlich daran zu erinnern -- warum nicht?
Gruss Peter
P.S. Aktuelle Zusammenfassung der 'Benchmark' auf:
http://www.luschny.de/math/factorial/approx/SimpleCases.html
HP
>
> Ich versuche zu einer *ersten vorläufigen* Bewertung zu kommen.
> Dabei setzte ich voraus, dass Cantrell einen Konvergenzbeweis
> für sein Verfahren hat.
>
> "Cantrells Formel sind die wichtigsten seit den
> Stirlingschen und machen diese gleichzeitig obsolet."
>
> Im anderen Extremfall sind die konvergenten Formeln so gut, dass
> die Verwendung asymptotischer keinen Sinn macht. Dies scheint
> mir in Bezug auf n! jetzt durch die Formeln von Cantrell erreicht
> worden zu sein. Cantrells Formeln liefern bei vergleichbarem
> Aufwand gleich gute oder bessere Näherungen als die Stirlingschen.
>
> Also die Nachricht wäre:
>
> Die Stirlingschen Formeln haben das Zeitliche gesegnet
> und sind nur mehr von historischem Interesse.
>
> Ich sage 'wäre', weil dies ausdrücklich eine erste vorläufige
> Bewertung ist, es kann leicht sein, dass ich wesentliche
> Einschränkungen der Formeln übersehen habe.
>
> Wirklich klären wird man dies erst können, wenn Cantrell sein
> Verfahren und seinen Konvergenzbeweis veröffentlicht hat.
> Ihn da freundlich daran zu erinnern -- warum nicht?
>
> Gruss Peter
>
> P.S. Aktuelle Zusammenfassung der 'Benchmark' auf:
> http://www.luschny.de/math/factorial/approx/SimpleCases.html
Da ich im Augenblick hinsichtlich Netzkommunikation ziemlich
eingeschraenkt bin, werd ich erst in ein paar Tagen dazu kommen, mit
David W. Cantrell Kontakt aufzunehmen. Vorher geb ich aber gerne mal
meine ganz subjektive Bewertung dazu ab, was man so fuer den taeglichen
Gebrauch anstelle der angegrauten altehrwuerdigen Stirling-Formel
hernehmen koennte. Das beste Preis-Leistungsverhaeltnis gibt Gosper0,
denn das (+1/6) in der Wurzel ist ja wirklich billig und gut zu merken.
Wenn's genauer sein soll, dann ist mein Favorit Wehmeier1, denn das zum
Gosperschen +1/6 dazukommende +1/(72*n) find ich am leichtesten zu
merken, bei praktisch identischer Genauigkeit mit Cantrell.
Deshalb die Empfehlung: Gosper+Wehmeier statt Stirling in die
Formelsammlung!
Hugo Pfoertner
Stefan Wehmeier
so allmählich entsteht ja eine richtige kleine Sammlung ...
Peter Luschny wrote:
> >> Kommt Cantrell. Will genau das in Lanczos Formeln reparieren.
>
> > Mein letzter Informationsstand ist naemlich, dass er nie
> > das zugehoerige Verfahren veroeffentlicht hat.
>
> Ich habe auch danach gesucht und nichts gefunden.
>
> > Vielleicht sollten wir mal versuchen, David ueber die aktuelle
> > Diskussionsfolge zu informieren und ihn dann nochmal auf die Herleitung
> > anzusprechen.
>
Das Verfahren ist an sich klar, man entwickelt kern(n+1/2)/n!
(mit kern(n) = sqrt(2*Pi)*n^n*exp(-n) wie auf deiner Webseite)
per Stirling-Formel in eine Taylorreihe in 1/(n + 1/2). Von dieser Reihe
macht man dann "einfach" eine Kettenbruchentwicklung, d.h. man
subtrahiert den ganzrationalen Teil a*(n+1/2)+b, invertiert, subtrahiert den
ganzrationalen Teil , ....
> Also die Nachricht wäre:
>
> Die Stirlingschen Formeln haben das Zeitliche gesegnet
> und sind nur mehr von historischem Interesse.
>
> Ich sage 'wäre', weil dies ausdrücklich eine erste vorläufige
> Bewertung ist, es kann leicht sein, dass ich wesentliche
> Einschränkungen der Formeln übersehen habe.
ich habe kein genaues Gefühl für den Aufwand, d.h. das Koeffizientenwachstum
in den verschiedenen Reihen.
>
> P.S. Aktuelle Zusammenfassung der 'Benchmark' auf:
> http://www.luschny.de/math/factorial/approx/SimpleCases.html
bei dem Verfahren, das du nach mir benannt hast, waere es vielleicht
stilechter, jeweils einen zusätzlichen Koeffizienten im Argument der Wurzel
hinzuzufügen. D.h. Du machst
series(sqrt(n + a_0 + a_1/n + a_2/n^2 + ...)/sqrt(n), n=infinity)
setzt die entstandenen Taylorreihe in 1/n mit der Stirlingschen
für n!/kern(n)/n^(1/2) gleich und erhältst den Koeffizientenvergleich
a_0/2 = 1/6, also a_0 = 1/6
a_1/2 - a_0^2/8 = 1/288, also a_1/2 - 1/288 = 1/288, also a_1 = 0
1/16*a_0^3 - 1/4*a_1*a_0 + 1/2*a_2 = -139/51840, also a_2/2 + 1/3456 = -
139/51840, und daher a_2 = -77/12960,
und so weiter.
Jedenfalls freue ich mich schon auf dein zweibändiges Werk
"1000 Methoden, n! zu berechnen" :-)
Grüße,
--
Stefan Wehmeier
ste...@math.uni-paderborn.de
Peter Luschny
> Jedenfalls freue ich mich schon auf dein zweibändiges Werk
> "1000 Methoden, n! zu berechnen" :-)
Du irrst, das Paper wird genau 2 Seiten lang und heißt:
"Die beste jemals erdachte Methode 1*2*3*... zu berechnen."
Daneben stell ich dann einen großen Papierkorb für die
dann obsolet gewordenen 1000 anderen Methoden. :-)
Gruss Peter
Peter Luschny
> bei dem Verfahren, das du nach mir benannt hast, waere es vielleicht
> stilechter, jeweils einen zusätzlichen Koeffizienten im Argument der Wurzel
> hinzuzufügen. D.h. Du machst
> series(sqrt(n + a_0 + a_1/n + a_2/n^2 + ...)/sqrt(n), n=infinity)
> setzt die entstandenen Taylorreihe in 1/n mit der Stirlingschen
> für n!/kern(n)/n^(1/2) gleich und erhältst den Koeffizientenvergleich
> a_0/2 = 1/6, also a_0 = 1/6
> a_1/2 - a_0^2/8 = 1/288, also a_1/2 - 1/288 = 1/288, also a_1 = 0
> 1/16*a_0^3 - 1/4*a_1*a_0 + 1/2*a_2 = -139/51840, also a_2/2 + 1/3456 = -
> 139/51840, und daher a_2 = -77/12960,
> und so weiter.
So, jetzt habe ich mir das mal angeschaut.
Ich bekomme aber andere Werte als du.
> restart;
> solve(a/2=1/12,a):a:=%;
> solve(b/2-a^2/8=1/288,b): b:=%;
> solve(c/2-a*b/4+a^3/16=-139/51840,c): c:=%;
> solve(d/2-a*c/4-b^2/8+(3/16)*a^2*b-(5/128)*a^4=-571/2488320,d): d:=%;
a=1/6; b=1/72; c=-31/6480; d=-139/155520;
Damit werden die Näherungen zu:
> kern := proc(n) sqrt(2*Pi)*n^n*exp(-n) end:
> wehmeier0 := proc(n) kern(n)*sqrt(n+1/6) end:
> wehmeier1 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n) end:
> wehmeier2 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n-(31/6480)/n^2) end:
> wehmeier3 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n-(31/6480)/n^2-(139/155520)/n^3) end:
Ist es das, was du willst? In der Tabelle stünde dann
Wehmeier 00010! 4.2 5.6 7.4 8.1
Wehmeier 00100! 6.2 8.6 11.4 13.1
Wehmeier 01000! 8.2 11.6 15.4 18.1
Gruss Peter
Christopher Creutzig
> Stefan Wehmeier schrieb:
>
>
>>Jedenfalls freue ich mich schon auf dein zweibändiges Werk
>>"1000 Methoden, n! zu berechnen" :-)
>
>
> Du irrst, das Paper wird genau 2 Seiten lang und heißt:
> "Die beste jemals erdachte Methode 1*2*3*... zu berechnen."
Och, wie langweilig. Ich fand “Nineteen Dubious Ways to Compute the
Exponential of a Matrix” (und den Nachfolgeartikel mit dem Zusatz “25
years later”) immer schon alleine vom Title her genial. :-)
> Daneben stell ich dann einen großen Papierkorb für die
> dann obsolet gewordenen 1000 anderen Methoden. :-)
Und die Generatinen nach Dir sollen nichts aus den von Dir verworfenen
Fehlern lernen? Das ist ja gerade so, als würde ein Werk über oublic
key-Verschlüsselung die Ansätze mit knapsacks gar nicht erst erwähnen. :)
Gruß,
Christopher
Hugo Pfoertner
Ich hab auch noch mal etwas nachgeforscht, ob ich eine Quelle zu Gospers
Formel finden kann, aber leider bisher ohne Erfolg. David W. Cantrell
hat in einem sci.math-Beitrag vom Januar 2002 auch schon mal nach mehr
Information dazu gefragt, ohne eine Reaktion darauf zu bekommen. Siehe
http://groups.google.com/group/sci.math/msg/c6d2aa58ea7b7648
In diesem Beitrag gibt Cantrell auch eine umgeformte Darstellung von
"Cantrell1" an, die sich auch recht gut merken laesst:
n! is approx. Sqrt(2*Pi)*((n+1/2)/e)^(n+1/2)*(48*n+23)/(48*n+25)
Uebrigens an alle hier der Hinweis: Peter hat seine Tabelle um weitere
Naehrungsformeln (de Moivre, Henrici) erweitert:
http://www.luschny.de/math/factorial/approx/SimpleCases.html
Hugo
Peter Luschny
>>>Boudewijn Moonen schrieb:
>>> Der Link hat bei mir nicht funktionoert.
.. dann muss da wohl mal der Ingenioer hoer...
> Ich hab auch noch mal etwas nachgeforscht, ob ich eine Quelle zu Gospers
> Formel finden kann, aber leider bisher ohne Erfolg.
> In diesem Beitrag gibt Cantrell auch eine umgeformte Darstellung von
> "Cantrell1" an, die sich auch recht gut merken laesst:
>
> n! is approx. Sqrt(2*Pi)*((n+1/2)/e)^(n+1/2)*(48*n+23)/(48*n+25)
Nice :-)
Allerdings werde ich das in meine Systematik umschreiben, was
einfach notwendig ist, um den Überblick zu behalten und um in den
Verzweigungspunkten des virtuellen Formelbaums, indem die alle
stehen, die dahinter liegenden Ideen besser zu sehen.
cantrell1' :: N=n+1/2; kern(N)*(48*N-1)/(48*N+1)
> Uebrigens an alle hier der Hinweis: Peter hat seine Tabelle um weitere
> Naehrungsformeln (de Moivre, Henrici) erweitert:
> http://www.luschny.de/math/factorial/approx/SimpleCases.html
... und im Übrigen versuche ich ein bisschen den Advocatus Diaboli
zu spielen und die numerische Brauchbarkeit von Cantrells Verfahren
abzuklopfen.
Dazu habe ich mir mal eine Weltklasse-Bibliothek gekrallt, ISML,
http://www.vni.com/ und deren ausgereifte Implementierung gegen
meine quick&dirty&naive Implementierung von Cantrell laufen lassen.
Hier eine kleine Stichprobe:
n -> True
n -> Cantrell
n -> Isml
10 -> 3628800
10 -> 3628799,99999999 rel: -2,44249065417534E-15
10 -> 3628800,00000001 rel: 2,33146835171283E-15
20 -> 2,43290200817664E+18
20 -> 2,43290200817664E+18 rel: 8,88178419700125E-16
20 -> 2,43290200817664E+18 rel: 0
30 -> 2,65252859812191E+32
30 -> 2,65252859812193E+32 rel: 6,21724893790088E-15
30 -> 2,65252859812189E+32 rel: -6,21724893790088E-15
40 -> 8,15915283247898E+47
40 -> 8,15915283247896E+47 rel: -1,55431223447522E-15
40 -> 8,15915283247894E+47 rel: -4,44089209850063E-15
60 -> 8,32098711274139E+81
60 -> 8,32098711274136E+81 rel: -3,33066907387547E-15
60 -> 8,32098711274168E+81 rel: 3,50830475781549E-14
80 -> 7,15694570462638E+118
80 -> 7,15694570462656E+118 rel: 2,50910403565285E-14
80 -> 7,15694570462668E+118 rel: 4,14113188185183E-14
100 -> 9,33262154439442E+157
100 -> 9,33262154439457E+157 rel: 1,65423230669148E-14
100 -> 9,33262154439422E+157 rel: -2,04281036531029E-14
120 -> 6,68950291344913E+198
120 -> 6,68950291344869E+198 rel: -6,50590692430342E-14
120 -> 6,68950291344939E+198 rel: 3,95239396766556E-14
130 -> 6,46685548922047E+219
130 -> 6,46685548922050E+219 rel: 4,66293670342566E-15
130 -> 6,46685548922021E+219 rel: -4,06341627012807E-14
140 -> 1,34620124757175E+241
140 -> 1,34620124757188E+241 rel: 9,54791801177635E-14
140 -> 1,34620124757187E+241 rel: 8,75965966429249E-14
170 -> 7,25741561530800E+306
170 -> 7,25741561530810E+306 rel: 1,38777878078145E-14
170 -> 7,25741561530888E+306 rel: 1,21680443498917E-13
Resultat: Die bewegen sich auf gleicher Augenhöhe!
Den größten relativen Fehler bekomme ich bei n=140
und der ist kleiner als der von ISML bei n=170.
Wenn man sich dann aber auch noch den Quellcode
der ISML-Routine anschaut, dann weiß man, dass
der *viel* aufwändiger ist. Da wird eine Tschebyschew
Reihe ausgewertet, die locker 30 und mehr Terme hat!
Gruss Peter
Stefan Wehmeier
> Stefan Wehmeier schrieb:
>
>> bei dem Verfahren, das du nach mir benannt hast, waere es vielleicht
>> stilechter, jeweils einen zusätzlichen Koeffizienten im Argument der
>> Wurzel hinzuzufügen. D.h. Du machst
>> series(sqrt(n + a_0 + a_1/n + a_2/n^2 + ...)/sqrt(n), n=infinity)
>> setzt die entstandenen Taylorreihe in 1/n mit der Stirlingschen
>> für n!/kern(n)/n^(1/2) gleich und erhältst den Koeffizientenvergleich
>> a_0/2 = 1/6, also a_0 = 1/6
>> a_1/2 - a_0^2/8 = 1/288, also a_1/2 - 1/288 = 1/288, also a_1 = 0
oops .. peinlich, das kommt davon, wenn man doch einmal von Hand rechnet :-)
>> 1/16*a_0^3 - 1/4*a_1*a_0 + 1/2*a_2 = -139/51840, also a_2/2 + 1/3456 = -
>> 139/51840, und daher a_2 = -77/12960,
>> und so weiter.
>
> So, jetzt habe ich mir das mal angeschaut.
> Ich bekomme aber andere Werte als du.
>
>
> a=1/6; b=1/72; c=-31/6480; d=-139/155520;
>
ja, ist ok.
> Damit werden die Näherungen zu:
>
> > kern := proc(n) sqrt(2*Pi)*n^n*exp(-n) end:
> > wehmeier0 := proc(n) kern(n)*sqrt(n+1/6) end:
> > wehmeier1 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n) end:
> > wehmeier2 := proc(n) kern(n)*sqrt(n+1/6+(1/72)/n-(31/6480)/n^2) end:
> > wehmeier3 := proc(n)
> > kern(n)*sqrt(n+1/6+(1/72)/n-(31/6480)/n^2-(139/155520)/n^3) end:
>
> Ist es das, was du willst? In der Tabelle stünde dann
>
> Wehmeier 00010! 4.2 5.6 7.4 8.1
> Wehmeier 00100! 6.2 8.6 11.4 13.1
> Wehmeier 01000! 8.2 11.6 15.4 18.1
>
einverstanden !
Gruß,
Peter Luschny
> oops .. peinlich, das kommt davon, wenn man doch einmal von Hand rechnet :-)
Deshalb rechne ich alle meine Übungsaufgaben mit MuPad!
[Psst, eure Werbeabteilung darf mich gerne kontaktieren, gegen ein
kleines Honorar lass ich diesen Satz hier noch öfter fallen :-) ]
> einverstanden !
Das beruhigt mich sehr. Denn jetzt kann ich den Nachzünder
starten: Deine Formel passt auf die von Cantrell wie die
Faust auf's Auge.
peter1(n) :: (wehmeier1(n)+cantrell1(n))/2
Und hier das Resultat:
Stirling 00010! 2.1 4.5 5.6 7.8
De Moivre 00010! 2.4 5.0 5.7 7.8
Gosper 00010! 4.2 5.5 7.7
Wehmeier 00010! 4.2 5.6 7.4 8.1
Henrici 00010! 2.1 5.5 8.4 10.5
StirSeries 00010! 2.1 5.7 8.3 10.5
Cantrell 00010! 2.4 5.7 8.3 10.6
Peter 00010! 6.8
Stirling 00100! 3.1 6.5 8.6 11.7
De Moivre 00100! 3.4 7.1 8.6 12.0
Gosper 00100! 6.2 8.5 11.9
Wehmeier 00100! 6.2 8.6 11.4 13.1
Henrici 00100! 3.1 8.4 13.2 17.2
StirSeries 00100! 3.1 8.6 13.1 17.3
Cantrell 00100! 3.4 8.6 13.2 17.5
Peter 00100! 11.4
Stirling 01000! 4.1 8.5 11.6 15.6
De Moivre 01000! 4.4 9.1 11.6 16.0
Gosper 01000! 8.2 11.4 15.9
Wehmeier 01000! 8.2 11.6 15.4 18.1
Henrici 01000! 4.1 11.4 18.2 24.2
StirSeries 01000! 4.1 11.6 18.1 24.2
Cantrell 01000! 4.4 11.6 18.2 24.5
Peter 01000! 13.7
Stirling 10000! 5.1 10.5 14.6 19.6
De Moivre 10000! 5.4 11.1 14.6 20.0
Gosper 10000! 10.2 14.4 19.9
Wehmeier 10000! 10.2 14.6 19.3 23.1
Henrici 10000! 5.1 14.4 23.2 31.2
StirSeries 10000! 5.1 14.6 23.1 31.2
Cantrell 10000! 5.4 14.6 23.2 31.5
Peter 10000! 16.7
Gruss Peter
Hugo Pfoertner
[...]
Da ich an solider Numerik immer interessiert bin, geb ich hier auch mal
meine eigenen Ergebnisse eines Vergleichs des Cantrell-Verfahrens mit
einer "klassischen" Implementierung der Gammafunktion zum Besten. Als
Vergleichsreferenz nehm ich den Microsoft Rechner, der nach meiner
Erfahrung allemal fuer 30 richtige Stellen beim Berechnen von n!
brauchbar ist. Da ich keinen Zugriff auf einen IMSL-Quellcode hab, hab
ich wieder, wie schon in der damaligen Diskussionsfolge, auf den
Cody/Stoltz-Code aus
ALGORITHM 765, COLLECTED ALGORITHMS FROM ACM.
TRANSACTIONS ON MATHEMATICAL SOFTWARE,
VOL. 23, NO. 1, March, 1997, P. 81--90.
http://www.netlib.org/toms/765
zurueckgegriffen, aus dem ich mir dann die Gammafunktion rausgezogen
hab, Quelle auf
http://www.enginemonitoring.com/gamma/gamcw.for
Dieses Mal hab ich allerdings in Double precision gerechnet, um
Vergleichbares mit den oben von Peter dargestellten Ergebnissen zu
bekommen. Das Cantrell'sche Verfahren hab ich in der etwas
umgeschriebenen Form verwendet, wie sie im wesentlichen in
http://www.enginemonitoring.com/gamma/gamtry.for
dargestellt ist. Wiederum mit den erforderlichen einfachen Anpassungen,
um konsistent in Double Precision zu rechnen. Die wesentlichste
Aenderung zum Original ist wohl die Umrechnung des Kettenbruchs in eine
einfache rationale Funktion:
QCF(Z) = Z* ( 1.4034722222222222222*Z^2 + 2.1949405419284926536 ) /
(Z^2 * (Z^2 + 1.8214351694378563721) + 0.260008182821928029307 ),
die dann zweimal in der "HALF-SHIFT GAMMA EXPANSION"
HSG(Z) = sqrt(2*Pi)*(Z/e)^Z *
(48*Z + 2*QCF(Z) - 1) / (48*Z + 2*QCF(Z) + 1) verwendet wird.
Die Ergebnisse mit der Cody/Stolz-Funktion, in der im wesentlichen
rationale Approximationen verwendet werden (vermutlich nicht ganz so
bizarr wie die Tschebyschew-Reihe mit 30++ Termen):
n Naehrung fuer n! relativer Fehler
10 3628800.000000000 0.0
20 2.4329020081766400E+18 0.0
30 2.6525285981219316E+32 7.922115e-15
40 8.1591528324789411E+47 -4.442061e-15
60 8.3209871127416814E+81 3.500254e-14
80 7.1569457046266771+118 4.148006e-14
100 9.3326215443947553+157 3.643476e-14
120 6.6895029134497719+198 9.639616e-14
130 6.4668554892202120+219 -4.046364e-14
140 1.3462012475718704+241 8.760905e-14
170 7.2574156153080564+306 7.913644e-15
Zun Vergleich meine Implementierung von Cantrell:
10 3628800.000000000 0.0
20 2.4329020081766467E+18 2.753913e-15
30 2.6525285981219121E+32 5.706392e-16
40 8.1591528324790044E+47 3.316097e-15
60 8.3209871127413832E+81 -8.345496e-16
80 7.1569457046264137+118 4.676648e-15
100 9.3326215443944236+157 8.927642e-16
120 6.6895029134491675+198 6.045653e-15
130 6.4668554892204691+219 -7.070681e-16
140 1.3462012475717680+241 1.154316e-14
170 Ueberlauf eines Zwischenergebnisses; keine Lust zum Rumbasteln ohne
Auftrag und Bezahlung.
Da bei einigen der Cody/Stoltz-Ergebnisse fast genau die gleiche
Abweichung wie bei Peters IMSL-Ergebnissen rauskommt, liegt die
Vermutung nahe, dass die in IMSL verwendeten Approximationen eng mit
denen im TOMS 765-Programm verwendeten "verwandt" sind. Die
IMSL-Dokumentation schweigt sich bei der Gammafunktion leider ueber die
Herkunft des Verfahrens aus; zumindesten hab ich nichts gefunden.
Cantrell ist bei mir noch um eine Ecke genauer als in Peters
Implementierung. Allerdings hat das wohl eher mit
Operationsreihenfolgen, Rundung und aehnlichem Hexenwerk rund um die
hintersten Bits als mit der Guete der Approximationsfunktion zu tun. Mit
einigem Aufwand liesse sich da vermutlich noch mal eine halbe Stelle an
Genauigkeit rausquetschen.
Hugo
Hugo Pfoertner
>
[...]
>
> Das beruhigt mich sehr. Denn jetzt kann ich den Nachzünder
> starten: Deine Formel passt auf die von Cantrell wie die
> Faust auf's Auge.
>
> peter1(n) :: (wehmeier1(n)+cantrell1(n))/2
Da nehm ich aber doch lieber "wehmeier2(n)". Da brauch ich das "teure"
kern(x) nur fuer ein Argument auszurechnen und "w..2" ist immer gleich
gut oder besser als "peter1".
Hugo
Peter Luschny
>> Peter Luschny schrieb:
>> Das beruhigt mich sehr. Denn jetzt kann ich den Nachzünder
>> starten: Deine Formel passt auf die von Cantrell wie die
>> Faust auf's Auge.
>> peter1(n) :: (wehmeier1(n)+cantrell1(n))/2
> Da nehm ich aber doch lieber "wehmeier2(n)". Da brauch ich das "teure"
> kern(x) nur fuer ein Argument auszurechnen und ..
Aber mit diesem Argument kannst du natürlich meine
ganze kleine Käfersammlung totschlagen ;-) Das billigste
ist in jedem Fall einfach den Grad der Näherung zu erhöhen.
Mein nächstes Ziel, das sich hier schon andeutet, ist es,
einfache Inklusionen zu finden, um damit eine
Fehlerabschätzung mitberechnen zu können.
Gruss Peter
Peter Luschny
> Mit
> einigem Aufwand liesse sich da vermutlich noch mal eine halbe Stelle an
> Genauigkeit rausquetschen.
Aber locker, zumindest bei meiner.
Den Knackpunkt der ganzen Sache, den weder ich
noch du dabei angegangen sind, liegt in der
exp-log-Hintereinanderschaltung. Der relative
Fehler gleicht hier nämlich ungefähr dem absoluten
Fehler in log. Das heißt jeder Rundungsfehler auf
Maschinenebene wird hier dicke vergrößert. Will man
hier absolut sauber rechnen (und dabei übrigens
IMSL weit hinter sich lassen), müßte man an dieser
Stelle die Rechnung aufreißen und temporär mit
erweiterter Genauigkeit rechnen.
Gruss Peter
P.S. Dazu kann man etwa lesen Clark, Cody [Clark:1969:SCE],
wie in dieser Bibliographie zu finden.
http://www.math.utah.edu/pub/bibnet/authors/c/cody-william-j.ps.gz
Peter Luschny
> Mein nächstes Ziel, das sich hier schon andeutet, ist es,
> einfache Inklusionen zu finden, um damit eine
> Fehlerabschätzung mitberechnen zu können.
Done. Ich habe gerade eine weitere sehr gute Näherungsformel
hinzugefügt. In Ermangelung eines Namens habe ich sie
'MacLaurin series' genannt. Ob sie MacLaurin kannte, weiß ich
nicht, aber mit seiner "Treatise of Fluxions", 1742,
hatte er alle Mittel in der Hand um sie zu entwickeln.
Vielleicht kennt sie jemand hier und kann Genaueres sagen?
Mein Ziel der einfachen Inklusion ist damit auch erreicht,
und zwar systematisch (und nicht so wie gestern eklektisch).
Gruss Peter
http://www.luschny.de/math/factorial/approx/SimpleCases.html
HP
> Hugo Pfoertner schrieb:
>
> > Mit
> > einigem Aufwand liesse sich da vermutlich noch mal eine halbe Stelle an
> > Genauigkeit rausquetschen.
>
> Aber locker, zumindest bei meiner.
>
> Den Knackpunkt der ganzen Sache, den weder ich
> noch du dabei angegangen sind, liegt in der
> exp-log-Hintereinanderschaltung.
In meiner Implementierung verwende ich die in Fortran definierte
allgemeine Exponentiation x**y. Wie die genau in den einzelnen
Compilern bzw. den zugehoerigen Laufzeitbibliotheken implementiert ist,
weiss ich nicht. Nachdem sich ein erheblicher Teil der ehemals
legendaeren DigitalEquipment-Compiler-Entwickler inzwischen auf Umwegen
beim Riesen Intel befindet, der ja wohl langfristig alleine die
Massstaebe in Sachen Rechnerarithmetik bestimmt, waers sicher mal eine
nette Aufgabe, sowas wie die allgemeine x**y-Funktion hinsichtlich
Genauigkeit gegen eine hochpraezise Referenz (irgendein
Multiprecision-Paket) laufen zu lassen. Das gibts wahrscheinlich auch
schon laengst fix und fertig, aber aktuell hab ichs nicht auf dem
Schreibtisch.
> Der relative
> Fehler gleicht hier nämlich ungefähr dem absoluten
> Fehler in log. Das heißt jeder Rundungsfehler auf
> Maschinenebene wird hier dicke vergrößert. Will man
> hier absolut sauber rechnen (und dabei übrigens
> IMSL weit hinter sich lassen), müßte man an dieser
> Stelle die Rechnung aufreißen und temporär mit
> erweiterter Genauigkeit rechnen.
>
> Gruss Peter
>
> P.S. Dazu kann man etwa lesen Clark, Cody [Clark:1969:SCE],
> wie in dieser Bibliographie zu finden.
> http://www.math.utah.edu/pub/bibnet/authors/c/cody-william-j.ps.gz
Was die Cantrell'sche Naeherung tatsaechlich hergaebe, wenn nicht als
begrenzender Faktor die Rechengenauigkeit mit den ueblichen 64bit
Gleitkommazahlen auftraete, kann man sehen, wenn man mal probehalber
statt mit 64bit mit 128bit Gleitkomma rechnet. Ich hab mal mein
gamdwctry.for entsprechend umgeschrieben, wobei ich allerdings die
Konstanten mit 20 Stellen stehen gelassen habe.
n Cantrells Naeherung fuer n! relativer Fehler
10 3628799.99999999984620316249143284 -4.238229e-17
20 2432902008176640000.29759642492623 1.223216e-19
30 265252859812191058747943720061056. 4.208635e-19
40 8.159152832478977347919374182522262E+047 5.470251e-19
60 8.320987112741390150855630869585221E+081 7.906862e-19
80 7.156945704626380236881755096474495E+118 1.034045e-18
100 9.332621544394415280091396469829124E+157 1.277398e-18
120 6.689502913449127067761187387660069E+198 1.520751e-18
130 6.466855489220473683128647031844115E+219 1.642428e-18
140 1.346201247571752462962446729321380E+241 1.764104e-18
170 7.257415615307998982848737668696705E+306 2.129134e-18
Ob jetzt hier die erreichte Genauigkeit durch meine 20-stelligen
Koeffizienten oder durch die Approximationsguete der
Kettenbruchentwicklung bestimmt wird, hab ich nicht nachgeprueft.
P.S.: Gratulation zur wunderschoenen Kaefersammlung; ich trachte keinem
einzigen nach dem Leben. Allerdings ist es vermutlich ein schwieriges
Unterfangen, den tuechtigsten Exemplaren aus dem Luschnyschen Herbarium
einen Weg in die freie Wildbahn zu bahnen. Auch in 10 Jahren wird wohl
gebetsmuehlenhaft "Stirlingformel" empfohlen werden, wenns um eine
Naeherung fuer n! geht. Und in Bibliotheken wie der IMSL sind fast
alle Verfahren auf einem Stand so um die 1980 eingefroren, weil man ja
hauptsaechlich mit dem Verkauf Geld verdienen und dieses nicht an
irgendwelche weltfremden Mathematiker oder Numeriker weitergeben will,
die den ohnehin in der Praxis fast nie benoetigten letzten
Dezimalstellen nachjagen :-(.
Hugo Pfoertner
HP
>
> Was die Cantrell'sche Naeherung tatsaechlich hergaebe, wenn nicht als
> begrenzender Faktor die Rechengenauigkeit mit den ueblichen 64bit
> Gleitkommazahlen auftraete, kann man sehen, wenn man mal probehalber
> statt mit 64bit mit 128bit Gleitkomma rechnet. Ich hab mal mein
> gamdwctry.for entsprechend umgeschrieben, wobei ich allerdings die
> Konstanten mit 20 Stellen stehen gelassen habe.
>
> n Cantrells Naeherung fuer n! relativer Fehler
>
> 10 3628799.99999999984620316249143284 -4.238229e-17
> 20 2432902008176640000.29759642492623 1.223216e-19
> 30 265252859812191058747943720061056. 4.208635e-19
> 40 8.159152832478977347919374182522262E+047 5.470251e-19
> 60 8.320987112741390150855630869585221E+081 7.906862e-19
> 80 7.156945704626380236881755096474495E+118 1.034045e-18
> 100 9.332621544394415280091396469829124E+157 1.277398e-18
> 120 6.689502913449127067761187387660069E+198 1.520751e-18
> 130 6.466855489220473683128647031844115E+219 1.642428e-18
> 140 1.346201247571752462962446729321380E+241 1.764104e-18
> 170 7.257415615307998982848737668696705E+306 2.129134e-18
>
> Ob jetzt hier die erreichte Genauigkeit durch meine 20-stelligen
> Koeffizienten oder durch die Approximationsguete der
> Kettenbruchentwicklung bestimmt wird, hab ich nicht nachgeprueft.
>
http://groups.google.com/group/sci.math.num-analysis/msg/214bf7476f29e405
anstelle meiner 20stelligen (die fuer 64bit-Rechnung aber ausreichend
sind), bekomme ich folgende Tabelle:
n Cantrells Naeherung fuer n! relativer Fehler
10 3628799.99999999984554159889666405 -4.256459e-17
20 2432902008176639999.55802890305028 -1.816642e-19
30 265252859812191058635035581482148. -4.798812e-21
40 8.159152832478977343453553030560443E+0047 -3.137171e-22
60 8.320987112741390144276293649190500E+0081 -5.712547e-24
80 7.156945704626380229481151190086317E+0118 -3.049112e-25
100 9.332621544394415268169923603670575E+0157 -3.021195e-26
120 6.689502913449127057588118024143272E+0198 -4.476771e-27
130 6.466855489220473672507304383078972E+0219 -1.926514e-27
140 1.346201247571752460587607384709171E+0241 -8.803747e-28
170 7.257415615307998967396728210316661E+0306 -1.120596e-28
Da sieht man jetzt wirklich schoen, dass die Cantrell-Entwicklung fuer
grosse n immer besser wird. Wenn es in ferner Zukunft mal eine
Bibliothek fuer spezielle Funktionen mit 128bit Genauigkeit geben
sollte, dann waere Cantrells Entwicklung fuer die Gammafunktion dafuer
wirklich die allererste Wahl. Allerdings muesste man noch ein
Pochhammer-Polynom hoeheren Grades bei den kleinen Argumenten
verwenden, damit auch dieser Bereich genuegend genau behandelt wird.
Peter hat ja in seinem C#-Programm schon ein bisschen in diese Richtung
vorgehalten.
Der Vollstaendigkeit haeng ich ans Ende mal meine Fortranquelle rein,
damit sich das auch nach der Kompostierung aller persoenlichen
Webseiten noch nachvollziehen laesst.
Hugo Pfoertner
-----------------------------
REAL*16 FUNCTION GAMDWC ( X )
C
C APPROXIMATION OF THE GAMMA FUNCTION USING DAVID W. CANTRELL'S
C CONVERGENT EXPANSION POSTED IN NEWSGROUP SCI.MATH 2001-11-05
C "A new convergent expansion for the gamma function"
C
C Change history:
C 08.03.06 Quadruple precision (REAL*16)
C 07.11.01 ALL CONSTANTS WITH 20 DECIMAL DIGITS TO ENABLE
C CONVERSION TO DOUBLE PRECISION
C exponentiation now performed in double precision
C 06.11.01 INITIAL VERSION
C Author: Hugo Pfoertner
C Contact Info: http://www.pfoertner.org/
C
REAL*16 X
C
C LOCAL VARIABLES
REAL*16 SQPP, C1, C2, C3, C4, CF, HSG, Z, PP5, XT, DEINV
C
C CONSTANTS (higher precision required, see below)
C 1/e
C*** PARAMETER ( DEINV = 0.36787944117144232160D+0 )
C SQRT(2*PI)
C*** PARAMETER ( SQPP = 2.5066282746310005024E+0 )
C
C COEFFICIENTS FOR CONTINUED FRACTION
PARAMETER ( C1 = 1440Q0/2021Q0,
& C2 = 686186088Q0/125896643Q0,
& C3 = 1521596612992267104Q0/4596084813365743279Q0,
& C4 = 6.1441227298035761673076437188243880Q0/
& 2.0539143739435534417826656817767471Q0)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C STATEMENT FUNCTIONS
C
C CONTINUED FRACTION
CF(Z) = 1.0Q0 /
& ( C1 * Z + 1.0Q0 /
& ( C2 * Z + 1.0Q0 /
& ( C3 * Z + 1.0Q0 / ( C4 * Z ) ) ) )
C
C HALF-SHIFT GAMMA EXPANSION,
C EXPONENTIATION PERFORMED IN DOUBLE PRECISION
HSG (Z) = SQPP * (Z * DEINV )**Z /
& ( 1.0Q0 + 1.0Q0 /
& ( 24.0Q0*Z - 0.5Q0 + CF(Z) ) )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C
C Quadruple precision constants (computed to get full precision,
C should be moved to PARAMETER section)
SQPP = QSQRT(8.0Q0*QATAN(1.0Q0))
DEINV = 1.0Q0/QEXP(1.0Q0)
C
C POCHHAMMER POLYNOMIAL
PP5 = X * (X+1.0Q0) * (X+2.0Q0) * (X+3.0Q0) * (X+4.0Q0)
C
C RECURRENCE FORMULA
XT = X + 4.5Q0
GAMDWC = HSG ( XT ) / PP5
C
RETURN
C END OF FUNCTION GAMDWC
END
C**************************************
C Test program
REAL*16 GAMDWC
write(*,*) 10, GAMDWC(11.0Q0)
write(*,*) 20, GAMDWC(21.0Q0)
write(*,*) 30, GAMDWC(31.0Q0)
write(*,*) 40, GAMDWC(41.0Q0)
write(*,*) 60, GAMDWC(61.0Q0)
write(*,*) 80, GAMDWC(81.0Q0)
write(*,*) 100, GAMDWC(101.0Q0)
write(*,*) 120, GAMDWC(121.0Q0)
write(*,*) 130, GAMDWC(131.0Q0)
write(*,*) 140, GAMDWC(141.0Q0)
write(*,*) 170, GAMDWC(171.0Q0)
end
Peter Luschny
So, ich habe jetzt meine Käfersammlung mit der eigentlichen 'Serviceseite'
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm verlinkt.
> Allerdings ist es vermutlich ein schwieriges
> Unterfangen, den tuechtigsten Exemplaren aus dem Luschnyschen Herbarium
> einen Weg in die freie Wildbahn zu bahnen. Auch in 10 Jahren wird wohl
> gebetsmuehlenhaft "Stirlingformel" empfohlen werden, wenns um eine
> Naeherung fuer n! geht.
Sicher. 'Aber steter Tropfen...' und ... 'Das Bessere ist... '.
Was ich echt toll finde ist, dass ich von den drei Formeln, die ich heute
auf meiner Website empfehle, zwei noch nicht kannte, als ich diesen
Thread los trat, wobei eine davon wohl noch gar nicht existierte :-)
Außerdem habe ich das Gefühl gewonnen, dass man noch bessere als
Cantrell's Formeln finden kann. Aber das in einer anderen Runde.
Danke an alle die mitgemacht haben, besonders an Hugo und Stefan.
Gruss Peter
HP
Vielen Dank fuer die schoene Zusammenstellung; mir hat's Spass gemacht,
mich wieder mal ein bisschen mit dem Thema zu beschaeftigen. Vielleicht
noch eine kleine Ergaenzung zur Herkunft der "Windschitl"-Formel. Man
findet sie und einiges weitere zu Naeherungen der Gamma-Funktion auf
der Webseite "Calculators and the Gamma Function"
http://www.rskey.org/gamma.htm
Ein kleiner Ergaenzungswunsch zu Deiner Formelsammlung waere vielleicht
noch ein Hinweis darauf, dass man mit den meisten (allen?) n!-Formeln
fuer nicht ganzzahlige Argumente auch Naeherungen der Gammafunktion
ausrechnen kann. Das ist fuer mich sogar im realen Leben interessant,
weil naemlich der Mittelwert fuer die durch die
Wahrscheinlichkeitsdichtefunktion f(x)=s*x^(s-1)*exp(-x^s) gegebene
Normalform der Weibullverteilung gleich Gamma((s+1)/s) ist. Mit
Weibullverteilungen lassen sich Ueberlebensdauern von hoch belasteten
Bauteilen, aber auch z.B.die Schaedigungsverteilungen einer grossen
Zahl von Belastungsereignissen recht ordentlich wiedergeben. Die Werte
der Gammafunktion braucht man dann aber eher fuer Argumentwerte von 1
bis 5 als fuer 100 oder 1000.
Hugo Pfoertner
Peter Luschny
> Ein kleiner Ergaenzungswunsch zu Deiner Formelsammlung waere vielleicht
> noch ein Hinweis darauf [...]
Ich habe einen Abschnitt 'Comments' eingefügt und deinen Kommentar
leicht gekürzt dort eingefügt. Die Seite steht auch für weitere
Kommentare offen.
Gruss Peter
Hugo Pfoertner
So wie mein deutscher "Comment" auf der ansonsten englischsprachigen
Webseite jetzt dasteht, ist das natuerlich eine Falle, in die ein
Newcomer auf dem Gebiet leicht reintappen koennte. Fuer niedrige
Argumentwerte taugen die empfohlenen Naeherungen naemlich nicht viel.
Wenn man z.B. Gamma(1.5) (also sozusagen "0.5!") wissen moechte, dann
kommt bei der Gosper-Wehmeierschen Naeherung
k = sqrt(2 * Pi) * exp(z * (log( z ) - 1))
w = sqrt(z + 1 / 6 + 1 / (72 * z))
n! ~= k*w
ca.
2.506628274631*exp(0.5*(log(0.5)-1))*sqrt(0.5+1/6+1/(72*0.5))=0.89587
raus, was angesichts des Sollergebnisses von 0.8862270 (siehe
http://www.enginemonitoring.net/gamma/gtd.txt ) nur bescheidene
Genauigkeitswuensche erfuellt. Es ist aber nicht alles verloren. So wie
in den Programmen mit der Cantrell-Naeherung muss man mit einer
verschobenen Funktion arbeiten, die man dann durch Division mit einem
Korrekturpolynom z*(z-1)*(z-2)..... wieder zuruecktransformiert. Fuer
die Genauigkeiten, die ich fuer kleine Argumente so brauche, langt eine
Verschiebung um 2 oder 3.
Beispiel: Gamma(1.5) ~= 2.506628274631*exp(3.5*(log(3.5)-1))*
sqrt(3.5+1/6+1/(72*3.5))/(1.5*2.5*3.5)=0.886275
Hugo
Peter Luschny
> Peter Luschny schrieb:
>> HP schrieb:
>>> Ein kleiner Ergaenzungswunsch zu Deiner Formelsammlung waere vielleicht
>>> noch ein Hinweis darauf [...]
>> Ich habe einen Abschnitt 'Comments' eingefügt und deinen Kommentar
>> leicht gekürzt dort eingefügt. Die Seite steht auch für weitere
>> Kommentare offen.
> So wie mein deutscher "Comment" auf der ansonsten englischsprachigen
> Webseite jetzt dasteht, ist das natuerlich eine Falle, in die ein
> Newcomer auf dem Gebiet leicht reintappen koennte.
Angebot: Schicke mir den Text, den du haben möchtest, nicht länger
als 10 Zeilen, in welcher Sprache auch immer, per Email zu und
ich stelle das dort ein.
> Fuer niedrige
> Argumentwerte taugen die empfohlenen Naeherungen naemlich nicht viel.
> Wenn man z.B. Gamma(1.5) (also sozusagen "0.5!") wissen moechte, dann
> kommt bei der Gosper-Wehmeierschen Naeherung
> k = sqrt(2 * Pi) * exp(z * (log( z ) - 1))
> w = sqrt(z + 1 / 6 + 1 / (72 * z))
> n! ~= k*w
> ca.
> 2.506628274631*exp(0.5*(log(0.5)-1))*sqrt(0.5+1/6+1/(72*0.5))=0.89587
> raus, was angesichts des Sollergebnisses von 0.8862270 (siehe
> http://www.enginemonitoring.net/gamma/gtd.txt ) nur bescheidene
> Genauigkeitswuensche erfuellt. Es ist aber nicht alles verloren. So wie
> in den Programmen mit der Cantrell-Naeherung muss man mit einer
> verschobenen Funktion arbeiten, die man dann durch Division mit einem
> Korrekturpolynom z*(z-1)*(z-2)..... wieder zuruecktransformiert. Fuer
> die Genauigkeiten, die ich fuer kleine Argumente so brauche, langt eine
> Verschiebung um 2 oder 3.
>
> Beispiel: Gamma(1.5) ~= 2.506628274631*exp(3.5*(log(3.5)-1))*
> sqrt(3.5+1/6+1/(72*3.5))/(1.5*2.5*3.5)=0.886275
Ja, aber da gäbe es noch soviel Weiteres zu sagen - was ich alles
nicht will. Ich wollte nur meine Käfersammlung präsentieren.
Deinen hilfreichen Kommentar (sozusagen: Gamma fuer Dummies)
nehme ich aber gerne noch auf.
Gruss Peter
Hugo Pfoertner
>
[...]
Wenn man jetzt fleissig (und ein guter C-Programmier) waere, koennte man
getreu der Empfehlung "Hic Rhodus, hic saltus" bei der Implementierung
der math-Bibliothek von clc-wiki mithelfen und dort gleich der
Cantrell-Naeherung zu einiger Popularitaet verhelfen. Eine Aufforderung
dazu gibt's von einem Herrn Pietsch ganz aktuell in der sci.math-NG
http://groups.google.com/group/sci.math/msg/922e94d6d1ec634d
Die augenblickliche Quelle von (Log)Gamma steht auf
http://www.clc-wiki.net/wiki/Lgammaf.c
was natuerlich nur ein winziger Teil der benoetigten Funktionen ist:
http://www.clc-wiki.net/wiki/Category:Implementation_of_math.h
Hat jemand Lust, da mitzumachen?
Hugo
Hugo Pfoertner
>
[...]
> Wenn man jetzt fleissig (und ein guter C-Programmier) waere, koennte man
> getreu der Empfehlung "Hic Rhodus, hic saltus" bei der Implementierung
> der math-Bibliothek von clc-wiki mithelfen und dort gleich der
> Cantrell-Naeherung zu einiger Popularitaet verhelfen. Eine Aufforderung
> dazu gibt's von einem Herrn Pietsch ganz aktuell in der sci.math-NG
> http://groups.google.com/group/sci.math/msg/922e94d6d1ec634d
>
> Die augenblickliche Quelle von (Log)Gamma steht auf
> http://www.clc-wiki.net/wiki/Lgammaf.c
> was natuerlich nur ein winziger Teil der benoetigten Funktionen ist:
> http://www.clc-wiki.net/wiki/Category:Implementation_of_math.h
>
> Hat jemand Lust, da mitzumachen?
>
> Hugo
Da hat sich inzwischen was ziemlich Erschreckendes rausgestellt. Der
gute Gregory Pietsch hat sich vermutlich bei einem nicht genau bekannten
Teil der "math"-Funktionen allzu eng an die kommerzielle "The Standard C
Library" der Firma Dinkumware "angelehnt". Siehe
http://www.dinkumware.com/libdual.html , und dort steht in Fettdruck:
<<Beginn Zitat>>
Please note also that the Dinkum Unabridged Library is not a part of
Project Gnu. It is not covered by the Gnu Copyleft Agreement and it is
not shareware. It contains no Gnu code. Equally, Comeau Computing is not
involved in the development or distribution of this library. The Dinkum
Unabridged Library is a commercial product that you license from
Dinkumware, Ltd. It is protected by copyright and is subject to our
usual licensing restrictions.
<<Ende Zitat>>
Auf einen entsprechenden Hinweis
http://groups.google.com/group/sci.math/msg/fe09d19978e3499e von P.J.
Plauger, der zum Urgestein im Fundament der Sprache C zaehlt, siehe
http://www.plauger.com/books.html , hat jetzt der Administrator der
c-wiki den kompletten Inhalt der oben zitierten Seite gesperrt, siehe
http://groups.google.com/group/sci.math/msg/c4a3960050c7cb9c
Als jemand, der auch davon lebt, dass Arbeit, die man in Software
gesteckt hat, auch gelegentlich bezahlt wird, kann ich diesen
Schreckschuss nur begruessen. Jetzt koennten die clc-wiki-Leute aber
wahrscheinlich erst recht bessere Implementierungen der Gamma-Funktion
(und noch vieles mehr) brauchen ;-)
Hugo Pfoertner
Peter Luschny
>> Wenn man jetzt fleissig (und ein guter C-Programmier) waere, [...]
>> Hat jemand Lust, da mitzumachen?
'C'? >:-(( Unter Niveau.
> Jetzt koennten die clc-wiki-Leute aber
> wahrscheinlich erst recht bessere Implementierungen der Gamma-Funktion
> (und noch vieles mehr) brauchen ;-)
Ich glaube nicht, dass da Leute arbeiten, die sich vorher fragen,
was die Anzahl der Potenzreihenglieder ist, die man zur Bestimmung
einer vorgegebenen Stellenzahl braucht. Wildes 'hacken' hilft da
nicht weiter. Aber ohne solches Wissen ist die Arbeit eh für den Müll...
Gruss Peter
Peter Luschny
> Allerdings ist es vermutlich ein schwieriges
> Unterfangen, den tuechtigsten Exemplaren aus dem Luschnyschen Herbarium
> einen Weg in die freie Wildbahn zu bahnen. Auch in 10 Jahren wird wohl
> gebetsmuehlenhaft "Stirlingformel" empfohlen werden, wenns um eine
> Naeherung fuer n! geht.
Man muss anfangen! Und eine solche Empfehlung hat schon ein
bisschen Gewicht:
http://www.nist.gov/dads/HTML/stirlingappx.html
Gruss Peter
Peter Luschny
> So wie mein deutscher "Comment" auf der ansonsten englischsprachigen
> Webseite jetzt dasteht, ist das natuerlich eine Falle, in die ein
Ok, ich habe das jetzt wieder herausgenommen, aber deine anderen
Hinweise mit eingearbeitet.
> Newcomer auf dem Gebiet leicht reintappen koennte. Fuer niedrige
> Argumentwerte taugen die empfohlenen Naeherungen naemlich nicht viel.
> Wenn man z.B. Gamma(1.5) (also sozusagen "0.5!") wissen moechte, dann
> kommt bei der Gosper-Wehmeierschen Naeherung
> k = sqrt(2 * Pi) * exp(z * (log( z ) - 1))
g = sqrt(n + 1 / 6)
> w = sqrt(z + 1 / 6 + 1 / (72 * z))
> n! ~= k*w
> 2.506628274631*exp(0.5*(log(0.5)-1))*sqrt(0.5+1/6+1/(72*0.5))=0.89587
> raus, was angesichts des Sollergebnisses von 0.8862270 (siehe
> http://www.enginemonitoring.net/gamma/gtd.txt ) nur bescheidene
> Genauigkeitswuensche erfuellt.
Das ist überhaupt nicht schlimm, weil eine Approximation ein
Paar ist, bestehend aus einer Näherung an den gesuchten Wert
und einer Näherung an den Fehler.
Und so habe ich auch die Gosper-Wehmeier Näherung dargestellt.
Man schaut sich also in diesem Fall die Fehlerschätzung an
err = k * (w - g) / 2 = 0.009
relerr = (w - g) / (2 * w) = 0.01
und entscheidet dann, ob einem diese Genauigkeit genügt oder nicht.
> Es ist aber nicht alles verloren. So wie
> in den Programmen mit der Cantrell-Naeherung muss man mit einer
> verschobenen Funktion arbeiten, die man dann durch Division mit einem
> Korrekturpolynom z*(z-1)*(z-2)..... wieder zuruecktransformiert. Fuer
> die Genauigkeiten, die ich fuer kleine Argumente so brauche, langt eine
> Verschiebung um 2 oder 3.
>
> Beispiel: Gamma(1.5) ~= 2.506628274631*exp(3.5*(log(3.5)-1))*
> sqrt(3.5+1/6+1/(72*3.5))/(1.5*2.5*3.5)=0.886275
Yes, das ist eine von vielen Möglichkeiten um zu genaueren Näherungen
zu kommen. Eine andere wäre einfach auf eine Formel mit höherer
Approximationsordnung umzusteigen.
Aber ob das Resultat nun 'langt' für meine Genauigkeitsanforderung?
Ohne Fehlerschätzung kann ich es nicht wissen. Und mit Fehlerschätzung
ist es kein Problem, ob für kleine oder große Argumente.
Gruss Peter