# A320920 - a(n) is the smallest number m such that binomial(m,n) is nonzero and is divisible by n!. 
# Python program by Bert Dobbelaere

# factorint(n!)
def fact_factors(n):
	flst={}
	for k in range(1,n+1):
		x=k
		while x>1:
			p=2
			while x%p != 0:
				p+=1
			x//=p
			if not p in flst:
				flst[p]=0
			flst[p]+=1
	return flst

# For each position with window of size n, find minimal exponent of prime power divisor p^e
# so that the sum of all p powers over the window is at least totalpwr 
def scan_minpwr(n, p, totalpwr):
	l=[]
	for k in range(n):
		t=0
		for q in range(n):
			if q!=k:
				d=q-k
				qq=0
				while (d%p)==0:
					t+=1
					d//=p
		l.append(totalpwr-t)
	return l

# For each position within window of size n, find primes p and their minimal powers for each location
def findminpwrs(n):
	flst=fact_factors(n)
	mp={}
	for p in flst:
		mp[p]=scan_minpwr(n,p, 2*flst[p])
	return mp			

# List of prime powers and associated remainder lists for valid a(n)
def findresidu_spec(n):
	mp=findminpwrs(n)
	spec={}
	for p in mp:
		maxpwr=max(mp[p])
		lst=[]
		for k in range(n):
			res=[k]
			pwr=mp[p][k]
			while pwr<maxpwr:
				res=[x+m*p**pwr for m in range(p) for x in res]
				pwr+=1
			lst+=res
		spec[p**maxpwr]=sorted(lst)
	return spec

# Extended Euclidean algorithm
def bezout_coeffs(a, b):
    s = 0; old_s = 1
    t = 1; old_t = 0
    r = b; old_r = a  
    while r != 0 :
        quotient = old_r // r
        (old_r, r) = (r, old_r - quotient * r)
        (old_s, s) = (s, old_s - quotient * s)
        (old_t, t) = (t, old_t - quotient * t)
    assert old_r==1 # gcd, no shared primes in arguments
    return (old_s, old_t)

# Chinese remainder theorem
# gcd(m1, m2) = 1, 0<=r1<m1, 0<=r2<m2
def chinese_rt(r1, m1, r2, m2):
	a,b=bezout_coeffs(m1,m2)
	x= r1 * m2 * b + r2 * m1 * a
	return x%(m1*m2)

# Max number of prime powers to use in our "wheel" (performance tradeof)
MAXFIRSTSPECS=3

for N in range(2,51):
	spec=findresidu_spec(N)
	numspecs=len(spec)
	ppwrs=sorted(spec.keys(), reverse=True)
	speclst=[spec[ppwrs[k]] for k in range(numspecs)]
	# Build the wheel
	m1=ppwrs[0] ; r1_lst=spec[m1]
	firstspecs=min(MAXFIRSTSPECS,numspecs)
	for k in range(1,firstspecs):
		m2=ppwrs[k] ; r2_lst=spec[m2]
		m12=m1*m2; r12_lst=sorted([chinese_rt(r1,m1,r2,m2) for r1 in r1_lst for r2 in r2_lst])
		m1=m12 ; r1_lst=r12_lst
	print("  Info: wheel size = %d / #residues  = %d"%(m1,len(r1_lst)))
	ofs=0
	found=False
	while not found:
		for r in r1_lst:
			m=ofs+r
			ok=True
			if m<N:
				ok=False
			for k in range(firstspecs, numspecs):
				pp=ppwrs[k]
				s=speclst[k]
				if not m%pp in s:
					ok=False
					break
			if ok:
				print("A320920(%d)=%d"%(N,m))
				found=True
				break
		ofs+=m1