Is there an efficient algorithm to enumerate the factors of a number n, in ascending order, without sorting? By “efficient” I mean:
The algorithm avoids a brute-force search for divisors by starting with the prime-power factorization of n.
The runtime complexity of the algorithm is O(d log₂ d) or better, where d is the divisor count of n.
The spatial complexity of the algorithm is O(d).
The algorithm avoids a sort operation. That is, the factors are produced in order rather than produced out of order and sorted afterward. Although enumerating using a simple recursive approach and then sorting is O(d log₂ d), there is a very ugly cost for all the memory accesses involved in sorting.
A trivial example is n = 360 = 2³ × 3² × 5, which has d=24 factors: { 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 18, 20, 24, 30, 36, 40, 45, 60, 72, 90, 120, 180, 360 }.
A more serious example is n = 278282512406132373381723386382308832000 = 2⁸ × 3⁴ × 5³ × 7² × 11² × 13² × 17 × 19 × 23 × 29 × 31 × 37 × 41 × 43 × 47 × 53 × 59 × 61 × 67 × 71 × 73 × 79, which has d=318504960 factors (obviously too many to list here!). This number, incidentally, has the greatest number of factors of any number up to 2^128.
I could swear that I saw a description of such algorithm a few weeks ago, with sample code, but now I can't seem to find it anywhere. It used some magic trick of maintaining a list of progenitor indexes in the output list for each prime factor. (Update: I was confusing factor generation with Hamming numbers, which operate similarly.)
I ended up using a solution which is O(d) in runtime, has extremely low memory overhead creating the O(d) output in-place, and is significantly faster than any other method I am aware of. I've posted this solution as answer, with C source code. It is a heavily optimized, simplified version of a beautiful algorithm presented here in Haskell by another contributor, Will Ness. I've selected Will's answer as the accepted answer, as it provided a very elegant solution that matched all the requirements as originally stated.
We can merge streams of multiples, produced so there are no duplicates in the first place.
Starting with the list [1], for each unique prime factor p, we multiply the list by p iteratively k times (where k is the multiplicity of p) to get k new lists, and merge them together with that incoming list.
In Haskell,
ordfactors n = -- <----------------------<---<---<-----\
foldr g [1] -- g (19,1) (g (7,1) (g (3,2) (g (2,3) [1])))
. reverse -- [(19,1),(7,1),(3,2),(2,3)] ^
. primePowers -- [(2,3),(3,2),(7,1),(19,1)] |
$ n -- 9576 --->--->--->----/
where
g (p,k) xs = mergeAsTree
. take (k+1) -- take 3 [a,b,c,d..] = [a,b,c]
. iterate (map (p*)) -- iterate f x = [x, y, z,...]
$ xs -- where { y=f x; z=f y; ... }
{- g (2,3) [1] = let a0 = [1]
a1 = map (2*) a0 -- [2]
a2 = map (2*) a1 -- [4]
a3 = map (2*) a2 -- [8]
in mergeAsTree [ a0, a1, a2, a3 ] -- xs2 = [1,2,4,8]
g (3,2) xs2 = let b0 = xs2 -- [1,2,4,8]
b1 = map (3*) b0 -- [3,6,12,24]
b2 = map (3*) b1 -- [9,18,36,72]
in mergeAsTree [ b0, b1, b2 ] -- xs3 = [1,2,3,4,6,8,9,12,...]
g (7,1) xs3 = mergeAsTree [ xs3, map (7*) xs3 ] -- xs4 = [1,2,3,4,6,7,8,9,...]
g (19,1) xs4 = mergeAsTree [ xs4, map (19*) xs4 ]
-}
mergeAsTree xs = foldi merge [] xs -- [a,b] --> merge a b
-- [a,b,c,d] --> merge (merge a b) (merge c d)
foldi f z [] = z
foldi f z (x:xs) = f x (foldi f z (pairs f xs))
pairs f (x:y:t) = f x y : pairs f t
pairs _ t = t
foldi arranges the binary merges (which presuppose that there's no duplicates in the streams being merged, for streamlined operation) in a tree-like fashion for speed; whereas foldr works in linear fashion.
primePowers n | n > 1 = -- 9576 => [(2,3),(3,2),(7,1),(19,1)]
map (\xs -> (head xs,length xs)) -- ^
. group -- [[2,2,2],[3,3],[7],[19]] |
. factorize -- [2,2,2,3,3,7,19] |
$ n -- 9576 --->--->--->---/
group is a standard function that groups adjacent equals in a list together, and factorize is a well-known trial-division algorithm that produces prime factors of a number in non-decreasing order:
factorize n | n > 1 = go n (2:[3,5..]) -- 9576 = 2*2*2*3*3*7*19
where -- produce prime factors only
go n (d:t)
| d*d > n = [n]
| r == 0 = d : go q (d:t)
| otherwise = go n t
where
(q,r) = quotRem n d
factorize _ = []
(.) is a functional composition operator, chaining the output of its argument function on its right as input into the one on its left. It (and $) can be read aloud as "of".