Back to Questions
3.5k
views6
comments
Sieve of Atkin
f#

Hello all,

I have a question on the Sieve of Atkins. I have implemented the Sieve of Eratosthenes, and have confirmed that it is working without any problems. Now, I wanted to implement SoA to save on time when calculating primes.

My implementation is based on the pseudo-code found at [link:en.wikipedia.org]. I know I probably shouldn't trust wikipedia too much, but the code seems to make sense. However, there is a problem with generating primes in that the code below doesn't weed out all composite numbers.

If I change the multiples I am looking for (function mults to search for all multiples of a number vs. all multiples of its square, it eliminates all composites. However, from reading discussions around the web, the consensus seems to be that eliminating multiples of a prime numbers square *should* be enough. Obviously, I am missing a step. However, I'm not sure what that step is.

EDIT: 65, 85 and 91 are showing up when running soa 100I, and they shouldn't (as they are composites). I have worked it out by hand, and it makes sense for the function to let those numbers through as primes. So, my question is, am I misunderstanding the pseudocode or is the pseudocode incorrect?

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

// Sieve of Atkin **********DOESN'T WORK YET*********************
let soa limit =
  let sqlmt = BigInt.to_float limit |> sqrt |> ceil |> Float.to_string
              |> BigInt.of_string
  let cand = [ for x in 1I .. sqlmt
               for y in 1I .. sqlmt
               let a = 4I*x*x + y*y
               let b = 3I*x*x + y*y
               let c = 3I*x*x - y*y
               if a <= limit && (a % 12I = 1I || a % 12I = 5I) then yield a
               if b <= limit && b % 12I = 7I then yield b
               if x > y && c <= limit && c % 12I = 11I then yield c ]
             |> Set.of_list
  let mults p = p :: [ (p*p) .. (p*p) .. limit ] |> Set.of_list
  let prune y = Set.min_elt y, Set.min_elt y |> mults |> Set.diff y
  let rec nsoe' (p, pset) ans =
    match pset with
    | x when x = Set.empty || p > sqlmt -> Set.add p ans |> Set.union pset
    | _ -> nsoe' (prune pset) (Set.add p ans)
  nsoe' (prune cand) Set.empty
  |> Set.add 2I |> Set.add 3I |> Set.to_list

Any help is greatly appreciated.

Thanks and regards,

z.

.

The problem is in the definition of cand. Numbers don't get into the cand list by satisfying one of the if-conditions. Instead they get flipped in and out of the list each time they satisfy an if-condition. Quoting the pseudo-code comment:

<quote> <i>candidate primes [are] integers which have an odd number of representations by certain quadratic forms<i></quote> One way to fix this is <code lang=fsharp>cand |> List.sort compare |> group |> List.choose (fun xs -> if List.length xs % 2 = 1 then Some (List.hd xs) else None)</code> where <code lang=fsharp>let group xs = let rec span' p xs acc = match xs with | x::xs' when p x -> span' p xs' (x::acc) | _ -> List.rev acc, xs let rec group' xs acc = match xs with | x :: _ -> let ys, zs = span' (fun y -> x=y) xs [] group' zs (ys :: acc) | [] -> List.rev acc group' xs []</code>

By on 1/10/2008 10:43 PM ()

Acutally since you're using sets not lists anyway, it might be better to write it as

1
2
3

cand
|> List.map Set.singelton
|> List.fold_left xor Set.empty

where

1

let xor a b = Set.diff (Set.union a b) (Set.inter a b)
By on 1/10/2008 11:09 PM ()

Thanks for the very informative post. I had read the line that said "integers which have an odd number of representations by certain quadratic forms", but for some reason I just didn't connect it to the bit/bool-flipping implementation method.

After reading your explanation on how integers are put into / taken out of cand, I came to pretty much the same conclusion you did on using Set union/intersection.

EDIT: Well, after some very basic testing, it appears that working on eliminating duplicates through Sets is VERY slow (4 times slower than Sieve of Eratosthenes on primes < 10000). So, instead, I am using the function dup (below) to find duplicates.

1
2
3
4
5

  let rec dup ans lst=
    match ans, lst with
    | _, [] -> List.rev ans
    | [], h::t -> dup [ h ] t
    | a::b, h::t -> if a = h then dup b t else dup (h :: ans) t

Thanks again for your explanation,

z.

By on 1/10/2008 11:19 PM ()

Hi zakaluka,

If you can and do not mind, please post your implementations of SoE and SoA. I was going to write them, but you could spare me some time. :)

Kind regards,

Carl

By on 1/14/2008 1:40 AM ()

Sure, I've posted the code below. However, a few caveats:

  1. I have not extensively tested the code for correctness, save when using them for Project Euler problems.
  2. They are not optimized in the least. My focus was on readability, correctness, functional approach and KISS, not speed. I am sure you can get more performance by introducing mutability, especially with SoA.

Good luck,

z.

1
2
3
4
5
6
7
8

let soe limit =
  let mults p = p :: [ (p*p) .. p .. limit ] |> Set.of_list
  let prune y = Set.min_elt y, Set.min_elt y |> mults |> Set.diff y
  let rec nsoe' (p, pset) ans =
    match pset with
    | x when x = Set.empty -> Set.add p ans
    | _ -> nsoe' (prune pset) (Set.add p ans)
  nsoe' ([ 2I .. limit ] |> Set.of_list |> prune) Set.empty |> Set.to_list
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

let soa limit =
  let sqlmt = BigInt.to_float limit |> sqrt |> ceil |> Float.to_string
              |> BigInt.of_string
  let rec dup ans lst=
    match ans, lst with
    | _, [] -> List.rev ans
    | [], h::t -> dup [ h ] t
    | a::b, h::t -> if a = h then dup b t else dup (h::ans) t
  let cand = [ for x in 1I .. sqlmt
               for y in 1I .. sqlmt
               let a = 4I*x*x + y*y
               let b = 3I*x*x + y*y
               let c = 3I*x*x - y*y
               if a <= limit && (a % 12I = 1I || a % 12I = 5I) then yield a
               if b <= limit && b % 12I = 7I then yield b
               if x > y && c <= limit && c % 12I = 11I then yield c ]
             |> List.sort BigInt.compare |> dup [] |> Set.of_list
  let mults p = p :: [ p*p .. p*p .. limit ] |> Set.of_list
  let prune y = Set.min_elt y, Set.min_elt y |> mults |> Set.diff y
  let rec nsoe' (p, pset) ans =
    match pset with
    | x when x = Set.empty || p > sqlmt -> Set.add p ans |> Set.union pset
    | _ -> nsoe' (prune pset) (Set.add p ans)
  nsoe' (prune cand) Set.empty
  |> Set.add 2I |> Set.add 3I |> Set.to_list
By on 1/14/2008 11:39 AM ()

Thanks zakaluka,

I will let you know if I get it better optimized, but I am still learning F#.

By on 1/16/2008 3:02 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