Sieve of Pritchard

Sieve of Pritchard
You are encouraged to solve this task according to the task description, using any language you may know.

The Sieve of Pritchard is a modern algorithm for finding prime numbers. It takes many fewer operations than the Sieve of Eratosthenes (better time complexity), at the cost of greater storage requirements (worse space complexity).

Conceptually, it works by constructing a series of "wheels" marked along their circumference with the pattern of primes up to the value of successive primorial numbers (where the Nth primorial is the product of the first N primes). Those wheels are then rolled along the number line, and only the numbers touched by the marks are considered as candidate primes, in contrast to Eratosthenes' sieve in which all the integers in the range start out as candidates. (The Sieve of Pritchard is an example of the "wheel-based optimizations" mentioned in the Eratosthenes task.)

For example, the second-order wheel has size 6 (the product of the first two primes, 2 and 3) and is marked only at the numbers between 1 and 6 that are not multiples of 2 or 3, namely 1 and 5. As this wheel is rolled along the number line, it will pick up only numbers of the form 6k+1 or 6k+5 (that is, n where n mod 6 is in {1,5}). By the time it stops at 30 (2x3x5) it has added only 8 of the numbers between 6 and 30 as candidates for primality, only one of which is actually composite and must be removed (25). In the process it has constructed the next wheel, which will add only nine out of every 30 numbers as it rolls up to 210.

This YouTube video tells a story to help motivate the algorithm's design;this one presents the execution of the algorithm for N=150 in a format that permits single-stepping forward and backward through the run. In that implementation, the list of primes is populated into a sparse global array s such that s[p] contains the next prime after p iff p is itself a prime in the target range; this allows numbers to be removed from consideration quickly without any the copying/shifting that would be required from a normally-packed array.

Task

Write a program/subprogram that uses the Sieve of Pritchard algorithm to find all primes up to a specified limit. Show the result of running it with a limit of 150.

Related tasks

•   Sieve of Eratosthenes
•   Emirp primes
•   count in factors
•   prime decomposition
•   factors of an integer
•   extensible prime generator
•   primality by trial division
•   factors of a Mersenne number
•   trial factoring of a Mersenne number
•   partition an integer X into N primes
•   sequence of primes by Trial Division

C#

Loosely based on the Python version. I cut a couple of things out and it still worked. Not too crazy about having to create temporary lists to add or remove from the SortedSet, seems inefficient. But that is the work-around I employed, since SortedSets can't be accessed by indexing, and are non-mutable in a foreach loop. I haven't yet directly tested this against a Sieve of Eratosthenes to compare performance. The Wikipedia article suggests using a doubly linked list, so this C# incarnation is a kludge at best.

Compared to the prototype algorithm, it appears there isn't any code to do the follow-up end-of-wheel additions when necessary. But the main loop limit has been changed to go to the next prime, and the existing code handles the additions.

using System;
using System.Collections.Generic;

class Program {

    // Returns list of primes up to limit using Pritchard (wheel) sieve
    static List<int> PrimesUpTo(int limit) {
        var members = new SortedSet<int>{ 1 };
        int stp = 1, prime = 2, n, nxtpr, rtlim = 1 + (int)Math.Sqrt(limit), nl = 2;
        var primes = new List<int>();
        while (prime <= rtlim) {
            if (stp < limit) {
                var nu = new List<int>(); 
                foreach (var w in members)
                    for (n = w + stp; n <= nl; n += stp) nu.Add(n);
                members.UnionWith(nu);
            }
            stp = nl; // update wheel size to wheel limit
            nxtpr = 0; // for obtaining the next prime
            var wb = new List<int>();
            foreach (var w in members) {
                if (nxtpr == 0 && w > prime) nxtpr = w;
                if (members.Contains(n = prime * w)) wb.Add(n);
            }
            foreach (var itm in wb) members.Remove(itm);
            primes.Add(prime);
            prime = prime == 2 ? 3 : nxtpr;
            nl = Math.Min(limit, prime * stp);
        }
        members.Remove(1);
        primes.AddRange(members);
        return primes;
    }

    static void Main(string[] args) {
        Console.WriteLine("[{0}]", string.Join(", ", PrimesUpTo(150)));
    }
}

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149]

J

Implementation:

pritchard=: {{
  spokes=. $.6$4{.1
  primes=. 2, p=.3
  while. y > #spokes do.
    primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime
    rim=. #spokes NB. "length" of "circumference" of wheel
    spokes=. (y<.p*rim)$spokes NB. roll next larger wheel
    spokes=. 8 $.0 ((#~ y>])_1+p*1+i.rim)} spokes NB. remove newly recognized prime from wheel
  end.
  while. y > p*p do.
    primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime
    spokes=. 0 ((#~ y>])_1+p*1+i.rim)} spokes NB. scrub it out of wheel
  end.
  primes,1+}.,I.spokes
}}

Task example:

   pritchard 150
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149

Julia

Translation of: Raku

Added add/remove statistics. "Removed" figure is a combination of composites and primes under sqrt of limit.

Running on Julia 1.0

""" Rosetta Code task rosettacode.org/wiki/Sieve_of_Pritchard """

""" Pritchard sieve of primes up to limit, uses limit for type of the primes """

""" using Julia 1.0 (on Tio.run) """

function pritchard(limit::T) where T <: Integer
    members = Set(one(T))
    steplength = 1 # wheel size 
    prime = T(2)
    primes = T[]
    nlimit = prime * steplength # wheel limit
    ac = 2 # added count, since adding 1 & 2 during initialization
    rc = 1 # removed count, since 1 will be removed at the end
    rtlim = T(floor(sqrt(limit))) # this allows the main loop to go
    while prime <= rtlim # one extra time, eliminating the follow-up for
                         # the last partial wheel (if present)
        if steplength < limit
            for w in members
                n = w + steplength
                while n <= nlimit
                    push!(members, n)
                    ac += 1
                    n += steplength
                end
            end
            steplength = nlimit # advance wheel size
        end
        np = 5
        for w in sort!(collect(members))
            if np == 5 && w > prime np = w end
            n = prime * w
            if n > nlimit break end 
            rc += 1
            delete!(members, n)
        end
        if np < prime break end
        push!(primes, prime)
        # this works, but doubles the runtime for limit = 1000000
        # prime = prime == 2 ? 3 : sort!(collect(members))[2]
        # this is faster
        prime = prime == 2 ? 3 : np
        
        nlimit = min(steplength * prime, limit) # advance wheel limit 
    end
    delete!(members, 1)
    println("up to ",limit, ", added:", ac, " removed:", rc, " prime count:", length(primes) + length(members) )
    return sort!(append!(primes, members))
end

println(pritchard(150))

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149]

Phix

with javascript_semantics
function pritchard(integer limit)
    sequence members = {1}, primes = {}
    integer steplength = 1, prime = 2
    while prime * prime <= limit do
        if steplength < limit then
            integer mpsll = min(prime * steplength, limit)
            for w in members do
                integer n = w + steplength
                while n <= mpsll do
                    members &= n
                    n += steplength
                end while
            end for
            steplength = mpsll
        end if
        members = sort(filter(members,"out",sq_mul(members,prime)))
        primes &= prime
        prime = iff(prime=2?3:members[2])
    end while
    primes &= members[2..$]
    return primes
end function

printf(1,"%s\n",{join_by(pritchard(150),1,7," ",fmt:="%3d")})

Output:

  2   3   5   7  11  13  17
 19  23  29  31  37  41  43
 47  53  59  61  67  71  73
 79  83  89  97 101 103 107
109 113 127 131 137 139 149

Python

Translation of: Julia

""" Rosetta Code task rosettacode.org/wiki/Sieve_of_Pritchard """

def pritchard(limit):
    """ Pritchard sieve of primes up to limit """
    members = set([1])
    steplength, prime = 1, 2
    primes = []
    while prime * prime <= limit:
        if steplength < limit:
            nlimit = min(prime * steplength, limit)
            for w in list(members):
                n = w + steplength
                while n <= nlimit:
                    members.add(n)
                    n += steplength
            steplength = nlimit
        for w in sorted(members):
            n = prime * w
            if n > steplength: break # no use trying to remove items that can't even be there
            members.remove(n) # no checking necessary now
        primes.append(prime)
        prime = 3 if prime == 2 else min(m for m in members if m > 1)
    if steplength < limit:
        for w in sorted(members):
            n = w + steplength
            if n > limit: break # no use etc...
            while n <= limit:
                members.add(n)
                n += steplength 
    return primes + sorted(members)[1:]

print(pritchard(150))

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149]

Raku

First, a direct translation of the implementation in the YouTube video:

unit sub MAIN($limit = 150);

my $maxS = 1;
my $length = 2;
my $p = 3;
my @s = ();

while $p*$p <= $limit {
  if $length < $limit {
    extend-to [$p*$length, $limit].min;
  }
  delete-multiples-of($p);
  $p = next(1);
}
if $length < $limit {
  extend-to $limit;
}

# Done, build the list of actual primes from the array
$p = 3;
my @primes = 2, |gather while $p <= $limit {
  take $p;
  $p = next($p);
};
say @primes;

exit;

sub extend-to($n) {
  my $w = 1;
  my $x = $length + 1;
  while $x <= $n {
     append $x;
     $w = next($w);
     $x = $length + $w;
  }
  $length = $n;
  if $length == $limit {
    append $limit+2;
  }
}

sub delete-multiples-of($p) {
  my $f = $p;
  while $p*$f <= $length {
    $f = next($f);
  }
  while $f > 1 {
    $f = prev($f);
    delete($p*$f);
  }
}

sub append($w) {
  @s[$maxS-1] = $w;
  @s[$w-2] = $maxS;
  $maxS = $w;
}

sub next($w) { @s[$w-1]; }
sub prev($w) { @s[$w-2]; }

sub delete($pf) {
  my $temp1 = @s[$pf-2];
  my $temp2 = @s[$pf-1];
  @s[$temp1-1] = $temp2;
  @s[($temp2-2)%@s] = $temp1;
}

Output:

[2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149

Then a slightly more Raku-ish implementation based on the description in the Wikipedia article:

unit sub MAIN($limit = 150);

class Wheel {
  has $.members is rw;
  has $.length is rw;
  method extend(*@limits) {
    my @members = $.members.keys;
    for @members -> $w {
      my $n = $w + $.length;
      while $n <= @limits.all {
        $.members.set($n);
        $n += $.length;
      }
    }
    $.length = @limits.min;
  }
}

# start with W₀=({1},1)
my $wheel = Wheel.new: :members(SetHash(1)), :length(1);
my $prime = 2;
my @primes = ();

while $prime * $prime <= $limit {
  if $wheel.length < $limit {
    $wheel.extend($prime*$wheel.length, $limit);
  }
  for $wheel.members.keys.sort(-*) -> $w {
    $wheel.members.unset($prime * $w);
  }
  @primes.push: $prime;
  $prime = $prime == 2 ?? 3 !! $wheel.members.keys.grep(*>1).sort[0];
}

if $wheel.length < $limit {
  $wheel.extend($limit);
}
@primes.append:  $wheel.members.keys.grep: * != 1;
say @primes.sort;

The only difference in the output is that the result of `.sort` is a list rather than an array, so it's printed in parentheses instead of square brackets:

Output:

(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149)

Wren

Library: Wren-sort

Library: Wren-fmt

import "./sort" for SortedList
import "./fmt" for Fmt

var extend = Fn.new { |W, length, n|
    var w = 1
    var x = length + 1
    while (x <= n) {
        W.add(x)
        var ix = W.indexOf(w)
        w = W[ix+1]
        x = length + w
    }
}

var deleteMultiples = Fn.new { |W, length, p|
    var w = p
    while (p * w <= length) {
        var ix = W.indexOf(w)
        w = W[ix+1]
    }
    while (w > 1) {
        var ix = W.indexOf(w)
        w = W[ix-1]
        W.remove(p*w)
    }
}

var sieveOfPritchard = Fn.new { |N|
    if (N < 2) return []
    var W  = SortedList.fromOne(1)
    var Pr = SortedList.fromOne(2)
    var k = 1
    var length = 2
    var p = 3
    while (p * p <= N) {
        if (length < N) {
            var n = N.min(p*length)
            extend.call(W, length, n)
            length = n
        }
        deleteMultiples.call(W, length, p)
        Pr.add(p)
        k = k + 1
        p = W[1]
    }
    if (length < N) extend.call(W, length, N)
    return (Pr + W)[1..-1]
}

Fmt.tprint("$3d", sieveOfPritchard.call(150), 7)

Output:

  2   3   5   7  11  13  17 
 19  23  29  31  37  41  43 
 47  53  59  61  67  71  73 
 79  83  89  97 101 103 107 
109 113 127 131 137 139 149