#!/usr/bin/python

import sys
import gmpy
import time
import psyco
import itertools
import miller_rabin

psyco.full()

def odds_from(n):
	while True:
		yield n
		n += 2

likely_prime_file = open('likely-prime.txt', 'w')
def write_likely_prime(n):
	likely_prime_file.write('%s\n' % n.digits())
	likely_prime_file.flush()

def main():
	# for real life:
	superexponent = gmpy.mpz(8)
	# for testing:
	#superexponent = gmpy.mpz(4)

	subexponent = gmpy.mpz(10)
	base = gmpy.mpz(10)

	sys.stderr.write('Exponent is %d^%d\n' % (subexponent, superexponent))
	sys.stderr.write('Computing first number to test\n')
	first_candidate = pow(base, pow(subexponent, superexponent) - 1) + 1
	digits = first_candidate.digits()
	len_digits = len(digits)
	sys.stderr.write('Number of digits is %d\n' % len_digits)
	if len_digits < 100:
		sys.stderr.write('%s\n' % first_candidate.digits())

	prev_time = time.time()

	sys.stderr.write('Reading small primes\n')
	small_primes_file = open('primes-below-100000000', 'r')
	small_primes = [ gmpy.mpz(int(small_prime)) for small_prime in small_primes_file.readlines() ]
	small_primes.remove(gmpy.mpz(2))
	small_primes_file.close()
	candidate_no = 0
	sys.stderr.write('Starting tests at %s\n' % time.ctime())
	t0 = time.time()
	for candidate_prime in odds_from(first_candidate):
		t1 = time.time()
		difference = t1 - t0
		if difference != 0 and candidate_no != 0:
			sys.stderr.write('%f candidates per day (%s)\n' % (24 * 60 * 60 * candidate_no / (t1 - t0), time.ctime()))
		candidate_no += 1
		digits_of_interest = 30
		sys.stderr.write('Checking candidate # %d' % candidate_no)
		sys.stderr.write(' (candidate is %d ^ (%d ^ %d - 1) + %d)\n' % (base, subexponent, superexponent, (candidate_no-1) * 2 + 1))
		if len_digits < 1000000:
			sys.stderr.write('(last %d digits are %s)\n' % (digits_of_interest, candidate_prime.digits()[-digits_of_interest:]))
		# this is a quick test that can eliminate the need for a lot of miller rabin'ing
		definitely_composite = miller_rabin.composite_by_small_primes(small_primes, candidate_prime, verbose=True)
		if definitely_composite:
			print '  Composite based on small prime divisibility'
			continue
		# miller-rabin
		definitely_composite = miller_rabin.composite_by_miller_rabin(candidate_prime)
		if definitely_composite:
			print '  Composite based on Miller-Rabin'
			continue
		# eventually we should have AKS or something here - but for now, we just write out the first likely prime we find
		cur_time = time.time()
		sys.stderr.write('%s: candidate # %d is likely prime\n' % (sys.argv[0], candidate_no))
		write_likely_prime(candidate_prime)
		prev_time = cur_time
	likely_prime_file.close()

main()