Question

Euler published the remarkable quadratic formula:

$\displaystyle 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.

$\displaystyle 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.

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   0m8.044s
user   0m7.948s
sys    0m0.040s

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 0m29.248s user 0m29.068s sys 0m0.012s 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   0m10.635s
user   0m10.324s
sys    0m0.240s

Questions? Comments? Send me an email: z@chdenton.com