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:

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.

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.

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.

Back to Riemann, Part 14


Finally, the culmination of all my messing around – some animations.

Direct link, b from 1 to 100

Direct link, b from 1 to 500

Direct link, b from 5,000 to 5,100

Back to Riemann, Part 13


Ok, so I’m now ready to sit down and write up an Excel spreadsheet to graph the first few zeta zeroes. As I’m preparing for this, I noticed that I made a mistake in part 8 of the Riemann posts, a couple hours after I uploaded it to wordpress, so I had to go back and fix that. I’m also using part 8 for getting the math I need for the spreadsheet.

First, the zeta function is:
zeta(x) = 1 + 1/2^x + 1/3^x + 1/4^x + 1/5^x…

1/2^x translates to: (1/(2^a)) * (cos(b * ln(2)) – i*sin(b * ln(2)))

The generalized form is:
1/n^x = (1/(n^a)) * (cos(b * ln(n)) – i*sin(b * ln(n)))

Second, the eta function is:
eta(x) = 1 – 1/2^x + 1/3^x – 1/4^x + 1/5^x – 1/6^x…

Writing the generalized translated form of the eta function just means having alternating minus signs for the even terms.

Where things got trickier than I was expecting was with the conversion factor to go from the eta(x) value to the zeta(x) value.

zeta(x) = eta(x) / (1 – 2^(1-x))

Remembering that
2^x = 2^(a + bi) = (2^a) * (2^bi) = (2^a) * e^(bi * ln(2))
e^bi = cos(b) + i*sin(b)
(2^a) * e^(bi * ln(2)) = (2^a) * (cos(b * ln(2)) + i*sin(b * ln(2)))

This gives:
1/(1 – 2^(1-x)) = 1/(1 – ((2^(1-a)) * (cos(b * ln(2)) + i*sin(b * ln(2)))))

Or, since a = 1/2,
1/(1 – sqrt(2) * cos(b * ln(2)) + i * sqrt(2) *sin(b * ln(2)))

Simplified:
1/(j + ik)

And, yeah, I didn’t remember how to take the reciprocal of a complex number, and this doesn’t seem to have been covered in the der Veen and de Craats book. Fortunately, this is one of the rarer cases where I could find the answer quickly on the net.

Step 1, find the conjugate of the complex number (i.e., j – ik).
Step 2, multiply the numerator and denominator by the conjugate.
Step 3, simplify.

If j = 1 – sqrt(2) * cos(b * ln(2))
and k = sqrt(2) * sin(b * ln(2))

(j – ik)/(j + ik)(j – ik)
(j – ik)/(j^2 + ijk – ijk + k^2)
(j – ik)/(j^2 + k^2)

Writing all this out here in the blog would be too messy. But, doing this in Excel isn’t that bad because I can treat “j” and “k” as cell names.

Getting the zeta function results is even messier, but in simplified form:
If eta(x) = m + in

zeta(x) = (j*(m+n) – ik*(m+n)) / (j^2 + k^2)
Believe me, it makes more sense when you put everything into Excel and refer to the cell names this way.

This then gives me two choices as to how to proceed next. First is to just do a simple spreadsheet for one zero, and write a little VBScript to update the spreadsheet to make an animated dot running around the graph. Second is to copy the first set of columns multiple times to draw a more detailed line segment showing the zeta going to zero.

The second approach has the advantage of overcoming some of the rounding errors introduced by Excel when adding up the zeta function terms. The results of the Excel sheet aren’t going to be all that accurate, so it’s not going to go to |zeta(x)| = 0, it’s just going to get close and then move away again. But that’s ok.

So, which approach did I pick?
#2

Without the VBscript animation part, this is what I’ve got for the first zero.

This graph shows the real and imaginary components of zeta(x) for a=0.5 and b approaching the first zero at 14.134

Here we’ve got the magnitude of zeta(x), and it’s pretty obvious it’s not going to zero. I blame rounding errors.

And then, the zeta(x) function from a + ib = 0.5 + i*14.0 to 0.5 + i*14.5 in 0.05 steps on the complex plane. This is the line segment I’d animate with a VBScript driving the Excel graphs. If my math were perfect, the segment would be crossing the axes at 0,0. As it is, I have the eta(x) calculated out to 300 terms. Not sure if it’s worth the effort of going to 1,000 terms or not. Depends on how bored I get, I guess.

All that’s left now is to write the VBScript (I can modify the one I wrote for the first Riemann prose blog series), and to decide if I want to figure out a VBScript version of the Gamma approximation.

Back to Riemann, Part 12


One of the issues I hadn’t covered in much detail in the last post was what happens if Riemann’s hypothesis fails somewhere. That is, what if one of the zeta zeros is not on the critical line x = 1/2.

Remember that in the explicit formula for the prime counting function, we’re summing up the effects of ALL of the zeros for a specific value of x (x being the range of counting numbers that contains the prime numbers we want). To make that easier to do, we pair up each zero with its mirror version: u = v + wi with mirror(u) = v – wi.

The equation for that zero-pair then becomes:

-(2/|u|) * (x^v) * cos(w*ln(x) – a)

If the zero lies on the critical line, we get:
-(2/sqrt((1/2)^2 + w^2) * sqrt(x) * cos(w * ln(x) – atan(w))

If it doesn’t, then we end up with two pairs of zeros:
v + wi, v – wi, 1 – v + wi and 1 – v – wi

And we get two equations:

-(2/sqrt((v)^2 + w^2) * x^v * cos(w * ln(x) – atan(w/v))
-(2/sqrt((1-v)^2 + w^2) * x^(1-v) * cos(w * ln(x) – atan(w/(1-v)))

Right away, it should be obvious that if v is really close to 0.5, that we’re going to have two results that are very similar, in effect doubling the impact of this zero on the total for the psi(x) summation. And, although size of w for the zero is going to be much larger than v is (say we use the 410th zero: 693.1769701, and v is 0.51), meaning that the variation in v won’t affect the absolute value of v + wi very much, the variations will become much more noticeable in the x^v and x^(1-v) parts as we calculate psi(x) for larger values of x (i.e. -> x = 100,000^(0.1) or 100,000^(0.9)).

Then, what would that look like? Well, first, it depends on how far from v = 1/2 the zeros are. We know that there are no zeros for v < 0 or v > 0, and it turns out that there are no zeros directly on the boundaries v = 0 or v = 1. In fact, the zeros are going to be closer to v = 1/2, anyway. But let’s assume that the rogue zero is somewhere between 0.1 and 0.9. If we have one zero at 0.9 + w, then we’ll get matching zeros at 0.9 – w, and 0.1 + w and 0.1 + w.

Now, let’s arbitrarily pick one of the zeros for w, say the 13th one, at 59.347044003.

Then, the plots for 2/sqrt((v)^2 + w^2) and 2/sqrt((1-v)^2 + w^2), with v running from 0.1 to 0.9 are above.
This is pretty flat.
If we add them together, remembering that we’re not counting v = 0.5 twice…

They are flat, except that, as mentioned before, having 4 zeros instead of 2 is doubling the magnitude of these pairs.

Now, the place where we really get a lot more contribution from these zeros is in the parts x^v and x^(1-v). If v is 0.9 and 1 – v is 0.1, then if x is 1000 and v is 0.9, (2/(|u|)) * x^v = 8.4, compared with 0.533 when v=1/2.

The reason for using x=1000 is simply that the contribution for (2/(|u|)) * x^0.5 is so small if x is small. Besides, looking at psi(1000) starts getting us up into somewhat larger values of x, and the entire point is what happens when x tends towards infinity.

Adding the results from the above graph together, we get:

where the v-notch represents the smallest value for the contribution for v = 0.5 and w = 59.347044003. If the zeros here are off the critical line, they really throw off the results of the psi(x) approximation.


(Correct approximation for psi(100.)

Compare the correct graph above, for psi(100), with the 11th zero set so that v=0.9 and 1-v=0.1, below.


(Summations of the first 11 zeros, comparing the correct values to v=0.1/v=0.9 for the 11th zero.)

der Veen and de Craats talk about the “music of the primes.” We saw this a little bit in the last blog entry. The cosine component for each non-trivial zero of the explicit formula creates a growing, slowing cosine result as x gets bigger. These ringing bell waveshapes get added together for ALL zeros, spanning the entire range of x, in the psi(x) equation. This is why the Riemann hypothesis is so important. If even one zero is off the critical line, you get the 2 pairs of zeros adding together instead of the one pair, throwing off the step-shape of the approximation. And that invalidates the prime number theorem.

The way der Veen and de Craats put it, each non-trivial zero, in the explicit formula, acts like the next octave up on the “musical scale”. These are smooth, tonal increments from one zero to the next. If any of the non-trivial zeros is not on the critical line, the amplitude of that doubled-up zero pair becomes an ugly discordant spike.

I’m showing the first 10 zero pairs from the explicit formula to show what I mean. I introduced the first and tenth zeros in the last blog entry. Adding them all together gives the approximation for psi(100) in the above graph.

Also, I just want to show the magnitude contributions for the first ten zeroes, as x goes to 6,250. The relative spacing between one zero and the next remains constant, independent of x.

But it’s not a linear, smooth or consistent spacing.


(Left: Relative spacing of the contributions of the first 10 zeros for the explicit formula. Right: For comparison’s sake, the relative spacing for the first 10 prime numbers. Note that spacing of the zeros doesn’t reflect the spacing of the primes themselves.)

Anyway, the objective for this blog entry was just to consider the impact of having one of the non-trivial zeros not being on the critical line. The answer is that it really depends on how far off from the line the two matching zero pairs (v and 1 – v) are. If they’re very close to 0.5, then they just double the magnitude of the -(2/|u|) * (x^v) * cos(w*ln(x) – a) contribution for that zero.

If they’re farther apart, such as being at v = 0.9 and 1 – v = 0.1, then the results can be up to 8 times greater. It all depends on how big x is where you’re looking.

Having a proof showing that the Riemann hypothesis is true, or even that it is false, would remove the uncertainty from the question. After 158 years, we’re still waiting for someone to figure out what that proof is.

Back to Riemann, Part 11


The entire point to this discussion is, “How do we estimate the distribution of the prime numbers?”

At one stage, I’d gotten into the idea of the prime counting functions, psi(x) and pi(x). pi(x) physically counts the primes less than x, and is approximated by pi(x) = x/ln(x). And psi(x) is a weighted summation of the powers of the primes. psi(x) = T2*ln(2) + T3*ln(3) + T5*ln(5) + T7*ln(7)… I.e. – psi(10) = 3*ln(2) + 2*ln(3) + ln(5) + ln(7) = 7.830 (there are three powers of 2 – 2, 4, 8; two powers of 3 – 3, 9; and one each of 5 and 7). As x gets larger, psi(x) approaches x, and pi(x) approaches x/ln(x). Therefore, psi(x) approaches pi(x)*ln(x).

der Veen and de Craats give a more accurate equation, called the “Explicit Formula for psi(x),” without going into the details for how Riemann derived it, or how von-Mangoldt simplified it. The details can be found in the wiki article.

Although psi(x) gets close to x as x goes to infinity, there is some error that can be accounted for through a correction factor. In this case, ln(2*pi) just shifts the diagonal line y = x down by 1.83787… If x is greater than 1,000, then this term can effectively be ignored.

The third term is a summation of all the zeros of the zeta function, both the trivial and non-trivial ones. If we separate them out, as follows,

as the summation of all the trivial zeros (x raised to the trivial zero divided by that zero; i.e. x^(-2)/(-2) + x^(-4)/(-4) + x^(-6)/(-6)…), and the summation of all the non-trivial zeros.

If we go back to the series expansion for ln(1 + y),
ln(1 + y) = y – (y^2)/2 + (y^3)/3 – (y^4)/4 + (y^5)/5…

plugging in y = -1/x^2 and multiplying both sides by 1/2 gives:
1/2 * ln(1  – 1/x^2) = 1/2 * ( -1/(x^2) – 1/(2*x^4) – 1/(3*x^6) – 1/(4x^8)…
= (x^-2)/(-2) + (x^-4)/(-4) + (x^-6)/(-6) + (x^8)/(-8)…

recognizing that this last form is the summation of the trivial zeros,
the trivial summation can be simplified as:

1/2 * ln(1 – 1/x^2).

If x = 50, then x^2 = 2,500, and 1 – 1/2500 = 0.9996.
-1/2*ln(0.9996) = 0.0002

Which is very small, and goes to 0 fast as x gets bigger. Therefore, the ln(2 * pi) and the 1/2 * ln(1 – 1/x^2) terms can be ignored as x goes to infinity (they total less than 2). They are needed to get the exact value of psi(x), but it’s pretty obvious that the main contribution to psi(x) is from the non-trivial zeros.

-1.837877… + 0.0002 = -1.837677…
If x = 50, and psi(x) approximates x, then the contribution from -ln(2*pi) – 1/2*ln(1 – 1/x^2) is only 3.67%


(The sum of the trivial zeros goes to zero pretty fast.)

Since what we’re really interested in is the distribution of prime numbers over 100 digits long, we can ignore the correction factors and just focus on the term with the non-trivial zeros.

The next step is referred to as pairing up the non-trivial zeros. Because the explicit formula is a sum of ALL the non-trivial zeros, we need to look at the ones in the negative half of the complex plane (b*i < 0) as well as those in the positive half (b*i > 0). We can do this by pairing them up and looking at the angle they form. (Note, I’m using mirror() here to represent the mirror value of a complex number.)

If u = v + iw is a zero, then so is mirror(u) = v – iw

The summation of the non-trivial zeros is in the form of -(x^p)/p, so adding our zero “u” and mirror(u) is going to look like:

eq. 1) -(x^u)/u – (x^mirror(u))/mirror(u) = -(2/|u|) * (x^v) * cos(w*ln(x) – a)

And we do this for all of the non-trivial zeros.
(|u| is the absolute value of u.)

Now, the thing that’s interesting about this pairing is it’s independent of the opposite mirror: 1 – v. That is, while we’re pairing v + bw up with v – bw, we need to make a separate pairing of 1 – v + wi with 1 – v – wi.

Why is this important? Because it gives us 4 zeros, not 2.
But, Riemann’s hypothesis says that all of the zeros are on the critical line, so we should only have 2 matching zeros, not 4, and v would = 1/2 in all cases. This simplifies the math, as well as minimizing the impact of the non-trivial zeros on the explicit formula for calculating psi(x). If the hypothesis fails, this double pair is going to have twice as much weight as the zeros that are on the critical line.

Ok, so what do these paired non-trivial zeros look like?
The easiest way to find out is to break eq. 1 down into individual components.
u = v + iw
|u| = sqrt(v^2 + w^2)
a = atan(2*w/(2*v)) = atan(2*w) (for v = 1/2)

Then, what are v and w?
They’re the zeros for the zeta function. I’ve already stated that all the zeros that have been found so far lie on the critical line, so v = 1/2. I’ve calculated the stats for the first 10 zeros in Excel. As the zeros get bigger, “a” gets closer to pi/2 = 1.570796

Then, what I did was to set up a spreadsheet where x runs from 1 to 100 in 0.2 increments in the F column. vk is in the A column, wk in the B column, |uk| is in the C column, and ak (atan(2*wk)) is in the D column.
G column is then the paired non-trivial zero equation:
=(2/$C$2)*POWER($F2,$A$2)*COS($B$2*LN($F2) – $D$2)

And, for the first non-trivial zero (0.5 + i*14.13457), the graph for x = 1 to 100 is:

The 10th zero is at 0.5 + 49.77383, and the graph is:

Summing the results of the first 10 zeros and subtracting from x, for x = 1 to 100, gives the following graph, which is plotted along with psi(x). Remember that psi(x) counts the powers of each prime (i.e. -> 2, 4, 8, 16, 32, 64) from 1 to x, and sums the weighted primes as:
psi(x) = T2*ln(2) + T3*ln(3) + T5*ln(5) + T7*ln(7)
where, if x=10, T2=3, T3=2, T5=1 and T7=1

This is just for the first 10 zeros. Adding all the other known zeros to infinity creates a step function that is a near-perfect match for psi(x).

Where does this leave us? Well, we know that the Riemann zeta function can be extended to the entire complex plane, with the exception of x = 1. We know that the zeros for the zeta function lie in the critical strip 0 < x < 1, that the zeros can be paired up in the explicit formula for psi(x), and that the known non-trivial zeros lie on the critical line x = 1/2, which is covered by the eta(x) function. If we plug the known zeros into the explicit formula and subtract this from x we get psi(x), which is a prime counting function that physically does count the prime powers and is then a weighted sum of those powers from 2 to x.

Actually, there are two things left I want to deal with. Creating an Excel file for confirming the first few known zeros, and addressing why having a zero that is not on the critical line messes everything up.

Back to Riemann, Part 10


Ok, back to the der Veen and de Craats version of the Riemann functional equation.

z( -x ) = ((-2 * x!) / (2*pi)^(x+1)) * sin(pi*x/2) * z(x + 1)

What this states is that if you have the zeta(x) value for any complex x (x = a + bi) then you can get zeta(x + 1), and from there apply the functional equation to define zeta(-x).

This extends the zeta function to the left hand side of the complex plane (a < 0). Combining the functional equation with eta(x) and the original version of zeta(x) for x>1, we now have complete coverage of the complex plane everywhere except for the pole point x = a = 1.

One interesting feature here is that the mirror relationship between -z and z + 1 gives a center of rotation at x = 0.5.
c = B + (A – B)/2 = -z + (z + 1 – (-z))/2
c = -z + (2z + 1)/2
c = -z + 2z/2 + 1/2
c = 0.5, for any z.

I’ve mentioned before that the zeta(x) function converges to 1 for large x, meaning that there are no zeros for x>1.

And, the functional equation relates -z to z + 1, so again, there are no zeroes for x < 0, except for the trivial zeroes at x = -2, -4, -6, etc, because of the sin() component. (If there were any other zeroes to the left of x = 0, then there would have to be matching zeroes for x > 1, and there aren’t.)

That just leaves what is called the “critical strip,” for x = a + bi, 0 <= a < 1, which is covered by the eta(x) function.

Given the nature of the complex plane, if we know the value for x = a + bi, then we also have the values for a – bi, -a + bi and -a – bi. Which means that if zeta(a + bi) = 0, then zeta(a – bi) = zeta(-a + bi) = zeta(-a – bi) = 0.

That is, if there is a zero within the critical strip, it will have 3 more matching zeroes to go with it.

What Riemann did with his hypothesis was minimize the number of possible zeroes, by saying that all of the non-trivial zeroes lie on the critical line x = 0.5. Or, instead of having 4 zeroes for any given a + bi combination, there are only 2 (a + bi and a – bi).

So far, no one has been able to prove OR disprove this statement. Billions of zeroes have been found by using computers, and they’re all on the critical line. But, that’s not a proof. If even just one zero exists that is not on the critical line, then the entire hypothesis collapses, as do all the other theorems that are based on it.

Ok, so that’s the basic concept of the zeta function. We have the regular version,
z = 1 + 1/2^x + 1/3^x + 1/4^x…

For x > 1 across the complex plane.

We have the eta(x) function,
eta(x) = 1 – 1/2^x + 1/3^x – 1/4^x…

followed by,
zeta(x) = eta(x) / (1 – 2^(1-x))

for 0 < x < 1 along the complex plane.

And Riemann’s functional equation above, for x < 0.

If we want to find the non-trivial zeroes, we need to use the eta(x) plus the conversion formula.

In order to follow the examples in the book, I’ve been playing around with Wolfram Alpha. It’s a pretty decent online math graphing tool, but there are some inconsistencies between the version online now and what der Veen and de Craats apparently used (dVdC say that the “|” character is used to calculate absolute values, i.e. “plot |x| for x = -10 to 10”. But that doesn’t work for me. I have to use “plot abs(x) x = -10 to 10). Most of the inconsistencies are minor, but it still takes time to figure them out.

The main thing is that because Wolfram Alpha already has the zeta function implemented, it’s very easy to see right away what it looks like near the origin. If we plot zeta(x) for (-30 -30i) to (30 + 30i), we get two graphs, one of the real component of the results and the other for the imaginary component.


(plot zeta[x + y*I] x = -30 to 30 y = -30 to 30)

Note that the results for both components get clipped at z = 1000. The important part of the complex plane is that strip between 0 and 1, where the results weave around z = 0, but that’s getting swamped out.

Plotting the absolute value of the above formula, we get the combined magnitudes of the real and imaginary components, but clipping is now occurring at 10^9. No zeroes visible from this altitude.


(plot abs(zeta[x + y*I]) x = -30 to 30 y = -30 to 30)

The next step is to concentrate only on the critical strip, and to try to include at least one of the zeroes. Remember, the zeroes themselves are on the critical line, x = 0.5.


(plot zeta[x + y*I] x = 0 to 1 y = -30 to 30)

And the absolute value of the plot as well. Note the pole at x = 1.


(plot abs(zeta[x + y*I]) x = 0 to 1 y = -30 to 30)

We can try plotting just the critical line, and to find the zeroes we need to see where both the real and imaginary components cross the x-axis at the same time.


(plot zeta[0.5 + y*I] x = -30 to 30)

Or, we can look at the absolute value of the graph, where the zeroes are the points where the line touches down on the axis. Regardless, the mirror nature of the complex plane is visible here. So, we only need to find the zeroes for y > 0 and we automatically get the matching zeroes for y < 0 at the same time.


(plot abs(zeta [1/2 + x*I]) x = -30 to 30)

This just leaves the graphs of the trivial zeroes.


(plot zeta [x + y*I] x = -30 to 0 y = -0.5 to 0.5)


(plot abs(zeta [x + y*I]) x = -30 to 0 y = -0.5 to 0.5)

As x goes more negative, the results get much bigger faster, swamping the y-axis crossings. So, I’m only showing the first few trivial zeroes, for x = -10 to 0.


(plot zeta [x] x = -10 to 0)


(plot abs(zeta [x]) x = -10 to 0)

While I’m at it, I might as well include a little of the gamma function. I’ll just use the part closer to the origin to show the best detail.


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


(plot abs(gamma[x + I*y ]) x = -10 to 0 y = -10 to 10)


(plot gamma[x ] x = -10 to 0. Note the zero crossings at -2, -4, -6, etc.)

Next time: Tying the non-trivial zeroes to the Prime Counting Function

Pi*r^2


Pi*r*not^2

Pi*r*2^2

Happy Pi*r^2 Day.