Using Functional Programming to compute prime numbers efficiently

Question

I'm getting acquainted with F# by going over Project Euler and solving some of the problems. Many of the early problems consist of prime numbers. After looking around I came up with the following solution:

let primesL =
    let rec prim n sofar = 
        seq { if (sofar |> List.forall (fun i->n%i <>0L)) then
                  yield n
                  yield! prim (n+1L) (n::sofar)
              else
                  yield! prim (n+1L) sofar  }
    prim 2L []

This works well, but then I generate all the prime numbers up to 2000000:

let smallPrimes = primesL |> Seq.takeWhile (fun n->n<=2000000)

This takes ages. It's quite obvious something is done in O(N^2) or worst.

I know I can write an imperative version and implement a sieve of some sort, but I want to stick to functional code. If I wanted imperative, I would have stayed with C#.

What am I missing?

You're missing _encapsulation_. Write a sieve, and give it a functional _interface_ -- no one needs to know the implementation details. — ildjarn, Mar 19 '12 at 07:59
I should have been clear about the purpose of learning F# - I'm considering it for *some* parts of a big project. The other parts will be C#. If I need to resort to imperative code for most of the calculations, I'll just keep everything in C#. — zmbq, Mar 19 '12 at 08:51
you can implement Sieve of Eratosthenes algorithm using purely functional approach, [here's one in Haskell](http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Basic_list-based_sieve). — jfs, Mar 19 '12 at 10:03

score 8 · Answer 1 · answered Mar 19 '12 at 17:32

8

Rather than write a long answer here, I refer you to Melissa O'Neill's great paper on the sieve of Eratosthenes.

answered Mar 19 '12 at 17:32

Norman Ramsey

198,648
61
360
533

2

Also F# implementation of Melissa O'Neill's approach can be found in [another great post of Dustin Campbell](http://diditwith.net/2009/01/20/YAPESProblemSevenPart2.aspx). – Gene Belitski Mar 19 '12 at 19:22

Gene Belitski · Accepted Answer · 2012-03-19T15:10:43.640

6

You may want to compare your approach with my variant of Problem Euler 10 solution

let rec primes = 
    Seq.cache <| seq { yield 2; yield! Seq.unfold nextPrime 3 }
and nextPrime n =
    if isPrime n then Some(n, n + 2) else nextPrime(n + 2)
and isPrime n =
    if n >= 2 then
        primes 
        |> Seq.tryFind (fun x -> n % x = 0 || x * x > n)
        |> fun x -> x.Value * x.Value > n
    else false

It is purely functional, uses sequence cashing, optimized for primality check; also it yields very useful isPrime n function as a co-result.

And being applied to the original problem

let problem010 () =
    primes
    |> Seq.takeWhile ((>) 2000000)
    |> (Seq.map int64 >> Seq.sum)

it completes in decent 2.5 s. This is not blasting fast, but was good enough to use this primes sequence in handful of my other Project Euler solutions (27, 35, 37, 50, 58, 69, 70, 77 so far).

As to what you're missing in your solution - from your code I believe you're building a brand new sequence on each internal call to prim, i.e. for each natural while my approach uses a single sequence for already found primes and only enumerates its cached instance when producing each next prime.

edited Mar 19 '12 at 15:10

answered Mar 19 '12 at 14:29

Gene Belitski

10,270
1
34
54

This seems both quick *and* concise. I'll give it a try, I just have to make sure I understand this first. – zmbq Mar 19 '12 at 18:35
1

@zmbq: It may help you to look at a [great post](http://diditwith.net/2008/09/17/YAPESProblemSevenPart1.aspx) by Dustin Campbell on the subject; he came up with pretty similar approach much ahead of me )c: – Gene Belitski Mar 19 '12 at 19:17
Cool. It has one optimization you didn't use - instead of using n+2 in nextPrime, I used `let inline next n = if n%6=1 then n+4 else n+2` It saved another 10% of the execution time. – zmbq Mar 19 '12 at 20:26
The trick here is that when isPrime is called for a candidate C, it checks primes up to sqrt(C), and all the primes up to C have already been saved in the cache? That's why it works fast? – zmbq Mar 19 '12 at 20:28

score 1 · Answer 3 · answered Mar 19 '12 at 09:20

1

First, it is O(n^2) - remember that you use List.forall on each iteration.

Second, if you use the generator a lot, you should cache the results (so that each prime number is only calculated once):

let primesL =
    let rec prim n sofar = 
        seq { if (sofar |> List.forall (fun i -> n % i <> 0UL)) then
                  yield n
                  yield! prim (n + 1UL) (n::sofar)
              else
                  yield! prim (n + 1UL) sofar }
    prim 2UL []
    |> Seq.cache

answered Mar 19 '12 at 09:20

Ramon Snir

7,520
3
43
61

I do not think placing `Seq.cache` where you placed it may help to speedup things as the cause of slowness is in implementation of `prim`. – Gene Belitski Mar 19 '12 at 15:38
Seq.cache definitely didn't help prim, but it helps if someone ever uses primesL more than once. – zmbq Mar 19 '12 at 18:33
I didn't say it would make the computation go faster, I only said it would help on second or later use. – Ramon Snir Mar 19 '12 at 19:33

score 1 · Answer 4 · edited May 23 '17 at 12:32

1

What you really want here is a sieve - I have written a pretty fast F# sieve before here:

F# parallelizing issue when calculating perfect numbers?

edited May 23 '17 at 12:32

Community

1
1

answered Mar 19 '12 at 09:43

John Palmer

25,356
3
48
67

score -1 · Answer 5 · answered Mar 19 '12 at 09:37

-1

Use "Miller–Rabin primality test" for more than some large prime number

answered Mar 19 '12 at 09:37

BLUEPIXY

39,699
7
33
70

This is not a good solution if you want to test the primality of a large amount of numbers - a sieve is faster – John Palmer Mar 19 '12 at 09:44
In large numbers, it is pulling his leg seems to be like List.forall so obvious, Be replaced by high-speed determination that part, (enough to downvote?) so bad I did not think to... – BLUEPIXY Mar 20 '12 at 12:25

Using Functional Programming to compute prime numbers efficiently

5 Answers5

Linked