Efficiently determine every factor of a given natural number

Question

It's not hard to find all of the factors of a given natural number. Below is a pretty fast python function that can do just that:
(Modified version of Steinar Lima's function, https://stackoverflow.com/a/19578818/5946921)

import itertools as it

def factor(n):
    step = 2 if n % 2 else 1
    factors = ([i, n//i] for i in range(1, int(sqrt(n))+1, step) if not n % i)
    return sorted(set(it.chain.from_iterable(factors)))

I am looking to make a function that can do this more efficiently for very large numbers. I have come up with the following function that does just that:

import itertools as it

def large_factor(n):
    primes = []
    ranges = []
    for prime, factors in it.groupby(primeFactor(n), lambda x: x):
        primes.append(prime)
        ranges.append(range(len(tuple(factors)) + 1))
    primes = primes[::-1]
    ranges = ranges[::-1]
    z = it.product(*(x for x in ranges))
    return sorted(product(a ** b for a, b in zip(primes, q)) for q in z)

I am looking for ways to improve this function. I want it to return the factors in ascending order, without duplicates. Originally, I was hoping I'd be able to have the entire thing be a lazy function (no eager evaluation at all). Unfortunately, I haven't been able to figure out a good way to do this. My approach right now is to get the prime factors from n, and a list of ranges [0, #primes factors + 1). These lists are then reversed, as this makes sorting easier later. The list of ranges is used to generate tuples that contain the powers that are applied to the prime numbers in the last line.

The following two functions are used in the large_factor() function:

import math

def primeFactor(n):
    ''' Generates prime factors of n'''
    if n <= 1: return
    prime = next((x for x in range(2, math.ceil(math.sqrt(n))+1) if n%x == 0), n)
    yield prime
    yield from primeFactor(n//prime)

---

from functools import reduce
from operator import mul

def product(iterable):
    ''' The product of the items in n '''
    return reduce(mul, iterable)

I've tested these functions with timeit, and found that factor() was significantly faster on small values on n, but on large values of n, large_factor() was significantly faster. I'm looking for any kind of optimization that can make large_factor() faster on smaller, or (especially) larger numbers, as well as any other suggestions you might have.

Answer 1

def primeFactor(n):
    ''' Generates prime factors of n'''
    if n <= 1: return
    prime = next((x for x in range(2, math.ceil(math.sqrt(n))+1) if n%x == 0), n)
    yield prime
    yield from primeFactor(n//prime)

For a start, consider that after yielding prime this can never yield a number less than prime, but the recursive call still starts again from 2. So one optimisation would be

def primeFactor(n):
    ''' Generates prime factors of n'''
    return primeFactorInner(n, 2)

def primeFactorInner(n, min):
    if n < min: return
    prime = next((x for x in range(min, math.ceil(math.sqrt(n))+1) if n%x == 0), n)
    yield prime
    yield from primeFactorInner(n//prime, prime)

---

Next, that range with an implicit step of 1 is a bit overkill. Once you've finished returning 2s you could just step through odd numbers. At this point I'm tempted to say that it would be convenient to change the return value to already do the grouping, and have an auxiliary method to get the power of a prime in n.

def primePower(n, p): return 0 if n%p > 0 else 1 + primePower(n//p, p)

Then we can write

def primeFactorInner(n, min):
    if n < min: return
    if min == 2:
        a = primePower(n, 2)
        if a > 0:
            yield (2, a)
            yield from primeFactorInner(n // (2 ** a), 3)
            return
    prime = next((x for x in range(min, math.ceil(math.sqrt(n))+1, 2) if n%x == 0), n)
    a = primePower(n, prime)
    yield (prime, a)
    yield from primeFactorInner(n // (prime ** a), prime + 2)

---

Now, if we carry that argument to extremes we would only want to loop over primes. So we could split into two phases:

for p in enumeratePrimes(math.floor(sqrt(n)) + 1):
    a = primePower(n, p)
    if (a > 0):
        yield (p, a)
        n = n // (p ** a)
        if n == 1:
            # Early abort
            return
if n > 1:
    yield (n, 1)

enumeratePrimes would probably best be implemented with a simple sieve.

---

I want it to return the factors in ascending order, without duplicates. Originally, I was hoping I'd be able to have the entire thing be a lazy function (no eager evaluation at all). Unfortunately, I haven't been able to figure out a good way to do this.

The short answer is a priority queue. Initially you insert 1. When you pop and yield a factor, you push certain multiples of the factor. The mildly interesting part is ensuring that you only push each number once. I've found the Java code where I did a variant on this, and ported it to Python:

import heapq

def divisors(n):
    primes = [(1, 1)] + list(primeFactors(n))
    q = [(1, 0, 1)]
    while len(q) > 0:
        # d is the divisor
        # i is the index of its largest "prime" in primes
        # a is the exponent of that "prime"
        (d, i, a) = heapq.heappop(q)
        yield d
        if a < primes[i][1]:
            heapq.heappush(q, (d * primes[i][0], i, a + 1))
        if i + 1 < len(primes):
            heapq.heappush(q, (d * primes[i + 1][0], i + 1, 1))
            # The condition i > 0 is to avoid duplicates arising because
            # d == d // primes[0][0]
            if i > 0 and a == 1:
                heapq.heappush(q, (d // primes[i][0] * primes[i + 1][0], i + 1, 1))

I'm not really a Pythonista, so although I suspect that you can make this even lazier by avoiding the list(), I'll leave that to someone who knows the language better than me.