Preliminaries

Suppose x is a variable representing some equal division of the octave. For example, if x = 80, x reflects 80edo with a step size of 15 cents and with pure octaves. Suppose that x can also be continuous, so that it can also represent fractional or "nonoctave" divisions as well. The Bohlen-Pierce scale, 13 equal divisions of 3/1, is approximately 8.202 equal divisions of the "octave" (although the octave itself does not appear in this tuning), and would hence be represented by a value of x = 8.202.

Now suppose that ||x|| denotes the difference between x and the integer nearest to x. For example, ||8.202|| would be .202, since it's the difference between 8.202 and the nearest integer, which is 8. ||7.95|| would be .05, which is the difference between 7.95 and the nearest integer, which is 8. Mathematically speaking, ||x|| denotes the function |x - floor(x+1/2)|.

For any value of x, we can construct a p-limit generalized patent val. We do so by rounding log2(q)*x to the nearest integer for each prime q up to p. Now consider the function

$\xi(x) = \sum_2^p (\frac{||x \log_2 q||}{\log_2 q})^2$

This function has local minima, corresponding to associated generalized patent vals. The minima occur for values of x which are the Tenney-Euclidean tunings of the octaves of the associated vals, while ξ for these minima is the square of the Tenney-Euclidean relative error of the val - equal to the TE error times the TE complexity, and sometimes known as "TE simple badness."

Now suppose we don't want a formula for any specific prime limit, but which applies to all primes. We can't take the above sum to infinity, since it doesn't converge. However, we could change the weighting factor to a power so that it does converge:

$\sum_2^\infty \frac{||x \log_2 q||^2}{q^s}$

If s is greater than one, this does converge. However, we might want to make a few adjustments. For one thing, if the error is low enough that the tuning is consistent, then the error of the square of a prime is twice that of the prime, of the cube tripled, and so forth until the error becomes inconsistent. When the weighting uses logarithms and error measures are consistent, then the logarithmic weighting cancels this effect out, so we might consider that prime powers were implicitly included in the Tenney-Euclidean measure. We can go ahead and include them by adding a factor of 1/n for each prime power p^n. A somewhat peculiar but useful way to write the result of doing this is in terms of the Von Mangoldt function, an arithmetic function on positive integers which is equal to ln p on prime powers p^n, and is zero elsewhere. This is written using a capital lambda, as Λ(n), and in terms of it we can include prime powers in our error function as

$\sum_1^\infty \frac{\Lambda(n)}{\ln n} \frac{||x \log_2 n||^2}{n^s}$

where the summation is taken formally over all positive integers, though only the primes and prime powers make a nonzero contribution.

Another consequence of the above definition which might be objected to is that it results in a function with a discontinuous derivative, whereas a smooth function be preferred. The function ||x||^2 is quadratically increasing near integer values of x, and is periodic with period 1. Another function with these same properties is 1 - cos(2πx), which is a smooth and in fact an entire function. Let us therefore now define for any s > 1

$E_s(x) = \sum_1^\infty \frac{\Lambda(n)}{\ln n} \frac{1 - \cos(2 \pi x \log_2 n)}{n^s}$

For any fixed s > 1 this gives a real analytic function defined for all x, and hence with all the smoothness properties we could desire.

We can clean up this definition to get essentially the same function:

$F_s(x) = \sum_1^\infty \frac{\Lambda(n)}{\ln n} \frac{\cos(2 \pi x \log_2 n)}{n^s}$

This new function has the property that $F_s(x) = F_s(0) - E_s(x)$, so that all we have done is flip the sign of $E_s(x)$ and offset it vertically. This now increases to a maximum value for low errors, rather than declining to a minimum. Of more interest is the fact that it is a known mathematical function, which can be expressed in terms of the real part of the logarithm of the Riemann zeta function:

$F_s(x) = \Re \ln \zeta(s + 2 \pi i x/\ln 2)$

If we take exponentials of both sides, then

$\exp(F_s(x)) = |\zeta(s + 2 \pi i x/\ln 2)|$

so that we see that the absolute value of the zeta function serves to measure the relative error of an equal division.

Into the critical strip

So long as s is greater than or equal to one, the absolute value of the zeta function can be seen as a relative error measurement. However, the rationale for that view of things departs when s is less than one, particularly in the critical strip, when s lies between zero and one. As s approaches the value s=1/2 of the critical line, the information content, so to speak, of the zeta function concerning higher primes increases and it behaves increasingly like a badness measure (or more correctly, since we have inverted it, like a goodness measure.) The quasi-symmetric functional equation of the zeta function tells us that past the critical line the information content starts to decrease again, with 1-s and s having the same information content. Hence it is the zeta function between s=1/2 and s=1, and especially the zeta function along the critical line s=1/2, which is of the most interest.

As s>1 gets larger, the Dirichlet series for the zeta function is increasingly dominated by the 2 term, getting ever closer to simply 1 + 2^(-z), which approaches 1 as s = Re(z) becomes larger. When s >> 1 and x is an integer, the real part of zeta is approximately 1 + 2^(-s), and the imaginary part is approximately zero; that is, zeta is approximately real. Starting from s = +∞ with x an integer, we can trace a line back towards the critical strip on which zeta is real. Since when s >> 1 the derivative is approximately -ln(2)/2^s, it is negative on this line of real values for zeta, meaning that the real value for zeta increases as s decreases. The zeta function approaches 1 uniformly as s increases to infinity, so as s decreases, the real-valued zeta function along this line of real values continues to increase though all real values from 1 to infinity monotonically. When it crosses the critical line where s=1/2, it produces a real value of zeta on the critical line. Points on the critical line where ζ(1/2 + i g) are real are called "Gram points", after Jørgen Pedersen Gram. We thus have associated pure-octave edos, where x is an integer, to a value near to the pure octave, at the special sorts of Gram points which corresponds to edos.

Because the value of zeta increased continuously as it made its way from +∞ to the critical line, we might expect the values of zeta at these special Gram points to be relatively large. This would be especially true if -ζ'(z) is getting a boost from other small primes as it travels toward the Gram point. A complex formula due to Bernhard Riemann which he failed to publish because it was so nasty becomes a bit simpler when used at a Gram point. It is named the Riemann-Siegel formula since Carl Ludwig Siegel went looking for it and was able to reconstruct it after rooting industriously around in Riemann's unpublished papers. From this formula, it is apparent that when x corresponds to a good edo, the value of ζ(1/2 + i g) at the corresponding Gram point should be especially large.

The Z function

The absolute value ζ(1/2 + i g) at a Gram point corresponding to an edo is near to a local maximum, but not actually at one. At the local maximum, of course, the partial derivative of ζ(1/2 + i t) with respect to t will be zero; however this does not mean its derivative there will be zero. In fact, the Riemann hypothesis is equivalent to the claim that all zeros of ζ'(s + i t) occur when s > 1/2, which is where all known zeros lie. These do not have values of t corresponding to good edos. For this and other reasons, it is helpful to have a function which is real for values on the critical line but whose absolute value is the same as that of zeta. This is provided by the Z function.

In order to define the Z function, we need first to define the Riemann-Siegel theta function, and in order to do that, we first need to define the Log Gamma function. This is not defined as the natural log of the Gamma function since that has a more complicated branch cut structure; instead, the principal branch of the Log Gamma function is defined as having a branch cut along the negative real axis, and is given by the series

$\Upsilon(z) = -\gamma z - \ln z + \sum_{k=1}^\infty \frac{z}{k} - \ln(1 + \frac{z}{k})$

where the lower-case gamma is Euler's constant. We now may define the Riemann-Siegel theta function as

$\theta(z) = (\Upsilon((1 + 2 i z)/4) - \Upsilon((1 - 2 i z)/4))/(2 i) - \ln(\pi) z/2$

Another approach is to substitute z = (1 + 2i t)/4 into the series for Log Gamma and take the imaginary part, this yields

$\theta(t) = -\frac{\gamma + \log \pi}{2}t - \arctan 2t + \sum_{n=1}^\infty \left(\frac{t}{2n} - \arctan\left(\frac{2t}{4n+1}\right)\right)$

Since the arctangent function is holomorphic in the strip with imaginary part between -1 and 1, it follows from the above formula, or arguing from the previous one, that theta is holomorphic in the strip with imaginary part between -1/2 and 1/2. It may be described for real arguments as an odd real analytic function of x, increasing when |x| > 6.29. Plots of it may be studied by use of the Wolfram online function plotter.

Using the theta and zeta functions, we define the Z function as

$Z(t) = \exp(i \theta(t)) \zeta(1/2 + it)$

Since theta is holomorphic on the strip with imaginary part between -1/2 and 1/2, so is Z. Since the exponential function has no zeros, the zeros of Z in this strip correspond one to one with the zeros of zeta in the critical strip. Since the exponential of an imaginary argument has absolute value 1, the absolute value of Z along the real axis is the same as the absolute value of zeta at the corresponding place on the critical line. And since theta was defined so as to give precisely this property, Z is a real even function of the real variable t.

Using the online plotter we can plot Z in the regions corresponding to scale divisions, using the conversion factor t = 2π/ln(2) x, for x a number near or at an edo number. Hence, for instance, to plot 12 plot around 108.777, to plot 31 plot around 281.006, and so forth. An alternative plotter is the applet here.

If you have access to Mathematica, which has Z, zeta and theta as a part of its suite of initially defined functions, you can do even better. Below is a Mathematicia-generated plot of Z(2πx/ln(2)) in the region around 12edo:

The peak around 12 is both higher and wider than the local maximums above 11 and 13, indicating its superiority as an edo. Note also that the peak occurs at a point slightly larger than 12; this indicates the octave is slightly compressed in the zeta tuning for 12. The size of a step in octaves is 1/x, and hence the size of the octave in the zeta peak value tuning for Nedo is N/x; if x is slightly larger than N as here with N=12, the size of the zeta tuned octave will be slightly less than a pure octave. Similarly, when the peak occurs with x less than N, we have stretched octaves.

For larger edos, the width of the peak narrows, but for strong edos the height more than compensates, measured in terms of the area under the peak (the absolute value of the integral of Z between two zeros.) Note how 270 completely dominates its neighbors:

Note that for one of its neighbors, 271, it isn't entirely clear which peak value corresponds to the line of real values from +∞. This can be determined by looking at the absolute value of zeta along other s values, such as s=1 or s=3/4, and in this case the local minimum at 271.069 is the value in question. However, other peak values are not without their interest; the local maximum at 270.941, for instance, is associated to a different mapping for 3.

Zeta EDO lists

If we examine the increasingly larger peak values of |Z(x)|, we find they occur with values of x such that Z'(x) = 0 near to integers, so that there is a sequence of edos 1, 2, 3, 4, 5, 7, 10, 12, 19, 22, 27, 31, 41, 53, 72, 99, 118, 130, 152, 171, 217, 224, 270, 342, 422, 441, 494, 742, 764, 935, 954, 1012, 1106, 1178, 1236, 1395, 1448, 1578, 2460, 2684, 3395, 5585, 6079, 7033, 8269, 8539, 11664, 14348, 16808, 28742, 34691 ... of zeta peak edos. This is listed in the On-Line Encyclopedia of Integer Sequences as OEIS Sequence A117536.

Similarly, if we take the integral of |Z(x)| between successive zeros, and use this to define a sequence of increasing values for this integral, these again occur near integers and define an edo. This sequence, the zeta integral edos, goes 2, 5, 7, 12, 19, 31, 41, 53, 72, 130, 171, 224, 270, 764, 954, 1178, 1395, 1578, 2684, 3395, 7033, 8269, 8539, 14348, 16808, 36269, 58973 ... This is listed in the OEIS as OEIS Sequence A117538. The zeta integral edos seem to be, on the whole, the best of the zeta function sequences, but the other two should not be discounted; the peak values seem to give more weight to the lower primes, and the zeta gap sequence discussed below to the higher primes.

Finally, taking the midpoints of the successively larger normalized gaps between the zeros of Z leads to a list of zeta gap edos. These are 2, 3, 5, 7, 12, 19, 31, 46, 53, 72, 270, 311, 954, 1178, 1308, 1395, 1578, 3395, 4190 ... Since the density of the zeros increases logarithmically, the normalization is to divide through by the log of the midpoint. These edos are listed in the OEIS as OEIS Sequence A117537. The zeta gap edos seem to weight higher primes more heavily and have the advantage of being easy to compute from a table of zeros on the critical line.

Optimal Octave Stretch

Another use for the Riemann zeta function is to determine the optimal tuning for an EDO, meaning the optimal octave stretch. This is because the zeta peaks are typically not integers. The fractional part can give us the degree to which the generator diverges from what you would need to have the octave be a perfect 1200 cents. Here is a list of successively higher zeta peaks, taken to five decimal places:

0.00000

1.12657

1.97277

3.05976

3.90445

5.03448

6.95669

10.00846

12.02318

18.94809

22.02515

27.08661

30.97838

40.98808

52.99683

71.95061

99.04733

117.96951

130.00391

152.05285

170.99589

217.02470

224.00255

270.01779

341.97485

422.05570

441.01827

494.01377

742.01093

764.01938

935.03297

953.94128

1012.02423

1105.99972

1177.96567

1236.02355

1394.98350

1447.97300

1577.98315

2459.98488

2683.99168

3395.02659

5585.00172

6079.01642

7032.96529

8268.98378

8539.00834

11664.01488

14347.99444

16807.99325

28742.01019

34691.00191

Removing primes

The Euler product for the Riemann zeta function is

$\zeta(s) = \prod_p (1 - p^{-s})^{-1}$

where the product is over all primes p. The product converges for values of s with real part greater than or equal to one, except for s=1 where it diverges to infinity. We may remove a finite list of primes from consideration by multiplying ζ(s) by the corresponding factors (1-p^(-s)) for each prime p we wish to remove. After we have done this, the smallest prime remaining will dominate peak values for s with large real part, and as before we can track these peaks backwards and, by analytical continuation, into the critical strip. In particular if we remove the prime 2, (1-2^(-s))ζ(s) is now dominated by 3, and the large peak values occur near equal divisions of the "tritave", ie 3.

Along the critical line, |1 - p^(-1/2-i t)| may be written

$\sqrt{1 + \frac{1}{p} - \frac{2 \cos(\ln p\ t)}{\sqrt{p}}}$

Multiplying the Z-function by this factor of adjustment gives a Z-function with the prime p removed from consideration. Zeta peak and zeta integral tunings may then be found as before.

Removing 2 leads to increasing adjusted peak values corresponding to the division of 3 (the "tritave") into 4, 7, 9, 13, 15, 17, 26, 32, 39, 45, 52, 56, 71, 75, 88, 131, 245, 316 ... parts. A striking feature of this list is the appearance not only of 13edt, the Bohlen-Pierce division of the tritave, but the multiples 26, 39 and 52 also.

The Black Magic Formulas

When Gene Smith discovered these formulas in the 70s, he thought of them as "black magic" formulas not because of any aura of evil, but because they seemed mysteriously to give you something for next to nothing. They are based on Gram points and the Riemann-Siegel theta function θ(t). Recall that a Gram point is a point on the critical line where ζ(1/2 + ig) is real. This implies that exp(iθ(g)) is real, so that θ(g)/π is an integer. Theta has an asymptotic expansion

$\theta(t) \sim \frac{t}{2}\log \frac{t}{2\pi} - \frac{t}{2} - \frac{\pi}{8}+\frac{1}{48t}+ \frac{7}{5760t^3}+\cdots$

From this we may deduce that θ(t)/π ≈ r ln(r) - r - 1/8, where r = t/2π = x/ln(2); hence while x is the number of equal steps to an octave, r is the number of equal steps to an "e-tave", meaning the interval of e, 1200/ln(2) = 1731.234 cents.

Recall that Gram points near to pure-octave edos, where x is an integer, can be expected to correspond to peak values of |ζ| = |Z|. We can find these Gram points by Newton's method applied to the above formula. If r = x/ln(2), and if n = floor(r ln(r) - r + 3/8) is the nearest integer to θ(2πr)/π, then we may set r⁺ = (r + n + 1/8)/ln(r). This is the first iteration of Newton's method, which we may repeat if we like, but in fact no more than one iteration is really required. This is the first black magic formula, giving an adjusted "Gram" tuning from the orginal one.

For an example, consider x = 12, so that r = 12/ln(2) = 17.312. Then r ln(r) - r - 1/8 = 31.927, which rounded to the nearest integer is 32, so n = 32. Then (r + n + 1/8)/ln(r) = 17.338, corresponding to x = 12.0176, which means a single step is 99.853 cents and the octave is tempered to twelve of these, which is 1198.238 cents.

The fact that x is slightly greater than 12 means 12 has an overall sharp quality. We may also find this out by looking at the value we computed for θ(2πr)/π, which was 31.927. Then 32 - 31.927 = 0.0726, which is positive but not too large; this is the second black magic formula, evaluating the nature of an edo x by computing floor(r ln(r) - r + 3/8) - r ln(r) + r + 1/8, where r = x/ln(2). This works more often than not on the clearcut cases, but when x is extreme it may not; 49 is very sharp in tendency, for example, but this method calls it as flat; similarly it counts 45 as sharp.

Computing zeta

There are various approaches to the question of computing the zeta function, but perhaps the simplest is the use of the Dirichlet eta function which was introduced to mathematics by Johann Peter Gustav Lejeune Dirichlet, who despite his name was a German and the brother-in-law of Felix Mendelssohn.

The zeta function has a simple pole at z=1 which forms a barrier against continuing it with its Euler product or Dirichlet series representation. We could subtract off the pole, or multiply by a factor of (z-1), but at the expense of losing the character of a Dirichlet series or Euler product. A better method is to multiply by a factor of (1-2^(1-z)), leading to the eta function:

$\eta(z) = (1-2^{1-z})\zeta(z) = \sum_{n=1}^\infty (-1)^{n-1} n^{-z} = \frac{1}{1^z} - \frac{1}{2^z} + \frac{1}{3^z} - \frac{1}{4^z} + \cdots$

The Dirichlet series for the zeta function is absolutely convergent when s>1, justifying the rearrangement of terms leading to the alternating series for eta, which converges conditionally in the critical strip. The extra factor introduces zeros of the eta function at the points 1 + 2πix/ln(2) corresponding to pure octave divisions along the line s=1, but no other zeros, and in particular none in the critical strip and along the critical line. The convergence of the alternating series can be greatly accelerated by applying Euler summation.

Relationship to Harmonic Entropy

The expression

$|\zeta(0.5+it)|^2 \cdot \overline {\phi(t)}$

is, up to a flip in sign, the Fourier transform of the unnormalized Harmonic Shannon Entropy for $N=\infty$, where $\phi(t)$ is the characteristic function (aka Fourier transform) of the spreading distribution and $\overline {\phi(t)}$ denotes complex conjugation.

Note that in the most common case where the spreading distribution is symmetric (as in the case of the Gaussian and Laplace distributions), the characteristic function is purely real and hence the conjugate is unnecessary. In particular, when the spreading distribution is a Gaussian, the characteristic function is also a Gaussian.

More can be found at the page on Harmonic Entropy, including a generalization to Renyi entropy for arbitrary $a$.