fastest way to produce a list of all divisors of a number-CodePudding

Given a number n, I want to generate a sorted list of all the unique divisors of n (with no duplications). Solving this problem is really straight forward, but what I'm interested in is the efficiency. What is the fastest way to do it?

Here is one way, with pure python:

def get_divisors(n):
    """
    :param n: positive integer.
    :return: list of all different divisors of n.
    """
    if n <= 0:
        return []
    divisors = [1, n]
    for div in range(1, int(n ** 0.5   1)):
        if n % div == 0:
            divisors.extend([n // div, div])
    return sorted(list(set(divisors)))

Any suggestions on how to optimize this? Numpy and other optimized libraries are welcome.

CodePudding user response：

You already have the square root optimization. Next would be to leverage numpy's parallelism:

import numpy as np

def npDivs(N):
    divs = np.arange(1,int(N**0.5) 1)  # potential divisors up to √N
    divs = divs[N % divs==0]           # divisors
    comp = N//divs[::-1]               # complement quotients
    return np.concatenate((divs,comp[divs[-1]==comp[0]:])) # combined
    
print(getDivs(1001**2))
[      1       7      11      13      49      77      91     121     143
     169     539     637     847    1001    1183    1573    1859    5929
    7007    8281   11011   13013   20449   77077   91091  143143 1002001]

comp[divs[-1]==comp[0]:] avoids repeating the square root if it is an integer.

An even faster approach, would be to get prime factors and combine them in a resulting set:

def getDivs(N):
    factors = {1}
    maxP  = int(N**0.5)
    p,inc = 2,1
    while p <= maxP:
        while N%p==0:
            factors.update([f*p for f in factors])
            N //= p
            maxP = int(N**0.5)
        p,inc = p inc,2
    if N>1:
        factors.update([f*N for f in factors])
    return sorted(factors)

Benchmarks:

from timeit import timeit
N = 1010101**2
print(timeit(lambda:getDivs(N),number=100))      # 0.0015
print(timeit(lambda:npDivs(N),number=100))       # 0.9753
print(timeit(lambda:get_divisors(N),number=100)) # 8.5605

CodePudding user response：

I expect the fastest way is to build the divisors "ground up" from a full prime factorization of the input. But you're not going to get an answer here for "the fastest way" to get a prime factorization to begin with: that's a huge topic all on its own, and still a very active area of research for large integers.

The way below simply uses trial division, up through the square root of what remains to be factored, but skipping multiples of 2 and 3 (except for 2 and 3 themselves). This is reasonably zippy through 32-bit ints.

Given a prime factorization

n == p0**e0 * p1**e1 * p2**e2 ...

Then the divisors of n are all and only the integers of that form with exponents less than or equal equal to the e0, e1, e2, .... itertools.product() can straightforwardly generate all such tuples of exponents, and then it's just a matter of doing the powers and multiplying the results.

def get_divisors(n):
    from math import isqrt, prod
    from itertools import accumulate, chain, cycle, product
    if n <= 1:
        return [n]
    ps = []
    exps = []
    limit = isqrt(n)
    for cand in chain([2, 3], accumulate(cycle([2, 4]),
                                               initial=5)):
        if cand > limit:
            break
        if n % cand == 0:
            count = 1
            n //= cand
            while n % cand == 0:
                count  = 1
                n //= cand
            ps.append(cand)
            exps.append(count)
            limit = isqrt(n)
    if n > 1:
        ps.append(n)
        exps.append(1)
    result = []
    for pows in product(*(range(exp   1) for exp in exps)):
        result.append(prod(p**e for p, e in zip(ps, pows)))
    return sorted(result)

>>> for i in range(22):
...     print(i, get_divisors(i))
0 [0]
1 [1]
2 [1, 2]
3 [1, 3]
4 [1, 2, 4]
5 [1, 5]
6 [1, 2, 3, 6]
7 [1, 7]
8 [1, 2, 4, 8]
9 [1, 3, 9]
10 [1, 2, 5, 10]
11 [1, 11]
12 [1, 2, 3, 4, 6, 12]
13 [1, 13]
14 [1, 2, 7, 14]
15 [1, 3, 5, 15]
16 [1, 2, 4, 8, 16]
17 [1, 17]
18 [1, 2, 3, 6, 9, 18]
19 [1, 19]
20 [1, 2, 4, 5, 10, 20]
21 [1, 3, 7, 21]