jzakiya/twinprimes_ssoz.go

## twinprimes_ssoz.go
// This Go 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.

// Code originally 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: $ go build twinprimes_ssoz.go
// To reduce binary size do: $ strip twinprimes_ssoz
// Single val: $ echo val | ./twinprimes_ssoz
// Range vals: $ echo val1 val2 | ./twinprimes_ssoz

// 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/fbc77b8fdd12b0581a0ff7c2476373d9

// 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-2023; Jabari Zakiya -- jzakiya at gmail dot com
// Last update: 2023/08/15

package main

import (
  "fmt"
  "time"
  "sort"
  "sync"
  "math"
  "math/bits"
  "runtime"
  "sync/atomic"
)

// Customized gcd for prime generators; n > m; m odd
func gcd(m, n int) int {
  for m | 1 != 1 { m, n = n % m, m }
  return m
}

// Compute modular inverse a^-1 to base m, e.g. a*(a^-1) mod m = 1
func modinv(a0, m0 int) int {
  if m0 == 1 { return 1 }
  a, m := a0, m0
  x0, inv := 0, 1
  for a > 1 {
    inv -= (a / m) * x0
    a, m = m, a % m
    x0, inv = inv, x0
  }
  if inv < 0 { inv += m0 }
  return inv
}

func gen_pg_parameters(prime int) (uint, uint, int, []int, []int) {
  // Create prime generator parameters for given Pn
  fmt.Printf("using Prime Generator parameters for P%d\n", prime)
  primes := []int{2, 3, 5, 7, 11, 13, 17, 19, 23}
  modpg, res_0 := 1, 0                   // compute Pn's modulus and res_0 value
  for _, prm := range primes { res_0 = prm; if prm > prime { break }; modpg *= prm }

  restwins := []int{}                    // save upper twinpair residues here
  inverses := make([]int, modpg + 2)     // save Pn's residues inverses here
  pc, inc, res := 5, 2, 0                // use P3's PGS to generate pcs
  for pc < (modpg >> 1) {                // find PG's 1st half residues
    if gcd(pc, modpg) == 1 {             // if pc a residue
      var mc = modpg - pc                // create its modular complement
      inverses[pc] = modinv(pc, modpg)   // save pc and mc inverses
      inverses[mc] = modinv(mc, modpg)   // if in twinspair save both hi residues
      if res + 2 == pc { restwins = append(restwins, pc, mc + 2) }
      res = pc                           // save current found residue
    }
    pc += inc; inc ^= 0b110              // create next P3 sequence pc: 5 7 11 13 17 19 ...
  }
  sort.Ints(restwins); restwins = append(restwins, modpg + 1) // last residue is last hi_tp
  inverses[modpg + 1] = 1; inverses[modpg - 1] = modpg - 1    // last 2 residues are self inverses
  return uint(modpg), uint(res_0), len(restwins) , restwins, inverses
}

func set_sieve_parameters(start_num, end_num uint) (uint, uint, uint, uint, uint, uint, int, []int, []int) {
  // Select at runtime best PG and segment size parameters for input values.
  // These are good estimates derived from PG data profiling. Can be improved.
  trange := end_num - start_num
  var bn, pg = 0, 3
  if end_num < 49 {
    bn, pg = 1, 3
  } else if trange < 77_000_000 {
    bn, pg = 16, 5
  } else if trange <  1_100_000_000 {
    bn, pg = 32, 7
  } else if trange < 35_500_000_000 {
    bn, pg = 64, 11
  } else if trange < 14_000_000_000_000 {
    pg = 13
    if        trange > 7_000_000_000_000 { bn = 384
    } else if trange > 2_500_000_000_000 { bn = 320
    } else if trange >   250_000_000_000 { bn = 196
    } else { bn = 128 }
  } else {
    bn, pg = 384, 17
  }
  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
  var ks, n = uint(0), 0
  if krange < 37_500_000_000_000 { n = 4 } else if krange < 975_000_000_000_000 { n = 6 } else { n = 8}
  b := uint(bn * 1024 * n)               // set seg size to optimize for selected PG
  if krange < b { ks = krange } else { ks = b } // segments resgroups size

  fmt.Printf("segment size = %d resgroups; seg array is [1 x %d] 64-bits\n", ks, (((ks-1) >> 6) + 1))
  maxpairs := krange * uint(pairscnt)    // maximum number of twinprime pcs
  fmt.Printf("twinprime candidates = %d; resgroups = %d\n", maxpairs, krange)
  return modpg, res_0, ks, kmin, kmax, krange, pairscnt, restwins, resinvrs
}

func sozp5(val, res_0, start_num, end_num uint) []uint {
  // Return the primes r0..sqrt(end_num) within range (start_num...end_num)
  md, rescnt := 30, 8                    // P5's modulus and residues count
  res  := []int{7,11,13,17,19,23,29,31}  // P5's residues
  bitn := []int{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) / uint(md) + 1       // number of resgroups upto input value
  prms := make([]uint8, kmax)            // byte array of prime candidates, init '0'
  sqrt_n := int(math.Sqrt(float64(val))) // compute integer sqrt of val
  modk, r, k := 0, -1, 0                 // initialize resgroup parameters

  for true {                             // for r0..sqrtN primes mark their multiples
    r++; if r == rescnt { r = 0; modk += md; k += 1 } // update resgroup parameters
    if (prms[k] & (1 << r)) != 0 { continue }         // skip pc if not prime
    var prm_r = res[r]                   // if prime save its residue value
    var prime = modk + prm_r             // numerate the prime value
    if prime > sqrt_n { break }          // stop sieving when true
    for _, ri := range res {             // mark prime's multiples in prms
      prod  := prm_r * ri - 2            // compute cross-product for prm_r|ri pair
      bit_r := uint8(bitn[prod % md])         // bit mask value for prod's residue
      var kpm = k * (prime + ri) + prod / md  // resgroup for 1st prime multiple
      for uint(kpm) < kmax { prms[kpm] |= bit_r; kpm += prime }
  } }
  // prms now contains the prime multiples positions marked for the pcs r0..N
  // now along each restrack, identify|numerate the primes in each resgroup
  // return only the primes with a multiple within range (start_num...end_num)
  primes := []uint{}                     // create empty dynamic array for primes
  for i, r_i := range res {              // along each restrack|row til end
     for k, resgroup := range prms  {    // for each resgroup along restrack
	if resgroup & (1 << i) == 0 {    // if bit location a prime
	  prime := uint(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
          rem := start_num % prime
  	  if (res_0 <= prime && prime <= val) && (prime-rem <= end_num-start_num || rem == 0) {primes = append(primes, prime)}
  } } }
  return primes                          // primes stored in restrack|row sequence order
}

func nextp_init(rhi, kmin, modpg uint, primes []uint, resinvrs []int) []uint {
  // 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 := make([]uint,len(primes)*2)    // 1st mults array for twinpair
  r_hi, r_lo := rhi, rhi - 2             // upper|lower twinpair residue values
  for j, prime := range primes {         // 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 := uint(resinvrs[r])           // and residue inverse
    rl := (r_lo * r_inv - 2) % modpg + 2 // compute r's ri for r_lo
    rh := (r_hi * r_inv - 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
    if kl < kmin { kl = (kmin - kl) % prime; if kl > 0 { kl = prime - kl } } else { kl -= kmin }
    if kh < kmin { kh = (kmin - kh) % prime; if kh > 0 { kh = prime - kh } } else { kh -= kmin }
    nextp[j * 2] = uint(kl)              // prime's 1st mult lo_tp resgroup val in range
    nextp[j * 2 | 1] = uint(kh)          // prime's 1st mult hi_tp resgroup val in range
  }
  return nextp
}

func twins_sieve(rhi int, kmin, kmax, ks, start_num, end_num, modpg uint, primes []uint, resinvrs []int) (uint, uint) {
  // 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.
  s := 6                                               // shift value for 64 bits
  bmask := uint((1 << s) - 1)                          // bitmask val for 64 bits
  sum, ki, kn, r_hi := 0, kmin-1, ks, uint(rhi)        // init these parameters
  hi_tp, k_max := uint(0), kmax                        // max twinprime|resgroup
  seg := make([]uint64, ((ks - 1) >> s) + 1)           // seg array of ks resgroups
  if r_hi - 2 < (start_num - 2) % modpg + 2 {ki += 1}  // ensure lo tp in range
  if r_hi > (end_num - 2) % modpg + 2 {k_max -= 1}     // ensure hi tp in range
  nextp := nextp_init(r_hi, ki, modpg, primes, resinvrs) // init nextp array
  for ki < k_max {                       // sieve ks size slices upto kmax
    if ks > k_max - ki {kn = k_max - ki} // adjust kn size for last seg
    for j, prime := range primes {       // for each prime r0..sqrt(N)
                                         // for lower residue track
      var k1 uint = nextp[j * 2]         // starting from this resgroup in seg
      for k1 < kn   {                    // mark primenth resgroup bits prime mults
        seg[k1 >> s] |= uint64(1) << (k1 & bmask)
        k1 += prime }                    // set resgroup for prime's next multiple
      nextp[j * 2] = k1 - kn             // save 1st resgroup in next eligible seg
                                         // for upper residue track
      var k2 uint = nextp[j * 2 | 1]     // starting from this resgroup in seg
      for k2 < kn   {                    // mark primenth resgroup bits prime mults
        seg[k2 >> s] |= uint64(1) << (k2 & bmask)
        k2 += prime }                    // set resgroup for prime's next multiple
      nextp[j * 2 | 1] = k2 - kn         // save 1st resgroup in next eligible seg
    }                                    // set as nonprime unused bits in last seg[n]
                                         // so fast, do for every seg[i]
    seg[(kn - 1) >> s] |= ^uint64(1) << ((kn - 1) & bmask)
    cnt := 0                             // count the twinprimes in the segment
    for i := 0; i <= int(kn - 1) >> s; i++ { cnt += bits.OnesCount64(^seg[i]) }
    if cnt > 0 {                         // if segment has twinprimes
      sum += cnt                         // add the segment count to total count
      var upk = kn - 1                   // from end of seg count back to largest tp
      for seg[upk >> s] & (uint64(1) << (upk & bmask)) != 0 { upk -= 1 }
      hi_tp = upk + ki                   // set its full range resgroup value
    }
    ki += ks                             // set 1st resgroup val of next seg slice
    if ki < k_max { for i := range seg { seg[i] = 0 } } // set seg bits to all primes
  }                                      // when sieve done, numerate largest twin in range
                                         // for small ranges w/o twins set largest to 1
  if r_hi > end_num || sum == 0 { hi_tp = 1 } else { hi_tp = hi_tp * modpg + r_hi }
  return hi_tp, uint(sum)                // return largest twinprime|twins count in range
}

func main() {
  var start_num, end_num = uint(3), uint(3)
  _, err := fmt.Scan(&end_num, &start_num)
  if end_num < 3 { end_num = 3 }
  if err != nil { start_num = 3 } else { if start_num < 3 { start_num = 3 } }
  if start_num > end_num { start_num, end_num = end_num, start_num }
  start_num |= 1                         // if start_num even increase by 1
  end_num = (end_num - 1) | 1            // if end_num even decrease by 1
  if end_num - start_num < 2 { start_num, end_num = 7, 7 }

  fmt.Printf("threads = %d\n", runtime.NumCPU())
  ts := time.Now()                       // 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
  var primes = []uint{}
  if end_num < 49 { primes = []uint{5}
  } else { primes = sozp5(uint(math.Sqrt(float64(end_num))), res_0, start_num, end_num) }

  fmt.Printf("each of %d threads has nextp[2 x %d] array\n", pairscnt, len(primes))

  twinscnt := uint(0)                    // init twinprimes range count
  lo_range := uint(restwins[0] - 3)      // lo_range = lo_tp - 1
  for _, tp := range []uint{3, 5, 11, 17} { // excluded low tp values PGs used
    if end_num == 3 { break }            // if 3 end of range, no twinprimes
    if start_num <= tp && tp < lo_range { twinscnt += 1 } // cnt any small tps
  }
  te := time.Now()
  fmt.Printf("setup time = %v \n", te.Sub(ts)) // display sieve setup time
  fmt.Println("perform twinprimes ssoz sieve")
  t1 := time.Now()                       // start timing ssoz sieve execution

  cnts := make([]uint, pairscnt)         // number of twinprimes found per thread
  lastwins := make([]uint, pairscnt)     // largest twinprime val for each thread
  var threadscnt uint64;                 // count of finished threads

  var wg sync.WaitGroup
  for i, r_hi := range restwins {        // sieve twinpair restracks
    wg.Add(1)
    go func(i, r_hi int) {
      defer wg.Done()
      l, c := twins_sieve(r_hi, kmin, kmax, ks, start_num, end_num, modpg, primes, resinvrs)
      lastwins[i] = l; cnts[i] = c
      fmt.Printf("\r%d of %d twinpairs done", atomic.AddUint64(&threadscnt, 1), pairscnt)
    }(i, r_hi)
  }
  wg.Wait()
  fmt.Printf("\r%d of %d twinpairs done", pairscnt, pairscnt)

  last_twin := uint(0)                   // init val for largest twinprime
  for i := 0; i < pairscnt; i++ {        // find largest twinprime|sum in range
    twinscnt += cnts[i]
    if last_twin < lastwins[i] { last_twin = lastwins[i] }
  }
  if end_num == 5 && twinscnt == 1 { last_twin = 5 }
  kn := krange % ks                      // set number of resgroups in last slice
  if kn == 0 { kn = ks }                 // if multiple of seg size set to seg size
  t2 := time.Now()                       // sieve execution time

  fmt.Printf("\nsieve time = %v", t2.Sub(t1)) // ssoz sieve time
  fmt.Printf("\ntotal time = %v", t2.Sub(ts)) // setup + sieve time
  fmt.Printf("\nlast segment = %d resgroups; segment slices = %d\n", kn, (krange - 1)/ks + 1)
  fmt.Printf("total twins = %d; last twin = %d+/-1", twinscnt, last_twin - 1)
}
	// This Go 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.

	// Code originally 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: $ go build twinprimes_ssoz.go
	// To reduce binary size do: $ strip twinprimes_ssoz
	// Single val: $ echo val \| ./twinprimes_ssoz
	// Range vals: $ echo val1 val2 \| ./twinprimes_ssoz

	// 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/fbc77b8fdd12b0581a0ff7c2476373d9

	// 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-2023; Jabari Zakiya -- jzakiya at gmail dot com
	// Last update: 2023/08/15

	package main

	import (
	"fmt"
	"time"
	"sort"
	"sync"
	"math"
	"math/bits"
	"runtime"
	"sync/atomic"
	)

	// Customized gcd for prime generators; n > m; m odd
	func gcd(m, n int) int {
	for m \| 1 != 1 { m, n = n % m, m }
	return m
	}

	// Compute modular inverse a^-1 to base m, e.g. a*(a^-1) mod m = 1
	func modinv(a0, m0 int) int {
	if m0 == 1 { return 1 }
	a, m := a0, m0
	x0, inv := 0, 1
	for a > 1 {
	inv -= (a / m) * x0
	a, m = m, a % m
	x0, inv = inv, x0
	}
	if inv < 0 { inv += m0 }
	return inv
	}

	func gen_pg_parameters(prime int) (uint, uint, int, []int, []int) {
	// Create prime generator parameters for given Pn
	fmt.Printf("using Prime Generator parameters for P%d\n", prime)
	primes := []int{2, 3, 5, 7, 11, 13, 17, 19, 23}
	modpg, res_0 := 1, 0 // compute Pn's modulus and res_0 value
	for _, prm := range primes { res_0 = prm; if prm > prime { break }; modpg *= prm }

	restwins := []int{} // save upper twinpair residues here
	inverses := make([]int, modpg + 2) // save Pn's residues inverses here
	pc, inc, res := 5, 2, 0 // use P3's PGS to generate pcs
	for pc < (modpg >> 1) { // find PG's 1st half residues
	if gcd(pc, modpg) == 1 { // if pc a residue
	var mc = modpg - pc // create its modular complement
	inverses[pc] = modinv(pc, modpg) // save pc and mc inverses
	inverses[mc] = modinv(mc, modpg) // if in twinspair save both hi residues
	if res + 2 == pc { restwins = append(restwins, pc, mc + 2) }
	res = pc // save current found residue
	}
	pc += inc; inc ^= 0b110 // create next P3 sequence pc: 5 7 11 13 17 19 ...
	}
	sort.Ints(restwins); restwins = append(restwins, modpg + 1) // last residue is last hi_tp
	inverses[modpg + 1] = 1; inverses[modpg - 1] = modpg - 1 // last 2 residues are self inverses
	return uint(modpg), uint(res_0), len(restwins) , restwins, inverses
	}

	func set_sieve_parameters(start_num, end_num uint) (uint, uint, uint, uint, uint, uint, int, []int, []int) {
	// Select at runtime best PG and segment size parameters for input values.
	// These are good estimates derived from PG data profiling. Can be improved.
	trange := end_num - start_num
	var bn, pg = 0, 3
	if end_num < 49 {
	bn, pg = 1, 3
	} else if trange < 77_000_000 {
	bn, pg = 16, 5
	} else if trange < 1_100_000_000 {
	bn, pg = 32, 7
	} else if trange < 35_500_000_000 {
	bn, pg = 64, 11
	} else if trange < 14_000_000_000_000 {
	pg = 13
	if trange > 7_000_000_000_000 { bn = 384
	} else if trange > 2_500_000_000_000 { bn = 320
	} else if trange > 250_000_000_000 { bn = 196
	} else { bn = 128 }
	} else {
	bn, pg = 384, 17
	}
	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
	var ks, n = uint(0), 0
	if krange < 37_500_000_000_000 { n = 4 } else if krange < 975_000_000_000_000 { n = 6 } else { n = 8}
	b := uint(bn * 1024 * n) // set seg size to optimize for selected PG
	if krange < b { ks = krange } else { ks = b } // segments resgroups size

	fmt.Printf("segment size = %d resgroups; seg array is [1 x %d] 64-bits\n", ks, (((ks-1) >> 6) + 1))
	maxpairs := krange * uint(pairscnt) // maximum number of twinprime pcs
	fmt.Printf("twinprime candidates = %d; resgroups = %d\n", maxpairs, krange)
	return modpg, res_0, ks, kmin, kmax, krange, pairscnt, restwins, resinvrs
	}

	func sozp5(val, res_0, start_num, end_num uint) []uint {
	// Return the primes r0..sqrt(end_num) within range (start_num...end_num)
	md, rescnt := 30, 8 // P5's modulus and residues count
	res := []int{7,11,13,17,19,23,29,31} // P5's residues
	bitn := []int{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) / uint(md) + 1 // number of resgroups upto input value
	prms := make([]uint8, kmax) // byte array of prime candidates, init '0'
	sqrt_n := int(math.Sqrt(float64(val))) // compute integer sqrt of val
	modk, r, k := 0, -1, 0 // initialize resgroup parameters

	for true { // for r0..sqrtN primes mark their multiples
	r++; if r == rescnt { r = 0; modk += md; k += 1 } // update resgroup parameters
	if (prms[k] & (1 << r)) != 0 { continue } // skip pc if not prime
	var prm_r = res[r] // if prime save its residue value
	var prime = modk + prm_r // numerate the prime value
	if prime > sqrt_n { break } // stop sieving when true
	for _, ri := range res { // mark prime's multiples in prms
	prod := prm_r * ri - 2 // compute cross-product for prm_r\|ri pair
	bit_r := uint8(bitn[prod % md]) // bit mask value for prod's residue
	var kpm = k * (prime + ri) + prod / md // resgroup for 1st prime multiple
	for uint(kpm) < kmax { prms[kpm] \|= bit_r; kpm += prime }
	} }
	// prms now contains the prime multiples positions marked for the pcs r0..N
	// now along each restrack, identify\|numerate the primes in each resgroup
	// return only the primes with a multiple within range (start_num...end_num)
	primes := []uint{} // create empty dynamic array for primes
	for i, r_i := range res { // along each restrack\|row til end
	for k, resgroup := range prms { // for each resgroup along restrack
	if resgroup & (1 << i) == 0 { // if bit location a prime
	prime := uint(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
	rem := start_num % prime
	if (res_0 <= prime && prime <= val) && (prime-rem <= end_num-start_num \|\| rem == 0) {primes = append(primes, prime)}
	} } }
	return primes // primes stored in restrack\|row sequence order
	}

	func nextp_init(rhi, kmin, modpg uint, primes []uint, resinvrs []int) []uint {
	// 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 := make([]uint,len(primes)*2) // 1st mults array for twinpair
	r_hi, r_lo := rhi, rhi - 2 // upper\|lower twinpair residue values
	for j, prime := range primes { // 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 := uint(resinvrs[r]) // and residue inverse
	rl := (r_lo * r_inv - 2) % modpg + 2 // compute r's ri for r_lo
	rh := (r_hi * r_inv - 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
	if kl < kmin { kl = (kmin - kl) % prime; if kl > 0 { kl = prime - kl } } else { kl -= kmin }
	if kh < kmin { kh = (kmin - kh) % prime; if kh > 0 { kh = prime - kh } } else { kh -= kmin }
	nextp[j * 2] = uint(kl) // prime's 1st mult lo_tp resgroup val in range
	nextp[j * 2 \| 1] = uint(kh) // prime's 1st mult hi_tp resgroup val in range
	}
	return nextp
	}

	func twins_sieve(rhi int, kmin, kmax, ks, start_num, end_num, modpg uint, primes []uint, resinvrs []int) (uint, uint) {
	// 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.
	s := 6 // shift value for 64 bits
	bmask := uint((1 << s) - 1) // bitmask val for 64 bits
	sum, ki, kn, r_hi := 0, kmin-1, ks, uint(rhi) // init these parameters
	hi_tp, k_max := uint(0), kmax // max twinprime\|resgroup
	seg := make([]uint64, ((ks - 1) >> s) + 1) // seg array of ks resgroups
	if r_hi - 2 < (start_num - 2) % modpg + 2 {ki += 1} // ensure lo tp in range
	if r_hi > (end_num - 2) % modpg + 2 {k_max -= 1} // ensure hi tp in range
	nextp := nextp_init(r_hi, ki, modpg, primes, resinvrs) // init nextp array
	for ki < k_max { // sieve ks size slices upto kmax
	if ks > k_max - ki {kn = k_max - ki} // adjust kn size for last seg
	for j, prime := range primes { // for each prime r0..sqrt(N)
	// for lower residue track
	var k1 uint = nextp[j * 2] // starting from this resgroup in seg
	for k1 < kn { // mark primenth resgroup bits prime mults
	seg[k1 >> s] \|= uint64(1) << (k1 & bmask)
	k1 += prime } // set resgroup for prime's next multiple
	nextp[j * 2] = k1 - kn // save 1st resgroup in next eligible seg
	// for upper residue track
	var k2 uint = nextp[j * 2 \| 1] // starting from this resgroup in seg
	for k2 < kn { // mark primenth resgroup bits prime mults
	seg[k2 >> s] \|= uint64(1) << (k2 & bmask)
	k2 += prime } // set resgroup for prime's next multiple
	nextp[j * 2 \| 1] = k2 - kn // save 1st resgroup in next eligible seg
	} // set as nonprime unused bits in last seg[n]
	// so fast, do for every seg[i]
	seg[(kn - 1) >> s] \|= ^uint64(1) << ((kn - 1) & bmask)
	cnt := 0 // count the twinprimes in the segment
	for i := 0; i <= int(kn - 1) >> s; i++ { cnt += bits.OnesCount64(^seg[i]) }
	if cnt > 0 { // if segment has twinprimes
	sum += cnt // add the segment count to total count
	var upk = kn - 1 // from end of seg count back to largest tp
	for seg[upk >> s] & (uint64(1) << (upk & bmask)) != 0 { upk -= 1 }
	hi_tp = upk + ki // set its full range resgroup value
	}
	ki += ks // set 1st resgroup val of next seg slice
	if ki < k_max { for i := range seg { seg[i] = 0 } } // set seg bits to all primes
	} // when sieve done, numerate largest twin in range
	// for small ranges w/o twins set largest to 1
	if r_hi > end_num \|\| sum == 0 { hi_tp = 1 } else { hi_tp = hi_tp * modpg + r_hi }
	return hi_tp, uint(sum) // return largest twinprime\|twins count in range
	}

	func main() {
	var start_num, end_num = uint(3), uint(3)
	_, err := fmt.Scan(&end_num, &start_num)
	if end_num < 3 { end_num = 3 }
	if err != nil { start_num = 3 } else { if start_num < 3 { start_num = 3 } }
	if start_num > end_num { start_num, end_num = end_num, start_num }
	start_num \|= 1 // if start_num even increase by 1
	end_num = (end_num - 1) \| 1 // if end_num even decrease by 1
	if end_num - start_num < 2 { start_num, end_num = 7, 7 }

	fmt.Printf("threads = %d\n", runtime.NumCPU())
	ts := time.Now() // 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
	var primes = []uint{}
	if end_num < 49 { primes = []uint{5}
	} else { primes = sozp5(uint(math.Sqrt(float64(end_num))), res_0, start_num, end_num) }

	fmt.Printf("each of %d threads has nextp[2 x %d] array\n", pairscnt, len(primes))

	twinscnt := uint(0) // init twinprimes range count
	lo_range := uint(restwins[0] - 3) // lo_range = lo_tp - 1
	for _, tp := range []uint{3, 5, 11, 17} { // excluded low tp values PGs used
	if end_num == 3 { break } // if 3 end of range, no twinprimes
	if start_num <= tp && tp < lo_range { twinscnt += 1 } // cnt any small tps
	}
	te := time.Now()
	fmt.Printf("setup time = %v \n", te.Sub(ts)) // display sieve setup time
	fmt.Println("perform twinprimes ssoz sieve")
	t1 := time.Now() // start timing ssoz sieve execution

	cnts := make([]uint, pairscnt) // number of twinprimes found per thread
	lastwins := make([]uint, pairscnt) // largest twinprime val for each thread
	var threadscnt uint64; // count of finished threads

	var wg sync.WaitGroup
	for i, r_hi := range restwins { // sieve twinpair restracks
	wg.Add(1)
	go func(i, r_hi int) {
	defer wg.Done()
	l, c := twins_sieve(r_hi, kmin, kmax, ks, start_num, end_num, modpg, primes, resinvrs)
	lastwins[i] = l; cnts[i] = c
	fmt.Printf("\r%d of %d twinpairs done", atomic.AddUint64(&threadscnt, 1), pairscnt)
	}(i, r_hi)
	}
	wg.Wait()
	fmt.Printf("\r%d of %d twinpairs done", pairscnt, pairscnt)

	last_twin := uint(0) // init val for largest twinprime
	for i := 0; i < pairscnt; i++ { // find largest twinprime\|sum in range
	twinscnt += cnts[i]
	if last_twin < lastwins[i] { last_twin = lastwins[i] }
	}
	if end_num == 5 && twinscnt == 1 { last_twin = 5 }
	kn := krange % ks // set number of resgroups in last slice
	if kn == 0 { kn = ks } // if multiple of seg size set to seg size
	t2 := time.Now() // sieve execution time

	fmt.Printf("\nsieve time = %v", t2.Sub(t1)) // ssoz sieve time
	fmt.Printf("\ntotal time = %v", t2.Sub(ts)) // setup + sieve time
	fmt.Printf("\nlast segment = %d resgroups; segment slices = %d\n", kn, (krange - 1)/ks + 1)
	fmt.Printf("total twins = %d; last twin = %d+/-1", twinscnt, last_twin - 1)
	}