Sieve of Pritchard

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

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)