hubFS: THE place for 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 http://en.wikipedia.org/wiki/Sieve_of_Atkin.  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?

// 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:

candidate primes [are] integers which have an odd number of representations by certain quadratic forms

One way to fix this is

cand 
|> List.sort compare
|> group 
|> List.choose (fun xs -> if List.length xs % 2 = 1 then Some (List.hd xs) else None)

where

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 []

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

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

where

let xor a b = Set.diff (Set.union a b) (Set.inter a b)

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.

  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.

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

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.

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

Thanks zakaluka,

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