# Primes and SemiPrimes in Binary

I'm in the midst of re-learning python. My current effort is an implementation around prime factorization. I am attempting to reproduce something similar to the table found here (in binary): https://math.stackexchange.com/questions/586010/relative-size-of-most-factors-of-semiprimes-close

Performance seems to be rather dreadful. I've looked at profiling of the code in cProfile but I don't see any obvious improvements.

I'm considering: (a) Installing PyPy to see if that improves performance and (b) Leveraging Numpy further, although precisely how is a little vague to me. I'm working through the following to see if it can help: https://www.freecodecamp.org/news/if-you-have-slow-loops-in-python-you-can-fix-it-until-you-cant-3a39e03b6f35/

I've already refactored to eliminate the need for a list of semiprimes and that helped a bit, but I'm not sure if there are obvious inefficiencies that I'm not seeing.

miller_rabin.py is not mine. I include it here for completeness. My effort centers on the performance of prime_stats_bin.py. I recognize that primality testing is going to be a large component of the cost, I'm just not understanding why I can't even get to 32 bits within a reasonable amount of time. I'm running this on a 5 year old Windows 10 laptop (Core i5, 2.5GHz, 8GB RAM), under Python 3.7.3 (installed yesterday).

As for the goal of this effort, I'd like to know whether its reasonable to expect output for primes up to 32 bits from this program within a day.

With 32 bits in the prime expansion, the list of primes I generate is about 2xe⁸ which given my 32 bit installation of Python, I'm not sure it will maintain.

prime_stats_bin.py
import math
import random
import hashlib
import numpy
from miller_rabin import *
from time import process_time

time1 = process_time()

iexp = 2 #init exponent
poe = 2 #pieces of eight
mexp = 8*poe #max exponent
cexp = 2 #curr exponent
primes =  #list of primes
counts = [] #list of counts (num primes, num of semiprimes, num of close semiprimes)

f = open("prime_stats_bin.txt", "w", encoding="utf-8")
f.write("BinDigits\tPrimes\tSemiprimes\tClose semiPrimes\n")
for x in range(iexp, mexp+1): #initialize lists for all exponents
counts.append([1,1,1]) #since are loop looks only at odd primes, initialize with 1s in all counts
for x in range(3, pow(2,mexp), 2):
if numpy.log2(x) >= cexp: #output to file when we pass a power of 2
f.write(str(cexp)+":\t"+str(counts[cexp-2])+"("+str(round((counts[cexp-2]*numpy.log(pow(2,cexp)))/pow(2,cexp),4))+")\t"+
str(counts[cexp-2])+" ("+str(round((counts[cexp-2]*numpy.log(pow(2,cexp)))/(pow(2,cexp)*numpy.log(numpy.log(pow(2,cexp)))),4))+")\t"+
str(counts[cexp-2])+" ("+str(round((counts[cexp-2]*numpy.log(pow(2,cexp)))/pow(2,cexp),4))+";"+str(round((100*counts[cexp-2])/counts[cexp-2],4))+")\n")
if cexp%8==0:
f.write('time elapsed: %f seconds'%(process_time()-time1)+"\n")
cexp += 1
if is_prime(x):
primes.append(x)
for z in range(iexp, mexp+1):
if x <= pow(2,z):
counts[z-2] += 1 #increment count of primes less than z taken to the power of 2
for y in primes:
if y*x > pow(2,mexp): #only need to worry about semiprimes that are less than our upper bound
break
for w in range(iexp, mexp+1):
if y*x <= pow(2,w):
counts[w-2] += 1 #increment count of semiprimes less than w taken to the power of 2
if x <= pow(y,2):
counts[w-2] += 1 #increment count of close (p*q=N, where p<q and q<=p^2) semiprimes less than w taken to the power of 2
f.write(str(mexp)+":\t"+str(counts[mexp-2])+"("+str(round((counts[mexp-2]*numpy.log(pow(2,mexp)))/pow(2,mexp),4))+")\t"+
str(counts[mexp-2])+" ("+str(round((counts[mexp-2]*numpy.log(pow(2,mexp-1)))/(pow(2,mexp)*numpy.log(numpy.log(pow(2,mexp)))),4))+")\t"+
str(counts[mexp-2])+" ("+str(round((counts[mexp-2]*numpy.log(pow(2,mexp-1)))/pow(2,mexp),4))+";"+str(round((100*counts[mexp-2])/counts[mexp-2],4))+")\n")
f.write('time elapsed: %f seconds'%(process_time()-time1)+"\n")
f.close

miller_rabin.py
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])

_known_primes = [2, 3]
_known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]

output and timing
BinDigits   Primes      Semiprimes  Close semiPrimes
2:  2(0.6931)   1 (1.061)   1 (0.3466;100.0)
3:  4(1.0397)   2 (0.7101)  2 (0.5199;100.0)
4:  6(1.0397)   6 (1.0196)  4 (0.6931;66.6667)
5:  11(1.1913)  10 (0.8714) 6 (0.6498;60.0)
6:  18(1.1697)  22 (1.0031) 9 (0.5848;40.9091)
7:  31(1.1751)  42 (1.008)  17 (0.6444;40.4762)
8:  54(1.1697)  82 (1.0369) 28 (0.6065;34.1463)
time elapsed: 0.140625 seconds
9:  97(1.1819)  157 (1.0449)    47 (0.5727;29.9363)
10: 172(1.1643) 304 (1.0629)    89 (0.6024;29.2763)
11: 309(1.1504) 589 (1.0795)    171 (0.6366;29.0323)
12: 564(1.1453) 1124 (1.0775)   311 (0.6315;27.669)
13: 1028(1.1308)    2186 (1.0937)   584 (0.6424;26.7155)
14: 1900(1.1253)    4192 (1.0926)   1086 (0.6432;25.9065)
15: 3512(1.1143)    8110 (1.099)    2093 (0.6641;25.8076)
16: 6542(1.1071)    15658 (1.1013)  4023 (0.6808;25.6929)
time elapsed: 27.906250 seconds
17: 12251(1.1014)   30253 (1.1026)  7617 (0.6848;25.1777)
18: 23000(1.0947)   58546 (1.1041)  14597 (0.6947;24.9325)
19: 43390(1.0899)   113307 (1.1041) 27817 (0.6987;24.5501)
20: 82025(1.0844)   219759 (1.105)  53301 (0.7047;24.2543)
21: 155611(1.0801)  426180 (1.1046) 101532 (0.7047;23.8237)
22: 295947(1.076)   827702 (1.1045) 195376 (0.7103;23.6046)
24: 1077871(1.0688) 3129211 (1.1036)    721504 (0.7154;23.0571)
time elapsed: 292.140625 seconds

Update: 6/22/2019

Many iterations later...

• Found that my initial code wouldn't operate beyond 2²⁷ due to memory limitations so I installed Python 64 bit.
• With Python 64 bit, I still encountered that my code aborted due to memory limitations between 2³⁰ and 2³¹.
• I jettisoned the Miller_Rabin as it was counterproductive since I was already calculating primes, replacing it with a Sieve of Eratosthenes.

At this point, getting this thing to run up to 2³² was a matter of pride.

I found that trying to keep 200M primes in memory @8 bytes per prime was leading to interminable delay as it was swapped to disk.

Tried a bit-string solution which didn't function as maintaining a 2³²-bit long int didn't perform.

Finally, resolved to use a segmented bit-string within a list to minimize the memory footprint while simultaneously limiting the length of the bit shifts that would be entailed.

Success!

While it takes 8 hours to run on my machine, I can count precisely the number of close semiprimes less than or equal to 2³².

With respect to the good suggestions I received, I have:

• implemented black formatting
• used f-string formatting of output
• refactored code to create functions, making the code more understandable
• modified my code import strategy, eliminated unused module and making clear the source of the functions I was leveraging

## TODO:

• typing
• investigate Sieve of Atkin

Thanks. Updated code below.

primes_bit.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

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
# 3 registers included: num of primes, num of semiprimes, num of close semiprimes
def init_prime_cnts():
for p in range(power_of_two):
prime_cnts.append([0, 0, 0])

digits_hdr = "Limit(Bin)"
primes_hdr = "Primes"
sprimes_hdr = "Semiprimes"
close_sprimes_hdr = "Close semiprimes"
print(f"{digits_hdr:^10}{primes_hdr:^20}{sprimes_hdr:^20}{close_sprimes_hdr:^20}")

def output_body(idx):
digits = idx + 2
# count of primes <= 2**digits
primes = prime_cnts[idx]
# how close is (primes) to estimate
prime_est = (primes * log(2 ** digits)) / 2 ** digits
# count of semiprimes <= 2**digits
sprimes = prime_cnts[idx]
# 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]
# 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
print(
f"{digits:>10}{primes:>10} ({prime_est:.4f}){sprimes:>10} ({sprime_est:.4f}){close_sprimes:>10} ({csprimes_est1:.4f}; {csprimes_est2:.4f})"
)
prime_cnts[idx + 1] += prime_cnts[idx]
prime_cnts[idx + 1] += prime_cnts[idx]
prime_cnts[idx + 1] += prime_cnts[idx]
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 = format(
prime_bits[segmnt], "0" + str(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] += 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 = format(
prime_bits[segmnt_innr], "0" + str(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
if inner_loop_num <= o_loop_num ** 2:
prime_cnts[inner_loop_idx] += 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
power_of_two = 32
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()

• Change poe to 3 to see what I'm currently limited to. poe of 4 is what I'm attempting to get to. Jun 18, 2019 at 2:19
• I tested this varying the bit_per_segment value and came to the conclusion that, even though it has a wide and shallow trough, 2**7 seems to offer me the best results from a performance perspective (e.g. shortest total execution time) Jun 22, 2019 at 19:39

Some suggestions related to review performance:

• This code is really hard to read. As far as I can tell every single variable except for primes is abbreviated to the point where I need to hold the entire program in memory in order to reason about any part of it. Naming is really important for readability, and readability is really important for comprehensibility. Relevant quote:

Everyone knows that debugging is twice as hard as writing a program in the first place. So if you're as clever as you can be when you write it, how will you ever debug it?

• Pulling out functions should also make the code much easier to understand. I would recommend adding type hints and validating them using a strict mypy configuration
• black can format your code to be much more idiomatic.
• flake8 can give you some hints about writing idiomatic Python.
• There are at least three unused imports.
• I would recommend avoiding import * in general; in a dynamic language like Python it's harder to reason about the code in that case, and it's easier to run into naming conflicts.

Performance-related suggestions:

• It looks like mexp is only set once, but pow(2,mexp) is calculated in a bunch of places. This can be calculated once outside both loops. There are other duplicate calculations - they should all be run only once, and in the outermost context.
• That answer might be worth a bounty just for the quote on debugging ;-) Jun 18, 2019 at 5:47
• Unfortunately, I don't yet have a reputation that will allow me to grant a bounty. Jun 22, 2019 at 19:40
• XD No worries @WolfLarson
– l0b0
Jun 23, 2019 at 6:21

My current effort is an implementation around prime factorization. I am attempting to reproduce something similar to the table found here

Prime factorisation is the wrong approach to build that table. The efficient way to build it is to generate a list of primes using a sieve of some kind (a good implementation of Eratosthenes is definitely good enough up to 24 bits, but for 32 bits you might want to think about Atkins-Bernstein, and for 64 bits I would definitely prefer Atkins-Bernstein).

Then with the list of primes you can generate a list of semiprimes with a double-loop (making sure to break early from the inner loop when you pass the threshold!) There are about 600 million 32-bit semiprimes, so this should meet your performance requirements.

    f.write(str(cexp)+":\t"+str(counts[cexp-2])+"("+str(round((counts[cexp-2]*numpy.log(pow(2,cexp)))/pow(2,cexp),4))+")\t"+
str(counts[cexp-2])+" ("+str(round((counts[cexp-2]*numpy.log(pow(2,cexp)))/(pow(2,cexp)*numpy.log(numpy.log(pow(2,cexp)))),4))+")\t"+
str(counts[cexp-2])+" ("+str(round((counts[cexp-2]*numpy.log(pow(2,cexp)))/pow(2,cexp),4))+";"+str(round((100*counts[cexp-2])/counts[cexp-2],4))+")\n")


That is frankly illegible. Pull out some variables for the calculations, and then look at f strings. And rather than copy-pasting it, define a function to compose the string.

        for z in range(iexp, mexp+1):
if x <= pow(2,z):
counts[z-2] += 1 #increment count of primes less than z taken to the power of 2


This probably isn't a major bottleneck, but it's still inefficient. If you maintain a count of how many integers in the ranges [a, b), [b, c), [c, d) etc. have property P then you can produce the counts for integers in the ranges [a, b), [a, c), [a, d) with the property P with just two additions when you output the tables, rather than two additions for each such integer in [a, b) and one for each in [b, c).