Project Euler Problem 27 Solution

Question

Euler published the remarkable quadratic formula:

n2+n+41\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, 402+40+41=40(40+1)+4140^2 + 40 + 41 = 40(40 + 1) + 41 is divisible by 41, and certainly when n = 41, 412+41+4141^2 + 41 + 41 is clearly divisible by 41.

Using computers, the incredible formula n279n+1601n^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:

n2+an+b, where a<1000 and b<1000\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.

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