Back to Riemann, Part 16


Back to the gamma function. There are a few more VBScript functions I have to add to my file. First, complex exponentials. From Part 13 of the blog we have:

‘ n^x = (n^a)) * (cos(b * ln(n)) – i*sin(b * ln(n))
‘ for a scalar value to a complex power
Function complexSPower(cPn, c)
dim cplx(2)                         ‘ where n = cPn is real only

cplx(re) =  cPn^(c(re)) * cos(c(im) * log(cPn))
cplx(im) =  cPn^(c(re)) * sin(c(im) * log(cPn))
complexSPower = cplx
end Function

Then, reciprocals, again from part 13, multiplying the top and bottom by (a – bi):

‘ 1/(a + bi) = (a – bi)/(a^2 – b^2)
Function complexRecip(c)
dim cplx(2)

cplx(re) =    c(re) / (c(re)^2 + c(im)^2)
cplx(im) = – (c(im) / (c(re)^2 + c(im)^2))
complexRecip = cplx
end Function

And this is where I ran into a road block with taking exponential complex numbers (i.e. – (a + bi)^(c + di)). The above code for complexSPower() takes a complex exponent of a scalar, and that’s fine. But, there’s one part of the gamma approximation that requires the power of a complex number, and I had problems finding other ways to approach this.

The issue revolves around taking the log of a complex number (ln(a + bi)). Remembering that n^x = (n^a)) * (cos(b * ln(n)) – i*sin(b * ln(n)), substituting c + di for n sticks us with (c + di)^a, and ln(c + di). Yes, I can get this far, but…

The issue is with the fact that e^pi*b*i describes a unit circle. As b increases and represents multiple integer values of pi, the circle repeats. So, the inverse function, the natural logarithm, ln(x), has an infinite number of solutions. There are several approaches to solve for ln(), but I had difficulty understanding them, or translating them into VBScript functions.

Ignoring this for the moment, let’s look at the wiki article Lanczos approximation Python code in detail (which I’m learning from scratch as I go along…)

“def gamma(z):”

This is the start of a function definition. In VBScript, this is:

Function gamma(zG)
z = zG

gamma = result
End Function

Then we have:

“epsilon = 0.0000001
def withinepsilon(x):
return abs(x) <= epsilon”

This defines a constant, called epsilon. It’s followed by inline code defining a function for testing whether the absolute value of the imaginary part of z is less than epsilon. If it is, the imaginary part will be discarded later on. In VBScript, this is:

epsilon       = 0.0000001
Function withinepsilon(we)
withinepsilon = abs(we) <= epsilon
end Function”

Then we get the import statement:

“from cmath import sin,sqrt,pi,exp”

which doesn’t have a VBScript equivalent.

And,

“p = [ 676.5203681218851,   -1259.1392167224028,  771.32342877765313,
-176.61502916214059,     12.507343278686905, -0.13857109526572012,
9.9843695780195716e-6, 1.5056327351493116e-7]
z = complex(z)”

Here, we’re creating an array of the pre-created constants for Ag(z), and a re-cast that forces the z passed to the gamma(z) function to a complex number. My version in VBScript is:

p = Array( 676.5203681218851,        -1259.1392167224028,          771.32342877765313, _
-176.61502916214059,          12.507343278686905,         -0.13857109526572012, _
0.0000099843695780195716,     0.00000015056327351493116)
z     = complex(0.6, 0.4)

Where I’m just picking 0.6 and 0.4 arbitrarily as test numbers for debugging the program.

Now, the reflection formula is tied to the fact that the Lanczos approximation only holds for z >= 0.5. If z < 0.5, which is true for the entire left side of the complex plane for z(re) < 0, then we can apply a conversion formula,

pi / (sin(pi*z) * gamma(1-z))

Python supports recursion, where you can call a function from within itself. VBScript can’t do this. Instead, I have to break the code up into two parts, and set a flag:

isReflect = False
pi        = 3.1415926535897932384

if complexReal(z) < 0.5 then
isReflect = True
z         = complexProd(z, -1)
z         = complexSum(1, z)
end if

What I’m doing is using my complex product and summation functions to calculate z = 1 – z.
(I could also do z = complexSum(1, complexProd(-1, z)).)
I then apply the following gamma code to the new value. Afterward, I use something like:

if( Not isReflect ) then
gamma     = result3
else
reflSin   = complexSin(complexProd(pi, zG))
reflect1  = complexProd(reflSin, result3)
reflRecip = complexRecip(reflect1)
reflect2  = complexProd(pi, reflRecip)
gamma     = reflect2
end if

Or,

gamma = complexProd(pi, complexRecip(complexProd(complexSin(complexProd(pi, zG)), result3)))

Where result3 is the value calculated for gamma(1 – z).
I’ll discuss this in the next blog entry.

The rest of the Python code continues in a fairly straightforward manner:

“z -= 1
x = 0.99999999999980993”

Ok, fine. In VBScript,

z = complexSum(z, -1)
x = complex(0.99999999999980993, 0)

Then comes a kind of tricky part:

“for (i, pval) in enumerate(p):
x += pval/(z+i+1)”

Python is stepping through the p list one element at a time to build up the value for Ag(z), adding each new result to x. I don’t really like the Python use of “i” as a variable in the for-loop because I keep confusing it with i as sqrt(-1). I changed it to cntr. I can do this as:

for cntr = 0 to ubound(p)
tSum     = complexSum(z, cntr + 1)
tRecip   = complexRecip(tSum)
tProduct = complexProd(p(cntr), tRecip)
x        = complexSum(x, tProduct)
next

The reason for breaking each part of the equation out and giving them unique names is that it makes things easier for debugging, and printing out the values of each step if needed.
I.e. – (within the for-loop)

if(1 = 0) then
wscript.echo “cntr                  = ” & cntr
wscript.echo “(z + cntr + 1) (re)   = ” & tSum(re)
wscript.echo “(z + cntr + 1) (im)   = ” & tSum(im)
wscript.echo “1/tSum (re)           = ” & tRecip(re)
wscript.echo “1/tSum (im)           = ” & tRecip(im)
wscript.echo “p(cntr)               = ” & p(cntr)
wscript.echo “tRecip * p(cntr) (re) = ” & tProduct(re)
wscript.echo “tRecip * p(cntr) (im) = ” & tProduct(im)
wscript.echo “x += tTemp (re)       = ” & x(re)
wscript.echo “x += tTemp (im)       = ” & x(im)
end if

Or, in one line, I’d use:

for cntr = 0 to ubound(p)
x      = complexSum(x, complexProd(p(cntr), complexRecip(complexSum(z, cntr + 1))))
next

Now, another part of the Python code that I’m not happy with is the use of len(p) on the list to get the total number of elements inside it.

“t = z + len(p) – 0.5”

This introduces a new glitch, because the VBscript equivalent, ubound(), returns the highest subscript of the array, not the total number of elements in it. ubound(p) = 7, and the Python len(p) = 8. I created a new function to address this.

Function lenAry(lA)
lenAry = ubound(lA) + 1
End Function

Which gives me:
t = complexSum(z, lenAry(p) – 0.5)

Eventually, I’m going to need -t, so I do it this way:

minusT  = complexProd(-1, t)

The real guts of the Lanczos approximation is this part here:

“result = sqrt(2*pi) * t**(z+0.5) * exp(-t) * x”

Let me skip over that for a second and look at the last part of the Python code:

“if withinepsilon(result.imag):
return result.real
return result”

All this says is that if the imaginary part of the calculation results is really small, punt and say that gamma is just comprised of the real part. Otherwise, return the calculation result. My version is:

if( withinepsilon(complexImg(result)) ) then
result(im) = 0
end if

gamma = result
End Function

Ok, the gamma guts.

“result = sqrt(2*pi) * t**(z+0.5) * exp(-t) * x”

I need to break this down into pieces.

From what I’ve learned about Python, exp(-t) is the same as e^(-t).
VBscript doesn’t know what “PI” and “e” are, so I have to create them as constants.
Also, since sqrt(2*pi) doesn’t change, I store this as a constant, too.

pi = 3.1415926535897932384
e  = exp(1)
sqrt2PI   = (2*pi)^0.5

And I can build up the entire product with:

t       = complexSum(z, lenAry(p) – 0.5)
t2z05   = complexPPower(t, complexSum(z, 0.5))

minusT  = complexProd(-1, t)
e2nT   = complexSPower(e, minusT)

result1 = complexProd(t2z05, sqrt2PI)
result2 = complexProd(result1, x)
result3 = complexProd(result2, complexSPower(e, minusT))

Which leaves me with “t**(z+0.5)”, which is “t^(z+0.5)”, and would be the same as:

temp = complexSum(z, 0.5)
t2z05 = t^temp

And this is where I bog down because t is a complex number, and I don’t know how to take ln(t).

As I said, the code is getting clunky already, and adding a function for taking the powers of a complex number is only going to make it worse.

The other thing is that I’m not that interested in the gamma function by itself. All I really cared about was being able to graph the zeta function along the critical line, and I’ve done that. I can graph gamma on the wolframalpha.com site, and the graph doesn’t look all that interesting from an animation viewpoint.  This isn’t sour grapes, per se. I’m just not all that compelled to spend another 2-3 days dealing with complex logarithms. But, who knows. I’m keeping the code I have now, and if I get an answer I like, I may finish it some day.

Advertisements
Leave a comment

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: