# A145294 - Smallest x >= 0 such that the Euler polynomial x^2 + x + 41 has a prime divisor of multiplicity n. 
# Python program by Bert Dobbelaere

# We need the nextprime() function, factorint() is just for sanity checking the end result
from sympy import nextprime,factorint

# The Euler polynomial
def f(x):
	return x**2+x+41

# Composes a list of values for which f(x) = 0 (mod p^k) with 0 <= x < min(p^k, cutoff)
# The cutoff value is for performance reasons: we stop the search if the best we can do
# is worse than what we already found before.
def modCandidates(p,k,cutoff):
	l=[]
	if k==1:
		# For k=1, the list is trivial.
		for x in range(p):
			if (f(x)%p) == 0 and (cutoff==None or x<cutoff):
				l.append(x)
	else:
		# We know f(x+q*m)=f(x) (mod m) for q>=0 and m>0.
		# Any candidate x for which f(x)=0 (mod p^k) must also satisfy f(x)=0 (mod p^(k-1))
		# Therefore, x must be of the form s + q * p^(k-1) with s an integer for which f(s)=0 mod (p^(k-1))  
		l_in=modCandidates(p,k-1,cutoff)
		pk1=p**(k-1)
		pk=p**k
		for j in range(p):
			for m in l_in:
				x= j*pk1 + m # solutions listed in sorted order: smallest is found first.
				if (cutoff==None) or (x<cutoff):
					if f(x)%pk==0:
						l.append(x)
	return l

# Computing a sequence term
def a145294(n):
	# Although the smallest primes have no chance of success, let's start from 2 for completeness
	p=2
	minx=None # No smallest solution yet
	while (minx==None) or (p**n < f(minx)): # Don't evaluate primes that can't lead to a smaller solution
		l=modCandidates(p,n,minx)
		while len(l)>0 and (f(l[0])%(p**(n+1)))==0: # Reject factors with a multiplicity that is too high.
			l=l[1:]
		if len(l)>0:
			minx=l[0] # We found a solution for f(x)=0 mod (p^k) that is smaller than any previous one
		p=nextprime(p)
	return minx
	
for n in range(1,101):
	a_n=a145294(n)
	factorlist=factorint(f(a_n),limit=100000)
	maxexp=max(factorlist.values())
	assert maxexp==n # just verify that f(a(n)) contains a small prime with multiplicity >= n
	print("%d %d"%(n, a_n))