Project Euler Solution 27: Quadratic primes

Problem 27: Quadratic primes of the Project Euler series is again about prime numbers.

Euler discovered the remarkable quadratic formula: $$ n^2 + n + 41,. $$

It turns out that the formula will produce 40 primes for the consecutive integer values $0 \lt n \lt 39$. However, when $n = 40$, $40^2 + 40 + 41$ is divisible by 41, and certainly when $n = 41$, $41^2 + 41 + 41$ is clearly divisible by 41.

The incredible formula $n^2 - 79 n + 1601$ was discovered, which produces 80 primes for the consecutive values $0 \lt n \lt 79$. The product of the coefficients, −79 and 1601, is −126479.

Considering quadratics of the form $n^2 + an + b$, where $|a| < 1000$ and $|b| \lt 1000$ where $|n|$ is the modulus/absolute value of $n$ e.g. $|11| = 11$ and $|-4| = 4$.

Find the product of the coefficients, $a$ and $b$, for the quadratic expression that produces the maximum number of primes for consecutive values of $n$, starting with $n = 0$.

The lattice spanned by $a$ and $b$ just has a million entries, so it should be possible to just test all combinations.

We take the prime generator from Solution 3: Largest Prime Factor and build a prime number set from that. With that we can easily test whether a given number is a prime number. If the set doesn't contain a prime number higher than the candidate, it will fill up the set further and then do the check. This way we don't need to specify an upper limit beforehand, but also only compute as many primes as we really need.

class PrimeSet:
    def __init__(self) -> None:
        self._primes: set[int] = []
        self._largest: int = 0
        self._prime_iterator = prime_generator()

    def contains(self, candidate: int) -> bool:
        while self._largest < candidate:
            new_prime = next(self._prime_iterator)
            self._largest = new_prime
            self._primes.append(new_prime)
        return candidate in self._primes

To make sure that this part works, we can also write a little test:

def test_prime_set() -> None:
    prime_set = PrimeSet()
    assert prime_set.contains(7)
    assert not prime_set.contains(9)

With that we can just go through all $a$ and $b$ and then through the $n$. As soon as we encounter a number which is not a prime, we check whether that combination $(a, b)$ has yielded more primes than any combination before. If so, we record the new maximum $n$ and also the coefficients $a$ and $b$. In the end we return the product $ab$ of the best combination.

def solution() -> int:
    n_max = 0
    best = (None, None)
    prime_set = PrimeSet()
    for a in range(-999, 1000):
        for b in range(-1000, 1001):
            n_end = 0
            for n in itertools.count():
                candidate = n**2 + a * n + b
                if candidate <= 0:
                    break
                if not prime_set.contains(candidate):
                    break
                n_end += 1
            if n_end > n_max:
                best = (a, b)
                n_max = n_end
    return best[0] * best[1]

This takes a bit longer to compute, 9.6 ms. That is still within one minute, so it should be okay.

One improvement is to start with $b \ge 2$ because for $n = 0$, we only produce $b$. But if $b < 2$, then we produce a negative number, which isn't a prime, and producing 0 or 1 isn't a prime number either. This doesn't make much of a difference, it still takes 8874.338 ms. This isn't surprising because these loops would be ended pretty quickly anyway.

Perhaps there are additional ways to make it faster, I don't see them, though.