Back to Riemann, Part 17


Where did I leave off? Oh yeah. Despair, despair. Woe is me, I can’t figure out complex logs because the code is too complex.

Yeah, well, whatever. Turns out that if you dig enough, and use the right magic search words, you can eventually find anything on the net (maybe. I haven’t tried ALL the magic search words yet.)

What I did find pointed me back to one specific part of the wiki page on complex logs, which had this little bit of info:

log(a + bi) = ln(sqrt(a^2 + b^2)) + i*atan2(b, a)

Ok, this is good, this is good. Shame that VBScript doesn’t support atan2. atn, yes, but not atan2.

So, what is atan2 now?

Turns out that when you take the arc tangent of a complex number, say (-1 – i*1), the result is the same as for 1 + i*1. This is a problem. It’s caused by the calculation:

b/a

which wipes out the signs of the numerator and denominator if they are both negative. Plus you get a divide by zero error if a is 0.

atan2 is a function implemented in many computer languages to “correctly” calculate atan.

I found an implementation for atan2 for vbscript at www.paulsadowski.com/wsh/suntimes.htm, which in turn was lifted from a techtips page. Unfortunately, this bit of code added the result of atn() to pi, and in some cases the final answer was off by that value of pi. If you plot the results on a graph, the extra rotation gets you into the correct quadrant, but when feeding the numbers into the function for complexLog(), it screwed things up. Also, the answers weren’t corresponding to the numbers from wolfram alpha. So, I made a couple tweaks to the code to align it with wolfram.

Function atan2(ys, xs)
‘ from the Draw Arrows tip
‘ @ http://www.devx.com/upload/free/features/VBPJ/techtips/techtips11.pdf
‘ by A-Jim Deutch, Syracuse, New York
‘ I got this from http://www.paulsadowski.com/wsh/suntimes.htm – CHH

Dim theta
If xs <> 0 Then
theta = Atn(ys / xs)
If xs < 0 and ys > 0 Then
theta = theta + pi
Elseif xs < 0 and ys < 0 Then
theta = theta – pi
End If
Else
If ys < 0 Then
theta =  (-1) * pi / 2  ‘ 270
Else
theta = pi / 2          ‘ 90
End If
End If
atan2 = theta

End Function

Then, complexLog() is simply:

Function complexLog(c)
dim cplx(2)

cplx(re) = log(complexMag(c))
cplx(im) = atan2(c(im), c(re))
complexLog = cplx
End Function

From where I left off,

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

“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)
minusT  = complexProd(-1, t)
t2z05   = complexPPower(t, complexSum(z, 0.5))
result1 = complexProd(t2z05, sqrt2PI)
result2 = complexProd(result1, x)
result3 = complexProd(result2, complexSPower(e, minusT))

Which leaves me with “t**(z+0.5)”, which is the Python version of “t^(z+0.5),” which I could handle as:

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

Now, I expect that pretty much everyone has noticed that I give away part of the answer in the above snippet, where I refer to the new VBscript function complexPPower:

t2z05   = complexPPower(t, complexSum(z, 0.5))

From the mathfacility site,
f(z) = z^c = exp(c * log(z))

‘ p1^p2 = exp(p2 * log(p1))
Function complexPPower(p1, p2)
cPPT1   = complexLog(p1)
cPPT2   = complexProd(p2, cPPT1)
complexPPower = complexSPower(e, cPPT2)
end Function

And, yeah, that’s all there is to it. So, I have all this code together, and the time comes to run it. I have some syntax errors, but I fix those eventually. The code runs, and I get an answer out. I compare it to wolframalpha.com, for gamma(0.6 + i*0.4), and the numbers don’t match.

I break the code apart into individual components. I print the components out with wscript.echo. I tweak a couple lines where I made more mistakes. The numbers still don’t come out right. I put the code away for 2 weeks and play Legend of Zelda. I finish Zelda, and break the code out again. I start verifying the values generated from the for-loop on each of the values from the p[] array, and they look fine. I sum up the values from the for-loop, and they look ok. I get to the part where I apply the complex log function on “t”, and suddenly I notice something.

The code I lifted for calculating complex logs has a minus sign in front of the part that calculates the imaginary component of the result. And my values for ln(z + 0.5) have result(im) coming out negative, as they should, but the value from wolfram alpha is positive.

I change my code to remove the minus sign, and lo and behold, my results for Gamma(0.6, 0.4) are an EXACT MATCH for wolfram alpha’s!!!

I run a few more test numbers, EXACT MATCHES.

I clean up my code, add more comments, and add the code for reflection at the end when z(re) < 0.5. I run more tests with gamma(-0.413, 0.62), and the like. EXACT MATCHES.

So, yeah, I’m getting giddy now. But, why not? My code works. It may not be readable to anyone else. It may not be useful to anyone else either, including me. But, this was supposed to be a learning exercise, and I learned a lot. I can’t ask for anything more than that. And the CODE WORKS. Heh heh heh. I think I may try doing an Excel animation anyway, if I can find part of the graph that looks interesting enough.

Here’s the VBscript code, minus the comments and extra complex functions.

Dim p, epsilon, pi, e, isReflect, re, im, sqrt2PI

re        = 0
im        = 1
epsilon   = 0.0000001
pi        = 3.1415926535897932384
sqrt2PI   = (2*pi)^0.5
e         = exp(1)

p  = Array( 676.5203681218851,        -1259.1392167224028,          771.32342877765313, _
-176.61502916214059,          12.507343278686905,         -0.13857109526572012, _
9.9843695780195716e-6,      1.5056327351493116e-7)

gammaVal = gamma(complex(-0.6, 0.4))

wscript.echo “gamma (re) = ” & gammaVal(re)
wscript.echo “gamma (im) = ” & gammaVal(im)

Function gamma(zG)
isReflect = False
z         = zG

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

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

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

‘  result = sqrt(2*pi) * t**(z+0.5) * exp(-t) * x
t       = complexSum(z, lenAry(p) – 0.5)
minusT  = complexProd(-1, t)
t2z05   = complexPPower(t, complexSum(z, 0.5))
result1 = complexProd(t2z05, sqrt2PI)
result2 = complexProd(result1, x)
result3 = complexProd(result2, complexSPower(e, minusT))

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

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

gamma = result

End Function

One last comment. The reason for storing the first value of gamma to result3, and then doing the silly bit in the if-statement to change it to result, and then later storing result to gamma is because of a limitation in VBscript. The way to return a value from a function to the caller is shown in the last line of the function.

gamma = result

This is because the name of the function is “gamma()”. Now, after storing the result to the variable gamma, any further attempts to manipulate that variable will result in syntax errors. This means that I can’t do the withinepsilon() test on the variable “gamma”. Instead, I have to do the silly bit of storing the value to “result” first if there’s no reflection, or calculate the reflection and store that to “result”, then do the withinepsilon() test of “result”. After zeroing result(im), I can finally do:

gamma = result

And not get those pesky syntax errors.

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: