# Project Euler Problem 27 Solution

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

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