Below is the Nim direct translation of the C++ code in my paper. I have
verified it on Nim 0.17 on Linux. In comparison on this same machine the Nim
code runs a bit faster for some (larger >= 10^9) test values (not
comprehensibly benchmarked yet).
I'll post it to a gist later to make it easier to download. I'm now
re-architecting it to make it parallel friendlier (I hope). Feedback welcome on
style, and performance tips.
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 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/21
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
import strutils, typetraits
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
]
let residues = [1, 7, 11, 13, 17, 19, 23, 29, 31]
# Global parameters
const
modp5 = 30 # prime generator mod value
rescnt = 8 # number of residues for prime generator
var
pcnt = 0 # number of primes from r1..sqrt(N)
primecnt = 0'u64 # number of primes <= N
next: seq[uint64] # array of primes first multiples
primes: seq[int] # array of primes <= sqrt(N)
seg: seq[uint8] # seg[B] segment byte array
proc sozp7(val: int | int64) =
let md = 210
let rscnt = 48
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](210)
for i in 0 .. <rscnt: posn[res[i]] = i-1
let num = (val-1) or 1 # if value even number then subtract 1
var k = num div md # kth residue group
var modk = md * k # base num value
var r = 1 # starting with first residue
while num >= modk+res[r]: r += 1 # find last pc position <= num
let maxprms = k*rscnt + r - 1 # maximum number of prime candidates
var prms = newSeq[bool](maxprms) # array of prime candidates set False
let sqrtN = int(sqrt float64(num))
modk=0; r=0; k=0
# sieve to eliminate nonprimes from small primes prms array
for prm in prms:
r += 1; if r > rscnt: (r = 1; modk += md; k += 1)
if prm: continue
let res_r = res[r]
let prime = modk + res_r
if prime > sqrtN: break
let prmstep = prime * rscnt
for ri in res[1..rscnt]:
let prod = res_r * ri # residues cross-product
var nonprm = (k*(prime + ri) + prod div md)*rscnt + posn[prod mod md]
while nonprm < maxprms: prms[nonprm] = true; nonprm += prmstep
# the prms array now has all the positions for primes r1..N
# extract prime numbers and count from prms into primes array
primes = @[7] # r1 prime for P5
modk = 0; r = 0
for prm in prms:
r += 1; if r > rscnt: (r = 1; modk += md)
if not prm: primes.add(modk + res[r])
pcnt = len(primes)
proc next_init() =
var pos = newSeq[int](modp5)
for i in 0 .. <rescnt: pos[residues[i]] = i-1
pos[1] = rescnt-1
# for each prime store resgroup on each restrack for prime*(modk+ri)
for j, prime in primes:
let k = uint((prime-2)) div uint(modp5) # find the resgroup it's
in
let r = uint((prime-2)) mod uint(modp5) + 2 # its base residue value
for ri in residues[1 .. rescnt]: # for each residue value
let prod: int = int(r) * ri # create residues
cross-product r*ri
let row: int = pos[prod mod modp5] * pcnt # create residue track
address
next[row + j] = k*(uint(prime) + uint(ri)) + uint(prod-2) div
uint(modp5)
proc segsieve(Kn: int) =
# for Kn resgroups in segment
for b in 0 .. <Kn: # for every byte in the segment
seg[b] = 0 # set every byte bit to prime (0)
for r in 0 .. <rescnt: # for each ith (of 8) residues for P5
let biti = uint8(1 shl r) # set the ith residue track bit mask
let row = r*pcnt # set address to ith row in next[]
for j in 0 .. <pcnt: # for each prime <= sqrt(N) for
restrack
if next[row+j] < uint(Kn): # if 1st mult resgroup index <= seg
size
var k: int = int(next[row+j]) # get its resgroup value
let prime = primes[j] # get the prime
while k < Kn: # for each primenth byte < segment
size
seg[k] = seg[k] or biti # set ith residue in byte as nonprime
k += int(prime)
next[row+j] = uint(k - Kn) # 1st resgroup in next eligible
segment
else: next[row+j] -= uint(Kn) # if 1st mult resgroup index > seg
size
# count the primes in the segment
for byt in seg[0..<Kn]: # for the Kn resgroup bytes
primecnt += uint(pbits[byt]) # count the '0' bits as primes
proc printprms(Kn: int, Ki: uint64) =
# Extract and print the primes in each segment:
# recommend piping output to less: ./ssozpxxx | less
# can then use Home, End, Pgup, Pgdn keys to see primes
for k in 0 .. <Kn: # for Kn residues groups|bytes
for r in 0 .. <8: # for each residue|bit in a byte
if (int(seg[k]) and (1 shl r)) == 0: # if it's '0' it's a prime
write(stdout, " ", uint64(modp5)*(Ki+uint64(k)) +
uint64(residues[r+1]))
#echo "\n"
proc ssozp5() =
stdout.write "Enter integer number: "
let val = stdin.readline.parseBiggestUInt # find primes <= val (7..2^64-1)
let B = 262144 # L2D_CACHE_SIZE 256*1024 bytes
let KB = B # number of segment resgroups
seg = newSeq[uint8](B) # create segment array of B bytes size
writeLine(stdout, "segment has ", B, " bytes and ", KB, " residues
groups")
let num = uint64((val-1) or 1) # if val even subtract 1
var k: uint64 = num div modp5 # resgroup of num
var modk = uint64(modp5) * k # base value of num
var r = 1 # starting with first residue
while num >= modk+uint64(residues[r]): r += 1 # find last pc position <=
num
let maxpcs = k * uint64(rescnt) + uint64(r-1) # maximum number of prime
candidates
let Kmax = uint64(num-2) div uint64(modp5) + 1 # maximum number of
resgroups for val
writeLine(stdout, "prime candidates = ", maxpcs, "; resgroups = ", Kmax)
let sqrtN = int(sqrt float64(num))
sozP7(sqrtN) # get pcnt and primes <= sqrt(nun)
writeLine(stdout, "create next[", rescnt, "x", pcnt, "] array")
next = newSeq[uint64](rescnt*pcnt) # create the next[] array
next_init() # load with first nonprimes 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 each
segment
echo "perform segmented SoZ\n";
while Ki < Kmax: # for KB size resgroup slices up to Kmax
if Kmax-Ki < uint64(KB): Kn = int(Kmax-Ki) # set sieve resgroups size
for last segment
segsieve(Kn) # sieve primes for current segment
#printprms(Kn,Ki) # print primes for the segment (optional)
Ki += uint64(KB) # next segment slice
var lprime = 0'u64 # get last prime and primecnt <= val
modk = uint64(modp5)*(Kmax-1) # mod for last resgroup in last segment
var b = Kn-1 # num bytes to last resgroup in segment
r = rescnt-1 # from last restrack in resgroup
while true: # repeat until last prime determined
if (int(seg[b]) and (1 shl r)) == 0:
lprime = modk+uint64(residues[r+1]) # determine the prime value
if lprime <= num: break # if <= num exit from loop with lprime
primecnt -= 1 # else reduce total primecnt
# reduce restrack, setup next resgroup
r -= 1; if r < 0: (r = rescnt-1; modk -= modp5; b -= 1) # if necessary
writeLine(stdout, "last segment = ", Kn, " resgroups; segments iterations
= ", Ki div uint64(KB))
writeLine(stdout, "last prime (of ", primecnt, ") is ", lprime)
ssozp5()
