#!/usr/bin/python import sys import math import time import psyco import list_bigno import miller_rabin psyco.full() def prime_via_trial_division(top, bn, candidate_prime): largest_possible_divisor = int(math.sqrt(candidate_prime)) + 1 candidate_is_prime = True for trial_divisor in bn.iterator(): if trial_divisor > largest_possible_divisor: break # this is a reasonable trial divisor if candidate_prime % trial_divisor == 0: # this is composite candidate_is_prime = False break return candidate_is_prime def generate(t0, top, bn): # basic trial division, except we add in a Miller-Rabin test prev_time = time.time() number_of_primes_found = 1 for candidate_prime in xrange(3L, top, 2L): #sys.stderr.write('Testing %d\n' % candidate_prime) definitely_composite = miller_rabin.composite_by_miller_rabin(candidate_prime, verbose=False) if definitely_composite: #sys.stderr.write('Composite by Miller-Rabin\n') continue if candidate_prime < 341550071728321: # http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test # Miller-Rabin is accurate up to 341550071728321 with the constants I've chosen candidate_is_prime = True else: #sys.stderr.write('sieving\n') # at this point, we're somewhat confident that candidate_prime is prime, so we do the expensive trial division candidate_is_prime = prime_via_trial_division(top, bn, candidate_prime) if candidate_is_prime: number_of_primes_found += 1 current_time = time.time() if current_time > prev_time + 1: prev_time = current_time sys.stderr.write('%dnth prime found, value %d, %d digits, %.1f primes found per second, %3.1f%% complete (assuming linear), %s\n' % \ (number_of_primes_found, \ candidate_prime, \ math.ceil(math.log10(candidate_prime)), \ number_of_primes_found / (time.time() - t0), \ float(candidate_prime) / top * 100.0, bn.fully_from_memory(candidate_prime))) bn.append(candidate_prime) def report(top, bn): print 2 for i in bn.iterator(): print i def main(): #maximum_prime_to_attempt = 10L**8 maximum_prime_to_attempt = int(sys.argv[1]) sys.stderr.write('Maximum prime to attempt: %d\n' % maximum_prime_to_attempt) #maximum_prime_to_attempt = 10**4 maximum_prime_stored_in_memory = 8*10L**7L # we look for primes up to top, satisfying primes with value up to top/2 from an in-memory set. Other primes are stored # on disk as a bitmap bn = list_bigno.list_bigno(maximum_prime_stored_in_memory) t0 = time.time() generate(t0, maximum_prime_to_attempt, bn) report(maximum_prime_to_attempt, bn) if True: main() else: import cProfile cProfile.run('main()')