Back to Questions
1.0k
views5
comments
Optimizing a solution for Project Euler Porblem 58
f#

Hello everyone!

This is my first post, so please be kind :).

As many people here do, I try to develop my skills in F# by solving Project Euler problems.

And I have a question about <a href=http://projecteuler.net/index.php?section=problems&id=58> Problem 58 </a> I wrote this code to solve it, and it successfully does the solving. But it runs about 2 full minutes on a modern computer. Could anyone please show me the slowest parts of the code? And maybe ways to improve them? Here is the actual code: <code lang=fsharp> #light let sqrt_int n = n |> float |> sqrt |> int64 let is_prime n = if n < 2L then false else [|2L .. sqrt_int n|] |> Array.exists (fun x -> n % x = 0L) |> not let diff2 (l1, l2) = ((l1 |> float) / (l2 |> float)) * 100. let numOfPrimesInSpiral ratio = let rec loop incBy acc count prev primes both currRatio = match currRatio with //We've finished the current loop of the spiral and current ratio is now below // target ratio -> finish. | _ when currRatio < ratio && count = 0L -> (incBy+1L, ratio) | _ -> match count with | nil when nil = 0L -> //start a new loop of spiral let i = prev + (incBy + 2L) if is_prime i then let prLen = primes |> List.length let bothLen = both |> List.length let r = ((prLen)+1, (bothLen)+1) |> diff2 loop (incBy+2L) (acc+i) 3L i (i::primes) (i::both) r else let prLen = primes |> List.length let bothLen = both |> List.length let r = ((prLen), (bothLen)+1) |> diff2 loop (incBy+2L) (acc+i) 3L i primes (i::both) r | _ -> //cotinue the current loop let i = prev + incBy if is_prime i then let prLen = primes |> List.length let bothLen = both |> List.length let r = ((prLen)+1, (bothLen)+1) |> diff2 loop incBy (acc+i) (count-1L) i (i::primes) (i::both) r else let prLen = primes |> List.length let bothLen = both |> List.length let r = ((prLen), (bothLen)+1) |> diff2 loop incBy (acc+i) (count-1L) i (primes) (i::both) r match ratio with | hundred when hundred = 100. -> (1L, 100.) | over when over > 100. -> failwith "ratio cannot be larger than 100%" | _ -> loop 2L 1L 4L 1L [] [1L] 0. let te = 10. printfn "test on %f is %A" te (numOfPrimesInSpiral te) </code> Thanks in advance! PS: I tried to get the indents right in the code block, but failed to do so even after reading <a href="/forums/post/12824.aspx">this helpful post</a>

.

Thank you for your help, everyone!

Took me about an hour to fully understand cfern's isPrime2 function, and have yet to understand the <a href = "http://cs.hubfs.net/forums/permalink/13544/13575/ShowThread.aspx#13575">Robin-Miller method</a>, <a href="/members/sharpdev.aspx">sharpdev</a> presented. But I'll definitely continue trying :). Thanks again for your help and this great community you make.

By on 4/2/2010 7:07 AM ()

I've solved this problem too, and the biggest speed gain is to be found in the primality test.

Compare these two prime tests, for instance:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

//ye olde-fashioned primality test
let isPrime n = {2L..int64 (sqrt (float n))} |> Seq.forall (fun d -> n%d <> 0L)


//A faster test, using the fact that primes>3 are of the form 6k +/- 1
//It's also faster because it's a (tail recursive) loop. It doesn't need
//to walk along a sequence iterator. 
//It's not as readable as the test above, but speed matters here.
let isPrime2 = function 
    | n when n<=1L -> false
    | 2L | 3L -> true
    | N when N &&& 1L = 0L || N%3L = 0L -> false
    | N -> let rec aux composite i num =
               if composite then false
                   elif i*i > num then true else
                       aux (num % i = 0L || num % (i+2L) = 0L) (i+6L) num
           aux false 5L N

I can bench these prime tests with the following problem 58 code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

let go primetest =
    let isPrimeN N = if primetest N then 1 else 0


    //numCPrimes N = how many corners of a size N spiral are prime. (corner N*N omitted for obvious reasons)
    let numCPrimes N =
        ((isPrimeN (N*N-(N-1L))) + (isPrimeN (N*N-(2L*(N-1L)))) + (isPrimeN (N*N-(3L*(N-1L)))))


    let cornerPrimes ratio = 
        let rec loop acc N =
            match acc with
            | acc when (( (float acc) / (-1. + 2.*(float N)) ) < ratio) -> (N,acc,( (float acc) / (-1. + 2.*(float N)) ))
            | _ -> loop (acc + numCPrimes (N+2L)) (N+2L)
        loop 3 3L


    let (euler058,acc,ratio) = cornerPrimes 0.1
    printfn "Ratio --> %f | Primes --> %i |Side length --> %i" ratio acc euler058

Results:

1
2
3
4
5
6
7
8

> go isPrime;;
Ratio --> 0.099998 | Primes --> 5248 |Side length --> 26241
Real: 00:00:09.394, CPU: 00:00:09.390, GC gen0: 2, gen1: 0, gen2: 0
val it : unit = ()
> go isPrime2;;
Ratio --> 0.099998 | Primes --> 5248 |Side length --> 26241
Real: 00:00:00.391, CPU: 00:00:00.390, GC gen0: 0, gen1: 0, gen2: 0
val it : unit = ()

This demonstrates that the efficiency of the primality test is important here.

By on 3/29/2010 12:43 AM ()

and the biggest speed gain is to be found in the primality test

You're right. With the following (simple) implementation of the Rabin Miller pseudo primality test speed is still increased.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

 

(* Rabin-Miller pseudo primality test*)
let decompose n = 
  let rec decompose' k r =
    if k % 2L = 0L then decompose' (k / 2L) (r + 1L)
    else (r,k)
  decompose' n 0L

let rec powerMod n a q = 
  if q = 0L then 1L
  else 
    let (q',r') = (q / 2L, q % 2L)
    let a' = powerMod n ((a * a) % n) q'
    if r' = 0L then a'
    else (a * a') % n

let oneOfEqualMinusOne n aq r = List.exists ((=) (n - 1L)) (List.scan (fun x _ -> (x * x) % n) aq [1L..(r - 1L)])

let rabinMiller a n =
  if n % 2L = 0L then false
  else 
    let n' = n - 1L
    let (r,q) = decompose n'
    let aq = powerMod n a q
    if aq = 1L then true
    else oneOfEqualMinusOne n aq r

// true primality if n < 2 152 302 898 747
let isPrime n =
  let l = [2L;3L;5L;7L;11L] 
  if List.exists ((=) n) l then true
  else List.forall (fun p -> rabinMiller p n) l

 

let rec pb58 n k = 
  let k' = k + (int64 (List.length (List.filter isPrime [n*n - n+1L;n*n-2L*n+2L;n*n-3L*n+3L])))
  if 10L*k' < 2L*n - 1L then n
  else pb58 (n+2L) k'

printfn "solution : %i" (pb58 3L 0L)
By on 3/30/2010 8:30 AM ()

use a sequence instead of an array in the is_prime function

for example

1
2
3
4
5
6
7
8
9
10
11
12

 

let is_prime n = seq{for i in 2..sqrt_int n -> i} |> Seq.exists (fun x -> n % x = 0) |> not


let rec pb58 n k  = 
  let k' = k + List.length (List.filter is_prime [n*n - n+1;n*n-2*n+2;n*n-3*n+3])
  if 10*k' < 2*n - 1 then n
  else pb58 (n+2) k'

printfn "solution : %i" (pb58 3 0) // solution = 26241
By on 3/28/2010 6:24 PM ()

Just some quickies:

- try "inline"ing the utiluty fn's like sqrt_int, diff2, etc.?
- calls to the O(n) List.length's stand out - use o(1) collection or pass in lengths?
- use higher-order fn's such as map/fold/reduce
- Looks like is_prime didn't get listed correctly - make sure it memoizes if the same # is likely to be queried again?

By on 3/28/2010 11:11 AM ()

    Topic tags

    IntelliFactory Offices Copyright (c) 2011-2012 IntelliFactory. All rights reserved.
    Home | Products | Consulting | Trainings | Blogs | Jobs | Contact Us | Terms of Use | Privacy Policy | Cookie Policy     		Built with WebSharper