Jump to content

Sieve of Pritchard

From Rosetta Code
Task
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 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: 1*5 and 5*5) 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.

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


AppleScript

on sieveOfPritchard(limit)
    if (limit < 2) then return {}
    script o
        property wheel : {2}
        property extension : missing value
    end script
    
    set {x, newCircumference, currentPrime, mv} to {0, 2, 1, missing value}
    repeat until (currentPrime * currentPrime > limit)
        -- Get the next confirmed prime from the wheel.
        set x to x + 1
        set currentPrime to o's wheel's item x
        -- Get an extension list nominally expanding the wheel to the lesser of
        -- its current circumference * currentPrime and the limit parameter.
        -- It'll be far longer than needed, but hey.
        set oldCircumference to newCircumference
        set newCircumference to oldCircumference * currentPrime
        if (newCircumference > limit) then set newCircumference to limit
        set o's extension to makeList(newCircumference - oldCircumference, mv)
        -- Insert numbers that are oldCircumference added to 1 and to each number currently in the
        -- unpartitioned part of the wheel, except where the results are multiples of currentPrime.
        set k to 0
        set listLen to (count o's wheel)
        repeat with augend from oldCircumference to (newCircumference - 1) by oldCircumference
            set n to augend + 1
            if (n mod currentPrime > 0) then
                set k to k + 1
                set o's extension's item k to n
            end if
            repeat with i from x to listLen
                set n to augend + (o's wheel's item i)
                if (n > newCircumference) then exit repeat
                if (n mod currentPrime > 0) then
                    set k to k + 1
                    set o's extension's item k to n
                end if
            end repeat
        end repeat
        -- Find and delete multiples of the current prime which occur in the old part of the wheel.
        set maxMultiple to oldCircumference div currentPrime
        set i to x
        repeat while ((o's wheel's item i) < maxMultiple)
            set i to i + 1
        end repeat
        repeat with i from i to x by -1
            set j to binarySearch((o's wheel's item i) * currentPrime, o's wheel, i, listLen)
            if (j > 0) then
                set o's wheel's item j to mv
                set listLen to j - 1
            end if
        end repeat
        -- Keep the undeleted numbers and any in the extension list.
        set o's wheel to o's wheel's numbers
        if (k > 0) then set o's wheel to o's wheel & o's extension's items 1 thru k
    end repeat
    
    return o's wheel
end sieveOfPritchard

on makeList(limit, filler)
    if (limit < 1) then return {}
    script o
        property lst : {filler}
    end script
    
    set counter to 1
    repeat until (counter + counter > limit)
        set o's lst to o's lst & o's lst
        set counter to counter + counter
    end repeat
    if (counter < limit) then set o's lst to o's lst & o's lst's items 1 thru (limit - counter)
    return o's lst
end makeList

on binarySearch(n, theList, l, r)
    script o
        property lst : theList
    end script
    repeat until (l = r)
        set m to (l + r) div 2
        if (o's lst's item m < n) then
            set l to m + 1
        else
            set r to m
        end if
    end repeat
    
    if (o's lst's item l is n) then return l
    return 0
end binarySearch

sieveOfPritchard(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}

BASIC

BASIC256

arraybase 1
call sieveOfPritchard(150, true)
call sieveOfPritchard(1e6, false)
end

function min(a, b)
  if a < b then return a else return b
end function

subroutine sieveOfPritchard(limit, imprime)
  dim members[limit + 1] fill false
  members[1] = true
  ub = members[?]
  stepLength = 1
  prime = 2
  rtlim = sqr(limit)
  nlimit = 2
  dim primes[1]
  cont = 0

  while prime <= rtlim
    if stepLength < limit then
      for w = 1 to ub
        if members[w] then
          dim n = w + stepLength
          while n <= nlimit
            members[n] = true
            n += stepLength
          end while
        end if
      next
      stepLength = nlimit
    end if

    np = 5
    dim mcpy[ub]
    for i = 1 to ub
      mcpy[i] = members[i]
    next

    for i = 1 to ub
      if mcpy[i] then
        if np = 5 and i > prime then np = i
        dim n = prime * i
        if n > limit then exit for
        members[n] = false
      end if
    next

    if np < prime then exit while
    cont += 1
    redim primes(cont)
    primes[cont] = prime
    if prime = 2 then prime = 3 else prime = np
    nlimit = min(stepLength * prime, limit)
  end while

  dim newPrimes(ub)
  for i = 2 to ub
    if members[i] then newPrimes[i] = i
  next

  cont = 0
  for i = 1 to primes[?]
    if imprime then print " "; primes[i];
    cont += 1
  next
  for i = 1 to ub
    if newPrimes[i] then
      cont += 1
      if imprime then print " "; i;
    end if
  next
  if not imprime then print : print "Number of primes up to "; limit; ": "; cont
end subroutine
Output:
Similar to FreeBASIC entry.

FreeBASIC

Translation of: Wren
#define min(a, b) iif((a) < (b), (a), (b))

Sub sieveOfPritchard(limit As Uinteger, imprime As Boolean)
    Dim As Boolean members(1 To limit + 1)
    members(1) = True
    Dim As Uinteger ub = Ubound(members)
    Dim As Uinteger stepLength = 1
    Dim As Uinteger prime = 2
    Dim As Uinteger rtlim = Sqr(limit)
    Dim As Uinteger nlimit = 2
    Dim As Integer primes()
    Dim As Integer i, cont = 0
    
    While prime <= rtlim
        If stepLength < limit Then
            For w As Integer = 1 To ub
                If members(w) Then
                    Dim As Integer n = w + stepLength
                    While n <= nlimit
                        members(n) = True
                        n += stepLength
                    Wend
                End If
            Next
            stepLength = nlimit
        End If
        
        Dim As Uinteger np = 5
        Dim As Boolean mcpy(ub)
        For i = 1 To ub
            mcpy(i) = members(i)
        Next
        
        For i = 1 To ub
            If mcpy(i) Then
                If np = 5 And i > prime Then np = i
                Dim As Uinteger n = prime * i
                If n > limit Then Exit For there.
                members(n) = False
            End If
        Next
        
        If np < prime Then Exit While
        cont += 1
        Redim Preserve primes(cont)
        primes(cont) = prime
        prime = Iif(prime = 2, 3, np)
        nlimit = min(stepLength * prime, limit)
    Wend
    
    Dim As Integer newPrimes(ub)
    For i = 2 To ub
        If members(i) Then newPrimes(i) = i
    Next
    
    cont = 0
    For i = 1 To Ubound(primes)
        If imprime Then Print primes(i); 
        cont += 1
    Next
    For i = 1 To ub
        If newPrimes(i) Then 
            cont += 1
            If imprime Then Print i; 
        End If
    Next
    If Not imprime Then Print !"\nNumber of primes up to "; limit; ":"; cont
End Sub

sieveOfPritchard(150, True)
sieveOfPritchard(1e6, False)

Sleep
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
Number of primes up to 1000000: 78498

True BASIC

FUNCTION min(a, b)
    IF a < b THEN LET min = a ELSE LET min = b
END FUNCTION

SUB sieveofpritchard (limit)
    DIM members(0)
    MAT REDIM members(limit+1)
    LET members(1) = 1
    LET ub = UBOUND(members)
    LET steplength = 1
    LET prime = 2
    LET rtlim = SQR(limit)
    LET nlimit = 2
    DIM primes(10)
    LET cnt = 0
    DIM mcpy(0)
    MAT REDIM mcpy(ub)

    DO WHILE prime <= rtlim
       IF steplength < limit THEN
          FOR w = 1 TO ub
              IF members(w)<>0 THEN
                 LET n = w+steplength
                 DO WHILE n <= nlimit
                    LET members(n) = 1
                    LET n = n+steplength
                 LOOP
              END IF
          NEXT w
          LET steplength = nlimit
       END IF

       LET np = 5
       FOR i = 1 TO ub
           LET mcpy(i) = members(i)
       NEXT i

       FOR i = 1 TO ub
           IF mcpy(i)<>0 THEN
              IF np = 5 AND i > prime THEN LET np = i
              LET n = prime*i
              IF n > limit THEN EXIT FOR
              LET members(n) = 0
           END IF
       NEXT i

       IF np < prime THEN EXIT DO
       LET cnt = cnt+1
       MAT REDIM primes(cnt)
       LET primes(cnt) = prime
       IF prime = 2 THEN LET prime = 3 ELSE LET prime = np
       LET nlimit = min(steplength*prime, limit)
    LOOP

    DIM newprimes(0)
    MAT REDIM newprimes(ub)
    FOR i = 2 TO ub
        IF members(i)<>0 THEN LET newprimes(i) = i
    NEXT i

    LET cnt = 0
    FOR i = 1 TO UBOUND(primes)
        PRINT primes(i);
        LET cnt = cnt+1
    NEXT i
    FOR i = 1 TO ub
        IF newprimes(i)<>0 THEN
           PRINT i;
           LET cnt = cnt+1
        END IF
    NEXT i
END SUB

CLEAR
CALL sieveofpritchard (150)
END
Output:
Similar to FreeBASIC entry.

Yabasic

sieveOfPritchard(150, true)
sieveOfPritchard(1e6, false)
end

sub sieveOfPritchard(limit, imprime)
  dim members(limit + 1)
  members(1) = true
  ub = arraysize(members(),1)
  stepLength = 1
  prime = 2
  rtlim = sqrt(limit)
  nlimit = 2
  dim primes(1)
  cont = 0

  while prime <= rtlim
    if stepLength < limit then
      for w = 1 to ub
        if members(w) then
          n = w + stepLength
          while n <= nlimit
            members(n) = true
            n = n + stepLength
          wend
        fi
      next
      stepLength = nlimit
    fi

    np = 5
    dim mcpy(ub)
    for i = 1 to ub
      mcpy(i) = members(i)
    next

    for i = 1 to ub
      if mcpy(i) then
        if np = 5 and i > prime  np = i
        n = prime * i
        if n > limit  break
        members(n) = false
      fi
    next

    if np < prime  break
    cont = cont + 1
    redim primes(cont)
    primes(cont) = prime
    if prime = 2 then prime = 3 else prime = np : fi
    nlimit = min(stepLength * prime, limit)
  wend

  dim newPrimes(ub)
  for i = 2 to ub
    if members(i)  newPrimes(i) = i
  next

  cont = 0
  for i = 1 to arraysize(primes(),1)
    if imprime  print " ", primes(i);
    cont = cont + 1
  next
  for i = 1 to ub
    if newPrimes(i) then
      cont = cont + 1
      if imprime  print " ", i;
    fi
  next
  if not imprime then print : print "Number of primes up to ", limit, ": ", cont : fi
end sub
Output:
Similar to FreeBASIC entry.


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.

Updated to include "numbers added / removed (to / from members)" and performance statistics. The "removed" figure includes both composite numbers and prime numbers less than the square root of limit. The Wikipedia article indicates only eight removals (for limit = 150) because it doesn't count the removed primes and the initial 1 that the members array is initialized with.

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, bool verbose = false) {
        var sw = System.Diagnostics.Stopwatch.StartNew();
        var members = new SortedSet<int>{ 1 };
        int stp = 1, prime = 2, n, nxtpr, rtlim = 1 + (int)Math.Sqrt(limit), nl, ac = 2, rc = 1;
        List<int> primes = new List<int>(), tl = new List<int>();
        while (prime < rtlim) {
            nl = Math.Min(prime * stp, limit);
            if (stp < limit) {
                tl.Clear(); 
                foreach (var w in members)
                    for (n = w + stp; n <= nl; n += stp) tl.Add(n);
                members.UnionWith(tl); ac += tl.Count;
            }
            stp = nl; // update wheel size to wheel limit
            nxtpr = 5; // for obtaining the next prime
            tl.Clear();
            foreach (var w in members) {
                if (nxtpr == 5 && w > prime) nxtpr = w;
                if ((n = prime * w) > nl) break; else tl.Add(n);
            }
            foreach (var itm in tl) members.Remove(itm); rc += tl.Count;
            primes.Add(prime);
            prime = prime == 2 ? 3 : nxtpr;
        }
        members.Remove(1); primes.AddRange(members); sw.Stop();
        if (verbose) Console.WriteLine("Up to {0}, added:{1}, removed:{2}, primes counted:{3}, time:{4} ms", limit, ac, rc, primes.Count, sw.Elapsed.TotalMilliseconds);
        return primes;
    }

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

Timing from Tio.run

Up to 150, added:45, removed:14, primes counted:35, time:13.2842 ms
[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]
Up to 1000000, added:186825, removed:108494, primes counted:78498, time:139.4323 ms

C++

We obtain a simple but high-performance implementation.

The starting idea is to represent W as simply as possible, with an array w[] containing the members in order; i.e. w[i] is the i'th member (indexing from 0). When the current wheel is extended by rolling it, the code simply iterates through the array w, adding the length of the wheel to each member w[i] and appending the result. The other step is to delete the composites formed by multiplying the values in the current wheel by the current prime p. However, this presents problems, firstly because each multiple cannot be found in w in O(1) time. Accordingly, a bit array d[] (for "deleted") is introduced such that d[x] is initialized to false when a value x is appended to W, and is set to true should x be deleted as a multiple p*w[i] of p. Deletions are now fast, but the array is left containing deleted elements.

So if the new W will be extended in the next iteration, because its length < N, then the array w is compressed by eliminating the deleted values. But once the length reaches N (which happens very quickly), it would be way too costly to compress w at the end of each iteration. However, only the values in W up to N/p will be used as factors in the next lot of deletions. So it suffices to compress only this initial section of w. When the remaining primes are gathered on completion, it is necessary to skip zero and deleted values in w.

Each low-level operation in the resulting algorithm can be associated with an abstract operation so that each of the latter gets O(1) operations. So the resulting program still runs in O(N/log log N) time, and the implicit constant factor is quite small.

/* Sieve of Pritchard in C++, as described at https://en.wikipedia.org/wiki/Sieve_of_Pritchard */
/* Simple but high-performance implementation using a simple array of integers and a bit array (using bytes for speed). */
/* 2 <= N <= 1000000000 */
/* (like the standard Sieve of Eratosthenes, this algorithm is not suitable for very large N due to memory requirements) */

#include <cstring>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <ctime>

void Extend (uint32_t w[], uint32_t &w_end, uint32_t &length, uint32_t n, bool d[], uint32_t &w_end_max) {
    /* Rolls full wheel W up to n, and sets length=n */
    uint32_t i, j, x;
    i = 0; j = w_end;
    x = length + 1; /* length+w[0] */
    while (x <= n) {
        w[++j] = x; /* Append x to the ordered set W */
        d[x] = false;
        x = length + w[++i];
    }
    length = n; w_end = j;
    if (w_end > w_end_max) w_end_max = w_end;
}

void Delete (uint32_t w[], uint32_t length, uint32_t p, bool d[], uint32_t &imaxf) {
    /* Deletes multiples p*w[i] of p from W, and sets imaxf to last i for deletion */
    uint32_t i, x;
    i = 0;
    x = p; /* p*w[0]=p*1 */
    while (x <= length) {
        d[x] = true; /* Remove x from W; */
        x = p*w[++i];
    }
    imaxf = i-1;
}

void Compress(uint32_t w[], bool d[], uint32_t to, uint32_t &w_end) {
    /* Removes deleted values in w[0..to], and if to=w_end, updates w_end, otherwise pads with zeros on right */
    uint32_t i, j;
    j = 0;
    for (i=1; i <= to; i++) {
        if (!d[w[i]]) {
            w[++j] = w[i];
        }
    }
    if (to == w_end) {
        w_end = j;
    } else {
        for (uint32_t k=j+1; k <= to; k++) w[k] = 0;
    }
}

void Sift(uint32_t N, bool printPrimes, uint32_t &nrPrimes, uint32_t &vBound) {
    /* finds the nrPrimes primes up to N, printing them if printPrimes */
    uint32_t *w = new uint32_t[N/4+5];
    bool *d = new bool[N+1];
    uint32_t w_end, length;
    /* representation invariant (for the main loop): */
    /* if length < N (so W is a complete wheel), w[0..w_end] is the ordered set W; */
    /* otherwise, w[0..w_end], omitting zeros and values w with d[w] true, is the ordered set W, */
    /* and no values <= N/p are omitted */
    uint32_t w_end_max, p, imaxf;
    /* W,k,length = {1},1,2: */
    w_end = 0; w[0] = 1;
    w_end_max = 0;
    length = 2;
    /* Pr = {2}: */
    nrPrimes = 1;
    if (printPrimes) printf("%d", 2);
    p = 3;
    /* invariant: p = p_(k+1) and W = W_k inter {1,...,N} and length = min(P_k,N) and Pr = the first k primes */
    /* (where p_i denotes the i'th prime, W_i denotes the i'th wheel, P_i denotes the product of the first i primes) */
    while (p*p <= N) {
        /* Append p to Pr: */
        nrPrimes++;
        if (printPrimes) printf(" %d", p);
        if (length < N) {
            /* Extend W with length to minimum of p*length and N: */
            Extend (w, w_end, length, std::min(p*length,N), d, w_end_max);
        }
        Delete(w, length, p, d, imaxf);
        Compress(w, d, (length < N ? w_end : imaxf), w_end);
        /* p = next(W, 1): */
        p = w[1];
        if (p == 0) break; /* next p is after zeroed section so is too big */
        /* k++ */
    }
    if (length < N) {
        /* Extend full wheel W,length to N: */
        Extend (w, w_end, length, N, d, w_end_max);
    }
    /* gather remaining primes: */
    for (uint32_t i=1; i <= w_end; i++) {
        if (w[i] == 0 || d[w[i]]) continue;
        if (printPrimes) printf(" %d", w[i]);
        nrPrimes++;
    }
    vBound = w_end_max+1;
}

int main (int argc, char *argw[]) {
    bool error = false;
    bool printPrimes = false;
    uint32_t N, nrPrimes, vBound;
    if (argc == 3) {
        if (strcmp(argw[2], "-p") == 0) {
            printPrimes = true;
            argc--;
        } else {
            error = true;
        }
    }
    if (argc == 2) {
        N = atoi(argw[1]);
        if (N < 2 || N > 1000000000) error = true;
    } else {
        error = true;
    }
    if (error) {
        printf("call with: %s N -p where 2 <= N <= 1000000000 and -p to print the primes is optional \n", argw[0]);
        exit(1);
    }
    int start_s = clock();
    Sift(N, printPrimes, nrPrimes, vBound);
    int stop_s=clock();
    printf("\n%d primes up to %lu found in %.3f ms using array w[%d]\n", nrPrimes,
      (unsigned long)N, (stop_s-start_s)*1E3/double(CLOCKS_PER_SEC), vBound);
}
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
35 primes up to 150 found in 0.038 ms using array w[40]

50847534 primes up to 1000000000 found in 2339.566 ms using array w[163588196]

EasyLang

Translation of: Julia
proc pritchard limit . primes[] .
   len members[] limit
   members[1] = 1
   steplength = 1
   prime = 2
   rtlimit = sqrt limit
   nlimit = 2
   primes[] = [ ]
   while prime <= rtlimit
      if steplength < limit
         for w to len members[]
            if members[w] = 1
               n = w + steplength
               while n <= nlimit
                  members[n] = 1
                  n += steplength
               .
            .
         .
         steplength = nlimit
      .
      np = 5
      mcpy[] = members[]
      for w to nlimit
         if mcpy[w] = 1
            if np = 5 and w > prime
               np = w
            .
            n = prime * w
            if n > nlimit
               break 1
            .
            members[n] = 0
         .
      .
      if np < prime
         break 1
      .
      primes[] &= prime
      if prime = 2
         prime = 3
      else
         prime = np
      .
      nlimit = lower (steplength * prime) limit
   .
   for i = 2 to len members[]
      if members[i] = 1
         primes[] &= i
      .
   .
.
pritchard 150 p[]
print p[]
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 model:

pritchard=: {{N=. y
  root=. >.@%: N
  spokes=. 1
  primes=. ''
  p=. 0
  while. p<:root do.
    primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime
    rim=. #spokes NB. "length" of "circumference" of wheel
    spokes=. (N<.p*rim)$spokes NB. roll next larger wheel
    NB. remove multiples of this next prime:
    spokes=. 0 ((#spokes) (>#]) _1+p*1+i.rim)} spokes NB. remove newly recognized prime from wheel
  end.
  N (>:#]) 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

However, this approach exhibits performance problems when N is large.

A faster approach recognizes when the wheel is large enough and treats all subsequent "next primes" specially:

pr =: {{N=.y
  root=. <.%:N  NB. performance optimzation
  circumference=. 1
  spokes=. ,1
  primes=. ''
  while. N > L=. circumference do.
    primes=. primes, p =. 1{ spokes,L+1  NB. next prime from sieve
    circumference=. N <. p * L           NB. next larger wheel:
    spokes=. circumference (>:#]), spokes +/~ L * i.circumference >.@% L
    NB. remove multiples of this next prime:
    spokes=. spokes -. p * spokes ( [{.~ >:@:(I.-(-.@e.)~))circumference<.@%p
  end.
  NB. set up for optimized version of above code
  comb=. root (>:#]) }. spokes NB. candidate next primes to consider
  discardp=. discard=. '' NB. what we'll be eliminating
  for_p. comb do.
    if. p e. comb =. comb (-. }.) discardp do.
      NB. remove multiples of this next prime:
      discardp=. p * spokes ( [{.~ >:@:(I.-(-.@e.)~))circumference<.@%p
      discard =. discard, discardp
    end.
  end.
  primes,comb,}.spokes-.discard
}}

Here, pr 150 gives the same result as pritchard 150 but pr 1e7 takes well under a second.

Java

import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;

public final class SieveOfPritchard {
 
	public static void main(String[] args) {		
		System.out.println(sieveOfPritchard(150) + System.lineSeparator());
		System.out.println("Number of primes up to 1,000,000 is " + sieveOfPritchard(1_000_000).size() + ".");
		System.out.println();
		
		final long start = System.currentTimeMillis();		
		System.out.print("Number of primes up to 100,000,000 is " + sieveOfPritchard(100_000_000).size());
		final long finish = System.currentTimeMillis();
		System.out.println(". Obtained in a time of " + ( (double) finish - start ) / 1_000 + " seconds.");
	}
	
	private static List<Integer> sieveOfPritchard(int limit) {		
		List<Integer> primes = new ArrayList<Integer>();
		BitSet members = new BitSet(limit + 1);
	    members.set(1);
	    List<Integer> deletions = new ArrayList<Integer>();
	    final int rootLimit = (int) Math.sqrt(limit);
	    int nLimit = 2;
	    int stepLength = 1;
	    int prime = 2;	  
	    	    		
	    while ( prime <= rootLimit ) {
	    	if ( stepLength < limit ) { 
	    		for ( int w = 1; w >= 0; w = members.nextSetBit(w + 1) ) {
                    int n = w + stepLength;	                    
                    while ( n <= nLimit ) {
                        members.set(n);
                        n += stepLength;
		            }
	    		}
	            stepLength = nLimit;
	    	} 
	    	
	    	deletions.clear();
	    	int nextPrime = 5;
	        for ( int w = 1; w < nLimit; w = members.nextSetBit(w + 1) ) {
                if ( nextPrime == 5 && w > prime ) {
                    nextPrime = w;
                }
                final int n = prime * w;
                if ( n > nLimit ) {
                	break;
                }
                deletions.add(n);
	        }
	        for ( int deletion : deletions ) {
	        	members.clear(deletion);
	        }
	        
	        if ( nextPrime < prime ) {
	            break;
	        }
	        
	        primes.add(prime);
	        prime = ( prime == 2 ) ? 3 : nextPrime;	        
	        nLimit = (int) Math.min((long) stepLength * prime, limit);	            
	    }
	    
	    if ( stepLength < limit ) {
	    	for ( int w = 1; w >= 0; w = members.nextSetBit(w + 1) ) {
                int n = w + stepLength;	                    
                while ( n <= limit ) {
                    members.set(n);
                    n += stepLength;
	            }
    		}
	    }
	    
	    members.clear(1);
	    for ( int i = members.nextSetBit(0); i >= 0; i = members.nextSetBit(i + 1) ) {
	    	primes.add(i);
	    }; 
	    
	    return primes;	
	}

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

Number of primes up to 1,000,000 is 78498.

Number of primes up to 100,000,000 is 5761455. Obtained in a time of 0.648 seconds.

Julia

Added add/remove statistics. "Removed" figure is a combination of composites and primes under sqrt of limit. Getting a >20x speedup from using a BitArray instead of a Set.

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

""" Pritchard sieve of primes up to limit. Uses type of `limit` arg for type of primes """
function pritchard(limit::T, verbose=false) where {T<:Integer}
    members = falses(limit)
    members[1] = true
    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(isqrt(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 1:steplength
                if members[w]
                    n = w + steplength
                    while n <= nlimit
                        members[n] = true
                        ac += 1
                        n += steplength
                    end
                end
            end
            steplength = nlimit # advance wheel size
        end
        np = 5
        mcopy = copy(members)
        for w in 1:nlimit
            if mcopy[w]
                np == 5 && w > prime && (np = w)
                n = prime * w
                n > nlimit && break
                rc += 1
                members[n] = false
            end
        end
        np < prime && break
        push!(primes, prime)
        prime = prime == 2 ? 3 : np
        nlimit = min(steplength * prime, limit) # advance wheel limit
    end
    members[1] = false
    newprimes = [i for i in eachindex(members) if members[i]]
    verbose && println(
        "up to $limit, added $ac, removed $rc, prime count ",
        length(primes) + length(newprimes),
    )
    return append!(primes, newprimes)
end

println(pritchard(150))
pritchard(1000000, true)
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]

up to 1000000, added 186825, removed 108494, prime count 78498

Nim

Translation of: Python
import std/[algorithm, math, sugar]

proc pritchard(limit: Natural): seq[int] =
  ## Pritchard sieve of primes up to "limit".
  var members = newSeq[bool](limit + 1)
  members[1] = true
  var
    stepLength = 1
    prime = 2
    rtlim = sqrt(limit.toFloat).int
    nlimit = 2
    primes: seq[int]

  while prime <= rtlim:
    if stepLength < limit:
      for w in 1..members.high:
        if members[w]:
          var n = w + stepLength
          while n <= nlimit:
            members[n] = true
            inc n, stepLength
      stepLength = nlimit

    var np = 5
    var mcpy = members
    for w in 1..members.high:
      if mcpy[w]:
        if np == 5 and w > prime:
          np = w
        let n = prime * w
        if n > limit:
          break  # No use trying to remove items that can't even be there.
        members[n] = false  # No checking necessary now.

    if np < prime:
      break
    primes.add prime
    prime = if prime == 2: 3 else: np
    nlimit = min(stepLength * prime, limit)   # Advance wheel limit.

  let newPrimes = collect:
                    for i in 2..members.high:
                      if members[i]: i
  result = sorted(primes & newPrimes)


echo pritchard(150)
echo "Number of primes up to 1_000_000: ", pritchard(1_000_000).len


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]
Number of primes up to 1_000_000: 78498

Pascal

Works with: Free Pascal

A console program written in Delphi 7. It follows the algorithm in the Wikipedia article "Sieve of Pritchard".

program Pritchard_console;

{$APPTYPE CONSOLE}

uses
  Math, SysUtils, Types;

// Function to return an array of all primes <= N
function PritchardSieve( const N : integer) : Types.TIntegerDynArray;
var
  j, j_max, k, len, nrPrimes, p : integer;
  marked : Types.TBooleanDynArray;
  smallPrimes : Types.TIntegerDynArray; // i.e. primes <= sqrt( N)
  spi : integer; // index into array smallPrimes
const
  SP_STEP = 16; // step when extending dynamic array smallPrimes
begin
  // Deal with trivial input
  result := nil;
  if (N <= 1) then exit;

  // Initialize
  SetLength( marked, N + 1); // 0..N for convenience; marked[0] is not used
  marked[1] := true; // no other initialization of "marked" is needed
  len := 1;
  p := 2;
  SetLength( smallPrimes, SP_STEP);
  spi := 0;

  while p*p <= N do begin
    // Roll the wheel
    if len < N then begin
      j_max := Math.Min( p*len, N);
      for j := len + 1 to j_max do marked[j] := marked[j - len];
      len := j_max;
    end;

    // Unmark multiples of p
    for k := len div p downto 1 do
      if marked[k] then marked[p*k] := false;

    // Store the prime p, extending the array if necessary
    if spi = Length( smallPrimes) then
      SetLength( smallPrimes, spi + SP_STEP);
    smallPrimes[spi] := p;
    inc(spi);

    // Find the next prime p
    if p = 2 then p := 3
    else repeat inc(p) until (p > N) or marked[p];
    // Condition p > N is a safety net; should always hit a marked value
    Assert(p <= N);
  end; // while

  // Final roll, if needed. It is not needed if N >= 49. This is because
  //  2 < 3^2, 2*3 < 5^2, 2*3*5 < 7^2, but thereafter 2*3*5*7 > 11^2, etc.
  if len < N then
    for j := len + 1 to N do marked[j] := marked[j - len];

  // Remove 1 and put the small primes back
  marked[1] := false;
  for k := 0 to spi - 1 do marked[smallPrimes[k]] := true;

  // Use the boolean array to return an array of prime integers
  nrPrimes := 0;
  for j := 2 to N do
    if marked[j] then inc( nrPrimes);
  SetLength( result, nrPrimes);
  k := 0;
  for j := 2 to N do
    if marked[j] then begin result[k] := j; inc(k); end;
end;

// Main routine. User types the program name,
// optionally followed by the limit N (defaults to 150)
var
  N, j : integer;
  primes : Types.TIntegerDynArray;
begin
  if ParamCount = 0 then N := 150
                    else N := SysUtils.StrToInt( ParamStr(1));
  primes := PritchardSieve(N);
  WriteLn( 'Number of primes = ', Length(primes));
  for j := 0 to Length(primes) - 1 do begin
    Write( ' ', primes[j]:4);
    if j mod 10 = 9 then WriteLn;
  end;
end.
Output:
Number of primes = 35
    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

Free Pascal

see Sieve_of_Eratosthenes#alternative_using_wheel

Perl

Translation of: Raku
use v5.36;
use List::Util 'min';

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

sub next_($w) { $s[$w-1] }
sub prev_($w) { $s[$w-2] }

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

sub delete_multiples_of($p) {
    my $f = $p;
    while ($p*$f <= $length) { $f = next_ $f                   }
    while (   $f >  1      ) { $f = prev_ $f; delete_pf($p*$f) }
}

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

sub extend_to($n) {
    my($w, $x) = (1, $length+1);
    while ($x <= $n) {
        append $x;
        $w = next_ $w;
        $x = $length + $w;
    }
    $length = $n;
    append $limit+2 if $length == $limit
}

do {
    extend_to min $p*$length, $limit if $length < $limit;
    delete_multiples_of $p;
    $p = next_ 1;
    extend_to $limit if $length < $limit
} until $p*$p > $limit;

my @primes = 2;
for (my $p = 3; $p <= $limit; $p = next_ $p) { push @primes, $p }
say "@primes";
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 """

from numpy import ndarray
from math import isqrt


def pritchard(limit):
    """ Pritchard sieve of primes up to limit """
    members = ndarray(limit + 1, dtype=bool)
    members.fill(False)
    members[1] = True
    steplength, prime, rtlim, nlimit = 1, 2, isqrt(limit), 2
    primes = []
    while prime <= rtlim:
        if steplength < limit:
            for w in range(1, len(members)):
                if members[w]:
                    n = w + steplength
                    while n <= nlimit:
                        members[n] = True
                        n += steplength
            steplength = nlimit

        np = 5
        mcpy = members.copy()
        for w in range(1, len(members)):
            if mcpy[w]:
                if np == 5 and w > prime:
                    np = w
                n = prime * w
                if n > nlimit:
                    break  # no use trying to remove items that can't even be there
                members[n] = False  # no checking necessary now

        if np < prime:
            break
        primes.append(prime)
        prime = 3 if prime == 2 else np
        nlimit = min(steplength * prime, limit)  # advance wheel limit

    newprimes = [i for i in range(2, len(members)) if members[i]]
    return sorted(primes + newprimes)


print(pritchard(150))
print('Number of primes up to 1,000,000:', len(pritchard(1000000)))
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]
Number of primes up to 1,000,000: 78498

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 
Cookies help us deliver our services. By using our services, you agree to our use of cookies.