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).
Any comments would be appreciated.
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 = [2] #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][0])+"("+str(round((counts[cexp-2][0]*numpy.log(pow(2,cexp)))/pow(2,cexp),4))+")\t"+
str(counts[cexp-2][1])+" ("+str(round((counts[cexp-2][1]*numpy.log(pow(2,cexp)))/(pow(2,cexp)*numpy.log(numpy.log(pow(2,cexp)))),4))+")\t"+
str(counts[cexp-2][2])+" ("+str(round((counts[cexp-2][2]*numpy.log(pow(2,cexp)))/pow(2,cexp),4))+";"+str(round((100*counts[cexp-2][2])/counts[cexp-2][1],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][0] += 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] += 1 #increment count of semiprimes less than w taken to the power of 2
if x <= pow(y,2):
counts[w-2][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][0])+"("+str(round((counts[mexp-2][0]*numpy.log(pow(2,mexp)))/pow(2,mexp),4))+")\t"+
str(counts[mexp-2][1])+" ("+str(round((counts[mexp-2][1]*numpy.log(pow(2,mexp-1)))/(pow(2,mexp)*numpy.log(numpy.log(pow(2,mexp)))),4))+")\t"+
str(counts[mexp-2][2])+" ("+str(round((counts[mexp-2][2]*numpy.log(pow(2,mexp-1)))/pow(2,mexp),4))+";"+str(round((100*counts[mexp-2][2])/counts[mexp-2][1],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
Feel free to make whatever additional comments you see fit.
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])
def output_header():
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][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
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][0] += prime_cnts[idx][0]
prime_cnts[idx + 1][1] += prime_cnts[idx][1]
prime_cnts[idx + 1][2] += prime_cnts[idx][2]
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][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 = 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] += 1
if inner_loop_num <= o_loop_num ** 2:
prime_cnts[inner_loop_idx][2] += 1
def main():
init_prime_bits()
pop_prime_bits()
output_header()
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()
