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

twinprimes_ssozp5.nim

#[
 This Nim source file is a single threaded implementation to perform an
 extremely fast Segmented Sieve of Zakiya (SSoZ) to find Twin 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}.
 For just Twin Primes, use generator: Pn = 30*k + {11,13,17,19,29,31}

 Adjust segment byte length parameter B (usually L1|l2 cache size)
 for optimum operational performance for cpu being used.
 
 Verified on Nim 0.17, using clang (smaller exec) or gcc

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

 Then run executable: $ ./twinprimes_ssozp5 <cr>, and enter value for N.
 As coded, input values cover the range: 13 -- 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/30

 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

# Convert segment byte values into number of twin primes in byte.
# For P5 r1=7 and r6=23 aren't used, so byte-mask is b_00100001
# A byte can identify 3 tps, '001xxxx1', 'xx100xx1', 'xx1xx001'
# Thus, convert byte values 0-255 into table of twin primes
# Ex: 33=b_00100001 -> 3; 35=b_00100011 -> 2; 43=b_00101011 -> 1
let pbits = [
 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,3,0,2,0,2,0,2,0,2,0,1,0,1,0,1,0,2,0,1,0,1,0,1,0,2,0,1,0,1,0,1
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,2,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,2,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,2,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,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|residue pair
      nextp[row + j] = uint(k*(prime + ri) + (prod-2) div modpg)

# This routine performs the twin prime sieve for a segment of Kn resgroups.
# Each 'seg' bit represents a restrack of Kn resgroups, but restracks for 
# ri 7 and 23 are set as nonprime. The 'nextp' resgroup vals mark the prime
# mults in 'seg' on each restrack, and are udpated for each prime for the
# next segment. The twin prime count for each 'seg' byte is then added.
proc segsieve(Kn: int) =             # for Kn resgroups|bytes in segment
  for b in 0 .. <Kn: seg[b] = 0x21   # set all bits to prime except ri 7|23
  for r in 1 .. <rescnt:             # for restracks, starting at ri=11, r=1
    if r == 0x20: continue           # skip sieve for ri = 23, r = 0x20
    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
                                     # count the twin primes in the segment
  for byt in seg[0..<Kn]:            # for every resgroup|byte in the segment
    primecnt += uint(pbits[byt])     # count the '0' bit pairs as twin primes

# Print even mid values 'n' between twin primes, e.g. n-1|n+1 are the primes. 
# Can store as: $ ./ssozp5 {| less} or {> <primesfile>} for easier viewing
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:                 # then for each of Kn seg resgroups|bytes
    if (int(not seg[k]) and 0x06)==0x06: # if restracks ri = 11|13 are prime
      write(stdout, modk+12, " ")    # print even value between twin primes
    if (int(not seg[k]) and 0x18)==0x18: # if restracks ri = 17|19 are prime
      write(stdout, modk+18, " ")    # print even value between twin primes
    if (int(not seg[k]) and 0xc0)==0xc0: # if restracks ri = 29|31 are prime
      write(stdout, modk+30, " ")    # print even value between twin primes
    modk += modpg                    # incr base num to next resgroup in seg 
  echo("\n")

proc twinprimes() =                # Use P5 prime generator for SSoZ
  stdout.write "Enter integer number: "
  let val = stdin.readline.parseBiggestUInt # find primes <= val (13..2^64-1)

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

  echo("segment has ", KB, " bytes and 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                        # starting with first residue
  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
  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 = 2                     # for first two twin primes (3,5)|(5,7)
  var Kn = KB                      # set sieve resgroups size to segment size
  var Ki = 0'u64                   # starting resgroup index for each segment
  let kb = uint(KB)                # to make code shorter and less noisy

  echo("perform Twin Prime Segmented SoZ")

  while Ki < Kmax:                 # for KB size resgroup 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 middle tp vals for the seg (opt)
    Ki += kb                       # first resgroup index for next seg slice

  var lprime = 0'u64               # get last twin 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
  while true:                      # step backwards from end of last resgroup
    if (int(not seg[b]) and 0xc0)==0xc0: # if ri = 29|31 restracks both prime
      lprime = modk + 31           # numerate|save larger twin prime value
      if lprime <= num: break      # if its <= num its the last prime, so exit
      primecnt -= 1                # else reduce primecnt, keep backtracking
    if (int(not seg[b]) and 0x18)==0x18: # if ri = 17|19 restracks both prime
      lprime = modk + 19           # numerate|save larger twin prime value
      if lprime <= num: break      # if its <= num its the last prime, so exit
      primecnt -= 1                # else reduce primecnt, keep backtracking
    if (int(not seg[b]) and 0x06)==0x06: # if ri = 11|13 restracks both prime
      lprime = modk + 13           # numerate|store larger twin prime value
      if lprime <= num: break      # if its <= num its the last prime, so exit
      primecnt -= 1                # else reduce primecnt, keep backtracking
                                   # reduce modk val for next smaller resgroup
    modk -= modpg; b -= 1          # if twin prime <= num not in this resgroup

  echo("last segment = ", Kn, " resgroups; segment slices = ", Ki div kb)
  echo("total twins = ", primecnt, "; last twin = ", lprime-1, "+/-1\n")

twinprimes()