Subways of the old Romans

This is very cool. Someone unearthed one of the subway maps used by the Romans. Now we know why they were able to run their Empire so efficiently. And a one-month rail pass cost only IV denarii.

Hachette 3D Puzzle Series, vol. 9

(All rights belong to their owners. Images used here for review purposes only.)

Well, it’s been a couple months since I first wrote about Hachette’s new serialized 3D wood puzzle magazine series. Of the 9 puzzles so far, I only had an interest in #1 (because it was the cheapest one, and I wanted to see what the magazine looked like), and #6. Unfortunately, for some reason #6 is the one that sold out right away and I haven’t seen a copy in the stores at all. These things do run $15 each, so it’s not like I really feel a desperate need to buy it over the net, but if I held the box in my hands in the store, I’d probably get it. (#6 is a 3x3x3 wooden cube slightly rotated and designed to set within a restraining dowel peg holder.)

I also liked the look of #9, so I bought the only copy still at the Junkudo bookstore here

(Ad for the display shelves.)

At first glance, “Mysterious Ball” is a puzzle in that there’s nothing saying what you’re supposed to do with it. I pushed and pulled the pieces sticking out of the ball with no effect. Then I checked out the magazine to see what was in it. One of the first pages gave a hint, and another showed the thing fully disassembled. I didn’t even bother looking all that closely at the photos. 30 seconds later, I had it apart, and back together again. I spent another 5 minutes really examining the MB and making sure I knew exactly why the parts moved sometimes and not others. And now, yeah, I can confidently claim that I’ve got it mastered. It’s got a difficulty level of 3 out of 5 stars, which is probably overstating things a bit. I’d say that it’s an “Easy.”

(Philidor article.)

The magazine feels kind of short this time. While the full thing is 20 pages, the outer sheet is detachable and is primarily advertising for the full series, with a page dedicated to the clear shelving structure you can buy for 5,000 yen ($48 USD) to display all the puzzles in your collection. That leaves 16 pages for the actual magazine itself, where you get the front cover, the inside front cover with publication information and table of contents, 2 pages describing the Mysterious Ball, 5 pages of text puzzles plus their solutions, 2 pages on Francois-Andre Danican Philidor, 1 page for the classic painting memorization game, and 4 pages detailing the solution for the Gyro puzzle from volume 8 (which I used to have a long time ago, and felt no need to buy again). If you’re keeping track, that’s just 7-8 pages of useful material (4 pages of text puzzles, the memorization painting and the piece on Philidor). (The solution for the Mysterious Ball will appear in vol. 10.)

(Look at this painting for 1 minute. Then look away, and answer the questions about it.)

Philidor (1726-1795) was a French composer, best known musically for his comic opera, Tom Jones (1765). However, he may have found a longer-lasting fame as a chess player. He’s considered the best player of his age, and had been friends with Benjamin Franklin. He wrote one of the most definitive books on chess, “Analyse du jeu des Echecs” (1749), with other editions in 1777 and 1790, and it went through another 70 editions and had been translated into English, Russian, German and Italian by 1871. He was the first player to recognize the importance of the pawn in the game. The Philidor Defense, Philidor’s Legacy and Philidor’s position are all named after him. In 1882, Amedee Dutacq (music) and Abraham Dreyfus (libretto) premiered a one act opera-comique entitled “Battez Philidor,” in which an impoverished musician has to beat Philidor at chess in order to win the hand of his sweetheart. Although Philidor agrees to throw the match, he gets distracted and wins by accident.

(Example page for the text puzzles.)

Summary: I like the Mysterious Ball puzzle and plan to challenge my students to solve it. I found the article on Philidor to be interesting, and I plan on tackling the remaining text puzzles at some point.

The next puzzle is going to be the cube snake, scheduled for a June 21 publication. This is the same elastic band wrap-around puzzle that I got from the capsule ball dispensers for $2, so I’m not going to buy volume 10. #11 is the Tower of Hanoi. #12 is the “Devil’s Star” (also from the capsule ball series). #13 is “Turn Table,” another scramble/descramble game in the vein of Rubik’s Cube. It’s pretty, so I might get it, although I’m useless with the Rubik’s Cube, so maybe I won’t. #14 is the White Wood Pig, one of the rare 3D burr puzzles that I kind of want to get. #15 is the Devil’s Box, one of only 2 puzzles rated 5 stars. It looks nice, so I might get that. #16 is a set of wooden disks with color markers that you’re supposed to align so the colored dots match up on all disks. It’s kind of a Penrose tile game, so I might get that, too. #17 is the Cross Piece Cube, which is similar to one of the other capsule ball puzzles. And #18 (the series might go past 18 volumes, but there’s no confirmation of that) is the Devil’s Molecule, which also looks pretty. At the moment, I’m mildly interested in #13 and #16, slightly more interested in #14 and #17, and will probably definitely get #15 and #18. But, these are coming out one every 2 weeks, and, again, they’re about $15 apiece. Even if I do get the ones listed here, the first one won’t hit the shelves until August. See you then.

Ch. 2 Solutions

Number of regions that can be made for a pancake of n straight cuts:
1/2 * n * (n-1)

Number of designs for n beads of 2 colors:
There’s no way to solve this with finite differences. It has to be approached with a recursive program.

Number of triangles for n straight lines:
You need 3 lines to get one triangle, so
0 1 2 3 4  5 : Lines
0 0 0 1 4 10
0 0 1 3 6
0 1 2 3
1 1 1

Assume that no two lines are parallel, and no more than two intersect at the same point. No three lines can form more than 1 triangle, and any set of three lines must form one triangle. The problem becomes similar to: How many different ways can n lines be taken three at a time. The answer can be found using Newton’s formula, as 1/6*n*(n-1)*(n-2).

Comments: Finite differences can be fun. Pick something you can measure at fixed intervals (car speeds or distances, average numbers of cars on the freeway during the day, prices of haircuts at different hairstylists) and see if you can find patterns for them as a factor of n.

Colossal Gardner, ch. 2

Chapter 2 is a continuation of the section on Arithmetic and Algebra.

This time, we get the Calculus of Finite Differences. This concept was first published between 1715 and 1717 by Brook Taylor, who developed the Taylor Theorem. Finite differences are an early form of integral calculus, and the idea is to look at the results of an equation in steps, effectively taking the differences of that equation at those steps to determine the underlying properties of an event. It’s not an infallible way of building up equations to explain specific behavior, but it’s a good starting point in many cases.

Gardner describes a “mind reading” party trick used by W. W. Sawyer, who taught math at Wesleyan University. You ask someone to think of a quadratic equation, with no powers greater than x^2 (to keep things simple).

For example, use 5*x^2 + 3*x – 7.

You ask the person to give you the answers for x = 0, 1 and 2. You can turn your back on them so you can’t see them working the math. In this case, the answers would be -7, 1 and 19. With practice, you could do the following steps in your head in a couple seconds.

The process is to write out the answers in a row:
-7 1 19

Subtract the number on the left from its neighbor to the right and put the results on the next row:

(1 – -7) = 8, and (19 – 1) = 18
8 18

And do the same thing to build up the third row.

18 – 8 = 10

-7 1 19
8 18

The original equation can be built up by the following rules:
The coefficient for x^2 is 1/2 of the bottom number.
The coefficient for x is the first number of the middle row minus half the bottom number.
The constant is the first number of the top row.

a*x^2 + b*x + c
a = 10/2 = 5
b = 8 – 10/2 = 3
c = -7

The calculus of finite differences is a study of these kinds of equations. It had been further developed by Leonhard Euler and George Boole (creator of boolean logic) but fell out of favor in the 1800’s. In the later 1900’s, it became useful again for use in statistics and the social sciences. You can determine the gravitational constant this way. Time the fall of a stone and record the distances at one second intervals:

0 16 64 144 256
16 48 80 112
32 32 32

distance = 32/2 * x^2 + (16 – 32/2) + 0 = 16*s^2

There’s no guarantee that this equation holds absolutely in all cases, but it’s one way to start approaching the development of a theory.

I’d done something like this myself many years ago when I was bored one day, looking at the differences of different strings of powers:

1 4 9 16 25 36
3 5 7  9  11
2 2  2  2

1 8 27 64 125 216
7 19 37 61 91
12 18 24 30
6  6  6

So, I’d discovered finite differences without having any idea what it was. Note that as you get higher powers of x, you get more lines on the pyramid.

Gardner gives a couple additional puzzles to play with – What’s the maximum number of pieces a pancake can be cut into by n straight cuts? How many different designs can you get when making circular necklaces of n beads of two colors (the beads are white or black). And, what’s the maximum number of triangles that can by made with n straight lines?


Power Pi

I was asked recently how people calculate PI. That is, what method is employed on a supercomputer to get the value of PI to a specific number of decimal places. I didn’t know off-hand, so I looked it up. Wolfram gives 32 power expansions for PI. For those of you unfamiliar with the math, a power expansion (or power series) is just another way to write out a number as the summation of individual arithmetic operations (i.e. – terms).

The simplest one to implement in VBScript is:
pi = 4 * sum [(-1)^k / (2*k + 1)] where k goes from 0 to infinity in steps of 1.

The (-1)^k part causes each iteration of the formula to oscillate between two points that gradually close in on the correct value. Then the 1/(2*k+1) is the portion of the power expansion that actually contributes to the final value. It’s easy to write, but it is slow. On my laptop, I need to run k up to 5,000,000 terms, and it takes 15 seconds to finish. Even then, the result is only accurate to 6 decimal places.

I tried to streamline the program a bit, so the part where I multiply the result by 4 isn’t until the final wscript.echo statement. Note also that there’s little point in calculating (-1)^k EVERY SINGLE TIME, so I just use sn2 to hold the current value of the sign flipping component.

for k=0 to 500
tmp = sn2/(2*k + 1)

is the same as

for k=1 to 1001 Step 2
tmp = sn2/k

but is maybe 25% faster to run.

res = 0
sn = -1
sn2 = 1
max = 2 * 5000000
for k = 1 to max Step 2
tmp = sn2 / k
res = res + tmp
sn2 = sn2 * sn

wscript.echo “res = ” & cstr(4 * res) & ” term ” & cstr(k) & ” = +/-” & cstr(2 * tmp)

The actual value of PI is: 3.14159265359
My value is 3.14159285359 At step 5,000,000
and 3.14159245359 At step 5,000,001
So, I’m at +/- 0.000000199999.

Not a really useful program, but it is fun to see the answer coming out close.


There’s a more direct method with just two terms (using Excel expressions):

pi = 4 * (4 * atan(1/5) – atan(1/239))

This of course now involves finding the irrational values of two separate numbers to whatever accuracy you’re seeking for pi… Maybe not an improvement.


Many of the formulas described at the Mathworld page are circular, in that to find pi, you have to use pi within the power series. Which to me seems to be kind of hard to do in software. However, I did like the version which comes from the special value of the Riemann zeta function for (pi^2)/6 = zeta(2n) for all positive integers n:

pi = sqrt(6*(1 + 1/2^2 + 1/3^2 + 1/4^2 + …))

In VBScript, this becomes:

t = 0
for n = 1 to 5000000
t = t + 1/n^2
wscript.echo “pi = ” & cstr((6 * t)^0.5)

pi = 3.14159246260382,

This is 29% faster than the first program above (5 seconds for 500 million terms, compared to 6.5 seconds for the first program), and they are both still accurate to only 6 places. (If both programs run for 500 million terms, the first program takes 357.6 seconds, and the second for 255.7 seconds, giving program #2 a 28.6% speed advantage, although #2 was only accurate to 7 places and #1 was already accurate to 8 places.)

After 500 million terms:
PI = 3.14159265359
#1 = 3.14159265158926
#2 = 3.14159264498239

Back to Riemann, Part 19

(plot[gamma(x + i*y) * gamma(x + i*y) x=-10 to 10 y=-10 to 10])

Right after writing up Part 18, I returned to reading the Martin Gardner math recreations book I’d received for Christmas (it’s a big book). There were a couple articles on fractals and fractal sound, and I got to thinking that I’d like to do something similar. The key element to the Mandelbrot fractals is that they’re based on reiterations of equations on the complex plane, and that’s exactly what the gamma function does, so, who knows, maybe performing gamma on itself in some permutation might be interesting.

After spending several days waiting for Excel to slowly grind through each of the animation approaches I decided to try, I gave up. The main problem is that the numbers start out really small and then get extremely big very fast, causing Excel to clip out most of the plots in half of the frames (even when treating the values as logs). So, from an animation viewpoint, gamma has its limitations. But, the Wolfram Alpha plots look interesting, so I’m posting them here.

(Looks like a Nazca drawing.)

(plot[gamma(x + i*y) * gamma(10-x + i*y)) x=-10 to 10 y=-10 to 10])

(plot[gamma(x + i*y) * gamma(x + i*(10-y)) x=-10 to 10 y=-10 to 10])

(plot[gamma(x + i*y) * gamma(10-x + i*(10-y)) x=-10 to 10 y=-10 to 10])

(plot[gamma(gamma(x + i*y))) x=-10 to 10 y=-10 to 10])

With gamma of gamma, Wolfram actually times out. Either it draws only the real component, the real and imaginary components, or just the real contour map. The graph takes too long for it to do all four plots.

(plot[gamma(1/gamma(x + i*y)) x=-10 to 10 y=-10 to 10])

And with this, I’m going to call it quits for gamma for a while.

Back to Riemann, Part 18

So, what does the gamma function look like in Wolfram Alpha?

From an Excel animation viewpoint, the really interesting part of the graph is when you get to about x=6. From the above plots, you can see that there are a couple small poles in the x< region, then the big one at x=0, which each provide a little activity leading up to x=6, but otherwise, the map is flat for the left hand half of the complex plane.

The animation also depends on which variable you treat as the line segment. In the first pass, the segment for x=-10 to 10 was calculated for varying values of y (for y=-10 to 10). In the second pass, it was the other way around. The second animation is much more interesting. I also zoomed in on the segment for the area from -0.5 to 1. Note that for the animation where the y segment is moved around, the segment disappears due to clipping as the results from gamma() get larger. The unraveling linked spirals are kind of hypnotic, but there aren’t quite enough data points between x=6 and x=8 to generate enough frames for looping them in movie maker.

Direct youtube link

And this is pretty much where I’m going to call it quits with gamma and zeta.

A final few words on the book itself. The Riemann Hypothesis, by der Veen and de Craats is a short book, at only 144 pages, but it is highly readable, and much more approachable than any of the sources I looked at. They include a lot of exercises that can help you link the different concepts together, as well as showing how to use Wolfram Alpha for doing the computer grunt work. I don’t know why it’s suddenly so expensive on Amazon, but if you can find a copy used for $15, it’s worth the price. Recommended if you want a decent introduction to Riemann and the zeta function. I’ve gotten about 3 month’s worth of activities from it, which isn’t bad for a little book like this.

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:


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, 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
‘ @
‘ by A-Jim Deutch, Syracuse, New York
‘ I got this from – 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
If ys < 0 Then
theta =  (-1) * pi / 2  ‘ 270
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, 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)

‘  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
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.

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.


“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
reflSin   = complexSin(complexProd(pi, zG))
reflect1  = complexProd(reflSin, result3)
reflRecip = complexRecip(reflect1)
reflect2  = complexProd(pi, reflRecip)
gamma     = reflect2
end if


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)

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))))

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 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.

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
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. 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, 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.