Problem Statement

The Sieve of Pritchard is an algorithm for finding the prime numbers up to a given limit N, published in 1981. It considers many fewer composite numbers than the Sieve of Eratosthenes (and has a better asymptotic time complexity). However, unlike the latter, it cannot be modified to greatly reduce its space requirement, making it unsuitable for very large limits. Conceptually, it works by building a wheel representing the repeating pattern of numbers not divisible by one of the first k primes, increasing k until the square of the k'th prime is at least N. Since wheels grow very quickly, the algorithm restricts attention to the initial portions of wheels up to N. (Small examples of the wheels constructed by the Sieve of Pritchard are used in the "wheel-based optimizations" mentioned in the Eratosthenes task.) For example, the second-order wheel has circumference 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. Those that are multiples of 5 (only 2: 15 and 55) are obtained by multiplying the members of the second-order wheel. Removing them gives the next wheel, and so on. This YouTube video 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 wheel is represented by a sparse global array s such that for each member w of the wheel, s[w] contains the next member of the wheel; along with a similar "previous member" value, this allows numbers to be removed in a constant number of operations. But the simple abstract algorithm is based on an ordered set, and there is plenty of scope for different implementations. 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.

Let's start with the solution:

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;
}


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;