Twin Primes counter in Nim using P5 prime generator residues

twinprimes_ssozp5a2a.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:gcc --d:release --gc:none twinprimes_ssozp5a2a.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: 2018/03/05

 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
import times, os              # for timing code execution

# 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 = [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 = [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)        # convert residue val -> residues[] index
  for i in 0..rscnt-1: posn[res[i]-2] = i

  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 = 0                         # from 1st residue in num's resgroup
  while num >= modk+res[r]: r += 1  # find last pc val|position <= num
  let maxpcs = k*rscnt + r          # max number of prime candidates <= num
  var prms = newSeq[bool](maxpcs)   # array of prime candidates set False

  let sqrtN =int(sqrt float64(num)) # compute integer sqrt of input num
  modk = 0; r = -1; 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 = 0; 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:                  # mark prime's multiples in prms
      let prod = prm_r * ri - 2     # 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 < maxpcs: 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 = -1
  for prm in prms:                  # numerate|store primes from pcs list
    r += 1; if r == rscnt: (r = 0; modk += md)
    if not prm: primes.add(modk + res[r]) # put prime in global 'primes' list
  pcnt = primes.len                       # 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-1: pos[residues[i]-2]= i # convert (x mod modpg) vals

  # 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:                       # for each prime|residue pair
      let prod = r * ri - 2                   # compute res cross-products
      var res_indx = pos[prod mod modpg] - 1  # determine restrack index value
      if res_indx == -1 or res_indx == 4: continue # skip if for ri = 7 or 23
      if res_indx > 4: res_indx -= 1
      let row = res_indx * pcnt
      # compute|store resgroup val of 1st prime mult for prime|residue pair
      nextp[row + j] = uint(k*(prime + ri) + prod 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 residue_sieve(Kn, r, indx: int) =
  let biti = uint8(1 shl r)        # set its residue track bit mask
  let row  = indx * pcnt           # set address to its restrack in 'nextp'
  for j, prime in primes:          # for each prime r1..sqrt(N)
    if nextp[row + j] < Kn.uint:   # 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] -= Kn.uint  # do if 1st mult resgroup not within seg     

# 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-1: seg[b] = 0x21    # set all bits to prime except ri 7|23
  for indx in 0..5:                  # for nextp row indexes, 
      var r = indx + 1               # indx|r: 0 -> 1, 1 -> 2, 2 -> 3, 3 -> 4
      if indx > 3: r += 1            # indx|r: 4 -> 6, 5 -> 7 restracks in seg
      residue_sieve(Kn,r,indx)       # sieve along just the required restracks
  for byt in 0..Kn-1:                # for every resgroup|byte in the segment
    primecnt += uint(pbits[seg[byt]])# count the '0' bit pairs as twin primess

# 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-1:                  # 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 ts = epochTime()             # start timing sieve setup execution

  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 = 0                        # starting with first residue
  while num >= modk+uint(residues[r]): r += 1 # find last pc position <= num
  let maxpcs = k*rescnt + r.uint   # 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-2, "x", pcnt, "] array")

  nextp = newSeq[uint64]((rescnt-2)*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 = KB.uint                 # to make code shorter and less noisy

  echo("perform Twin Prime Segmented SoZ")

  while Ki < Kmax:                 # for KB size resgroup slices up to Kmax
    Kn = min(KB, int(Kmax-Ki))     # set segment slice resgroup length
    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

  let t2 = (epochTime()-ts)        # total execution time

  echo("last segment = ", Kn, " resgroups; segment slices = ", Ki div kb)
  echo("total twins = ", primecnt, "; last twin = ", lprime-1, "+/-1")
  echo("total time = ", t2.formatFloat(ffDecimal, 3), " secs\n")

twinprimes()