The Gamma Function

[thornahawk]

I still need to write the documentation for this, so I thought some of you might want to test this for the time being:

[Weregoose]

I'm getting increment errors for any input for X not having a fraction part divisible by .5 or having an imaginary
part, and I'm also getting data type errors for lists. Some elaborate documentation would indeed be welcome.

Also, if I may: [font="courier new"]fnInt(T^(X-1)e^(‾T),T,0,E2

[thornahawk]

Goose, you may have mistyped. Download the attachment and see if that works. I've done (slightly) rigorous tests for random values in the complex plane, and I've no errors so far...

I'm familiar with your one-liner because 1) that's the defining integral due to Euler, and 2) it's in your site. No, the Lanczos approximation and the integral approximation are really different beasts, so I can't say anything about generalizing the integral, only that negative reals and complexes will not succumb to the treatment since fnInt( accepts only real inputs, and the integral is only valid for complex numbers with positive real part.

thornahawk

[Weregoose]

I downloaded the program, but ran into some errors again. The first one was because [font="courier new"]LGAMMA

[thornahawk]

Uh-oh... you need the list to be in the calculator as well. So that's why you were getting erroneous results. The GAMMA list and the program are both in the attached archive, and should not be separated . The GAMMA list included is the coefficient set I found to be most appropriate for TI calculators.

thornahawk

[Weregoose]

I see! I'm glad I read that article explaining the algorithm; it was a real eye-opener. Translating from the design of the Lanczos approximation, your program implements it quite nicely, and so there is no doubt that it works how it's supposed to. I'm starting to like this "hunt-and-kill" technique for developing math programs... I may have to join in the fun.

[thornahawk]

Let's just say I spend at most three weeks in research compared to at most four days in writing a math program. All that effort to make sure it works for all reasonable inputs. :D

thornahawk

P.S. For the interested, it's pronounced LAN-tsosh. Maybe the Hungarians here can confirm?

[Ray Kremer]

If you just want a gamma function (as opposed to wanting the experience of writing the program yourself):
http://www.technicalc.org/tiplist/en/files...tips/tip6_4.pdf

[thornahawk]

I must comment that the program Ray was pointing to uses the shifted Stirling approximation. Pretty accurate for large arguments, but not the best for small arguments. For about the same cost of storage, my professional opinion is that Lanczos has slightly higher accuracy than Stirling.

thornahawk

[Weregoose]

[font="courier new;font-size:9pt;line-height:100%;color:darkblue"](Ans√(Ans(e^(Ans‾¹)-e^(‾Ans‾¹))/2+Ans^‾6/810)/e)^Ans√(2[font="times new roman"]π

[thornahawk]

(repost of a reply made on June 2006, plus some additionals)

Windschitl's approximation was first mentioned in this page along with his inspired derivation for the approximation, FYI. However, there's a sinh() function available in the TI 83+ that could be used if one is to restrict him/herself to real inputs. Otherwise, it's fine the way Goose wrote it. :D

At any rate, the fact that this approximation is inaccurate for small arguments means that one might want to use the recursion Γ(x+1) == x*Γ(x) for a certain number of times before applying Windschitl's approximation. I have yet to experiment with the optimum number of recursions needed.

The Wolfram Functions Site provides for an online comparison as well:

http://functions.wolfram.com/webMathematic...i%20I&digits=10

which returns -0.1639600231 - 0.3810981009*i, close to the mark. ;)

(added 10/2007:

I came across this page by Peter Luschny showcasing a wide variety of simple approximations to the gamma function whenever only a few significant digits are required. One particularly striking approximation (at least to me) was the set of formulae found by Gergö Nemes; in particular, the formula referred to as "nemesgam(n)" is particularly compact and slightly more accurate than Windschitl's formula for real arguments.

I have yet to study the behavior of these approximations in the complex plane, and due to the fact that I am currently far away from computers that have Mathematica installed, I am not sure if the ~7 significant figure accuracy of "nemesgam(n)" (and more for large arguments; a bit like Stirling's approximation) carry over to the complex plane. Maybe somebody else interested in this stuff can look into it. ;)

)

thornahawk

[Weregoose]

Some test cases using Γ(Z) = √(2[font="times new roman"]π/Z)Z‾¹[font="times new roman"]^×

[thornahawk]

Hmm, looks good. As it stands, it, like most other simple approximations, does not do well for small magnitudes, but performs better for large magnitudes. What I think would be most illuminating would be a three-dimensional plot of |(Γ(x+iy)-nemesgam(x+iy-1)/Γ(x+iy)| for a variety of appropriate ranges. ;)

Anyway, for those who felt too lazy to stare at the web page I gave earlier, here is the formula Goose and I were looking at:



thornahawk

[DarkerLine]

My skills at illuminating plots in Mathematica clearly leave much to be desired. Nevertheless, here are two plots: one with the default range and box parameters, the other my best attempt at a closeup of the weird stuff near the origin.

[Weregoose]

I went at it again, and I must say that I am very happy this time with the results. The new program I made handles real-numbered inputs, both positive and negative, large and small, and returns answers quickly and with brilliant accuracy – the worst case has 8 correct digits; an input of n or greater, where Γ(n) = 100, suffices to get 12. In terms of its alacrity and precision, this is now virtually as good as just another function on the calculator! :lol:

See the Gamma Function here.

[spandiv]

That looks pretty nice and complicated In any case, I noticed you can simplify

:If 2≤abs(A
:Then
:sum(12‾¹/{Ans,‾30Ans³,105Ans^5,‾140Ans^7,99Ans^9
:Else
:2fnInt(tan‾¹(T/Ans)/(e^(2πT)-1),T,0,4
:End

to

:sum(12‾¹/{Ans,‾30Ans³,105Ans^5,‾140Ans^7,99Ans^9
:If 2>abs(A
:2fnInt(tan‾¹(T//(e^(2πT)-1),T,0,4

unless there is a reason for keeping the If-Then conditional.

[DarkerLine]

Well, there's the obvious speed reason - you don't want to do both calculations when you can only do one.

[Weregoose]

In fact, the entire program could be reduced to just the integral and the operations following it – speed was definitely the target here.

[Weregoose]

To compute Γ‾¹(x) for positive reals greater than one:

‾1+solve(Ans-([font="arial"]√(2[font="times new roman"]π

[DarkerLine]

If you wanted to, you could call it arcΓ.