Project Euler Solution 58: Spiral primes

Problem 58: Spiral primes is again about a number spiral.

Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.

37 36 35 34 33 32 31 38 17 16 15 14 13 30 39 18 5 4 3 12 29 40 19 6 1 2 11 28 41 20 7 8 9 10 27 42 21 22 23 24 25 26 43 44 45 46 47 48 49

It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.

If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?

This has similarities with Solution 28: Number spiral diagonals where one had to just add up the numbers on the diagonals for 500 windings.

We can use the same logic to produce an iterator over these windings. Let's start with a test for that:

def test_iter_diagonals() -> None:
    expected = [[3, 5, 7, 9], [13, 17, 21, 25], [31, 37, 43, 49]]
    actual = list(itertools.islice(iter_diagonals(), len(expected)))
    assert actual == expected

Writing the actual function then becomes pretty straightforward:

def iter_diagonals() -> Iterator[int]:
    number = 1
    for winding in itertools.count(1):
        elements = []
        for corner in range(4):
            number += 2 * winding
            elements.append(number)
        yield elements

We can then just go through all these diagonals, see whether the numbers are primes, and keep a ratio of primes:

def solution() -> int:
    num_primes = 0
    num_total = 1
    for winding, numbers in enumerate(iter_diagonals(), start=1):
        for number in numbers:
            if is_prime_accelerated(number):
                num_primes += 1
        num_total += 4
        if 10 * num_primes < num_total:
            return 2 * winding + 1

The crux is having a fast is_prime_accelerated here. Using the sieve like we usually do isn't sufficient. Also the function which checks for all prime factors up to the number is not good enough. What we rather want to do is only go up to the square root of the candidate. This makes it much faster.

def is_prime_accelerated(number: int) -> bool:
    for divisor in prime_generator():
        if number % divisor == 0:
            return False
        if divisor * divisor > number:
            return True

As always, we have used the prime_generator from Solution 3: Largest prime factor.

When we let it run, it takes 1.9 s to produce the solution. The fourth numbers in the last winding considered is 688,590,081. We would have to generate primes up to a billion in order to use the sieve here. And tthat takes too long. But checking with prime factors up until the square root means that we only have to generate primes up until 26,241; and that is pretty fast with the prime generator.