Disclaimer

I’m not a Mathematician by trade. Even though a lot of the stuff I work on is fairly mathematical (mostly DSP) it has a slightly different feel than pure numbers. So this is my first foray into Number Theory, and also my first encounter with $\LaTeX$, so kindly excuse any formatting issues :)

The setup

Is there a pattern in Prime Numbers? A notion that has always fascinated me but honestly never gave much attention to. I figured the search for a pattern has kept Mathematicians (much brighter than I) busy since the dawn of reason, so what more can be gleaned? Then I stumbled upon the beautiful Riemann hypothesis. Understanding the relationship of prime numbers and the Zeta function ($\zeta(s)=\sum _{n=1}^{\infty } n^{-s}$) has become a little bit of an obsession.

On the surface, primes seem pretty simple to grasp. It’s only when you start asking really big questions that you start to see how impossibly complex things get. Below is a simple prime number generator (using trial division) showing how simple the main concept is. This algorithm is fine as long as you’re only counting small, when looking for massive primes more sophisticated methods must be used.

is_prime(n) =
  i = 2
  while i * i <= n
    if n mod i == 0
      return 0
    i += 1
  return 1

primes = []
for n in 1..100
  if is_prime(n)
    primes << n

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, …

We care about primes because they’re critical to countless scientific fields, among them is cryptography. Currently most encryption techniques used online (namely RSA) completely rely on the simple fact that it’s hard to factor large semiprimes. If someone were to find a closed-form equation it would compromise online security as we know it. Though our Illuminati friends might have already done that to some degree. And no, that link isn’t a sleight of hand, the NSA actually owns the URL for Illuminatti spelled backwards :)

The first 100 primes, spy the pattern yet?

While looking into primes I started with traditional DSP techniques to gain some insight, things like Fourier analysis, Wavelets, etc… These are great tools for periodic signals, however, primes might be one of the most aperiodic sequences in nature. So in the end not much insight was gained by frequency analysis. It’s as if the prime sequence is just an echo of some other grand process.

Riemann’s hypothesis does hint at some regularity but in the most incredible way. First lets define some functions that are needed: Möbius function $\mu(n)$ which is used for it’s inversion properties. Next is the Prime-counting function $\pi(x)$ that exactly counts all primes $<=x$. Finally is the Logarithmic integral function $\text{li}(z)$ which is used as smooth the approximation to $\pi(x)$.

Basically Riemann claims that all the zeros (of his Zeta function $\zeta(s)$) along the critical-line, essentially ‘predict’ the distribution of primes. What this means is that by subtracting an error term $R(x)=\sum _n^{\infty} \frac{\mu(n) \text{li}\left(x^{1/n}\right)}{n}$ (wrt. the locations of the zeros) from $\text{li}(z)=\int _0^z\frac{dt}{\log}$ you can exactly calculate the prime counting function $\pi(x)=R(x)-\sum _p R(x)^p$. It’s an amazing, seemingly unrelated connection that Riemann came up with all the way back in 1859!

Unfortunately his hypothesis is still just that, unproven, though verified numerically for all practical purposes. If anyone can find a zero off of the critical-line ($\zeta(0.5 + i t)=0$), or conversely prove they all exist on the critical-line, they can win a million bucks via the Millennium Prize.

Some plots showing their relationships

Even with the Zeta function and all the research that has followed, we still haven’t gained the insight that we’re all looking for. This is because the location of the zeros are still seemingly random. Here’s the Zeta function and zeros:

And the first 100 locations of the zeros, notice how this function isn’t smooth because it trully encodes the actual error of $\text{li}(z)$:

The analysis

Another relative of primes is the Divisor function, denoted as $\sigma _k(n)$, which is the sum of the $k^\text{th}$ powers of $n$’s divisors. The special case where $k=0$ ie. $\sigma _0(n)=\sum _{i=1}^n \left\lfloor \frac{n}{i}\right\rfloor -\left\lfloor \frac{n-1}{i}\right\rfloor$, essentially counts all integer divisors of a given number $n$. Here’s the first 100 values of $\sigma _0(n)$:

1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6, 2, 4, 4, 5, 2, 6, 2, 6, 4, 4, 2, 8, 3, 4, 4, 6, 2, 8, 2, 6, 4, 4, 4, 9, 2, 4, 4, 8, 2, 8, 2, 6, 6, 4, 2, 10, 3, 6, 4, 6, 2, 8, 4, 8, 4, 4, 2, 12, 2, 4, 6, 7, 4, 8, 2, 6, 4, 8, 2, 12, 2, 4, 6, 6, 4, 8, 2, 10, 5, 4, 2, 12, 4, 4, 4, 8, 2, 12, 4, 6, 4, 4, 4, 12, 2, 6, 6, 9

There’s some neat properties that come out of this, one of the most useful for our purposes is the fact $\sigma _0(n) = 2$ when $n$ is a prime. Thus, every ‘2’ in the preceding sequence is a prime number by definition as it only has 2 divisors, 1 and itself. So this is a nice way to analytically understand what a prime number is.

Let’s take a closer look at a simpler form of $\sigma _0(n)$ where we ignore the second term in the sum leaving us with $\sum _{i=1}^n \left\lfloor \frac{n}{i}\right\rfloor$:

Here we show a listing of all the iterations of the sum. Each row represents the integer $n$ (the 1st column in red). The green column is whether or not $n$ is prime, magenta is the actual sum ($\sum _{i=1}^n \left\lfloor \frac{n}{i}\right\rfloor$), and blue the difference of the current and previous row’s sum. By taking the difference of sums we arrive at the same result as $\sigma _0(n)$, so every ‘2’ in the blue column represents a prime. It’s interesting to note that we can now see that primes share all divisors with the previous row except for the last which of course is itself, ie. $i=n$.

Next, lets see what happens if we apply a simple transform that counts all the similar terms on each row.

Here we have the same table but now with each number in the divisor’s location (black columns) representing the count of repeated divisors. For example let’s take $n=6$ which has the values {3, 1, 1, 0, 0, 1}. The first element ‘3’ represents how many 1s there were in the previous table at $n=6$, the next ‘1’ is how many 2s, then how many 3s, etc… Now, due to the fact that all the primes share the same divisors as the previous row, which still holds given our transformation if we ignore the 1st column (the 1s) and is easily verified visually. One thing to note though, is that when you sum all the numbers in a row the result is equal to $n$. This means that we can’t simply take the difference of sums (relative to the previous) row anymore, instead we must sum the absolute difference:

Which is what the blue column in the listing actually represents. ‘2’ still represents primes. It might not be obvious yet, but now there is a pattern in the rows. It’s easier to see if we rotate the table horizontally, along with some annotations:

Let $n, m$ represent the $n^\text{th}$ dot and $m^\text{th}$ row respectively, $P_{\text{0}}(n,m)$ equal the position of the green dots, $P_{\text{1}}(n,m)$ the position of the red dots, $L_{\text{0}}(m)$ the length of the yellow boxes, and $L_{\text{1}}(m)$ the length of the red boxes. Also, from this point onward lets assume that the origin starts at the first green dot in the upper-left with the value $(0, 0)$ but remember has the global value $(1, 2)$.

The first feature of interest is that no primes fall on any of the dots. This gives us a hint as to what the dots represent and that seems to be an overflow or carry. Remember that each row represents a divisor and the divisors are all related by their relative positions. So what we are seeing is modular arithmetic (of sorts) long-hand. The dots are where the divisors completely fit it seems.

So can we calculate all of these functions?

It turns out that you can and it’s pretty straightforward. Lets start with the simplest function which is $L_{\text{0}}(m)$. It’s easy to see that each row $L_{\text{0}}(m)$ increases by 2 and thus $L_{\text{0}}(m) = 2 (2 + m) + 1$. Note, I’m including the red dots (via the +1) as it makes the following definitions cleaner.

$L_{\text{1}}(m)$ is a bit more involved but it should be clear that it’s included twice, with the 2nd half being a reflection of the 1st, this means that we only need to focus on the first half. First lets define the sub-lengths (blue and magenta boxes) of $L_{\text{1}}(m)$ as $S(n,m)$. I printed out a much larger version of these tables for analysis and eventually came up with this pattern for the sub-lengths:

What’s going on here is that we have a very specific ordering such that, each pair (column-wise) always sum to a value greater than the row’s index (values in red). For example take row $n=10$ which has the pairs: {{1, 10}, {2, 9}, {3, 8}, {4, 7}, …}, if we sum each pair we always get a value > 10. Each row essentially represents all the permutations that cause an overflow. By looking at the data it’s clear that $S(n,m)$ is composed of 2 simpler equations:

We can now obtain $L_{\text{1}}(m)$ by a simple summation of $S(n,m)$ (and inclusion of the green dots). Each row increases the number of sub-lengths by 1 and thus:

Calculating $P_{\text{0}}(n,m)$ is now a simple matter of combining $L_0(i)$ and $L_1(i)$ as the total length $T(m)$:

Calculating $P_{\text{1}}(n,m)$ is trivial as it’s simply $P_{\text{0}}(n,m) + L_{\text{1}}(n,m)$.

Putting it all together

Now that we have closed-forms what can we do with them? Well one thing you could do is calculate prime numbers without division by keeping track of all the necessary rows needed for your current iteration (which turns out to be $\sqrt(n)$) and checking for the presence of dots. For large primes this method could really be an advantage because as you can see there are large sections that are the same, so you can take giant steps for bigger rows, and of course steps of 6 for the top row. Also you can take advantage of the fact that the divisors usually differ toward the upper rows, so you don’t have to go too far down to see if a number is prime or not. This is all by using simple addition and subtraction. You can also calculate a linear recurrence kernel for each row which turns out to be:

k0={            0, 1, 1, 0, -1}, x0={            1, 0, 1, 1, 1}
k1={      0, 0, 1, 1, 0, 0, -1}, x1={      1, 0, 0, 1, 1, 0, 1}
k2={0, 0, 0, 1, 1, 0, 0, 0, -1}, x2={1, 0, 0, 0, 1, 1, 0, 0, 1}
k3=...                         , x2=...

Though doing a convolution per row is definitely not the fastest method, it’s still interesting mathematically. Another idea is to only check if our current location includes a dot or not (no primes exist at these locations), and we have the equations so this should be pretty straightforward. First we need to create a periodic version of the equations so we can analytically calculate all dots. Lets see what happens if we fit sin waves to the dot locations for each row:

Here we have sin waves with zeros at the exact location of where the dots are (for the first several rows). In a way this helps show that primes aren’t random, there is definite structure it’s just on an unbelievably massive scale! Lets define this wave function as $W(x,m)$:

Now we can check for a dot at any number x for any row m. If we take the product of all the waves needed for a given number we get a function that can tell us if we’re on a prime or not, right?

Here we plot the differences of $Q(x)$ and the is_prime(x) function. At first glance it’s not exact, but if we look closer the errors happen at interesting locations! If we factor the outliers we get {5^2, 7^2, {5, 17}, {7, 13}} which shows that the differences we get are always semiprime (with 2 prime factors). So our function $Q(x)$ actually calculates all prime and semiprimes! The interesting thing is that as $x\to \infty$ the ratio of $primes / semiprimes \to 1$.

A more direct approach

What would happen if we directly took the product of differences from our dot functions instead of using sin waves? We still need 2 functions, one for each dot:

Note, subscript q represents the Pochhammer symbol. These functions do in fact give us the desired results but their values are massive. Lets try a more constraining form:

Where $H(x)$ is the Harmonic number and $\psi ^{(0)}$ is the Polygamma function both of which are related to the Zeta function. So maybe we’re on to something here?

Out of the entire equation, it turns out that we are mostly interested in $\psi ^{(0)}$, as it is the oscillatory component. Above is a plot of $\psi ^{(0)}\left(\frac{m-q}{T(m)}\right)$ for n=0..100, m=0..4. Interesting how it almost looks as if it’s a propagating wave/ripple starting opposite the origin… It’s clear that the asymptotes are located right on the zeros (or dots). We can transform the poles into zeros by taking the inverse and scaling by $e$ ie. $W_0(q,m)=\frac{e}{\psi ^{(0)}\left(\frac{m-q}{T(m)}\right)}$, we can now use this function in the same way as the sin based method:

Which gives us the same exact result when subtracted from is_prime(q):

Going further

The math is mostly here to help see what’s going on. If you wanted to implement this you probably wouldn’t want to be making tons of gamma library calls! Instead you would use some iterative technique based on the simpler, original functions. I’ve already proven this with a fairly easy implementation.

There’s another connection I made by using this perspective and it’s pretty amazing, to the point where I didn’t believe it myself at first!

I’m not going to go into complete detail as I’d like a chance to fully optimize and test everything, and be able to take credit for it! However, I’ll give some hints and maybe someone out there has already done this, but I’ve scoured the web looking for anything remotely related and have found nothing. So the first hint I’ll give is that we need to look at the prime numbers in binary form (along with the index):

{{2, 10}, {3, 11}, {5, 101}, {7, 111}, {11, 1011}, {13, 1101}, {17, 10001}, {19, 10011}, {23, 10111}, {29, 11101}, {31, 11111}, {37, 100101}, {41, 101001}, {43, 101011}, {47, 101111}, {53, 110101}, {59, 111011}, {61, 111101}, {67, 1000011}, {71, 1000111}, {73, 1001001}, {79, 1001111}, {83, 1010011}, {89, 1011001}, {97, 1100001}, {101, 1100101}, {103, 1100111}, {107, 1101011}, {109, 1101101}, {113, 1110001}, {127, 1111111}, {131, 10000011}, {137, 10001001}, {139, 10001011}, {149, 10010101}, {151, 10010111}, {157, 10011101}, {163, 10100011}, {167, 10100111}, {173, 10101101}, {179, 10110011}, {181, 10110101}, {191, 10111111}, {193, 11000001}, {197, 11000101}, {199, 11000111}, {211, 11010011}, {223, 11011111}, {227, 11100011}, {229, 11100101}}

The second hint is that instead of summing the absolute difference:

We take the product:

We arrive at a function that is very similar to the Möbius function $\mu(n)$. Now armed with that info, would it be possible to decompose a fixed set of primes ie. all the primes up to a certain $N$ where $N=2^k$, into a simple linear combination based on the binary coefficients of the index? It turns out that if we use 8 basis functions each a neighbor of the previous we can. I’m not saying why we need 8 just yet :)

l0(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc)   (1 - d) (1 - e) (1 - f)
l1(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc)   (1 - d) (1 - e) f
l2(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc)   (1 - d) e       (1 - f)
l3(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc)   (1 - d) e       f
l4(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc) d (1 - e)         (1 - f)
l5(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc) d (1 - e)         f
l6(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc) d         e       (1 - f)
l7(a, b, c, d, e, f) = (i + ja + kb + lab + mc + nac + obc) d         e       f

Where {i, j, k, l, m, n, o} are the unknowns and are unique per row. To calculate the result we simply sum all the basis functions, this of course could be done with matrices just as well:

prime(x) =
  sum = 0
  for i in 0..7
    sum += dec(l_i, x)
  return sum

Just a couple utility functions are needed, bit(x, k) returns the bit of x indexed by k, and dec(f, x) decomposes the integer x into 6 bits and calls the given function.

bit(x, k) = floor(2^-k x) mod 2
dec(f, x) = f(bit(x, 3), bit(x, 4), bit(x, 5), bit(x, 2), bit(x, 1), bit(x, 0))

If you’re curious what the coefficients might look like (and I know you are :) here they are:

Using those we can define the closed-form as:

l(a, b, c, d, e, f) =
  2 + 21a + 57b + 17ab + 135c + 21ac + 33bc + (1 + 5a + b - 3ab + c - 5ac - bc) z +
    y (3 + 5a + 5b - 7ab + 9c - 5 ac - 11bc + (1 - a + b - ab - c + ac + 3bc) z) +
    x (9 + 9a + 5b - 11ab + 11c - 11ac - 11bc + (1 - 5a + 3b + ab + 3c + ac + bc) z +
       y (3 - 5a - b + 15ab - 5c + 9ac + 13bc + (-1 + 5a - b - 3ab + c + 5ac - 7bc) z))

And here is the perfect reconstruction of primes in terms of the 8 basis functions.

The question I leave to you

Given $N=2^k$, is there currently any known algorithm(s) that can generate all the {i, j, k, l, m, n, o, …} unknowns, perhaps not in linear time, but close? Mainly I wanted to open this up for discussion to see what others thought, and to see if there’s currently any related work.

I plan to continue writing more on this topic so check back. And thanks for reading if you made it this far!

chris@jasuto.com