The Problem

Euler discovered the remarkable quadratic formula:

n²+n+41n²+n+41

It turns out that the formula will produce 40 primes for the consecutive integer values 0≤n≤39. However, when n=40,40²+40+41=40(40+1)+41 is divisible by 41, and certainly when n=41,41²+41+41 is clearly divisible by 41.

The incredible formula n²−79n+1601 was discovered, which produces 80 primes for the consecutive values 0≤n≤790≤n≤79. The product of the coefficients, −79 and 1601, is −126479.

Considering quadratics of the form:

n²+an+b, where |a|<1000 and |b|≤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 Solution

First, we'll define our boundaries.

[<Literal>]
let maxA = 999L

[<Literal>]
let maxB = 1000L

Next, we'll steal some code from another problem for determing if a 64-bit number is prime.

let int64sqrt (i : int64) = int64(sqrt(float i))

let isPrime64 (i : int64)  =
    if i <= 1L then false
    elif i = 2L then true
    elif (i &&& 1L) = 0L then false
    else
        let sr = int64sqrt i
        seq { 3L..2L..sr } |> Seq.forall (fun f -> i%f<>0L)

I can't think of a more efficient way to solve the problem other than testing all possible coefficients, so let's generate all possible pairs of a and b.

let coefficients = [|
    for a in -maxA..maxA do
        for b in -maxB..maxB do
            yield (a,b)
|] 

And we need a function for computing "quadratics of the form n² + an + b".

let quadratic (n : int64) (a : int64) (b : int64) = n*n + a*n + b

Now, for a given a and b we can find the number of primes produced. Here is a recursive function that keeps incrementing n until a non-prime occurs, then returns the count.

let numPrimesProduced (a : int64) (b : int64) =
    let rec numPrimes n =
        if isPrime64 (quadratic n a b) then numPrimes (n + 1L)
        else n

    numPrimes 0L

Finally, we can solve the problem:

#time
let a,b,n =
    coefficients
    |> Seq.map (fun (a,b) -> (a,b,(numPrimesProduced a b)))
    |> Seq.maxBy (fun (a,b,n) -> n)
#time
a*b

With all the coefficient pairs as inputs, we generate the number of primes produced, then find the pair with the maximum n. The answer is a*b.

Other Posts in This Series