from math import log limit = 10**6 primes = [2] for p in range(3,limit,2): ..primes.append(p) ..for q in primes: ....if not p%q: ......del primes[-1] ......break ....if q*q > p: break def exists(n,p): #does p^n appear in a factorial? ..cap = int(log(n*(p-1),p))+1 ..for m in range(cap,1,-1): ....if not (n+1)%((p**m-1)//(p-1)): return 0 ....n %= ((p**m-1)//(p-1)) ..return 1 arr = [] for n in range(limit): ..count = 0 ..for p in primes: ....if p > n: break ....count += 1 - exists(n,p) ..if not count in arr: ....arr.append(count) ....print(n,count)