# Back to Riemann, Part 15

There’s just one final thing I want to address in this series, and that’s how to implement the gamma function.

The Wiki article starts out with the following definition:
“Find a smooth curve that connects the points (x, y) given by y = (x – 1)! at the positive integer values for x.”

This is related to the factorial curve that I presented in part 9. While the answer can be found through integration and the use of limits, one of the more commonly used approaches is through the gamma function.

(Image from the wiki page.)

This formula isn’t very easy to use, so the next step is to look at the Lanczos approximation, published by Cornelius Lanczos in 1964.

(Image from the wiki page.)

where Ag(z) is:

(Image from the wiki page.)

then, g is a constant we choose such that the real part of z > 1/2. The pn coefficients can then be pre-calculated for a given g, using the formula in the wiki article. I’m not going to keep reproducing all these equations here, because then I’d have no room left for my part of the writing.

But, if you pre-calculate the constants, Ag(z) becomes:

(Image from the wiki page.)

Which leads me to the next part, which is the Python code example for calculating gamma. I’ve never used Python, and there are a couple parts that I had trouble understanding with just a cursory glance at the code. But, it seems that Python has a full implementation for complex number handling, which is nice. The thing is, I prefer VBScript because it’s (normally) easy to use, it comes with Windows, it’s free, the implementation doesn’t need to be downloaded from somewhere, I can run it directly from the command line, and I know how to tie it to Excel for drawing graphs. The downside is that VBScript has almost nothing in the way of math support, can’t handle complex numbers, and doesn’t support making classes.

In fact, VBScript is almost useless for calculating gamma, but I’m looking at this as a recreational math project, which could be an interesting learning experience. Plus, I’m almost done with the code now, so there’s little reason to quit at this stage. (Ha ha.)

First things first. We need a way to represent complex numbers. I’m going to use a simple 2-element array, and my own function naming conventions.

‘ Create two constants to act as array index constants,
‘ re for the real component, im for the imaginary one.
re = 0
im = 1

‘ Create a 2D array to hold the complex number: c(re) = a, c(im) = bi
Function complex(a, b)
dim cplx(2)

cplx(re) = a
cplx(im) = b
complex = cplx
end Function

‘ Return the real component of the complex number.
Function complexReal(c)
complexReal = c(re)
end Function

‘ Return the imaginary component of the complex number.
Function complexImg(c)
complexImg = c(im)
end Function

‘ Return the magnitude of the complex number (sqrt(a^2 + b^2))
Function complexMag(c)
complexMag = (c(re)^2 + c(im)^2)^0.5
end Function

Then, I want a few simple functions for handling multiplication and addition of complex numbers. Initially, I planned on having separate functions for pure complex numbers, and for mixed scalars with complex numbers. And, I was going to recast scalars as complex numbers in-line when needed. But that started getting really clunky and hard to read. So I figured I might as well put a test for whether a variable was an array within the function, and if it wasn’t, then I’d re-cast it. In my version of VBScript, vartype(array) returns 8204.

‘ Test if a variable is an array, and if not, return an array
Function makeComplex(mCc)
dim cplx(2)

if(vartype(mCc) <> 8204) then
cplx(re) = mCc
cplx(im) = 0
makeComplex = cplx
else
makeComplex = mCc
end if
end Function

‘ a1 + a2 + i*(b1 + b2)
Function complexSum(c1, c2)
dim cplx(2)

tc1 = makeComplex(c1)
tc2 = makeComplex(c2)

cplx(re) = tc1(re) + tc2(re)
cplx(im) = tc1(im) + tc2(im)
complexSum = cplx
end Function

‘ (a1 + b1*i) * (a2 + b2*i)
Function complexProd(c1, c2)
dim cplx(2)

tc1 = makeComplex(c1)
tc2 = makeComplex(c2)

cplx(re) = tc1(re) * tc2(re) – tc1(im) * tc2(im)
cplx(im) = tc1(re) * tc2(im) + tc1(im) * tc2(re)
complexProd = cplx
end Function

At some point, we’re going to need to take the sine of a complex number. While VBScript supports sin(), cos(), exp() (which is actually e^x) and log() (actually the natural log = ln()), they handle real numbers only. Fortunately, this is where the Internet steps in to save us. PlanetMath.org shows the first part.

sin(x+i*y) = sin(x)*cosh(y) + i*cos(x)*sinh(y)
cos(x+i*y) = cos(x)*cosh(y) – i*sin(x)*sinh(y)

But, VBScript doesn’t support hyperbolic sines and cosines. However, vbsedit.com has a page of Derived Math Functions, which is EXACTLY what I need here.

sinh(X) = (Exp(X) – Exp(-X)) / 2
cosh(X) = (Exp(X) + Exp(-X)) / 2

‘ sinh(a) = (e^a – e^(-a)) / 2
Function sinh(sa)
sinh = (exp(sa) – exp(-sa))/2
end Function

‘ cosh(a) = (e^a + e^(-a)) / 2
Function cosh(sa)
cosh = (exp(sa) + exp(-sa))/2
end Function

‘ sin(a+bi) = sin(a)*cosh(b) + i*cos(a)*sinh(b)
Function complexSin(c)
dim cplx(2)

cplx(re) = sin(c(re)) * cosh(c(im))
cplx(im) = cos(c(re)) * sinh(c(im))
complexSin = cplx
end Function

‘ cos(a+bi) = cos(a)*cosh(b) – i*sin(a)*sinh(b)
Function complexCos(c)
dim cplx(2)

cplx(re) =   cos(c(re)) * cosh(c(im))
cplx(im) = – sin(c(re)) * sinh(c(im))
complexCos = cplx
end Function

I’ll end here, because I’ll be getting into how Python handles for-loops on arrays, and there’s one bit I’m having trouble understanding, so I need to figure that out first.