# Distribution of Pairs of Semiprimes

I have a means of generating pairs of semiprimes of a given bit length. The process entails

• Step 1. Generate a semiprime P1 of bitlength N, with factors close in length (close as defined in the answer to the following [alpha=2]: relative size of most factors of semiprimes, close?)

• Step 2. Leaving intact the leading and trailing 1s, take the encapsulated bitlength of N-2 from P1 [The most significant and least significant bits are retained to assure that P2 in step 4 is the appropriate length and odd]

• Step 3. Apply a hash (MD5) to the N-2 bitlength number and truncate the digest in binary format to an N-2 bitlength

• Step 4. Produce P2, equivalent in length to P1, by concatenating 1+(truncated hash)+1

• Step 5. Test P2 to establish whether it is a semiprime with factors close in length

• Step 6. If not, discard P1 & P2 and start over

In reviewing the post I referenced previously, relative size of most factors of semiprimes, close? , I am attempting to calculate the number of paired semiprimes (P1 & P2) that exist up to a given bitlength.

From the estimate of $$\pi_{1/2;2}(x)$$ as $$x/log(x)$$

I calculated the probability that any pair of numbers up to a given length, $$x$$, related via the MD5 hash process above, satisfying the criteria that they are both close semiprimes at:

$$1/log(x)^2$$

To arrive at this I divided $$x/log(x)$$ by $$x$$ for the probability of a single semiprime and squared this to get the joint probability for both

I then multiplied this by $$x$$ to re-establish the population for a given bitlength,

I don't think I've got it right, however. When I generate the estimate $$\pi_{1/2;2;2}(x)$$ as $$x/log(x)^2$$, I get something which doesn't seem to jibe.

The number in parentheses below is calculated as

$$\frac{\pi_{1/2;2;2}(x) * log(x)^2}{x}$$

and it diverges from what I would expect. Mathematics is not my day job. I wonder where I've gone wrong and what an accurate estimate should look like.

List of bitstring primes up to 2**24 generated
time elapsed: 7.843750 seconds
LIMIT(Bin)    Primes          Semi        Close semi     Close semi pair
2              2              0              0              0 (0.0000)
3              4              2              2              0 (0.0000)
4              6              6              4              1 (0.4805)
5             11             10              6              2 (0.7507)
6             18             22              9              3 (0.8108)
7             31             42             17              5 (0.9196)
8             54             82             28             10 (1.2011)
time elapsed: 59.734375 seconds
9             97            157             47             17 (1.2922)
10            172            304             89             31 (1.4545)
11            309            589            171             68 (1.9303)
12            564           1124            311            113 (1.9087)
13           1028           2186            584            197 (1.9526)
14           1900           4192           1086            325 (1.8680)
15           3512           8110           2093            586 (1.9332)
16           6542          15658           4023           1062 (1.9931)
time elapsed: 364.593750 seconds
17          12251          30253           7617           1953 (2.0689)
18          23000          58546          14597           3712 (2.2043)
19          43390         113307          27817           6765 (2.2380)
20          82025         219759          53301          12672 (2.3225)
21         155611         426180         101532          23589 (2.3832)
22         295947         827702         195376          43998 (2.4393)
23         564163        1608668         374928          82250 (2.4920)
24        1077871        3129211         721504         153318 (2.5290)
time elapsed: 423.453125 seconds

Update (6/27/2019):

In reviewing my code for the output for the last column, I erred in at least two ways. I was comparing the factor bit length rather than the factor itself in determining 'closeness'. A second, smaller issue was that I was looking, post-hash, at only distinct semiprimes.

Once I corrected for these two issues, the estimate I generated previously looks inline with my results. This may be mere coincidence as I'm not convinced my logic is sound. If anyone would care to comment, feel free. I'm including my code (python) for anyone seeking to verify my results.

LIMIT: 16777216
List of bitstring primes up to 2**24 generated
time elapsed: 7.812500 seconds
LIMIT(Bin)    Primes          Semi        Close semi   Close semi pair
2              2              0              0              0 (0.0000)
3              4              2              2              0 (0.0000)
4              6              6              4              1 (0.4805)
5             11             10              6              1 (0.3754)
6             18             22              9              1 (0.2703)
7             31             42             17              1 (0.1839)
8             54             82             28              7 (0.8408)
time elapsed: 62.281250 seconds
9             97            157             47             10 (0.7601)
10            172            304             89             14 (0.6569)
11            309            589            171             27 (0.7664)
12            564           1124            311             44 (0.7432)
13           1028           2186            584             82 (0.8128)
14           1900           4192           1086            137 (0.7874)
15           3512           8110           2093            258 (0.8511)
16           6542          15658           4023            459 (0.8614)
time elapsed: 377.125000 seconds
17          12251          30253           7617            783 (0.8295)
18          23000          58546          14597           1454 (0.8634)
19          43390         113307          27817           2655 (0.8783)
20          82025         219759          53301           4940 (0.9054)
21         155611         426180         101532           9061 (0.9155)
22         295947         827702         195376          16908 (0.9374)
23         564163        1608668         374928          31200 (0.9453)
24        1077871        3129211         721504          57728 (0.9522)
time elapsed: 435.359375 seconds
primes_bit_close_factors.py
from time import process_time
from math import sqrt
from math import floor
from math import ceil
from numpy import log2
from numpy import log
from hashlib import md5
from miller_rabin import factors

time1 = process_time()

# initialize list of bitstring segments as all '1's [TRUE]
def init_prime_bits():
for n in range(LIMIT // NUM_PER_SEGMENT):
prime_bits.append(2 ** BIT_PER_SEGMENT - 1)

# populate list of prime bitstrings up to [LIMIT]
def pop_prime_bits():
print(f"LIMIT: {LIMIT}")
for x in range(3, SUB_LIMIT, 2):
segmnt = x // NUM_PER_SEGMENT
locatn = x % NUM_PER_SEGMENT // 2
if not (prime_bits[segmnt] >> locatn & 1):
continue
for y in range(x * x, LIMIT, x * 2):
segmnt = y // NUM_PER_SEGMENT
locatn = y % NUM_PER_SEGMENT // 2
prime_bits[segmnt] &= ~(1 << locatn)
print(f"List of bitstring primes up to 2**{POWER_OF_TWO} generated")
print(f"time elapsed: {process_time()-time1:.6f} seconds")

# initialize list of prime counts to '0' for each power of two being calculated
# 5 registers included: (1) num of primes, (2) num of semiprimes, (3) subset of [2] that are close,
#   (4) subset of [3] that are distinct, (5) subset of [4] that are of equal binary length
#   (6) subset of [5] whose bookended hash is also semiprime
def init_prime_cnts():
for p in range(POWER_OF_TWO):
prime_cnts.append([0, 0, 0, 0, 0, 0, 0])

digits_hdr = "LIMIT(Bin)"
primes_hdr = "Primes"
sprimes_hdr = "Semi"
csprimes_hdr = "Close semi"
cdsprimes_hdr = "Close dist"
cedsprimes_hdr = "Equal len"
hashsprimes_hdr = "Hash semi"
hasheqsprimes_hdr = "Hash eq semi"
hashcsprimes_hdr = "Close semi pair"
print(
f"{digits_hdr:^10}{primes_hdr:^15}{sprimes_hdr:^15}{csprimes_hdr:^15}{hashcsprimes_hdr:^15}"
)

def output_body(idx):
digits = idx + 2
# count of primes <= 2**digits
primes = prime_cnts[idx][0]
# how close is (primes) to estimate
prime_est = (primes * log(2 ** digits)) / 2 ** digits
# count of semiprimes <= 2**digits
sprimes = prime_cnts[idx][1]
# how close is (sprimes) to estimate
sprime_est = (sprimes * log(2 ** digits)) / (2 ** digits * log(log(2 ** digits)))
# count of close semiprimes
close_sprimes = prime_cnts[idx][2]
# how close is (close_sprimes) to estimate -- #1
csprimes_est1 = (close_sprimes * log(2 ** digits)) / 2 ** digits
# how close is (close_sprimes) to estimate -- #2
csprimes_est2 = 0 if sprimes == 0 else (100 * close_sprimes) / sprimes
# count of close distinct semiprimes
dist_sprimes = prime_cnts[idx][3]
# count of close equal length distinct semiprimes
equalen_sprimes = prime_cnts[idx][4]
# count of hash semiprimes
hash_sprimes = prime_cnts[idx][5]
# count of hash close length semiprimes
hash_close_sprimes = prime_cnts[idx][6]
# how close is (hash_close_sprimes) to estimate
hcsprimes = (hash_close_sprimes * log(2 ** digits) * log(2 ** digits)) / (2 ** digits)
print(
#        f"{digits:>10}{primes:>10} ({prime_est:.4f}){sprimes:>10} ({sprime_est:.4f}){close_sprimes:>10} ({csprimes_est1:.4f}; {csprimes_est2:.4f})"
f"{digits:>10}{primes:>15}{sprimes:>15}{close_sprimes:>15}{hash_close_sprimes:>15} ({hcsprimes:.4f})"
)
# update counters when a power of two is exceeded
prime_cnts[idx + 1][0] += prime_cnts[idx][0]
prime_cnts[idx + 1][1] += prime_cnts[idx][1]
prime_cnts[idx + 1][2] += prime_cnts[idx][2]
prime_cnts[idx + 1][3] += prime_cnts[idx][3]
prime_cnts[idx + 1][4] += prime_cnts[idx][4]
prime_cnts[idx + 1][5] += prime_cnts[idx][5]
prime_cnts[idx + 1][6] += prime_cnts[idx][6]
if digits % 8 == 0:
print(f"time elapsed: {process_time()-time1:.6f} seconds")

def outer_loop():
for n in range(0, LIMIT, 2):
segmnt = n // NUM_PER_SEGMENT
locatn = n % NUM_PER_SEGMENT // 2
if n % NUM_PER_SEGMENT == 0:
outer_loop_primes = f"{prime_bits[segmnt]:0{BIT_PER_SEGMENT}b}"[::-1]
if int(outer_loop_primes[locatn]):
outer_loop_num = n + 1
# this code implements a trick which labels the first bit in the bitstring, '1',
# as a prime and treats it for the purposes of the loops as '2'
if outer_loop_num == 1:
outer_loop_num = 2
outer_loop_idx = int(log2(outer_loop_num)) - 1
prime_cnts[outer_loop_idx][0] += 1
inner_loop(n, outer_loop_num, outer_loop_primes)
# print results when the power of two advances -- starting with 2 bits
if n != 0 and ceil(log2(n + 2)) == floor(log2(n + 2)):
output_body(int(log2(n + 2)) - 2)

def inner_loop(idx, o_loop_num, inner_loop_primes):
for p in range(idx, LIMIT, 2):
segmnt_innr = p // NUM_PER_SEGMENT
locatn_innr = p % NUM_PER_SEGMENT // 2
if p % NUM_PER_SEGMENT == 0:
inner_loop_primes = f"{prime_bits[segmnt_innr]:0{BIT_PER_SEGMENT}b}"[::-1]
if int(inner_loop_primes[locatn_innr]):
inner_loop_num = p + 1
# same trick as above, applied to the inner loop
if inner_loop_num == 1:
inner_loop_num = 2
inner_loop_prd = o_loop_num * inner_loop_num
if inner_loop_prd > LIMIT:
break
inner_loop_idx = int(log2(inner_loop_prd)) - 1
prime_cnts[inner_loop_idx][1] += 1
if inner_loop_num <= o_loop_num ** 2:
prime_cnts[inner_loop_idx][2] += 1
hash = md5(
bin(inner_loop_prd)[3:-1].encode()
).hexdigest()
hash_result = f"1{f'{int(hash,16):0128b}'[:inner_loop_idx+2]}1"
hash_sprime_factors = sorted(factors(int(hash_result, 2)))
if (len(hash_sprime_factors) == 3 and hash_sprime_factors[1] == sqrt(hash_sprime_factors[2])) or \
(len(hash_sprime_factors) == 4 and hash_sprime_factors[2] <= hash_sprime_factors[1]**2):
##                    print(
##                        f"{hash_sprime_factors[1]}:{len(bin(hash_sprime_factors[1]))-2}:{hash_sprime_factors[2]}:{len(bin(hash_sprime_factors[2]))-2}"
##                    )
prime_cnts[inner_loop_idx][6] += 1

def main():
init_prime_bits()
pop_prime_bits()
init_prime_cnts()
outer_loop()
print(f"time elapsed: {process_time()-time1:.6f} seconds")

# LIMIT is restricted to a power of 2 whose exponent is divisible by 4 (e.g. 2**12 is acceptable, 2**14 is not)
POWER_OF_TWO = 24
LIMIT = 2 ** POWER_OF_TWO
SUB_LIMIT = int(LIMIT ** 0.5)

# bitstring segment length within the list is a power of 2
BIT_PER_SEGMENT = 2 ** 7
# numbers per bitstring segment are doubled as we only store odd numbers
NUM_PER_SEGMENT = (BIT_PER_SEGMENT) * 2

# list of bitstring segments to maintain prime flags
prime_bits = []

# list of counts of primes by powers of two
prime_cnts = []

if __name__ == "__main__":
main()
miller_rabin.py
from functools import reduce

def _try_composite(a, d, n, s):
if pow(a, d, n) == 1:
return False
for i in range(s):
if pow(a, 2**i * d, n) == n-1:
return False
return True # n  is definitely composite

def is_prime(n, _precision_for_huge_n=16):
if n in _known_primes:
return True
if any((n % p) == 0 for p in _known_primes) or n in (0, 1):
return False
d, s = n - 1, 0
while not d % 2:
d, s = d >> 1, s + 1
# Returns exact according to http://primes.utm.edu/prove/prove2_3.html
if n < 1373653:
return not any(_try_composite(a, d, n, s) for a in (2, 3))
if n < 25326001:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))
if n < 118670087467:
if n == 3215031751:
return False
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7))
if n < 2152302898747:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
if n < 3474749660383:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
if n < 341550071728321:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
# otherwise
return not any(_try_composite(a, d, n, s)
for a in _known_primes[:_precision_for_huge_n])

def factors(n):