Problem

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.

Answer

-59231

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()

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]

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

Mathematica

CountConsecutivePrimes[a_, b_] := Block[{i}, For[i = 0, PrimeQ[i^2 + a * i + b], i++]; i]

maxConsecutive = 0;
maxProduct = 0;
For[a = -999, a < 1000, a++,
  For[b = -999, b < 1000, b++,
    temp = CountConsecutivePrimes[a, b];
    If[temp > maxConsecutive,
      maxConsecutive = temp;
      maxProduct = a * b]]]
maxProduct

Java

public final class p027 implements EulerSolution {
	
	public static void main(String[] args) {
		System.out.println(new p027().run());
	}
	
	
	public String run() {
		int bestNum = 0;
		int bestA = 0;
		int bestB = 0;
		for (int a = -1000; a <= 1000; a++) {
			for (int b = -1000; b <= 1000; b++) {
				int num = numberOfConsecutivePrimesGenerated(a, b);
				if (num > bestNum) {
					bestNum = num;
					bestA = a;
					bestB = b;
				}
			}
		}
		return Integer.toString(bestA * bestB);
	}
	
	
	private static int numberOfConsecutivePrimesGenerated(int a, int b) {
		for (int i = 0; ; i++) {
			int n = i * i + i * a + b;
			if (n < 0 || !Library.isPrime(n))
				return i;
		}
	}
	
}