Twinprimes generator, multi-threaded using bitarray, for SSoZ (Segmented Sieve of Zakiya), written in Crystal

twinprimes_ssoz_bitarray.cr

# This Crystal source file is a multiple threaded implementation to perform an
# extremely fast Segmented Sieve of Zakiya (SSoZ) to find Twin Primes <= N.

# Inputs are single values N, or ranges N1 and N2, of 64-bits, 0 -- 2^64 - 1.
# Output is the number of twin primes <= N, or in range N1 to N2; the last
# twin prime value for the range; and the total time of execution.

# This code was developed on a System76 laptop with an Intel I7 6700HQ cpu,
# 2.6-3.5 GHz clock, with 8 threads, and 16GB of memory. Parameter tuning
# probably needed to optimize for other hardware systems (ARM, PowerPC, etc).

# Compile as: $ crystal build twinprimes_ssozgistbary.cr -Dpreview_mt --release
# To reduce binary size do: $ strip twinprimes_ssoz
# Thread workers default to 4, set to system max for optimum performance.
# Single val: $ CRYSTAL_WORKERS=8 ./twinprimes_ssoz val1
# Range vals: $ CRYSTAL_WORKERS=8 ./twinprimes_ssoz val1 val2
# val1 and val2 can now be entered as either: 123456788 or 123_456_789

# Mathematical and technical basis for implementation are explained here:
# https://www.academia.edu/37952623/The_Use_of_Prime_Generators_to_Implement_Fast_
# Twin_Primes_Sieve_of_Zakiya_SoZ_Applications_to_Number_Theory_and_Implications_
# for_the_Riemann_Hypotheses
# https://www.academia.edu/7583194/The_Segmented_Sieve_of_Zakiya_SSoZ_
# https://www.academia.edu/19786419/PRIMES-UTILS_HANDBOOK
# https://www.academia.edu/81206391/Twin_Primes_Segmented_Sieve_of_Zakiya_SSoZ_Explained

# This source code, and its updates, can be found here:
# https://gist.github.com/jzakiya/2b65b609f091dcbb6f792f16c63a8ac4

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

# Copyright (c) 2017-2022; Jabari Zakiya -- jzakiya at gmail dot com
# Last update: 2023/02/10

require "bit_array"

# Customized gcd for prime generators; n > m; m odd
def gcd(m, n)
  while m|1 != 1; t = m; m = n % m; n = t end
  m
end

# Compute modular inverse a^-1 to base m, e.g. a*(a^-1) mod m = 1
def modinv(a0, m0)
  return 1 if m0 == 1
  a, m = a0, m0
  x0, inv = 0, 1
  while a > 1
    inv -= (a // m) * x0
    a, m = m, a % m
    x0, inv = inv, x0
  end
  inv += m0 if inv < 0
  inv
end

def gen_pg_parameters(prime)
  # Create prime generator parameters for given Pn
  puts "using Prime Generator parameters for P#{prime}"
  primes = [2, 3, 5, 7, 11, 13, 17, 19, 23]
  modpg, res_0 = 1, 0                    # compute Pn's modulus and res_0 value
  primes.each { |prm| res_0 = prm; break if prm > prime; modpg *= prm }

  restwins = [] of Int32                 # save upper twinpair residues here
  inverses = Array.new(modpg + 2, 0)     # save Pn's residues inverses here
  pc, inc, res = 5, 2, 0                 # use P3's PGS to generate pcs
  while pc < (modpg >> 1)                # find PG's 1st half residues
    if gcd(pc, modpg) == 1               # if pc a residue
      mc = modpg - pc                    # create its modular complement
      inverses[pc] = modinv(pc, modpg)   # save pc and mc inverses
      inverses[mc] = modinv(mc, modpg)   # if in twinpair save both hi residues
      restwins << pc << mc + 2 if res + 2 == pc
      res = pc                           # save current found residue
    end
    pc += inc; inc ^= 0b110              # create next P3 sequence pc: 5 7 11 13 17 19 ...
  end
  restwins.sort!;          restwins << (modpg + 1)         # last residue is last hi_tp
  inverses[modpg + 1] = 1; inverses[modpg - 1] = modpg - 1 # last 2 residues are self inverses
  {modpg, res_0, restwins.size, restwins, inverses}
end

def set_sieve_parameters(start_num, end_num)
  # Select at runtime best PG and segment size parameters for input values.
  # These are good estimates derived from PG data profiling. Can be improved.
  nrange = end_num - start_num
  bn, pg = 0, 3
  if end_num < 49
    bn = 1; pg = 3
  elsif nrange < 77_000_000
    bn = 16; pg = 5
  elsif nrange <  1_100_000_000
    bn = 32; pg = 7
  elsif nrange < 35_500_000_000
    bn = 64; pg = 11
  elsif nrange < 14_000_000_000_000
    pg = 13
    if    nrange > 7_000_000_000_000; bn = 384
    elsif nrange > 2_500_000_000_000; bn = 320
    elsif nrange >   250_000_000_000; bn = 196
    else  bn = 128
    end
  else
    bn = 384; pg = 17
  end
  modpg, res_0, pairscnt, restwins, resinvrs = gen_pg_parameters(pg)
  kmin = (start_num-2) // modpg + 1      # number of resgroups to start_num
  kmax = (end_num - 2) // modpg + 1      # number of resgroups to end_num
  krange = kmax - kmin + 1               # number of resgroups in range, at least 1
  n = krange < 37_500_000_000_000 ? 4 : (krange < 975_000_000_000_000 ? 6 : 8)
  b = bn * 1024 * n                      # set seg size to optimize for selected PG
  ks = krange < b ? krange : b           # segments resgroups size

  puts "segment size = #{ks} resgroups for seg bitarray"
  maxpairs = krange * pairscnt           # maximum number of twinprime pcs
  puts "twinprime candidates = #{maxpairs}; resgroups = #{krange}"
  {modpg, res_0, ks, kmin, kmax, krange, pairscnt, restwins, resinvrs}
end

def sozpg(val, res_0, start_num, end_num)
  # Compute the primes r0..sqrt(input_num) and store in 'primes' array.
  # Any algorithm (fast|small) is usable. Here the SoZ for P5 is used.
  md, rscnt = 30u64, 8                   # P5's modulus and residues count
  res  = [7,11,13,17,19,23,29,31]        # P5's residues
  bitn = [0,0,0,0,0,1,0,0,0,2,0,4,0,0,0,8,0,16,0,0,0,32,0,0,0,0,0,64,0,128]

  kmax = (val - 2) // md + 1             # number of resgroups upto input value
  prms = Array(UInt8).new(kmax, 0)       # byte array of prime candidates, init '0'
  modk, r, k = 0, -1, 0                  # initialize residue parameters

  loop do                                # for r0..sqrtN primes mark their multiples
    if (r += 1) == rscnt; r = 0; modk += md; k += 1 end # resgroup parameters
    next if prms[k] & (1 << r) != 0      # skip pc if not prime
    prm_r = res[r]                       # if prime save its residue value
    prime = modk + prm_r                 # numerate the prime value
    break if prime > Math.isqrt(val)     # exit loop when it's > sqrtN
    res.each do |ri|                     # mark prime's multiples in prms
      kn,rn = (prm_r * ri - 2).divmod md # cross-product resgroup|residue
      bit_r = bitn[rn]                   # bit mask for prod's residue
      kpm = k * (prime + ri) + kn        # resgroup for 1st prime mult
      while kpm < kmax; prms[kpm] |= bit_r; kpm += prime end
  end end
  # prms now contains the nonprime positions for the prime candidates r0..N
  # extract only primes that are in inputs range into array 'primes'
  primes = [] of UInt64                  # create empty dynamic array for primes
  prms.each_with_index do |resgroup, k|  # for each kth residue group
    res.each_with_index do |r_i, i|      # check for each ith residue in resgroup         
      if resgroup & (1 << i) == 0        # if bit location a prime
        prime = md * k + r_i             # numerate its value, store if in range
        # check if prime has multiple in range, if so keep it, if not don't
        n, rem = start_num.divmod prime  # if rem 0 then start_num is multiple of prime
        primes << prime if (res_0 <= prime <= val) && (prime * n <= end_num - prime || rem == 0)
  end end end
  primes
end

def nextp_init(rhi, kmin, modpg, primes, resinvrs)
  # Initialize 'nextp' array for twinpair upper residue rhi in 'restwins'.
  # Compute 1st prime multiple resgroups for each prime r0..sqrt(N) and
  # store consecutively as lo_tp|hi_tp pairs for their restracks.
  nextp = Slice(UInt64).new(primes.size*2) # 1st mults array for twinpair
  r_hi, r_lo = rhi, rhi - 2              # upper|lower twinpair residue values
  primes.each_with_index do |prime, j|   # for each prime r0..sqrt(N)
    k = (prime - 2) // modpg             # find the resgroup it's in
    r = (prime - 2) %  modpg + 2         # and its residue value
    r_inv = resinvrs[r].to_u64           # and residue inverse
    rl = (r_inv * r_lo - 2) % modpg + 2  # compute r's ri for r_lo
    rh = (r_inv * r_hi - 2) % modpg + 2  # compute r's ri for r_hi
    kl = k * (prime + rl) + (r * rl - 2) // modpg # kl 1st mult resgroup
    kh = k * (prime + rh) + (r * rh - 2) // modpg # kh 1st mult resgroup
    kl < kmin ? (kl = (kmin - kl) % prime; kl = prime - kl if kl > 0) : (kl -= kmin)
    kh < kmin ? (kh = (kmin - kh) % prime; kh = prime - kh if kh > 0) : (kh -= kmin)
    nextp[j * 2] = kl.to_u64             # prime's 1st mult lo_tp resgroup val in range 
    nextp[j * 2 | 1] = kh.to_u64         # prime's 1st mult hi_tp resgroup val in range
  end
  nextp
end

def twins_sieve(r_hi, kmin, kmax, ks, start_num, end_num, modpg, primes, resinvrs)
  # Perform in thread the ssoz for given twinpair residues for kmax resgroups.
  # First create|init 'nextp' array of 1st prime mults for given twinpair,
  # stored consequtively in 'nextp', and init seg array for ks resgroups.
  # For sieve, mark resgroup bits to '1' if either twinpair restrack is nonprime
  # for primes mults resgroups, and update 'nextp' restrack slices acccordingly.
  # Return the last twinprime|sum for the range for this twinpair residues.
  sum, ki, kn  = 0_u64, kmin-1, ks                     # init these parameters
  hi_tp, k_max = 0_u64, kmax                           # max twinprime|resgroup
  seg = BitArray.new(ks)                               # seg array of ks resgroups
  ki    += 1 if r_hi - 2 < (start_num - 2) % modpg + 2 # ensure lo tp in range
  k_max -= 1 if r_hi > (end_num - 2) % modpg + 2       # ensure hi tp in range
  nextp = nextp_init(r_hi, ki, modpg, primes,resinvrs) # init nextp array
  while ki < k_max                       # for ks size slices upto kmax
    kn = k_max - ki if ks > (k_max - ki) # adjust kn size for last seg
    primes.each_with_index do |prime, j| # for each prime r0..sqrt(N)
                                         # for lower twinpair residue track
      k = nextp.to_unsafe[j * 2]         # starting from this resgroup in seg
      while k < kn                       # until end of seg
        seg.unsafe_put(k, true)          # mark primenth resgroup bits prime mults
        k += prime  end                  # set resgroup for prime's next multiple
      nextp.to_unsafe[j * 2] = k - kn    # save 1st resgroup in next eligible seg
                                         # for upper twinpair residue track
      k = nextp.to_unsafe[j * 2 | 1]     # starting from this resgroup in seg
      while k < kn                       # until end of seg
        seg.unsafe_put(k, true)          # mark primenth resgroup bits prime mults
        k += prime  end                  # set resgroup for prime's next multiple
      nextp.to_unsafe[j * 2 | 1]= k - kn # save 1st resgroup in next eligible seg
    end                                   
    cnt = seg[...kn].count(false)        # count|store twinprimes in segment
    if cnt > 0                           # if segment has twinprimes
      sum += cnt                         # add segment count to total range count
      upk = kn - 1                       # from end of seg, count back to largest tp
      while seg.unsafe_fetch(upk); upk -= 1 end
      hi_tp = ki + upk                   # set its full range resgroup value
    end
    ki += ks                             # set 1st resgroup val of next seg slice
    seg.fill(false) if ki < k_max        # set next seg to all primes if in range
  end                                    # when sieve done, numerate largest twin 
                                         # for ranges w/o twins set largest to 1
  hi_tp = (r_hi > end_num || sum == 0) ? 1 : hi_tp * modpg + r_hi
  {hi_tp.to_u64, sum.to_u64}             # return largest twinprime|twins count
end

def twinprimes_ssoz() 
  end_num   = {(ARGV[0].to_u64 underscore: true), 3u64}.max
  start_num = ARGV.size > 1 ? {(ARGV[1].to_u64 underscore: true), 3u64}.max : 3u64
  start_num, end_num = end_num, start_num if start_num > end_num
  start_num |= 1                         # if start_num even increase by 1
  end_num = (end_num - 1) | 1            # if end_num even decrease by 1
  start_num = end_num = 7 if end_num - start_num < 2

  puts "threads = #{System.cpu_count}"
  ts = Time.monotonic                    # start timing sieve setup execution
                                         # select Pn, set sieving params for inputs
  modpg, res_0, ks, kmin, kmax, krange, pairscnt, restwins, resinvrs = set_sieve_parameters(start_num, end_num)

  # create sieve primes <= sqrt(end_num), only use those whose multiples within inputs range
  primes = end_num < 49 ? [5] : sozpg(Math.isqrt(end_num), res_0, start_num, end_num)

  puts "each of #{pairscnt} threads has nextp[2 x #{primes.size}] array"

  lo_range = restwins[0] - 3             # lo_range = lo_tp - 1
  twinscnt = 0_u64                       # determine count of 1st 4 twins if in range for used Pn
  twinscnt += [3, 5, 11, 17].select { |tp| start_num <= tp <= lo_range }.size unless end_num == 3

  te = (Time.monotonic - ts).total_seconds.round(6)
  puts "setup time = #{te} secs"         # display sieve setup time
  puts "perform twinprimes ssoz sieve"
  t1 = Time.monotonic                    # start timing ssoz sieve execution

  cnts = Array(UInt64).new(pairscnt, 0)  # number of twinprimes found per thread
  lastwins = Array(UInt64).new(pairscnt, 0) # largest twinprime val for each thread
  done = Channel(Nil).new(pairscnt)
  
  threadscnt = Atomic.new(0)             # count of finished threads
  restwins.each_with_index do |r_hi, i|  # sieve twinpair restracks
    spawn do
      lastwins[i], cnts[i] = twins_sieve(r_hi, kmin, kmax, ks, start_num, end_num, modpg, primes, resinvrs)
      print "\r#{threadscnt.add(1)} of #{pairscnt} twinpairs done"
      done.send(nil)
  end end
  pairscnt.times { done.receive }        # wait for all threads to finish
  print "\r#{pairscnt} of #{pairscnt} twinpairs done"

  last_twin = lastwins.max               # find largest hi_tp twinprime in range
  twinscnt += cnts.sum                   # compute number of twinprimes in range
  last_twin = 5 if end_num == 5 && twinscnt == 1
  kn = krange % ks                       # set number of resgroups in last slice
  kn = ks if kn == 0                     # if multiple of seg size set to seg size
  t2 = (Time.monotonic - t1).total_seconds       # sieve execution time

  puts "\nsieve time = #{t2.round(6)} secs"      # ssoz sieve time
  puts "total time = #{(t2 + te).round(6)} secs" # setup + sieve time
  puts "last segment = #{kn} resgroups; segment slices = #{(krange - 1)//ks + 1}"
  puts "total twins = #{twinscnt}; last twin = #{last_twin - 1}+/-1"
end

twinprimes_ssoz