My code works (as far as I can tell), so I'm looking mostly for readability and efficiency tips. Goal of my code: Generate a list of all numbers that have more divisors than any lower number up to a certain limit. In other words this sequence: https://oeis.org/A002182 Background: The number of divisors of \$n=p_1^{k1} \times p_2^{k2} \times p_3^{k3}\$ is just \$(k_1+1)\times(k_2+1)\times(k_3+1)\$... so the primes themselves don't matter, just the set of exponents. Thus to be the smallest possible number with that number of divisors, the exponents are required be in order from largest (on the smallest primes) to smallest (on the largest primes). My code breaks into 3 sections, the first figures out the largest prime needed. The second generates a list of candidate numbers which is a list of all numbers that fit the descending exponent requirement. Finally, the third part puts them in order and checks which of the candidates actually have more divisors than any previous numbers. Final Random question: **Is there a name for the middle algorithm?** An algorithm that finds a list of all lists that are under a certain threshold? A more generic version would probably set `pows[:marker] = [0] * marker` instead, but in my case I only want descending lists. from numpy import prod def highlydiv(n): """ Returns a list of all numbers with more divisors than any lower numbers, up to a limit n """ primeprod = 1 plist = [] for p in gen_prime(50): #gen_prime: generates primes up to n, can just use [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47] primeprod *= p if primeprod <= n: plist.append(p) else: break pows = [0 for x in plist] breakout = False candidate_list = [] calcsum = 1 while True: marker = 0 pows[marker] += 1 calcsum *= 2 while calcsum > n: marker += 1 if marker == len(pows): breakout = True break pows[marker] += 1 pows[:marker] = [pows[marker]] * marker calcsum = prod([p**k for p, k in zip(plist, pows)]) if breakout: break ndivs = prod([x+1 for x in pows]) candidate_list.append((calcsum,ndivs)) candidate_list.sort() maxdivs = 0 final_list = [] for candidate, ndivs in candidate_list: if ndivs > maxdivs: maxdivs = ndivs final_list.append(candidate) return final_list print(highlydiv(1000000))