Question

Euler published the remarkable quadratic formula:

$$n^2 + n + 41$$

It turns out that the formula will produce 40 primes for the consecutive values n = 0 to 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.

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

Considering quadratics of the form:

$$n^2 + an + b, \text{ where } |a| \lt 1000 \text{ and } |b| \lt 1000$$

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.

Haskell

import Data.Function (on)
import Data.List (maximumBy)

isPrime :: Int -> Bool
isPrime n | n < 1 = False
          | otherwise = not $ or [n `rem` x == 0 | x <- [2..floor $ sqrt $ fromIntegral n]]
 
coefficients :: [(Int, Int)]
coefficients = [(a, b) | a <- [-999..999], b <- filter isPrime [0..999]]

primesProduced :: (Int, Int) -> Int
primesProduced (a, b) = length $ takeWhile isPrime [n^2 + a*n + b | n <- [0..]]

main :: IO ()
main = print $ uncurry (*) $ maximumBy (compare `on` primesProduced) coefficients

$ ghc -O2 -o quadratic-primes quadratic-primes.hs
$ time ./quadratic-primes
real   0m2.612s
user   0m2.604s
sys    0m0.008s

Python

#!/usr/bin/env python
from collections import defaultdict
from itertools import product
from operator import mul
import math
from functools import reduce

def factorize(n):
    if n < 1:
        raise ValueError('fact() argument should be >= 1')
    if n == 1:
        return []  # special case
    res = []
    # iterate over all even numbers first.
    while n % 2 == 0:
        res.append(2)
        n //= 2
    # try odd numbers up to sqrt(n)
    limit = math.sqrt(n+1)
    i = 3
    while i <= limit:
        if n % i == 0:
            res.append(i)
            n //= i
            limit = math.sqrt(n+i)
        else:
            i += 2
    if n != 1:
        res.append(n)
    return res

def num_divisors(n):
    factors = sorted(factorize(n))
    histogram = defaultdict(int)
    for factor in factors:
        histogram[factor] += 1
    # number of divisors is equal to product of 
    # incremented exponents of prime factors
    from operator import mul
    try:
        return reduce(mul, [exponent + 1 for exponent in list(histogram.values())])
    except:
        return 1

def num_primes(formula):
    num = 0
    for n in range(1000):
        res = formula(n)
        if res < 1 or not is_prime(res):
            return num
        else:
            num += 1
 
def is_prime(num):
    if num_divisors(num) == 2 and num > 1:
        return True
    else:
        return False

def main():
    most = 0
    best = (0, 0)
    for a, b in product(list(range(-999,1000)), list(range(-999, 1000))):
        formula = lambda n: n**2 + a*n + b
        num = num_primes(formula) 
        if num > most:
            most = num
            best = (a, b)
    print(mul(*best))

if __name__ == "__main__":
    main()

$ time python3 quadratic-primes.py
real   0m18.794s
user   0m18.777s
sys    0m0.016s

Ruby

#!/usr/bin/env ruby
require 'mathn'
puts (-999..999).to_a.product((-999..999).to_a).map { |a, b|
  [(0..100).take_while { |n| (n**2 + a*n + b).prime? }.count, a * b]
}.max[1]

$ time ruby quadratic-primes.rb
real   0m5.754s
user   0m5.612s
sys    0m0.142s

Rust

fn is_prime(n: i64) -> bool {
    if n < 1 {
        return false;
    }
    let max = (n as f64).sqrt() as i64;
    for d in 2..max {
        if n % d == 0 {
            return false;
        }
    }
    true
}

fn main() {
    let mut best = (0, 0);
    let mut best_count = 0;
    let bs = (0..1000).filter(|&b| is_prime(b)).collect::<Vec<_>>();
    for a in -999..1000 {
        for &b in &bs {
            let count = (0..)
                .take_while(|n| is_prime(n * n + a * n + b))
                .map(|_| 1)
                .sum::<i64>();
            if count > best_count {
                best = (a, b);
                best_count = count;
            }
        }
    }
    println!("{}", best.0 * best.1);
}

$ rustc -C target-cpu=native -C opt-level=3 -o quadratic quadratic.rs
$ time ./quadratic
real   0m0.042s
user   0m0.042s
sys    0m0.000s