In todays post we will be implementing Willans’ Formula. This is an algorithm for finding the nth prime number in the set of prime numbers and was first described by Willans in his 1964 paper "On Formulae for the nth Prime Number" and the first version of the formula provided is as follows:

Don't panic at the sight of the formula, at its core it consists of basic arithmetic and trigonometric functions plus the numbers 1 and 2 and the constant PI. Pseudo code for the formula could be presented as follows:

fn X'(i, n) -> floor (Y'(i, n) ^ (1 / n))

fn Y'(i, n) -> n / sum ( floor ( cos ( ((j - 1)! + 1 / j) * PI ^ 2 ) ) )

n = 5
result = sum ( map i in range(1, power(2, n)) => 1 + X'(i, n) )

Sidenote 1:

This algorithm is not efficient, as per wikipedia:

In addition to the appearance of (j-1)!, it computes p ^ n by adding up p ^ n copies of 1; for example, p ^ 5 = 1+1+1+1+1+1+1+1+1+1+1+0+0...+0 = 11

Willans also presented us alternative formulas, we will be implementing the version shown above but be aware that due to its inefficiency and terrible runtime complexity, you may not be able to generate primes higher than those in a double digit order. Some of the other formulas presented to us by Willans are more efficient but for now we shall focus purely on formula number 1.

Sidenote 2:

For reference, you can see each version of Willans formula here.

Tests

For brevity we will test the numbers from 0 through 10, below are a collection of pytest tests covering these cases. Thus, in test_willans.py we have the following:

from main import willans
from pytest import raises

def test_0Throws():
    with raises(ZeroDivisionError):
       willans(0)

def test_1stPrimeIs2():
    assert willans(1) == 2

def test_2ndPrimeIs3():
    assert willans(2) == 3

def test_3rdPrimeIs5():
    assert willans(3) == 5

def test_4thPrimeIs7():
    assert willans(4) == 7

def test_5thPrimeIs11():
    assert willans(5) == 11

def test_6thPrimeIs13():
    assert willans(6) == 13

def test_7thPrimeIs17():
    assert willans(7) == 17

def test_8thPrimeIs19():
    assert willans(8) == 19

def test_9thPrimeIs23():
    assert willans(9) == 23

def test_10thPrimeIs29():
    assert willans(10) == 29

Implementation

In main.py we have the following implementation of Willans formula:

from math import cos, factorial, floor, pi


def X(i: int, n: int) -> int:
    y = Y(i, n)
    e = 1 / n
    r = pow(y, e)

    return floor(r)


def Y(i: int, n: int) -> float:
    def iterator(j: int) -> int:
        f = factorial(j - 1) + 1
        a = f / j
        b = a * pi
        c = cos(b)
        e = pow(c, 2)

        return floor(e)

    cs = [iterator(j) for j in range(1, i + 1)]

    return n / sum(cs)


def willans(n: int) -> int:
    xs = [X(i, n) for i in range(1, pow(2, n) + 1)]
    s = sum(xs)

    return 1 + floor(s)


lower_bound = 1
upper_bound = 6
result = [willans(n) for n in range(lower_bound, upper_bound + lower_bound)]

print(result)  # [2, 3, 5, 7, 11, 13]

Sidenote 3:

For the implementation of the algorithm I use pythons built in math module but this has issues with values over a certain size when conducting division operations on large integers and floats.

It would be worth using tools such as pandas or scikit learn if you wish to improve the efficiency and scale of the implementation.

Conclusions

I originally believed, as I am sure most people do, that there was no uniform algorithm or formula for generating the nth prime number but Willans does indeed stand up to scrutiny thus far. That said, it is so inefficient and exponentially slow as an algorithm that practically it is useless even if correct.

I came across this algorithm a few weeks ago when watching this fantastic video by Eric Rowland and I highly recommend giving it a watch.

After watching, I did further research as outlined in the links above and came up with the implementation of the formula seen in this article over that time. I find such algorithms and formulae to be utterly fascinating and perhaps in the future I will write more mathematics based algorithms in this series but for now, the article is at an end.

I hope that you found some value in todays post and if you have any questions, comments or suggestions, feel free to leave those in the comments area below the post!