List ALL prime-factorized natural numbers in ANY order

Question

This is a variant of List prime-factorized natural numbers up to N in ascending order, but the solutions can be very different.

Write a program that outputs prime factorizations of all natural numbers in any order. For example, the output could look like:

1:
4: 2^2
8: 2^3
2: 2^1
3: 3^1
5: 5^1
6: 2^1 3^1
21: 3^1 7^1
9: 3^2
...

Requirements:

1. You cannot just iterate over the numbers and factor each of them. This is too inefficient.
2. The output should be as in the example above: On each line the number and the list of its prime factors.
3. There is no given order in which the number appear in the list. However:
4. Each number must appear exactly once in the list (explain that your code satisfies this condition).

Edit - remark: Of course the program never finishes. But since (4) requires each number to appear in the list, for any n there is time t such that the list produced by the program within time t contains the factorization of n. (We don't care if t is comparable to the age of our universe or anything else.)

Note: This is not code-golf, it's code-challenge. This puzzle is about clever ideas, not code size. So please don't give the shortest solutions, give commented solutions that are easy to read and where the idea is easily understandable.

Answer 1

Generates the list by iterating over N. Each iteration generates all numbers which are the product of up-to-Nth powers of the first N primes. Starts with N=1 and keeps growing N. To avoid duplicates, each iteration makes sure not to generate anything generated in previous iterations. Those are exactly the numbers that don't include the Nth prime and which don't have at least one exponent equal to N.

import itertools

print '1 ='

n = 0
primes = []
while 1:
  n += 1

  # generate next prime                                                                                                                                    
  q = primes[-1] if primes else 1
  while 1:
    q += 1
    if all(q % p != 0 for p in primes):
      primes.append(q)
      break

  # generate all exponents of at most n of first n primes                                                                                                            
  for e in itertools.product(*(xrange(n+1) for i in xrange(n))):
    if max(e) == n or e[-1] > 0:  # we didn't do this in iteration n-1                                                                                  
      print reduce(lambda x, y: x * y, (primes[i]**e[i] for i in xrange(n))), '=', ' '.join('%d^%d' % (primes[i], e[i]) for i in xrange(n))

Answer 2

Haskell

import Data.List

merge :: [[a]] -> [[a]] -> [[a]]
merge [] ys = ys
merge (x:xs) ys = x : merge ys xs

atmostonefromeach :: [[a]]-> [[a]]
atmostonefromeach = foldr weave undefined

weave :: [a] -> [[a]] -> [[a]]
weave pwrs odds = merge (map (:[]) pwrs) $
foldr (\odd -> merge (odd: (map (:odd) (pwrs)))) undefined odds

primes = nubBy (\x y->gcd x y > 1) [2..]

shw factors = show (product (map (\(p,k)->p^k) factors)) ++ ":" ++ concatMap(\(p,k) -> ' ': show p ++ '^': show k) factors

powers p = map (\x->(p,x)) [1..]

main = mapM_ (putStrLn . shw) $ ([]:) $ atmostonefromeach $ map powers primes

output:

1:
2: 2^1
3: 3^1
4: 2^2
5: 5^1
8: 2^3
6: 2^1 3^1

...

The program constructs, from the list Q=[q0,q1,..] of composite numbers with all prime factors > p, the list of those with primes factors >= p, as merge(P, merge((map(q0*)P1), merge(map(q1*)P1),... where P=[p^1,p^2,...] and P1=[p^0,p^1,..].

The name odds in function weave describes the case p==2, where Q is the set of all odd numbers > 1.