views:

835

answers:

4

Hi all,

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.

+1  A: 

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.

Brian
Ah, that could be helpful... if there's no other obvious improvement to speed, I'll certainly give that a go. Thanks. I still need to find where the fault is first however...
Noldorin
At least use the sieve of Erastothenes.
Svante
@Harleqin: Unfortunately, I can't use the Sieve of Eratosthenes, as I don't know the upper bound of the prime numbers that may be generated (although I could calculate a safe limit, it would have to be done for each pair of coefficients). The algorithm now seems to be running fast enough anyway.
Noldorin
You can implement the sieve of Erastothenes as an infinite generator. The idea is, instead of marking off the multiples of one prime over a certain range, to just "put" this prime on its _next_ multiple. You keep a table of where the primes so far will be encountered next.
Svante
+5  A: 

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.

Dimitre Novatchev
Thanks - your first tip vastly improves the speed of the algorithm! Not sure how I missed something so obvious. Regarding your second point, I believe the code does in fact try n = 0, since the function Seq.init_infinite actually passes 0 as the parameter on the first iteration.
Noldorin
Have you tested the "is_prime" function to verify it really produces all primes belo some N?
Dimitre Novatchev
I've posted my complete solution below, but since two of your pointers have lead me in precisely the right direction, the answer is yours. Thank you again.
Noldorin
+1  A: 

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

Noldorin
A: 

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

LDomagala