This task has been clarified. Its programming examples are in need of review to ensure that they still fit the requirements of the task.


The Sieve of Eratosthenes is a simple algorithm that finds the prime numbers up to a given integer.



Task

Implement the   Sieve of Eratosthenes   algorithm, with the only allowed optimization that the outer loop can stop at the square root of the limit, and the inner loop may start at the square of the prime just found.

That means especially that you shouldn't optimize by using pre-computed wheels, i.e. don't assume you need only to cross out odd numbers (wheel based on 2), numbers equal to 1 or 5 modulo 6 (wheel based on 2 and 3), or similar wheels based on low primes.

If there's an easy way to add such a wheel based optimization, implement it as an alternative version.



Note

 It is important that the sieve algorithm be the actual algorithm used to find prime numbers for the task.


Related tasks

   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


 
 [Sieve of Eratosthenes]
 [EDSAC program. Initial Orders 2]
 
[Memory usage:
   56..87   library subroutine P6, for printing
   88..222  main program
  224..293  mask table: 35 long masks; each has 34 1's and a single 0
  294..1023 array of bits for integers 2, 3, 4, ...,
            where bit is changed from 1 to 0 when integer is crossed out.
  The address of the mask table must be even, and clear of the main program.
  To change it, just change the value after "T47K" below.
  The address of the bit array will then be changed automatically.]
 
[Subroutine M3, prints header, terminated by blank row of tape.
 It's an "interlude", which runs and then gets overwritten.]
      PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
      @&*SIEVE!OF!ERATOSTHENES!#2020
      @&*BASED!ON!CODE!BY!EIITI!WADA!#2001
      ..PZ
 
[Subroutine P6, prints strictly positive integer.
 32 locations; working locations 1, 4, 5.]
        T  56 K
  GKA3FT25@H29@VFT4DA3@TFH30@S6@T1FV4DU4DAFG26@TFTF
  O5FA4DF4FS4FL4FT4DA1FS3@G9@EFSFO31@E20@J995FJF!F
 
  [Store address of mask table in (say) location 47
  (chosen because its code letter M is first letter of "Mask").
  Address must be even and clear of main program.]
        T  47 K
        P 224 F  [<-------- address of mask table]
 
[Main program]
        T  88 K  [Define load address for main program.
                 Must be even, because of long values at start.]
        G     K [set @ (theta) for relative addressing]
 
[Long constants]
        T#Z PF TZ [clears sandwich digit between 0 and 1]
    [0] PD PF     [long value 1; also low word = short 1]
        T2#Z PF T2Z [clears sandwich digit between 2 and 3]
    [2] PF K4096F [long value 1000...000 binary;
                   also high word = teleprinter null]
 
 [Short constants
  The address in the following C order is the (exclusive) end of the bit table.
  Must be even: max = 1024, min = M + 72 where M is address of mask table set up above.
  Usually 1024, but may be reduced, e.g. to make the program run faster.]
    [4] C1024 D [or e.g. C 326 D to make it much faster]
    [5] U     F ['U' = 'T' - 'C']
    [6] K     F ['K' = 'S' - 'C']
    [7] H    #M [H order for start of mask table]
    [8] H  70#M [used to test for end of mask table]
    [9] P   2 F [constant4, or 2 in address field]
   [10] P  70 F [constant 140, or 70 in address field]
   [11] @     F [carriage return]
   [12] &     F [line feed]
 
 [Short variables]
   [13] P   1 F [p = number under test
                Let p = 35*q + r, where 0 <= r < 35]
   [14] P     F [4*q]
   [15] P   4 F [4*r]
 
 [Initial values of orders; required only for optional code below.]
   [16] C  70#M [initial value of a variable C order]
   [17] T    #M [initial value of a variable T order]
   [18] T  70#M [initial value of a variable T order]
 
   [19]
   [Enter with acc = 0]
 
  [Optional code to do some initializing at run time.
   This code allows the program to run again without being loaded again.]
         A  7 @ [initial values of variable orders]
         T 65 @
         A 16 @
         T 66 @
         A 17 @
         T 44 @
         A 18 @
         T 52 @
 
  [Initialize variables]
         A    @ [load 1 (short)]
         L    D [shift left 1]
         U 13 @ [p := 2]
         L  1 F [shift left 2]
         T 15 @ [4*r :=  8]
         T 14 @ [4*q :=  0]
  [End of optional code]
 
 [Make table of 35 masks 111...110,  111...101, ...,  011...111
  Treat the mask 011...111 separately to avoid accumulator overflow.
  Assume acc = 0 here.]
        S    #@ [acc all 1's]
        S  2 #@ [acc := 0111...111]
   [35] T 68 #M [store at high end of mask table]
        S    #@ [acc := -1]        
        L     D [make mask 111...1110]
        G 43  @ [jump to store it
 [Loop shifting the mask right and storing the result in the mask table.
  Uses first entry of bit array as temporary store.]
   [39] T     F [clear acc]
        A 70 #M [load previous mask]
        L     D [shift left]
        A    #@ [add 1]
   [43] U 70 #M [update current mask]
   [44] T    #M [store it in table (order changed at run time)]
        A 44  @ [load preceding T order]
        A  9  @ [inc address by 2]
        U 44  @ [store back]
        S 35  @ [reached high entry yet?]
        G 39  @ [loop back if not]
 [Mask table is now complete]
 
 [Initialize bit array: no numbers crossed out, so all bits are 1]
   [50] T     F [clear acc]
        S     #@ [subtract long 1, make top 35 bits all 1's]
   [52] T 70 #M [store as long value, both words all 1's  (order changed at run time)]
        A  52 @ [load preceding order]
        A   9 @ [add 2 to address field]
        U  52 @ [and store back]
        S   5 @ [convert to C order with same address (*)]
        S   4 @ [test for end of bit array]
        G  50 @ [loop until stored all 1's in bit table]
 [(*) Done so that end of bit table can be stored at one place only
      in list of constants, i.e. 'C m D' only, not 'T m D' as well.]
 
 [Start of main loop.]
 [Testing whether number has been crossed out]
   [59] T     F [acc := 0]
        A  66 @ [deriving S order from C order]
        A   6 @
        T  64 @
        S    #@ [acc := -1]
   [64] S     F [acc := 1's complement of bit-table entry (order changed at run time)]
   [65] H    #M [mult reg := start of mask array (order changed at run time)]
   [66] C  70#M [acc := -1 iff p (current number) is crossed out (order changed at run time)]
 [The next order is to avoid accumulator overflow if acc = max positive number]
        E  70 @ [if acc >= 0, jump to process new prime]
       A     #@ [if acc < 0, add 1 to test for -1]
       E 106  @ [if acc now >= 0 number is crossed out, jump to test next]
 [Here if new prime found.
  Send it to the teleprinter]
   [70] O  11 @ [print CR]
        O  12 @ [print LF]
        T     F [clear acc]
        A  13 @ [load prime]
        T     F [store in C(0) for print routine]
        A  75 @ [for subroutine return]
        G 56  F [print prime]
 
 [Cross out its multiples by setting corresponding bits to 0]
        A  65 @ [load H order above]
        T 102 @ [plant in crossing-out loop]
        A  66 @ [load C order above]
        T1 03 @ [plant in crossing-out loop]
 
 [Start of crossing-out loop. Here acc must = 0]
   [81] A 102 @ [load H order below]
        A  15 @ [inc address field by 2*r, where p = 35q + r]
        U 102 @ [update H order]
        S   8 @ [compare with 'H 70 #M']
        G  93 @ [skip if not gone beyond end of mask table]
        T     F [wrap mask address and inc address in bit tsble]
        A 102 @ [load H order below]
        S  10 @ [reduce mask address by 70]
        T 102 @ [update H order]
        A 103 @ [load C order below]
        A   9 @ [add 2 to address]
        T 103 @ [update C order]
   [93] T     F [clear acc]
        A 103 @ [load C order below]
        A  14 @ [inc address field by 2*q, where p = 35q + r]
        U 103 @ [update C order]
        S   4 @ [test for end of bit array]
        E 106 @ [if finished crossing out, loop to test next number]
        A  4  @ [restore C order]
        A  5  @ [make T order with same address]
        T 104 @ [store below]
 
 [Execute the crossing-out orders created above]
  [102] X     F [mult reg := mask (order created at run time)]
  [103] X     F [acc := logical and with bit-table entry (order created at run time)]
  [104] X     F [update entry (order created at run time)]
        E  81 @ [loop back with acc = 0]
 
  [106] T     F [clear acc]
        A  13 @ [load p = number under test]
        A     @ [add 1 (single)]
        T  13 @ [update]
        A  15 @ [load 4*r, where p = 35q + r]
        A   9 @ [add 4]
        U  15 @ [store back (r inc'd by 1)]
        S  10 @ [is 4*r now >= 140?]
        G 119 @ [no, skip]
        T  15 @ [yes, reduce 4*r by 140]
        A  14 @ [load 4*q]
        A   9 @ [add 4]
        T  14 @ [store back (q inc'd by 1)]
  [119] T     F [clear acc]
        A  65 @ [load 'H ... D' order, which refers to a mask]
        A   9 @ [inc mask address by 2]
        U  65 @ [update order]
        S   8 @ [over end of mask table?]
        G  59 @ [no, skip wrapround code]
        A   7 @ [yes, add constant to wrap round]
        T  65 @ [update H order]
        A  66 @ 
        A   9 @ [inc address by 2]
        U  66 @ [and store back]
        S   4 @ [test for end, as defined by C order at start]
        G  59 @ [loop back if not at end]
 
[Finished whole thing]
  [132] O   3 @ [output null to flush teleprinter buffer]
        Z     F [stop]
        E  19 Z [address to start execution]
        P     F [acc = 0 at start]
 