The sieve of Eratosthenes in F#

Question

I am interested in an implementation of the sieve of eratosthenes in purely functional F#. I am interested in an implementation of the actual sieve, not the naive functional implementation that isn't really the sieve, so not something like this:

let rec PseudoSieve list =     match list with     | hd::tl -> hd :: (PseudoSieve <| List.filter (fun x -> x % hd <> 0) tl)     | [] -> [] 

The second link above briefly describes an algorithm that would require the use of a multimap, which isn't available in F# as far as I know. The Haskell implementation given uses a map that supports an insertWith method, which I haven't seen available in the F# functional map.

Does anyone know a way to translate the given Haskell map code to F#, or perhaps knows of alternative implementation methods or sieving algorithms that are as efficient and better suited for a functional implementation or F#?

Answer 1

Reading that article I came up with an idea that doesn't require a multimap. It handles colliding map keys by moving the colliding key forward by its prime value again and again until it reaches a key that isn't in the map. Below primes is a map with keys of the next iterator value and values that are primes.

let primes =      let rec nextPrime n p primes =         if primes |> Map.containsKey n then             nextPrime (n + p) p primes         else             primes.Add(n, p)      let rec prime n primes =         seq {             if primes |> Map.containsKey n then                 let p = primes.Item n                 yield! prime (n + 1) (nextPrime (n + p) p (primes.Remove n))             else                 yield n                 yield! prime (n + 1) (primes.Add(n * n, n))         }      prime 2 Map.empty 

Here's the priority queue based algorithm from that paper without the square optimization. I placed the generic priority queue functions at the top. I used a tuple to represent the lazy list iterators.

let primes() =      // the priority queue functions     let insert = Heap.Insert     let findMin = Heap.Min     let insertDeleteMin = Heap.DeleteInsert      // skips primes 2, 3, 5, 7     let wheelData = [|2L;4L;2L;4L;6L;2L;6L;4L;2L;4L;6L;6L;2L;6L;4L;2L;6L;4L;6L;8L;4L;2L;4L;2L;4L;8L;6L;4L;6L;2L;4L;6L;2L;6L;6L;4L;2L;4L;6L;2L;6L;4L;2L;4L;2L;10L;2L;10L|]      // increments iterator     let wheel (composite, n, prime) =         composite + wheelData.[n % 48] * prime, n + 1, prime      let insertPrime prime n table =         insert (prime * prime, n, prime) table      let rec adjust x (table : Heap) =         let composite, n, prime = findMin table          if composite <= x then              table              |> insertDeleteMin (wheel (composite, n, prime))             |> adjust x         else             table      let rec sieve iterator table =         seq {             let x, n, _ = iterator             let composite, _, _ = findMin table              if composite <= x then                 yield! sieve (wheel iterator) (adjust x table)             else                 if x = 13L then                     yield! [2L; 3L; 5L; 7L; 11L]                  yield x                 yield! sieve (wheel iterator) (insertPrime x n table)         }      sieve (13L, 1, 1L) (insertPrime 11L 0 (Heap(0L, 0, 0L))) 

Here's the priority queue based algorithm with the square optimization. In order to facilitate lazy adding primes to the lookup table, the wheel offsets had to be returned along with prime values. This version of the algorithm has O(sqrt(n)) memory usage where the none optimized one is O(n).

let rec primes2() : seq<int64 * int> =      // the priority queue functions     let insert = Heap.Insert     let findMin = Heap.Min     let insertDeleteMin = Heap.DeleteInsert      // increments iterator     let wheel (composite, n, prime) =         composite + wheelData.[n % 48] * prime, n + 1, prime      let insertPrime enumerator composite table =         // lazy initialize the enumerator         let enumerator =             if enumerator = null then                 let enumerator = primes2().GetEnumerator()                 enumerator.MoveNext() |> ignore                 // skip primes that are a part of the wheel                 while fst enumerator.Current < 11L do                     enumerator.MoveNext() |> ignore                 enumerator             else                 enumerator          let prime = fst enumerator.Current         // Wait to insert primes until their square is less than the tables current min         if prime * prime < composite then             enumerator.MoveNext() |> ignore             let prime, n = enumerator.Current             enumerator, insert (prime * prime, n, prime) table         else             enumerator, table      let rec adjust x table =         let composite, n, prime = findMin table          if composite <= x then              table              |> insertDeleteMin (wheel (composite, n, prime))             |> adjust x         else             table      let rec sieve iterator (enumerator, table) =          seq {             let x, n, _ = iterator             let composite, _, _ = findMin table              if composite <= x then                 yield! sieve (wheel iterator) (enumerator, adjust x table)             else                 if x = 13L then                     yield! [2L, 0; 3L, 0; 5L, 0; 7L, 0; 11L, 0]                  yield x, n                 yield! sieve (wheel iterator) (insertPrime enumerator composite table)         }      sieve (13L, 1, 1L) (null, insert (11L * 11L, 0, 11L) (Heap(0L, 0, 0L))) 

Here's my test program.

type GenericHeap<'T when 'T : comparison>(defaultValue : 'T) =     let mutable capacity = 1     let mutable values = Array.create capacity defaultValue     let mutable size = 0      let swap i n =         let temp = values.[i]         values.[i] <- values.[n]         values.[n] <- temp      let rec rollUp i =         if i > 0 then             let parent = (i - 1) / 2             if values.[i] < values.[parent] then                 swap i parent                 rollUp parent      let rec rollDown i =         let left, right = 2 * i + 1, 2 * i + 2          if right < size then             if values.[left] < values.[i] then                 if values.[left] < values.[right] then                     swap left i                     rollDown left                 else                     swap right i                     rollDown right             elif values.[right] < values.[i] then                 swap right i                 rollDown right         elif left < size then             if values.[left] < values.[i] then                 swap left i      member this.insert (value : 'T) =         if size = capacity then             capacity <- capacity * 2             let newValues = Array.zeroCreate capacity             for i in 0 .. size - 1 do                 newValues.[i] <- values.[i]             values <- newValues          values.[size] <- value         size <- size + 1         rollUp (size - 1)      member this.delete () =         values.[0] <- values.[size]         size <- size - 1         rollDown 0      member this.deleteInsert (value : 'T) =         values.[0] <- value         rollDown 0      member this.min () =         values.[0]      static member Insert (value : 'T) (heap : GenericHeap<'T>) =         heap.insert value         heap          static member DeleteInsert (value : 'T) (heap : GenericHeap<'T>) =         heap.deleteInsert value         heap          static member Min (heap : GenericHeap<'T>) =         heap.min()  type Heap = GenericHeap<int64 * int * int64>  let wheelData = [|2L;4L;2L;4L;6L;2L;6L;4L;2L;4L;6L;6L;2L;6L;4L;2L;6L;4L;6L;8L;4L;2L;4L;2L;4L;8L;6L;4L;6L;2L;4L;6L;2L;6L;6L;4L;2L;4L;6L;2L;6L;4L;2L;4L;2L;10L;2L;10L|]  let primes() =      // the priority queue functions     let insert = Heap.Insert     let findMin = Heap.Min     let insertDeleteMin = Heap.DeleteInsert      // increments iterator     let wheel (composite, n, prime) =         composite + wheelData.[n % 48] * prime, n + 1, prime      let insertPrime prime n table =         insert (prime * prime, n, prime) table      let rec adjust x (table : Heap) =         let composite, n, prime = findMin table          if composite <= x then              table              |> insertDeleteMin (wheel (composite, n, prime))             |> adjust x         else             table      let rec sieve iterator table =         seq {             let x, n, _ = iterator             let composite, _, _ = findMin table              if composite <= x then                 yield! sieve (wheel iterator) (adjust x table)             else                 if x = 13L then                     yield! [2L; 3L; 5L; 7L; 11L]                  yield x                 yield! sieve (wheel iterator) (insertPrime x n table)         }      sieve (13L, 1, 1L) (insertPrime 11L 0 (Heap(0L, 0, 0L)))  let rec primes2() : seq<int64 * int> =      // the priority queue functions     let insert = Heap.Insert     let findMin = Heap.Min     let insertDeleteMin = Heap.DeleteInsert      // increments iterator     let wheel (composite, n, prime) =         composite + wheelData.[n % 48] * prime, n + 1, prime      let insertPrime enumerator composite table =         // lazy initialize the enumerator         let enumerator =             if enumerator = null then                 let enumerator = primes2().GetEnumerator()                 enumerator.MoveNext() |> ignore                 // skip primes that are a part of the wheel                 while fst enumerator.Current < 11L do                     enumerator.MoveNext() |> ignore                 enumerator             else                 enumerator          let prime = fst enumerator.Current         // Wait to insert primes until their square is less than the tables current min         if prime * prime < composite then             enumerator.MoveNext() |> ignore             let prime, n = enumerator.Current             enumerator, insert (prime * prime, n, prime) table         else             enumerator, table      let rec adjust x table =         let composite, n, prime = findMin table          if composite <= x then              table              |> insertDeleteMin (wheel (composite, n, prime))             |> adjust x         else             table      let rec sieve iterator (enumerator, table) =          seq {             let x, n, _ = iterator             let composite, _, _ = findMin table              if composite <= x then                 yield! sieve (wheel iterator) (enumerator, adjust x table)             else                 if x = 13L then                     yield! [2L, 0; 3L, 0; 5L, 0; 7L, 0; 11L, 0]                  yield x, n                 yield! sieve (wheel iterator) (insertPrime enumerator composite table)         }      sieve (13L, 1, 1L) (null, insert (11L * 11L, 0, 11L) (Heap(0L, 0, 0L)))   let mutable i = 0  let compare a b =     i <- i + 1     if a = b then         true     else         printfn "%A %A %A" a b i         false  Seq.forall2 compare (Seq.take 50000 (primes())) (Seq.take 50000 (primes2() |> Seq.map fst)) |> printfn "%A"  primes2() |> Seq.map fst |> Seq.take 10 |> Seq.toArray |> printfn "%A"  primes2() |> Seq.map fst |> Seq.skip 999999 |> Seq.take 10 |> Seq.toArray |> printfn "%A"  System.Console.ReadLine() |> ignore