Find primes <= N, using sequential Segmented Sieve of Zakiya (SSoZ) with P5 prime generator, written in Nim

ssozp5.nim

#[
 This Nim source file will compile to an executable program to
 perform the Segmented Sieve of Zakiya (SSoZ) to find primes <= N.
 It is based on the P5 Strictly Prime (SP) Prime Generator.

 Prime Genrators have the form:  mod*k + ri; ri -> {1,r1..mod-1}
 The residues ri are integers coprime to mod, i.e. gcd(ri,mod) = 1
 For P5, mod = 2*3*5 = 30 and the number of residues are
 rescnt = (2-1)(3-1)(5-1) = 8, which are {1,7,11,13,17,19,23,29}.

 Adjust segment byte length parameter B (usually L1|l2 cache size)
 for optimum operational performance for cpu|system being used.

 Verified on Nim 0.17, using clang (smaller exec) or gcc

 $ nim c --cc:[clang|gcc] --d:release  ssozp5.nim

 Then run executable: $ ./ssozp5 <cr>, and enter value for N.
 As coded, input values cover the range: 7 -- 2^64-1

 Related code, papers, and tutorials, are downloadable here:

 http://www.4shared.com/folder/TcMrUvTB/_online.html

 Use of this code is free subject to acknowledgment of copyright.
 Copyright (c) 2017 Jabari Zakiya -- jzakiya at gmail dot com
 Version Date: 2017/08/28

 This code is provided under the terms of the
 GNU General Public License Version 3, GPLv3, or greater.
 License copy/terms are here:  http://www.gnu.org/licenses/
]#

import math                   # for sqrt function
import strutils, typetraits   # for number input

# Global var used to count number of primes in 'seg' variable.
# Each value is number of '0' bits (primes) for values 0..255.
let pbits = [
 8,7,7,6,7,6,6,5,7,6,6,5,6,5,5,4,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3
,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2
,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2
,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1
,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2
,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1
,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1
,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,4,3,3,2,3,2,2,1,3,2,2,1,2,1,1,0
]

# The global residue values for the P5 prime generator.
let residues = [1, 7, 11, 13, 17, 19, 23, 29, 31]

# Global parameters
const
  modpg = 30           # mod value for P5 prime generator
  rescnt = 8           # number of residues for P5 prime generator
var
  pcnt = 0             # number of primes from r1..sqrt(N)
  primecnt = 0'u64     # number of primes <= N
  nextp: seq[uint64]   # table of regroups vals for primes multiples
  primes: seq[int]     # list of primes r1..sqrt(N)
  seg: seq[uint8]      # segment byte array to perform ssoz

# This routine is used to compute the list of primes r1..sqrt(input_num),
# stored in global var 'primes', and its count stored in global var 'pcnt'.
# Any algorithm (fast|small) can be used. Here the SoZ using P7 is used.
proc sozp7(val: int | int64) =      # Use P7 prmgen to finds primes upto val
  let md = 210                      # P7 mod value
  let rscnt = 48                    # P7 residue count and residues list
  let res = [1, 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,121,127,131,137,139
      ,143,149,151,157,163,167,169,173,179,181,187,191,193,197,199,209,211]

  var posn = newSeq[int](md)               # small array (or hash) to convert
  for i in 0 .. <rscnt: posn[res[i]] = i-1 # (x mod md) vals -> -1,0,1,2..46

  let num = (val-1) or 1            # if val even then subtract 1
  var k = num div md                # compute its residue group val
  var modk = md * k                 # compute its base num
  var r = 1                         # from 1st residue in num's resgroup
  while num >= modk+res[r]: r += 1  # find last pc val|position <= num 
  let maxprms = k*rscnt + r - 1     # max number of prime candidates <= num
  var prms = newSeq[bool](maxprms)  # array of prime candidates set False

  let sqrtN =int(sqrt float64(num)) # compute integer sqrt of input num
  modk = 0; r = 0; k = 0

  # process to mark the prime multiples in the list of pcs r1..sqrtN
  for prm in prms:                  # from list of pc positions in prms
    r += 1; if r > rscnt: (r = 1; modk += md; k += 1)
    if prm: continue                # if pc not prime go to next location
    let prm_r = res[r]              # if prime save its residue value
    let prime = modk + prm_r        # numerate the prime value
    if  prime > sqrtN: break        # we're finished if it's > sqrtN
    let prmstep = prime * rscnt     # prime's stepsize to mark its mults
    for ri in res[1..rscnt]:        # mark prime's multiples in prms
      let prod = prm_r * ri         # residues cross-product for this prime
      # compute resgroup val of 1st prime multiple for each restrack
      # starting there, mark all prime multiples on restrack upto end of prms
      var prm_mult = (k*(prime + ri) + prod div md)*rscnt + posn[prod mod md]
      while prm_mult < maxprms: prms[prm_mult] = true; prm_mult += prmstep
      
  # prms now contains the nonprime positions for the prime candidates r1..N
  # extract primes into global var 'primes' and count into global var 'pcnt'
  primes = @[7]                     # r1 prime for P5
  modk = 0; r=0
  for prm in prms:                  # numerate|store primes from pcs list
    r += 1; if r > rscnt: (r = 1; modk += md)
    if not prm: primes.add(modk + res[r]) # put prime in global 'primes' list
  pcnt = len(primes)                      # set global count of primes

# This routine initializes the [rescnt x pcnt] global var 'nextp' table with
# resgroup vals of the 1st prime multiples for each prime r1..sqrtN along the
# the residue tracks for the used prime generator.
proc nextp_init() =
  var pos = newSeq[int](modpg)                 # create small array to
  for i in 0..<rescnt: pos[residues[i]] = i-1  # convert (x mod modpg) vals
  pos[1] = rescnt-1                            # to restracks -> 0,1..rescnt-1

  # load 'nextp' with 1st prime multiples regroups vals on each residue track
  for j, prime in primes:                      # for each prime r1..sqrt(N)
    let k = (prime-2) div modpg                # find the resgroup it's in
    let r = (prime-2) mod modpg + 2            # and its residue value
    for ri in residues[1 .. rescnt]:           # for each prime|residue pair
      let prod = r * ri                        # compute res cross-products
      let row  = pos[prod mod modpg] * pcnt    # compute restrack address
      # compute|store resgroup val of 1st prime mult for prime|ri residue pair
      nextp[row + j] = uint(k*(prime + ri) + (prod-2) div modpg)

# This routine performs the prime sieve for a segment of Kn resgroups|bytes.
# Each 'seg' bit represents a restrack of Kn resgroups, which are processed
# sequentially. The 'nextp' resgroup values mark the prime multiples in 'seg'
# along each restrack, and are udpated for each prime for the next segment.
# Then the prime count for each 'seg' byte is added to global var 'primecnt'.
proc segsieve(Kn: int) =              # for Kn resgroups in segment
  for b in 0 .. <Kn: seg[b] = 0       # initialize seg bits to all prime '0'
  for r in 0 .. <rescnt:              # for each prime generator residue
    let biti = uint8(1 shl r)         # set its residue track bit mask
    let row  = r * pcnt               # set address to its restrack in 'nextp'
    for j, prime in primes:           # for each prime r1..sqrt(N)
      if nextp[row + j] < uint(Kn):   # if 1st mult resgroup is within 'seg'
        var k = int(nextp[row + j])   # starting from this resgroup in 'seg'
        while k < Kn:                 # for each primenth byte to end of 'seg'
          seg[k] = seg[k] or biti     # mark restrack bit in byte as nonprime
          k += prime                  # compute next prime multiple resgroup
        nextp[row + j] = uint(k - Kn) # save 1st resgroup in next eligible seg
      else: nextp[row+j] -= uint(Kn)  # do if 1st mult resgroup not within seg
                                      # now count the primes in the segment
  for byt in seg[0..<Kn]:             # for every resgroup|byte in the segment
    primecnt += uint(pbits[byt])      # count the '0' bits as primes

# Print segment primes. Can store as: $ ./ssozp5 {| less} or {> <primesfile>}
# to make viewing primes easier
proc printprms(Kn: int, Ki: uint64) = # starting at resgroup value Ki for seg
  var modk = Ki * modpg               # compute its base num
  for k in 0 .. <Kn:                  # for Kn residues groups|bytes in seg
    for r in 0 .. <rescnt:            # for each restrack bit in byte
      if (int(seg[k]) and (1 shl r)) == 0:  # if it's '0' it's a prime
        write(stdout, modk + uint(residues[r+1]), " ") # numerate|print value
    modk += modpg                     # incr base num to next resgroup in seg
  echo "\n"

proc ssozp5() =                    # Perform SSoZ using P5 prime generator
  stdout.write "Enter integer number: "
  let val = stdin.readline.parseBiggestUInt # input integer value (7..2^64-1)

  const B = 256 * 1024             # adjust size to optimize for system cache
  let KB = uint(B)                 # number of segment resgroups
  seg = newSeq[uint8](B)           # create segment array of B bytes size

  echo("segment has ", B, " bytes and ", KB, " residues groups")

  let num = uint64((val-1) or 1)   # if val even subtract 1
  var k = num div modpg            # compute its resgroup
  var modk = k * modpg             # compute its base num
  var r = 1                        # from 1st residue in num's resgroup
  while num >= modk+uint(residues[r]): r += 1 # find last pc position <= num
  let maxpcs = k*rescnt+ uint(r-1) # maximum number of prime candidates <= num
  let Kmax = (num-2) div modpg + 1 # maximum number of resgroups for num

  echo("prime candidates = ", maxpcs, "; resgroups = ", Kmax)

  let sqrtN=int(sqrt float64(num)) # compute integer sqrt of input num
  sozP7(sqrtN)                     # get pcnt and primes <= sqrt(num)

  echo("create nextp[", rescnt, "x", pcnt, "] array")

  nextp = newSeq[uint64](rescnt*pcnt) # create size of global 'nextp' table
  nextp_init()                     # load with first prime multiples resgroups

  primecnt = 3                     # 2,3,5 the P5 excluded primes count
  var Kn = int(KB)                 # set sieve resgroups size to segment size
  var Ki = 0'u64                   # starting resgroup index for first segment

  echo("perform segmented SoZ")

  while Ki < Kmax:                 # for KB resgroup size slices up to Kmax
    if Kmax-Ki < KB: Kn = int(Kmax-Ki) # this is to set last segment size
    segsieve(Kn)                   # sieve and count primes for this segment
    #printprms(Kn,Ki)              # print primes for the segment (optional)
    Ki += KB                       # first resgroup index for next seg slice

  var lprime = 0'u64               # get last prime and primecnt <= num
  modk = (Kmax-1) * modpg          # set mod for last resgroup in last segment
  var b = Kn-1                     # set val for last resgroup in last segment
  r = rescnt-1                     # starting from last restrack
  while true:                      # step backwards from end of last resgroup
    if (int(seg[b]) and (1 shl r)) == 0:  # if restrack bit indicates a prime
      lprime = modk + uint(residues[r+1]) # numerate the prime value
      if lprime <= num: break      # if its <= num its the last prime, so exit
      primecnt -= 1                # else reduce primecnt, keep backtracking
                                   # reduce restrack, next resgroup if needed
    r -= 1; if r < 0: (r = rescnt-1; modk -= modpg; b -= 1)

  echo("last segment = ", Kn, " resgroups; segment slices = ", Ki div KB)
  echo("total primes = ", primecnt, "; last prime = ", lprime)

ssozp5()