Project Euler Problem 27 in F#

Question

I've been trying to work my way through Problem 27 of Project Euler, but this one seems to be stumping me. Firstly, the code is taking far too long to run (a couple of minutes maybe, on my machine, but more importantly, it's returning the wrong answer though I really can't spot anything wrong with the algorithm after looking through it for a while.

Here is my current code for the solution.

/// Checks number for primality.
let is_prime n = 
    [|1 .. 2 .. sqrt_int n|] |> Array.for_all (fun x -> n % x <> 0)

/// Memoizes a function.
let memoize f = 
    let cache = Dictionary<_, _>()
    fun x -> 
        let found, res = cache.TryGetValue(x)
        if found then
            res
        else
            let res = f x
            cache.[x] <- res
            res

/// Problem 27
/// Find a quadratic formula that produces the maximum number of primes for consecutive values of n.
let problem27 n =
    let is_prime_mem = memoize is_prime
    let range = [|-(n - 1) .. n - 1|]
    let natural_nums = Seq.init_infinite (fun i -> i)
    range |> Array.map (fun a -> (range |> Array.map (fun b ->
        let formula n = n * n + a * n + b
        let num_conseq_primes = natural_nums |> Seq.map (fun n -> (n, formula n))
                                |> Seq.find (fun (n, f) -> not (is_prime_mem f)) |> fst
        (a * b, num_conseq_primes)) |> Array.max_by snd)) |> Array.max_by snd |> fst

printn_any (problem27 1000)

Any tips on how to a) get this algorithm actually returning the right answer (I think I'm at least taking a workable approach) and b) improve the performance, as it clearly exceeds the "one minute rule" set out in the Project Euler FAQ. I'm a bit of a newbie to functional programming, so any advice on how I might consider the problem with a more functional solution in mind would also be appreciated.

Answer 1

Two remarks:

1. You may take advantage of the fact that b must be prime. This follows from the fact that the problem asks for the longest sequence of primes for n = 0, 1, 2, ... So, formula(0) must be prime to begin with , but formula(0) = b, therefore, b must be prime.

2. I am not an F# programmer, but it seems to me that the code does not try n= 0 at all. This, of course, does not meet the problem's requirement that n must start from 0, therefore there are neglectable chances a correct answer could be produced.

Answer 2

Right, after a lot of checking that all the helper functions were doing what they should, I've finally reached a working (and reasonably efficient) solution.

Firstly, the is_prime function was completely wrong (thanks to Dimitre Novatchev for making me look at that). I'm not sure quite how I arrived at the function I posted in the original question, but I had assumed it was working since I'd used it in previous problems. (Most likely, I had just tweaked it and broken it since.) Anyway, the working version of this function (which crucially returns false for all integers less than 2) is this:

/// Checks number for primality.
let is_prime n = 
    if n < 2 then false
    else [|2 .. sqrt_int n|] |> Array.for_all (fun x -> n % x <> 0)

The main function was changed to the following:

/// Problem 27
/// Find a quadratic formula that produces the maximum number of primes for consecutive values of n.
let problem27 n =
    let is_prime_mem = memoize is_prime
    let set_b = primes (int64 (n - 1)) |> List.to_array |> Array.map int
    let set_a = [|-(n - 1) .. n - 1|]
    let set_n = Seq.init_infinite (fun i -> i)
    set_b |> Array.map (fun b -> (set_a |> Array.map (fun a ->
        let formula n = n * n + a * n + b
        let num_conseq_primes = set_n |> Seq.find (fun n -> not (is_prime_mem (formula n)))
        (a * b, num_conseq_primes))
    |> Array.max_by snd)) |> Array.max_by snd |> fst

The key here to increase speed was to only generate the set of primes between 1 and 1000 for the values of b (using the primes function, my implementation of the Sieve of Eratosthenes method). I also managed to make this code slightly more concise by eliminating the unnecessary Seq.map.

So, I'm pretty happy with the solution I have now (it takes just under a second), though of course any further suggestions would still be welcome...

Answer 3

You could speed up your "is_prime" function by using a probabilistic algorithm. One of the easiest quick algorithms for this is the Miller-Rabin algorithm.

Answer 4

to get rid of half your computations you could also make the array of possible a´s only contain odd numbers

Answer 5

my superfast python solution :P

flag = [0]*204
primes = []

def ifc(n): return flag[n>>6]&(1<<((n>>1)&31))

def isc(n): flag[n>>6]|=(1<<((n>>1)&31))

def sieve():
    for i in xrange(3, 114, 2):
        if ifc(i) == 0:
            for j in xrange(i*i, 12996, i<<1): isc(j)

def store():
    primes.append(2)
    for i in xrange(3, 1000, 2):
        if ifc(i) == 0: primes.append(i)

def isprime(n):
    if n < 2: return 0
    if n == 2: return 1
    if n & 1 == 0: return 0
    if ifc(n) == 0: return 1
    return 0    

def main():
    sieve()
    store()
    mmax, ret = 0, 0
    for b in primes:
        for a in xrange(-999, 1000, 2):
            n = 1
            while isprime(n*n + a*n + b): n += 1
            if n > mmax: mmax, ret = n, a * b
    print ret

main()