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

A new convergent expansion for the gamma function

129 views
Skip to first unread message

David W. Cantrell

unread,
Nov 5, 2001, 9:31:21 AM11/5/01
to
Summary: A convergent expansion for the gamma function -- different
from any expansion of Lanczos, and presumably new -- is presented.
It converges rapidly when Re(z) is large, and is thus well suited for
computation of the gamma function.

This article supersedes my articles from October 25 in the thread
<http://mathforum.org/epigone/sci.math/gayveephal>,
written naively, i.e., before I became aware of the work of C. Lanczos
in "A precision approximation of the gamma function", J. SIAM Numer.
Anal., Ser. B, Vol. 1 (1964), pp. 86-96. He gives convergent (not
merely asymptotic) expansions. The first such, given in his equation (21),
states that z! =

Sqrt(2*pi)*(z+1/2)^(z+1/2)*e^(-(z+1/2))*[1-(1/24)*(1/(z+1))
-(23/1152)*(1/((z+1)*(z+2)))-(11237/414720)*(1/((z+1)*(z+2)*(z+3)))-...]

which -- quite beautifully, I think -- takes into account contributions
from all the poles. But Lanczos does not stop there; since "the convergence
is too slow to be of any practical value", he proceeds to get more rapidly
convergent expansions, too messy to present here, which are frequently
used for computation of the gamma function nowadays.

Despite the slow convergence of his (21), I found it attractive, more so
than his other expansions, which have a feeling of "practical compromise"
about them. This fondness for his (21) led me to seek an expansion of
similar "purity" but with improved convergence. The result, the main one
of this article, is an expansion, convergent if Re(z) > 0, involving a
continued fraction: Gamma(z + 1/2), or "Half-Shift Gamma", is

(*) HSG(z) = Sqrt(2*pi)*(z/e)^z / [1 + 1/( 24*z - 1/2 + CF(z) )]

where CF(z) = 1/(c_1*z + 1/(c_2*z + 1/(c_3*z + ...))) with c_n rational:

c_1 = 1440/2021, c_2 = 686186088/125896643,
c_3 = 1521596612992267104/4596084813365743279, ...
or approximately
c_1 = 0.71251855517070757050, c_2 = 5.4503922554948506450,
c_3 = 0.33106364977586039145, c_4 = 2.9914210678640645176, ...

If the continued fraction is terminated at any point, the quantity in
the brackets above then becomes a rational function, which may be thought
of as a Padé approximant at infinity.

It is natural to ask now, esp. since only four coefficients have been
given (and their derivation is somewhat messy), what accuracy is
achieved when the continued fraction is replaced by one of its
convergents, and how this accuracy compares with that achieved by
some well known method of computation. As our standard for
comparison, the most accurate of the approximations of Lanczos (at
the top of his p. 95) will be used, for which |relative error| in
approximating Gamma(z) is claimed to be less than approx. 2*10^(-10)
when Re(z) > 1. Let AHSG_m(z) denote the approximation of HSG(z)
obtained by using 1/(c_1*z + 1/(c_2*z + 1/(c_3*z + ...+1/(c_m*z)))),
the (m+1)th convergent to CF(z), instead of the entire continued fraction.
It will be obvious from the few entries below, using convenient real
arguments, that my approximation is poor when x is small but excellent
when x is large. But this is to be expected: Lanczos designed his
approximations to work well throughout the half-plane Re(z) > 1,
whereas I designed mine to work well there only when |z| is large.

|Relative errors|
x Gamma(x) Lanczos' approx. AHSG_4(x - 1/2)
1 0! = 1 4*10^(-11) 10^(-3)
8 7! = 5040 2*10^(-10) 10^(-13)
64 63! 10^(-10) 8*10^(-24)
512 511! 2*10^(-11) 9*10^(-34)

[For those suspicious of the last entry, perhaps it should be noted that,
besides using a great number of digits of accuracy in a computer
algebra system, the values of the rational coefficients were specified
precisely.]

Using more of the continued fraction, better results can be obtained, of
course; nevertheless, when |z| is small, convergence is slow. Were we
to stop here, we would at least have the convergent expansion (*), from
which approximations for Gamma(z + 1/2) can be obtained easily, giving
exceptional accuracy in the right half-plane when |z| is large. But, by
simply employing the general recurrence formula
Gamma(z) = Gamma(z+n) / (z)_n, where (z)_n denotes the Pochhammer
polynomial z*(z+1)*(z+2)*...*(z+(n-1)), a trick sometimes used formerly
with Stirling's asymptotic expansion, we can obtain a result which
approximates Gamma(z) well throughout the half-plane Re(z) > 1.
[Only this half-plane need be considered since the reflection formula
Gamma(1-z)*Gamma(1+z) = pi*z/sin(pi*z)
allows Gamma(1-z) to be obtained from Gamma(1+z).]
In general, we may consider approximating Gamma(z) by

AG_(m,n)(z) = AHSG_m(z + n + 1/2) / (z)_n.

If we wish the accuracy when |z| is small to be comparable to that
afforded by the previously mentioned approximation of Lanczos, it is
reasonable to choose n = 5 (together with m = 4, as before).

|Relative errors|
x Gamma(x) Lanczos' approx. AG_(4,5)(x)
1 0! = 1 4*10^(-11) 3*10^(-12)
8 7! = 5040 2*10^(-10) 4*10^(-16)
64 63! 10^(-10) 3*10^(-24)
512 511! 2*10^(-11) 8*10^(-34)

[By the way, the reason it seems appropriate to use n = 5 when
comparing with this particular approximation of Lanczos is that the
expansions from which the approximations are derived then represent
the gamma function in the same region, namely, the half-plane
Re(z) > -9/2.]

Clearly we could have used somewhat smaller values of m or n, while
still having an accuracy comparable to that of Lanczos' approximation
when |z| is small and yet providing far greater accuracy when |z| is
large. It remains to be determined exactly which values of m and n
should be used to provide a specified level of accuracy at minimal
computational cost, but this is not a difficult task.

Three miscellaneous thoughts:

1) Perhaps it should be noted, although probably just as a curiosity,
that, based on (*), a representation of the gamma function can be given
which is valid throughout the complex plane:

Gamma(z) = limit_( n -> oo ) [ HSG(z + n + 1/2) / (z)_n ]

This is much like the last result in Lanczos' paper, about which he
comments: "The proper interpretation of this peculiar limit law ...
requires further investigation."

2) As I noted recently in the previously mentioned thread, the simple
approximations for N! based on Stirling's expansion, found in so many
references, should be replaced by better results. Instead of
N! ~ Sqrt(2*pi*N)*(N/e)^N, it would be better to give
N! ~ Sqrt(2*pi)*((N+1/2)/e)^(N+1/2). Similarly, instead of just
N! ~ Sqrt(2*pi*N)*(N/e)^N*(1+1/(12*x)), perhaps it would be nice to note
that, for N >= 3, N!/(Sqrt(2*pi)*((N+1/2)/e)^(N+1/2)) lies between
(48*N+23)/(48*N+25) and (24*N+12)/(24*N+13).
[By the way, it is very easy to use (*) to produce much tighter (but, of
course, more complicated) bounds than those above.]

3) 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),
so to speak. Nonetheless:
Curséd be those who would give serious thought to yet another
"normalization"! Things are messy enough already. Quoting from
Lanczos' opening paragraph: "... Gamma(n+1) = n!
The normalization of the gamma function to Gamma(n+1) instead of Gamma(n)
is due to Legendre and void of any rationality. This unfortunate
circumstance compels us to utilize the notation z! instead of Gamma(z+1)."

Assuming that expansion (*) is indeed new, the many details which
have not been supplied here will be given in a paper, yet to be written.
In conclusion: The important virtue of expansion (*) is its rapidity of
convergence in the right half-plane when |z| is large.

Questions and comments are welcome.

David W. Cantrell
DWCan...@sigmaxi.org

--
-------------------- http://NewsReader.Com/ --------------------
Usenet Newsgroup Service

Hugo Pfoertner

unread,
Nov 5, 2001, 5:02:29 PM11/5/01
to
David,

when I read your contribution on the Gamma function, an article by
John L. Spouge:
Computation of the Gamma, Digamma, and Trigamma Functions
SIAM J. Numer. Anal. Vol. 31, No. 3, pp. 931-944, June 1994
came to my mind. Among other results that might not be too far from
your's, 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.
>> end citation John L. Spouge.
I always wondered that the simple formula (13) seemed to stay widely
unknown. Maybe your "re-discovery" will contribute to more popularity
;-).

Best Regards
Hugo Pfoertner


"David W. Cantrell" wrote:
>
> Summary: A convergent expansion for the gamma function -- different
> from any expansion of Lanczos, and presumably new -- is presented.
> It converges rapidly when Re(z) is large, and is thus well suited for
> computation of the gamma function.

<<snip>>

>
> 2) As I noted recently in the previously mentioned thread, the simple
> approximations for N! based on Stirling's expansion, found in so many
> references, should be replaced by better results. Instead of
> N! ~ Sqrt(2*pi*N)*(N/e)^N, it would be better to give
> N! ~ Sqrt(2*pi)*((N+1/2)/e)^(N+1/2). Similarly, instead of just
> N! ~ Sqrt(2*pi*N)*(N/e)^N*(1+1/(12*x)), perhaps it would be nice to note
> that, for N >= 3, N!/(Sqrt(2*pi)*((N+1/2)/e)^(N+1/2)) lies between
> (48*N+23)/(48*N+25) and (24*N+12)/(24*N+13).
> [By the way, it is very easy to use (*) to produce much tighter (but, of
> course, more complicated) bounds than those above.]

<<snip>>

David W. Cantrell

unread,
Nov 5, 2001, 6:10:32 PM11/5/01
to
Hugo Pfoertner <hugo.pf...@talknet.de> wrote:
> David,
>
> when I read your contribution on the Gamma function, an article by
> John L. Spouge:
> Computation of the Gamma, Digamma, and Trigamma Functions
> SIAM J. Numer. Anal. Vol. 31, No. 3, pp. 931-944, June 1994
> came to my mind.

Hi Hugo,
Of course. I read his article the same day that I became aware of
Lanczos' work.

> Among other results that might not be too far from your's,

As far as I know, none of his results are comparable to (*), my
main result.

> 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.
> >> end citation John L. Spouge.
> I always wondered that the simple formula (13) seemed to stay widely
> unknown. Maybe your "re-discovery" will contribute to more popularity
> ;-).

Yes, it's hard to understand why this seems "to stay widely unknown".
After all, the corresponding approximation of Stirling, quite well
known, is much less accurate, while also being computationally slightly
more complex. My sci.math article of Oct. 25 at
<http://mathforum.org/epigone/sci.math/glonphyfle>
provides another reason why formula (13) is nice: it allows an
approximation for the principal branch of the inverse gamma function
to be obtained in terms of the Lambert W function.

Best regards,
David

Hugo Pfoertner

unread,
Nov 6, 2001, 8:55:47 PM11/6/01
to
David,

I have to admit that my comparison of your method with the Spouge
article was not careful enough - sorry. Since I am always interested in
numerical results, I have written a little Fortran program using a naive
implementation of your expansion (for real arguments):
http://www.enginemonitoring.com/gamma/gamdwc.for
http://www.enginemonitoring.com/gamma/gtest.for
I have run a comparison against the results of the Cody/Stoltz gamma
function implementation in ACM algo 765
http://www.netlib.org/toms/765 , (extracted gamma function for
convenience at
http://www.enginemonitoring.com/gamma/gamcw.for )
which uses rational approximations with a very high accuracy (20 decimal
digits). I have checked the 765-results against a high precision
calculation of gamma using "Reduce" and found that the results were
accurate to +-1 in the 8th decimal digit using single precision (Watcom
Fortran on PC). The comparison with my implementation of your method
gives the following results:
# x gamma 765 gamma DWC relative err
0.1000000 9.5135078 9.5135088 1.0024424E-007
0.5000000 1.7724539 1.7724540 6.7256636E-008
1.0000000 1.0000000 1.0000001 1.1920929E-007
1.5000000 0.8862270 0.8862271 1.3451327E-007
2.0000000 1.0000000 1.0000001 1.1920929E-007
2.5000000 1.3293403 1.3293406 1.7935105E-007
3.0000000 2.0000000 2.0000005 2.3841858E-007
3.5000000 3.3233509 3.3233516 2.1522125E-007
4.0000000 6.0000000 6.0000014 2.3841858E-007
4.5000000 11.6317282 11.6317310 2.4596713E-007
5.0000000 24.0000000 24.0000057 2.3841858E-007
5.5000000 52.3427773 52.3427925 2.9151661E-007
6.0000000 120.0000000 120.0000305 2.5431316E-007
6.5000000 287.8852844 287.8853455 2.1201207E-007
7.0000000 720.0000000 720.0001831 2.5431316E-007
7.5000000 1871.2542725 1871.2548828 3.2617243E-007
8.0000000 5040.0000000 5040.0014648 2.9064361E-007
8.5000000 14034.4072266 14034.4121094 3.4791725E-007
9.0000000 40320.0000000 40320.0117188 2.9064361E-007
9.5000000 119292.4609375 119292.5000000 3.2745154E-007
10.0000000 362880.0000000 362880.1250000 3.4446649E-007
11.0000000 3.6288000E+006 3.6288015E+006 4.1335980E-007
12.0000000 3.9916800E+007 3.9916816E+007 4.0083373E-007
13.0000000 4.7900160E+008 4.7900182E+008 4.6763935E-007
14.0000000 6.2270208E+009 6.2270239E+009 4.9333380E-007
15.0000000 8.7178289E+010 8.7178330E+010 4.6984175E-007
20.0000000 1.2164510E+017 1.2164517E+017 5.6491774E-007

The full table is at
http://www.enginemonitoring.com/gamma/gt.txt
plotted:
http://www.enginemonitoring.com/gamma/gt.pdf

Question: Where is the accuracy lost (according to your table the
relative error of AG_(4,5)(8) should be 4*10^(-16)? How should I
re-write my Fortran function to get a comparable accuracy to the
Cody/Stoltz gamma function (which also uses only single precision)?

BTW: My only application for the Gamma function is to determine the mean
value of the Weibull distribution from its shape parameter and 3 decimal
figures are usually sufficient :-)

Thanks and Best Regards,

Hugo Pfoertner

Charles R. Bond

unread,
Nov 6, 2001, 9:43:48 PM11/6/01
to
> Questions and comments are welcome.
>
> David W. Cantrell
> DWCan...@sigmaxi.org
>

David,

I am very interested in your results because I have been using Lanczos
method for some time. My application involves the gamma function
in an extended precision environment, so I wondered if you plan to
publish the algorithms for computing the coefficients for your
approximation. Any timeline for this?

Thanks..

Stuart L. Anderson

unread,
Nov 6, 2001, 11:16:38 PM11/6/01
to
In sci.math Hugo Pfoertner <hugo.pf...@talknet.de> wrote:

> Question: Where is the accuracy lost (according to your table the
> relative error of AG_(4,5)(8) should be 4*10^(-16)? How should I
> re-write my Fortran function to get a comparable accuracy to the
> Cody/Stoltz gamma function (which also uses only single precision)?

On most systems, the Fortran type "REAL" is 4-bytes or about 6 decimal
digits. Thus, replace "REAL" by "DOUBLE PRECISION" which is probably
about 15 digits. Your 20 digit constants will still be REAL unless you
append "D0" (without the quotes) to those that have no exponent. If
there is an "E" exponent, change it to "D". There may be other problems,
but this will do for a start.

If you are on a Cray or using a precision-doubling compiler switch, let me
know. That would change my interpretation of your code.

--Stu
___________________________________________________________________
Stu Anderson stua...@drizzle.com Renton, Washington, USA

David W. Cantrell

unread,
Nov 6, 2001, 11:53:34 PM11/6/01
to
Pre-script: I just noticed that Alan Miller has posted a response,
but have not had an opportunity to look at it yet. He may explain
matters better than I have below.

Hugo Pfoertner <hugo.pf...@talknet.de> wrote:
> I have to admit that my comparison of your method with the Spouge
> article was not careful enough - sorry.

No need to be sorry. I'm delighted you've taken a serious interest.

> Since I am always interested in
> numerical results, I have written a little Fortran program using a naive
> implementation of your expansion (for real arguments):
> http://www.enginemonitoring.com/gamma/gamdwc.for
> http://www.enginemonitoring.com/gamma/gtest.for
> I have run a comparison against the results of the Cody/Stoltz gamma
> function implementation in ACM algo 765
> http://www.netlib.org/toms/765 , (extracted gamma function for
> convenience at
> http://www.enginemonitoring.com/gamma/gamcw.for )
> which uses rational approximations with a very high accuracy (20 decimal
> digits). I have checked the 765-results against a high precision
> calculation of gamma using "Reduce" and found that the results were
> accurate to +-1 in the 8th decimal digit using single precision (Watcom
> Fortran on PC). The comparison with my implementation of your method
> gives the following results:
> # x gamma 765 gamma DWC relative err

[snip]
> 8.0000000 5040.0000000 5040.0014648 2.9064361E-007
[snip]


> Question: Where is the accuracy lost (according to your table the
> relative error of AG_(4,5)(8) should be 4*10^(-16)? How should I
> re-write my Fortran function to get a comparable accuracy to the
> Cody/Stoltz gamma function (which also uses only single precision)?

Yes, where is the accuracy lost? (It's at this point that I wish I had
some background in numerical analysis!) Clearly something is "going wrong",
but it seems just as clear that you have done a perfect naive
implementation of my AG_(4,5). I snipped all but one entry of your table;
if we can fix the problem there, it is to be hoped that it will then be
fixed elsewhere as well.

All of my numerical work had been done in computer algebra systems
which allow arbitrary precision. Paul Godfrey and I have discussed
matters a bit since yesterday; in light of those discussions, I will need
to make a few more comments -- eventually. (Beginning tomorrow
morning, I shall be computer incommunicado for at least two days.)
Anyway, things are quite different when dealing with finite precision
arithmetic. AG_(4,5) (8) should be about 5039.999999999997751254,
giving a _negative_ relative error, the absolute value of which I reported
correctly as 4*10^(-16). To check your implementation (and also force
finite precision on myself), I "walked through" it using my beloved old
HP15C (which has a nice gamma, by the way), and got 5040.000006
instead. I must wonder why this differs so much from your value above.
But it is obvious that something needs to be done so that adequate
accuracy can be obtained using finite precision arithmetic.

Thinking that the trouble lay perhaps in calculating (Z * EINV)**Z, I
simply rewrote things in terms of Ln(Gamma) first, and then used
Gamma = Exp(Ln(Gamma)) at the end. At least for my little
calculator, this had a good effect: It now gave 5039.999995. Much
better (but no cigar yet). So, Hugo, here's my naive rewrite for your
program. (I don't remember how to notate the natural log and exponential
functions in FORTRAN, so I'll just say LN and EXP here.)

Change your parameter SQPP for SQRT(2*PI) to a parameter P for
LN(2*PI)/2.

Replace HSG(Z) by
LHSG(Z) = P + Z*( LN(Z) -1 ) - LN(1+1/(24*Z-0.5+CF(Z)))

and then say
GAMDWC = EXP( LHSG(XT) - LN(PP5) ).

I'll be interested to see (in a few days, when I can get back to the
computer) if my suggestion -- still quite naive, of course -- has a good
effect for your implementation.

In the mean time,
Best wishes,

David W. Cantrell

unread,
Nov 6, 2001, 11:59:48 PM11/6/01
to
"Charles R. Bond" <cb...@ix.netcom.com> wrote:
> I am very interested in your results because I have been using Lanczos
> method for some time. My application involves the gamma function
> in an extended precision environment, so I wondered if you plan to
> publish the algorithms for computing the coefficients for your
> approximation.

Yes. It's a little involved, but I've automatted their precise rational
computation in Mathematica (using, in part, its ability to deal with the
Lambert W function).

> Any timeline for this?

Not yet on my timeline, sorry.

Regards,
David

HP

unread,
Nov 7, 2001, 6:08:47 AM11/7/01
to
Stuart,

thanks for your suggestion. Having worked on many types of machines
including some Crays for 32 years and 30 years of Fortran I am aware
of the relationship between using "E" and "D" in the exponent part of
constants and their representation. The reason for choosing the "E"
exponent in my program was to be compatible with algo 765 in its
single precision (4 byte real) implementation. The netlib source uses
CD or CS comments to mark the differences between single and double
precision and I un-commented the CS path. The algo 765 gamma function
does not use any double precision operations and produces an accuracy
near the theoretical limit given by the 4 byte constants (the
constants are the same for SP and DP with the exception of the
exponents E or D; 24 decimal digits covers all known architectures).
Of course double precision would produce a better accuracy, but my
intent was to understand the reason for the higher accuracy loss of my
implementation of David's method in comparison to the Cody/Stoltz
method using single precision (4 byte) calculation. I will try David's
suggestion when I am back from work this evening.

Best Regards,
Hugo Pfoertner
"Stuart L. Anderson" <stua...@drizzle.com> wrote in message news:<10051065...@yabetcha.sttl.drizzle.com>...

Hugo Pfoertner

unread,
Nov 7, 2001, 6:09:20 PM11/7/01
to
Alan,

thanks for your quick conversion of my program to F90 and the successful
test. David's method looks rather competitive. Being a little
old-fashioned I continued to try getting my single precision F77 fixed
and I eventually confirmed David's idea (the one decisive operation is
the exponentiation (z/e)^z ). See my reply to Davids post.

Thanks and Best Regards
Hugo Pfoertner

Alan Miller wrote:
>
> David,
> Thanks for making your code available.
> I suspect that you did not convert all of the parameter values to
> double precision.
> I found that the relative errors decreased steadily from 0.2E-10 at x = 0.1
> to 0.9E-14 at x = 5.0. Beyond that the relative errors start to oscillate
> as
> their true values fall below the rounding errors of the computer.
> Detailed results and program attached.
>
> --
> Alan Miller (Honorary Research Fellow, CSIRO Mathematical
> & Information Sciences)
> http://www.ozemail.com.au/~milleraj
> http://users.bigpond.net.au/amiller/
>
> Name: gamtest.out
> gamtest.out Type: Ohne Angabe (application/octet-stream)
> Encoding: x-uuencode
>
> Name: gamtest.f90
> gamtest.f90 Type: Ohne Angabe (application/octet-stream)
> Encoding: x-uuencode

Hugo Pfoertner

unread,
Nov 7, 2001, 6:14:09 PM11/7/01
to
David, Alan, Stuart,

there was no need to use double precision with one exception: David,
your idea that the exponentiation was the root cause of the problem was
100% ok. Instead of following your suggestion to use log and exp I
converted the calculation of (z/e)^z to double precision and left
everything else unchanged. See:
http://www.enginemonitoring.com/gamma/gamdwc2.for
The single precision results are now perfect, absolutely identical
results with your method and the Cody/Stoltz function. See:
http://www.enginemonitoring.com/gamma/gtd.txt
Before hitting this point I tried several modifications of your formulas
(e.g. to avoid the multiple divisions in the CF part). I have collected
all modifications in a separate Fortran routine:
http://www.enginemonitoring.com/gamma/gamtry.for
The results are identical with those of gamdwc2.for.
With some optimizations + safeguarding this function might be the basis
for a very efficient real gamma function implementation.
The coefficients in QCF(Z) were computed from your c_1,2,3,4's by
CF(z) := (z^3*c2*c3*c4 + z*(c2 + c4))/
(z^4*c1*c2*c3*c4 + z^2*(c1*c2 + c1*c4 + c3*c4) + 1)

Best Regards
Hugo Pfoertner

"David W. Cantrell" wrote:
>
> Pre-script: I just noticed that Alan Miller has posted a response,
> but have not had an opportunity to look at it yet. He may explain
> matters better than I have below.
>

Charles R. Bond

unread,
Nov 8, 2001, 9:23:29 AM11/8/01
to
Referring to my previous post, if you aren't certain about the
timeline for publishing your paper, would you mind posting
a few more (rational) coefficients in this newsgroup?

HP

unread,
Nov 11, 2001, 3:49:48 PM11/11/01
to
As Alan's post did not appear in the Google Groups archive (the post
is cited below) I have copied his two attachments to my
"enginemonitoring" page":
http://www.enginemonitoring.com/gamma/gamtest.f90 (Fortran 90 program)
and
http://www.enginemonitoring.com/gamma/gamtest.out (results of
comparison).

Regards
Hugo Pfoertner

>
> Alan Miller wrote:
> >
> > David,
> > Thanks for making your code available.
> > I suspect that you did not convert all of the parameter values to
> > double precision.
> > I found that the relative errors decreased steadily from 0.2E-10 at x = 0.1
> > to 0.9E-14 at x = 5.0. Beyond that the relative errors start to oscillate
> > as
> > their true values fall below the rounding errors of the computer.
> > Detailed results and program attached.
> >
> > --
> > Alan Miller (Honorary Research Fellow, CSIRO Mathematical
> > & Information Sciences)
> > http://www.ozemail.com.au/~milleraj
> > http://users.bigpond.net.au/amiller/
> >
> > Name: gamtest.out
> > gamtest.out Type: Ohne Angabe (application/octet-stream)
> > Encoding: x-uuencode

copy at: http://www.enginemonitoring.com/gamma/gamtest.out


> >
> > Name: gamtest.f90
> > gamtest.f90 Type: Ohne Angabe (application/octet-stream)
> > Encoding: x-uuencode

copy at: http://www.enginemonitoring.com/gamma/gamtest.f90

David W. Cantrell

unread,
Nov 13, 2001, 1:16:58 PM11/13/01
to
"Charles R. Bond" <cb...@ix.netcom.com> wrote:
> Referring to my previous post, if you aren't certain about the
> timeline for publishing your paper, would you mind posting
> a few more (rational) coefficients in this newsgroup?

Here's a list giving some coefficients c_n of CF(z) precisely.
It begins with c_1 and gives a few more than in my original article.

1440/2021,

686186088/125896643,

1521596612992267104/4596084813365743279,

61441227298035761673076437188243880/20539143739435534417826656817767471,

33216277034690456269201306591096663890958682442526052832/
154187684682287395130815676867766056654304274786409523983,

11032328226229460455855489633846011632863413789579121881600576134625626179
1133352536/
5375805591444238830052560265723765535361323652899078901434006830761123339
6794963869,

93031537512178442352094227069108760175265654144645574344023598901195682455
976853094770288729984214865280760719488767200/
58223518103369713001005282669819397550373262406557960675177252534536427896
5643722001124189437614592415486167547436141091,

33628368492593938606366221375475198529098028011657512026195575252177033834
63789811469717761279027851911640015172414646160598778507378894223951063478
4112163780424/
21567802201843997098160440415859375060513241460254317252146375741915711493
88439223194745212628922513941813781515823899145401025683605871576587831900
3377591156979,

crude decimal approximations of which are
0.712519, 5.45039, 0.331064, 2.99142, 0.215428, 2.05222, 0.159783,
and 1.55919.

I regret the delay in responding to your request, but I've had very
little time available for mathematics recently. I still plan to post
further general comments about my expansion.

Regards,
David Cantrell

Charles R. Bond

unread,
Nov 16, 2001, 12:00:59 PM11/16/01
to
David W. Cantrell wrote:

> I regret the delay in responding to your request, but I've had very
> little time available for mathematics recently. I still plan to post
> further general comments about my expansion.
>
> Regards,
> David Cantrell

David,

Thanks for taking the time to generate the numbers -- and
I would like to express my regret for taking so long to
post my appreciation.

0 new messages