# 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