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≤n≤39. However, when n = 40, 40^2 + 40 + 41 = 40(40 + 1) + 41 is divisible by 41, and certainly when n = 41, 41^2+ 41 + 41 is clearly divisible by 41.

The incredible formula n^2−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^2 + 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



In [24]:

    
let primeCalcLimit = 20000000

let compositesInRange fromNum toNum =
    [|2 .. toNum/fromNum|]
    |> Array.map (fun x -> fromNum * x)

// start off assuming every number from 0 .. n is prime
let isPrime' = Array.init (1+primeCalcLimit) (fun _ -> true)
   
// simple sieve of eratosthenes - update "primeness" of
// each number in the isPrime array
for x in 2 .. int(sqrt(double(primeCalcLimit))) do 
    if isPrime'.[x] then             
        compositesInRange x primeCalcLimit
        |> Array.map (fun i -> isPrime'.[i] <- false)
        |> ignore

let isPrime x = 
    if x <= 1 then false
    else isPrime'.[x]
    
let square x = x * x

let quadPrime a b n = 
    (square n) + (a * n) + b
    
let rec comb n l = 
    match n, l with
    | 0, _ -> [[]]
    | _, [] -> []
    | k, (x::xs) -> List.map ((@) [x]) (comb (k-1) xs) @ comb k xs

let infiniteIntegers = Seq.initInfinite (fun x -> 1+x)

let continuousPrimeCounts f = 
    infiniteIntegers 
    |> Seq.takeWhile (fun n -> (isPrime (f n)))
    |> Seq.length
    
let maxContPrime (maxLength, maxProduct) (currentProduct, currentFunc) =
    let currentLength = continuousPrimeCounts currentFunc
    if currentLength > maxLength then (currentLength, currentProduct)
    else (maxLength, maxProduct)

comb 2 [-1000..1000]
|> List.map (fun a -> ((a.[0] * a.[1]), (quadPrime a.[0] a.[1])))
|> List.fold (fun maxPair currentPair -> maxContPrime maxPair currentPair) (-1,-1)
|> snd









    Out[24]:





-59231



In [ ]: