Back to Questions
2.2k
views8
comments
Need some help with Seq.unfold
f#

Hello! I'm trying to compose an infinite list of prime numbers. I have the following code.

1
2
3
4
5
6
7
8

let primes =
   Seq.unfold
      (fun x ->
         if (is_prime x) then
            Some(x, x + 1)
         else
            Some(0, x + 1))
      2

is_prime is a function which returns a boolean value to whether the number is a prime or not. If you read the code, you'd see that I'm adding 0s to the sequence when the number is not prime. Later on I filter out the zeroes, and I have a sequence of only primes.

How else would I go about doing this?

PS: I'm still trying to wrap my head around FP, so I may not be understanding some concepts.

.

EFFICIENTLY generating a sequence of primes using functional techniques seems quite challenging.
An efficient algorithm would take advantage of two facts about primes not yet mentioned:

Fact 1. Testing only needs to be done against other PRIMES; so one should test against the list of primes built so far (stopping at sqrt(n) as observed). For example, any number that is divisible by 6 is also divisible by the factors of 6 (2, 3), so once you've tested against 2 & 3, it is a waste to also test against 6.

The natural and efficient way to build such a list (that is used while it is being

built) is to PREPEND each new prime to the front of the list, and then pass in the head of the growing list; e.g. when trying 12, the primes list we would have is [11; 7; 5; 3; 2]. But this won't work efficiently because of Fact 2 below.

Fact 2. It is more efficient to test against small primes first (many more numbers are divisible by 2 or 3 than are divisible by some large prime you just found.) That is, if you test against 2 first, you will quickly reject most non-primes, saving much time. So you should test against 2, then 3, then 5, ... stopping at sqrt(n).

So now we have a problem. We need a storage area that holds a growing list of primes, but we need to APPEND each new prime to the tail. This is inherently a changing structure; it does not lend well to applicative (functional) techniques.

The traditional imperative solution is to use a Resizable Array. I see that F# has a ResizeArray, but in keeping with its functional focus, I don't see a "push" operation to append to the tail. No doubt one could use a .Net System.Array since those are resizable, but I thought I'd throw out this challenge:

Any suggestions as to a functional prime-finding algorithm that makes use of the two facts above?
(WARNING: Don't spend too long trying to solve this if you have more important things to be doing!)

~TMSteve

By on 3/16/2008 9:52 PM ()

By making a table for the first 1000 primes (which you can easily get off the internet) and using that to do the division test, you can take both of your points into account for primes up to 62 million (since the top number in the table is 7919 and that squared is 62 million and some odd). I just append the computed sequence onto the prime number list. I suppose that if I thought I needed to go further I could truncate this sequence at 62 million and use it as my "prime number list" for numbers up to 7919 ^ 4. Dunno how efficient that would be when you reach the second sequence but if you're looking to recurse over primes up to 62 million, you've got a lot of computing to do regardless.

Darrell Plank

By on 6/6/2008 7:51 AM ()

By the way, this approach of using the table of the first 1000 primes solves the EulerProject 10 of adding all primes less than 2 million in 3 seconds on my machine.

Code (primeList is defined earlier as an array of the first 1000 primes):

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

// Should produce primes up to 62,710,561
let primes =
   let DivisibleByListPrime (n:int) =
      let maxTest = int(sqrt(float(n)))
      let lastTest = (primeList |> Seq.find (fun pTest -> pTest > maxTest || n % pTest = 0))
      lastTest <= maxTest
   Seq.append
      primeList
      (
         Seq.init_infinite (fun i -> 2 * i + 7921) |> // Last prime in prime list is 7919
         Seq.filter (fun n -> not (DivisibleByListPrime n))
      )
By on 6/9/2008 11:39 PM ()

Unfold works just fine once you define nextPrime:

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

let divides m n =
  m % n = 0L;;
 
let factors n =
  [ for x in 1L..n when (divides n x) -> x ];;

let prime n =   
  let max = int64 (sqrt (float n))
  not ({2L..max} |> Seq.filter(fun d -> divides n d) |> Seq.nonempty);;

let rec nextPrime n = 
  if (prime n) then
    n
  else
    nextPrime (n + 1L);;
 
let primes = 
  Seq.unfold 
    (fun n ->
      let next = (nextPrime n) in
      Some (next, next + 1L))
     2L;;

I'm fairly new to FP, and very new to F# so this is doubtless inept, but it's more or less translated from SICP.

The primality test could, I think, use primes itself as its input. i.e. instead of
{2L..max}, use takeWhile (fun x -> x <= max) primes. (for some definition of takeWhile) At least it could if I could figure out forward declarations, just one of the 10,000 things needlessly baffling me in F#.

Great to have a functional language I can use at work though!

By on 4/14/2008 12:16 PM ()

Another possibility is to use a sequence expression.

1
2
3
4
5
6

 let primes = seq { let x = ref 2
                   while true do
                     if is_prime !x then 
                       yield !x
                     do x := !x + 1 }
By on 3/3/2008 6:52 AM ()

Thanks guys. I appreciate it.

By on 3/3/2008 12:02 PM ()

Hi,

I think unfold is not the best function for this (you need a loop to find the next prime number). I suggest you to use comprehensions.

1
2
3
4
5
6

let is_prime x =
  x > 1 &&
  not <| Seq.exists (fun n -> x % n = 0) {2..x-1}

let primes =
  {for i in 2 .. Int32.max_int when is_prime i -> i}

Another solution:

1
2

let N = Seq.init_infinite (fun x -> x)
let primes = N |> Seq.filter is_prime

Seq.unfold is useful when you can immediately compute the next value. For instance, here are Fibonacci numbers:

1

let fibos = Seq.unfold (fun (x, y) -> Some (x, (y, x + y))) (1, 1)

I hope these examples will help.

Laurent.

By on 3/2/2008 4:16 PM ()
1
2
3

let is_prime x =
  x > 1 &&
  not <| Seq.exists (fun n -> x % n = 0) {2..x-1}

Note that when testing for primes, you only have to test up to the square root of the number you're testing. Thus, the last line could be

1

  not <| Seq.exists (fun n -> x % n = 0) {2..int(sqrt (float x))}
By on 3/12/2008 7:23 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