Finding integer with the most natural dividers

Problem description: You are given a segment [568023; 569230]. Find integer within given range that has the most natural dividers. If 2 or more integers are found, return the lowest integer among them.

Solution:

def allnums(k):
cnt = 0
for i in range(2, k // 2):
if k % i == 0:
cnt += 1
return cnt

MAX_NUM = 0
MAX = 0

for i in range(568023, 569231):
cnt = allnums(i)
if MAX < cnt:
MAX = cnt
MAX_NUM = i

print(MAX, ' ', MAX_NUM)

Output: 141 568260

Runtime: 17s

Question: Is there any way to make this code more concise, maybe also improving its runtime? Thank you in advance.

It's not clear what you mean with "natural dividers". Better specify that, show examples, and document allnums to say what it does. A better name would also help. For example allnums(6) returns only 1, despite having divisors 1, 2, 3 and 6. At the very least I'd expect answer 2 (for the divisors 2 and 3). So probably it's incorrect, but we can't even tell because you didn't say what exactly it's supposed to do.

Anyway, as usual it suffices to try divisors until the square root, not until half. If you find a divisor i of k, then don't just count that divisor but also count its partner divisor k/i. Unless that's the same.

And max can take an iterable and key.

So:

from math import isqrt

def divisors(k):
"""Compute the number of k's divisors 1 < divisor < k."""
cnt = 0
for i in range(2, isqrt(k) + 1):
if k % i == 0:
cnt += 2 - (k // i == i)
return cnt

print(max(range(568023, 569231), key=divisors))

• you need (k // i == i) in cnt += 2 - (k // i == i) to make it not count the square root of the integer twice, right? Dec 10 '20 at 12:48
• @JoshJohnson Yes, that's the "Unless that's the same" from the explanation. Dec 10 '20 at 16:09

As far as conciseness, the naive approach is probably about as good as you can get. A slightly simpler solution would use a generator expression + sum

def allnums_comprehension(k):
return sum(1 for i in range(2, k // 2) if k % i == 0)

In order to improve performance, we'll need to examine some traits of the numbers that do (or don't) meet our criteria.

A few observations:

1. If a number $$\k\$$ is not divisible by $$\x\$$, then no multiple of $$\x\$$ will be a natural divisor of $$\k\$$
2. Every natural divisor $$\x\$$ of $$\k\$$ has another number $$\y\$$ that is also a natural divisor of $$\k\$$ (e.g. $$\x * y = k\$$). (If $$\x = \sqrt k\$$ then $$\x=y\$$)

If we want to pursue the first observation, we might try to use a sieving approach, similar to one you might use for the Sieve of Eratosthenes. I'll save you the time of trying this one out - its way worse (~10x slower on my computer). The $$\O\$$ notation for this approach is left as an exercise for the reader.

The second observation is much nicer; if we look at the original algorithm, the cost is $$\O(n)\$$. Under the new algorithm we only need to try up to $$\\sqrt k\$$, giving us a complexity of $$\O(\sqrt n)\$$.

The new function is somewhat less readable, but still fairly simple overall:

def allnums_sqrt(k):
count = 1 if k % 2 == 0 else 0
return count + sum(2 if i * i != k else 1
for i in range(3, math.floor(math.sqrt(k)))
if k % i == 0)

Please note I've made some... questionable decisions here, because (as pointed out by superb rain, your definition of natural divisor is somewhat odd. The magic numbers here are a result of my attempt to replicate your results.

You're doing trial division. As for primality checking, something like the sieve of Eratosthenes is faster. Have a list representing the counts of the desired range, then go through candidate divisors and increase the counts of their multiples. Only use candidates up to the largest square root, and also count their partner divisors (as described in the other answers) unless they're among the candidates themselves and thus increase the counts themselves.

Execution time: 1.9 ms on my laptop (AJNeufeld's takes 4.0 ms there, didn't measure the others).

from math import isqrt

start, stop = 568023, 569231

size = stop - start
cnt =  * size
candidates = range(2, isqrt(stop-1) + 1)
for i in candidates:
for multiple in range(-start % i, size, i):
cnt[multiple] += 1
if (start + multiple) // i not in candidates:
cnt[multiple] += 1

print(start + cnt.index(max(cnt)))

(Note that multiple is the index of the multiple, the actual value of the multiple is start + multiple.)

Mixing mine with AJNeufeld's

Taking my Eratosthenes-like list approach, combined with AJNeufeld's calculation of numbers of divisors from prime powers. Instead of finding the prime divisor powers with trial division, I directly visit only multiples of the prime powers. So AJNeufeld's structure is like

for n in TheRange:
for p in PRIMES:
if p divides n:
for powers of p:
...

while I now do this:

for p in PRIMES:
for P in powers of p:
for multiple of P in TheRange:
...

Where AJNeufeld has the single current n and its number of divisors divisors, I have lists n and divisors representing the same, but for all numbers in the desired range. So I treat all numbers in the range in parallel, one prime at a time.

AJNeufeld determines the number k of exactly how often a number is divisible by a prime factor p. That then multiplies the number of divisors by k + 1. Here I simulate that by first going over all multiples of p and multiplying the number of divisors by 2. Then I go over all multiples of p**2 and instead multiply the number of divisors by 3 (by dividing by 2 and then multiplying by 3). And so on for larger powers of p. In the end, each number divisible by p exactly k times will have its number of divisors multiplied by k + 1.

Execution time: 0.76 ms on my laptop:

PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461,
463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563,
569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739,
743, 751, ]

start, stop = 568023, 569231

size = stop - start
n = list(range(start, stop))
divisors =  * size

for p in PRIMES:
P = p
k = 1
while P < stop:
K = k + 1
for multiple in range(-start % P, size, P):
divisors[multiple] = divisors[multiple] // k * K
n[multiple] //= p
P *= p
k += 1
for i in range(size):
if n[i] != 1:
divisors[i] *= 2

print(start + divisors.index(max(divisors)))
• Nice solution. I ran both on range(2, 10_000_000). Your solution is about 3 times faster. Eratosthenes 1, Math 0. (I'm disappointed that the math solution doesn't win.) Dec 10 '20 at 5:22
• Wow! Nice mashup of solutions! Try for i, n_i in enumerate(n): if n_i != 1: Dec 10 '20 at 6:45
• @AJNeufeld Didn't make it faster. That part doesn't take much time anyway. And I like it better the way it is, as it is a bit shorter and I want to keep treating n and divisors equally and keep it close to your style (I even changed my usual exponent name e to your k now :-). Dec 10 '20 at 6:58

Naming

PEP 8 - the Style Guide for Python Code list several conventions for Python coding.

First, function and variable names formed by multiple words should be snake_case. The allnum() function violates this; it should be named all_num(), but a better name would be num_divisors().

Second, only constants should be in UPPER_CASE. MAX_NUM and MAX are not constants. They should be lowercase. Perhaps max_divisors and number_with_max_divisors.

Efficiency

Concise and efficiency do not necessarily go hand-in-hand.

The number of natural divisors of an integer is a well studied area of mathematics.

If N is expressed in its prime factorization form:

$$N = \prod_{i} p_i^{k_i}$$

then the number of natural divisors of N can be directly computed as:

$$NumDivisors = \prod_{i} (k_i + 1)$$

For example: $$\200 = 2^3 \cdot 5^2\$$, so 200 has $$\(3 + 1)\cdot(2+1) = 12\$$ natural divisors (1, 2, 4, 5, 8, 10, 20, 25, 40, 50, 100, 200).

Since the largest number we're looking for the number of natural divisors of is 569230, the largest possible prime factor is $$\\lfloor\sqrt{569230}\rfloor = 754\$$. We only need prime numbers up to 754!

PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461,
463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563,
569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739,
743, 751, ]

def num_divisors(n):

divisors = 1
limit = isqrt(n)

for p in PRIMES:

if p > limit:
break

if n % p == 0:
k = 1
n //= p
while n % p == 0:
k += 1
n //= p

divisors *= k + 1
limit = isqrt(n)

if n != 1:
divisors *= 2

return divisors

As prime factors are found, the value n is reduced, which reduces the search space. Consider the last value in the range: 569230. Very quickly, both 2 and 5 are factored out, leaving n = 56923, which is a prime number, but we don't know that yet. As successive prime candidate factors are examined, we come to 239, which squared is larger than 56923. At this point, we can stop, ignoring all larger primes; none can be a factor.

Using this function, we can quickly (and concisely) determine the number of factors in each value of our range, and determine the number which produces the maximum:

print(max(range(568023, 569230 + 1), key=num_divisors))

Execution time: 7ms.

Conciseness

For a concise solution, look for a library which does all of the heavy lifting for you:

from sympy import divisor_count

print(max(range(568023, 569230 + 1), key=divisor_count))
• Tried a few more variations of this. Here are two, the first is my fastest and the second is my shortest. Dec 10 '20 at 5:38
• @KellyBundy I've added a few optimizations to mine, and it is down to 3.9ms -vs- 3.6ms for yours. I'm not going to post the not n % p improvement -- it does not improve readability. Dec 10 '20 at 5:40
• Ah yes, that's a bigger improvement. Though on my laptop, it's now 4.0 ms vs my 1.9 ms. I do find not n % p quite readable, I read it like the common expression "(leaves) no remainder". Dec 10 '20 at 5:50
• @KellyBundy Yes, it should be divisor_count. And I didn't brag about the time because it wasn't worth bragging about at 36ms. Dec 10 '20 at 15:45
• Woah, that seems surprisingly slow. Indeed not worth "bragging" about, rather worth warning about :-D Ok, your hardcoded just-long-enough list of primes is cheating a bit, but I wouldn't expect that sympy solution to be that much slower. Would've guessed anything from much faster to a bit slower. Dec 10 '20 at 16:42