Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Twinprimes generator, multi-threaded, using SSoZ (Segmented Sieve of Zakiya), written in C++.
/*
This C++ 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 hadware systems (ARM, PowerPC, etc).
On Linux use -O compiler flag to compile for optimum speed:
$ g++ -O3 -fopenmp twinprimes_ssoz.cpp -o twinprimes_ssoz
To reduce binary size do: $ strip twinprimes_ssoz
Single val: $ echo val1 | ./twinprimes_ssoz
Range vals: $ echo val1 val2 | ./twinprimes_ssoz
Range values can be provided in either order (larger or smaller first).
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
This C++ source code, and updates, can be found here:
https://gist.github.com/jzakiya/fa76c664c9072ddb51599983be175a3f
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/
Use of this code is free subject to acknowledgment of copyright.
Copyright (c) 2017-2022 Jabari Zakiya -- jzakiya at gmail dot com
Last update: 2022/05/22
*/
#include <cmath>
#include <algorithm>
#include <vector>
#include <array>
#include <tuple>
#include <cstdlib>
#include <iostream>
#include <stdint.h>
#include <chrono>
#include <omp.h>
#include <atomic>
using namespace std;
typedef uint64_t uint64;
// Compute modular inverse a^(-1) of a to base b, e.g. a*(a^-1) mod b = 1.
int modinv(int a0, int b0) {
int b = b0, a = a0, x = 0;
int inv = 1;
if (b == 1) return 1;
while (a > 1) {
int q = a / b;
int t = b; b = a % b; a = t;
t = x; x = inv - q * x; inv = t;
}
if (inv < 0) inv += b0;
return inv;
}
// Create prime generator parameters for given Pn
tuple<int, int, int, vector<int>, vector<int> > genPgParameters(int prime) {
cout << "using Prime Generator parameters for P" << prime << endl;
int primes[9] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
int modpg = 1; int res_0 = 0;
for (int prm : primes) { res_0 = prm; if (prm > prime) break; modpg *= prm; }
vector<int> restwins; // save upper twin pair residues here
vector<int> inverses(modpg + 2, 0); // save Pn's residues inverses here
int pc = 5; int inc = 2; int res = 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
int 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
if (res + 2 == pc) { restwins.push_back(pc); restwins.push_back(mc + 2); }
res = pc; // save current found residue
}
pc += inc; inc ^= 0b110; // create next P3 sequence pc: 5 7 11 13 17 19 ...
}
std::sort(restwins.begin(), restwins.end()); restwins.push_back(modpg + 1);
inverses[modpg + 1] = 1; inverses[modpg - 1] = modpg - 1;
return make_tuple(modpg, res_0, restwins.size(), restwins, inverses);
}
// Select at runtime best PG and segment size factor to use for input vals.
// These are good estimates derived from PG data profiling. Can be improved.
tuple<int, int, int, uint64, uint64, uint64, int, vector<int>,
vector<int>> setSieveParameters(uint64 start_num, uint64 end_num) {
uint64 nrange = end_num - start_num;
int Bn; int pg = 3;
if (end_num < 49) {
Bn = 1; pg = 3;
}
else if (nrange < 77000000) {
Bn = 16; pg = 5;
}
else if (nrange < 1100000000) {
Bn = 32; pg = 7;
}
else if (nrange < 35500000000) {
Bn = 64; pg = 11;
}
else if (nrange < 14000000000000) {
pg = 13;
if (nrange > 7000000000000) Bn = 384;
else if (nrange > 2500000000000) Bn = 320;
else if (nrange > 250000000000) Bn = 196;
else Bn = 128;
} else {
Bn = 384; pg = 17;
}
int modpg, res_0, pairscnt; vector<int> restwins, resinvrs;
tie(modpg, res_0, pairscnt, restwins, resinvrs) = genPgParameters(pg);
uint64 kmin = (start_num-2) / modpg + 1; // number of resgroups to start_num
uint64 kmax = (end_num - 2) / modpg + 1; // number of resgroups to end_num
uint64 krange = kmax - kmin + 1; // number of resgroups in range, at least 1
uint n = (krange < 37500000000000) ? 4 : (krange < 975000000000000) ? 6 : 8;
uint b = Bn * 1024 * n; // set seg size to optimize for selected PG
uint ks = (krange < b) ? krange : b; // segments resgroups size
cout << "segment size = " << ks << " resgroups; seg array is [1 x " << ((ks-1) >> 6) + 1 << "] 64-bits\n";
uint64 maxpairs = krange * pairscnt; // maximum number of twinprime pcs
cout << "twinprime candidates = " << maxpairs << "; resgroups = " << krange << endl;
return make_tuple(modpg, res_0, ks, kmin, kmax, krange, pairscnt, restwins, resinvrs);
}
// Compute primes r1..sqrt(input_num), and store in global 'primes' array.
// Any algorithm (fast|small) is usable. Here the SoZ for P5 is used.
vector<uint> sozpg(uint val, uint res_0, uint64 start_num, uint64 end_num) {
uint md = 30; // P5's modulus value
int rscnt = 8; // P5's residue count
int res[8] = {7,11,13,17,19,23,29,31}; // P5's residues list
int posn[30] = {0,0,0,0,0,0,0,0,0,1,0,2,0,0,0,3,0,4,0,0,0,5,0,0,0,0,0,6,0,7};
uint kmax = (val - 7) / md + 1; // number of resgroups upto input value
vector<uint8_t> prms(kmax, 0); // byte array of prime candidates init '0'
uint sqrtN = (uint) sqrt((double) val); // integer sqrt of val
uint modk = 0; int r = -1; uint k = 0; // init residue parameters
// mark the multiples of the primes r1..sqrtN in 'prms'
while (true) { // from list of pc positions in prms
++r; if (r == rscnt) {r = 0; modk += md; ++k;}
if ((prms[k] & (1 << r)) != 0) continue; // skip if not prime
uint prm_r = res[r]; // if prime save its residue value
uint prime = modk + prm_r; // numerate the prime value
if (prime > sqrtN) break; // we're finished if it's > sqrtN
for (int ri : res) { // mark prime's multiples in prms
uint prod = prm_r * ri - 2; // compute cross-product for prm_r|ri pair
uint8_t bit_r = 1 << posn[prod % md]; // bit mask for prod's residue
uint kpm = k * (prime + ri) + prod / md; // 1st resgroup for prime mult
for (; kpm < kmax; kpm += prime) prms[kpm] |= bit_r; // mark multiples
} }
// prms now contains the nonprime positions for the prime candidates r1..N
// extract only primes that are in inputs range into array 'primes'
vector<uint> primes; // dynamic array to store primes
for (uint k = 0; k < kmax; ++k) { // for each kth residue group
for (int r = 0; r < rscnt; ++r) { // check for each rth residue in resgroup
if ((prms[k] & (1 << r)) == 0){ // if bit location is prime
uint prime = md * k + res[r]; // numerate its value, store if in range
// check if prime has multiple in range, if so keep it, if not don't
uint64 n = start_num / prime; uint rem = start_num % prime;
if ((res_0 <= prime && prime <= val) && (prime * (n + 1) <= end_num || rem == 0)) primes.push_back(prime);
} } }
return primes;
}
// 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.
vector<uint64> nextp_init(uint rhi, uint64 kmin, uint modpg, vector<uint> primes, vector<int> resinvrs) {
vector<uint64> nextp(primes.size() * 2); // 1st mults array for twinpair
uint r_hi = rhi; // upper twinpair residue
uint r_lo = rhi - 2; // lower twinpair residue
for (int j = 0; j < primes.size(); ++j) { // for each prime r1..sqrt(N)
uint prime = primes[j]; // get the prime
uint k = (prime - 2) / modpg; // find the resgroup it's in
uint r = (prime - 2) % modpg + 2; // and its residue value
uint64 r_inv = resinvrs[r]; // and its residue inverse
uint64 rl = (r_inv * r_lo - 2) % modpg + 2; // compute the ri for r for lo_tp
uint64 rh = (r_inv * r_hi - 2) % modpg + 2; // compute the ri for r for hi_tp
uint64 kl = k * (prime + rl) + (r * rl - 2) / modpg; // kl 1st mult restroup
uint64 kh = k * (prime + rh) + (r * rh - 2) / modpg; // kh 1st mult restroup
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 << 1] = kl;
nextp[j << 1 | 1] = kh;
}
return nextp;
}
// 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
tuple<uint64, uint64> twins_sieve(uint r_hi, uint64 kmin, uint64 kmax, uint ks, uint64 start_num,
uint64 end_num, uint modpg, vector<uint> primes, vector<int>resinvrs) {
uint s = 6; // shift value for 64 bits
uint bmask = (1 << s) - 1; // bitmask val for 64 bits
uint64 sum = 0; // init tp cnt for this twin pair
uint64 ki = kmin - 1; // 1st resgroup seg val for each slice
uint kn = ks; // segment resgroup size
uint64 hi_tp = 0; // init value of largest tp in segment
uint64 k_max = kmax; // max number of resgroups for segment
vector<uint64> seg(((ks-1) >> s) + 1); // seg byte array for ks resgroups
if ((ki * modpg + r_hi - 2) < start_num) ki++; // ensure lo tp in range
if (((k_max - 1) * modpg + r_hi) > end_num) k_max--; // and hi tp in range
vector<uint64> nextp = nextp_init(r_hi, ki, modpg, primes, resinvrs); // 1st prime mults
while (ki < k_max) { // for kn resgroup size slices upto kmax
if (ks > (k_max - ki)) kn = (k_max - ki); // adjust kn for last seg size
for (int j = 0; j < primes.size(); ++j) { // for each prime index r1..sqrt(N)
uint prime = primes[j]; // for this prime
// for lower pair residue track
uint64 k = nextp[j << 1]; // starting from this lower resgroup in seg
for (; k < kn; k += prime) // mark primenth resgrouup bits prime mults
seg[k >> s] |= uint64(1) << (k & bmask);
nextp[j << 1] = k - kn; // save 1st lower resgroup in next eligible seg
// for upper pair residue track
k = nextp[j << 1 | 1]; // starting from this upper resgroup in seg
for (; k < kn; k += prime) // mark primenth resgrouup bits prime mults
seg[k >> s] |= uint64(1) << (k & bmask);
nextp[j << 1 | 1] = k - kn; // save 1st upper 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);
uint cnt = 0; // count the twinprimes in the segment
for (int k = 0; k <= (kn-1) >> s; ++k) cnt += __builtin_popcountl(~seg[k]); // tp count
if (cnt > 0) { // if segment has twin primes
sum += cnt; // add the segment count to total count
uint upk = kn - 1; // starting from end of segment, count backwards to largest tp
while ((seg[upk >> s] & (uint64(1) << (upk & bmask))) != 0) upk--;
hi_tp = ki + upk; // numerate its full resgroup value
}
ki += ks; // set 1st resgroup val of next seg slice
if (ki < k_max) {for (int b = 0; b <= ((kn-1) >> s); ++b) seg[b] = 0;} // set all seg bits prime
}
hi_tp = ((r_hi > end_num) || (sum == 0)) ? 1 : hi_tp * modpg + r_hi;
return make_tuple(hi_tp, sum);
}
int main() {
vector<uint64> x;
uint64 num;
while ( cin >> num ) x.push_back(num);
uint64 end_num = (x[0] > 3) ? x[0] : 3;
uint64 start_num = (x[1] > 3) ? x[1] : 3;
if (start_num > end_num) swap(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
if (end_num - start_num < 2) { start_num = 7; end_num = 7; }
cout << "threads = " << omp_get_max_threads() << endl;
auto ts = chrono::system_clock::now(); // start timing sieve setup execution
int modpg, res_0, pairscnt; vector<int> restwins, resinvrs;
uint64 ks, kmin, kmax, krange;
tie(modpg, res_0, ks, kmin, kmax, krange, pairscnt, restwins, resinvrs) = setSieveParameters(start_num, end_num);
// create sieve primes <= sqrt(end_num), only use those whose multiples within inputs range
vector<uint> primes;
(end_num < 49) ? primes = {5} : primes = sozpg((uint) sqrt((double) end_num), res_0, start_num, end_num);
cout << "each of "<< pairscnt << " threads has nextp["<< 2 << " x " <<primes.size() << "] array\n";
uint64 twinscnt = 0; // number of twin primes in range
int lo_range = restwins[0] - 3; // lo_range = lo_tp - 1
for (int tp : {3, 5, 11, 17}) { // lo tps for available Pn's
if (end_num == 3) break; // if 3 end of range, no twin primes
if (start_num <= tp && tp < lo_range) twinscnt++; // cnt small tps if any
}
auto te = chrono::system_clock::now() - ts; // sieve setup time
chrono::duration<double> seconds = te;
cout << "setup time = " << seconds.count() << " secs\n";
cout << "perform twinprimes ssoz sieve\n";
auto t1 = chrono::system_clock::now(); // start timing ssoz sieve execution
vector<uint64> cnts(pairscnt); // twins cnts for each twinpair
vector<uint64> lastwins(pairscnt); // last|largest twin per twinpair
std::atomic<uint> threadscnt(0); // count of finished threads
#pragma omp parallel for // process twinpairs in parrallel
for (int i = 0; i < pairscnt; i++) { // process each twinpair in own thread
uint64 l, c;
tie(l, c) = twins_sieve(restwins[i], kmin, kmax, ks, start_num, end_num, modpg, primes, resinvrs);
lastwins[i] = l; cnts[i] = c;
cout << "\r" << threadscnt++ << " of " << pairscnt << " twinpairs done";
}
cout << "\r" << pairscnt << " of " << pairscnt << " twinpairs done";
uint64 last_twin = 0; // init val for largest twinprime
for (int 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;
uint kn = krange % ks;
if (kn == 0) kn = ks;
auto t2 = chrono::system_clock::now() - t1;
chrono::duration<double> secs1 = t2;
chrono::duration<double> secs2 = t2 + te;
cout << "\nsieve time = " << secs1.count() << " secs\n";
cout << "total time = " << secs2.count() << " secs\n";
cout << "last segment = "<<kn<<" resgroups; segment slices = "<< (krange-1)/ks + 1<<"\n";
cout << "total twins = " << twinscnt << "; last twin = " << (last_twin - 1) << "+/-1"<< endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment