Is there an efficient algorithm for generating all numbers with n distinct prime factors in order?

Question

Bit of an x y problem here, so in full disclosure, I am attempting to find the next term of A152617, "Smallest number m such that m has exactly n distinct prime factors and sigma(m) has exactly n distinct prime factors."

But my main question is, how do you generate all numbers with n distinct prime factors in order? For small n, it's feasible to iterate through the natural numbers and factor each, but for the next term of A152617, getting to even the first number with 13 distinct prime factors (the 13th primorial) in this way is computationally intractable.

Looking at small cases, such as n=4, it seems to be some sort of combinatoric operation of primes and prime powers:

210 = 2 · 3 · 5 · 7
330 = 2 · 3 · 5 · 11
390 = 2 · 3 · 5 · 13
420 = 2^2 · 3 · 5 · 7
462 = 2 · 3 · 7 · 11
510 = 2 · 3 · 5 · 17
546 = 2 · 3 · 7 · 13
570 = 2 · 3 · 5 · 19
630 = 2 · 3^2 · 5 · 7
660 = 2^2 · 3 · 5 · 11
...

Normally, if there was a sporadically growing pattern with two degrees of freedom, to enumerate them in order, I would "batch" all of them in some range, and then sort. But with this problem, the prime powers replace primes, the prime powers can be any size, primes can become arbitrarily large, primes can be missing, so I'm unsure how to even go about choosing an order to enumerate the space due to all the degrees of freedom.

It seems there is an efficient PARI program on A033993 that might encode the algorithm I seek, but I have not yet been able to translate this code into a mental model of an algorithm, and found it more in my current ability to write this question. I've copied the code here for posterity:

(PARI) A246655(lim)=my(v=List(primes([2, lim\=1]))); for(e=2, logint(lim, 2), forprime(p=2, sqrtnint(lim, e), listput(v, p^e))); Set(v)
list(lim, pr=4)=if(pr==1, return(A246655(lim))); my(v=List(), pr1=pr-1, mx=prod(i=1, pr1, prime(i))); forprime(p=prime(pr), lim\mx, my(u=list(lim\p, pr1)); for(i=1, #u, listput(v, p*u[i]))); Set(v) \\ Charles R Greathouse IV, Feb 03 2023

Answer 1

I’d try it like this: Maintain a priority queue that you initialize with the $n$-th primorial, $p_n\#$. In each step, remove the least element from the queue, process it (in this case count the distinct prime factors of its divisor sum) and enter all numbers into the priority queue (if they’re not already there) that you can reach by increasing one of the prime factors to the next prime without changing the number of distinct prime factors. That is, you can increase a prime factor either if it has exponent $1$ and the next prime has exponent $0$ or if it has exponent $\gt1$ and the next prime has exponent $\gt0$. When you process a number of the form $2^kp_n\#$ (including the initial $p_n\#$), also enter $2^{k+1}p_n\#$ into the queue. (These are the smallest numbers with $n+k$ prime factors of which $n$ are distinct.)

I’m not sure exactly how quickly the queue will grow, but I think it should remain manageable up to quite large numbers – it seems plausible that you could reach the next number in that sequence without running out of space.

Answer 2

Here is an algorithm that works pretty well.

The idea is keeping a priority queue that stores the choices for the next smallest number.

Initially the priority queue contains only one entry that is the product of first $n$ primes.

Each time, we pop the priority queue to obtain the next smallest number, $num$. Then we push at most $2n$ numbers that are "barely larger" than $num$ into the priority queue as follows.

• for each prime factor $p$ of $num$:
  • push $num\times p$ into the priority queue.
  • If $p^2$ is not a factor of $num$ and the next prime after $p$, named $q$, is not a factor of $num$, push $\frac{num\times q}p$ into the priority queue.

---

Here is an implementation in Python.

For ease of implementation, we use a pair, (num, signature) to represent the number $num$, where signature tells the prime factorization of $num$.

import heapq
from math import isqrt, prod


def next_prime(k):
    """ return the smallest prime number that is larger than `k`"""
    p = k + 1
    while True:
        for i in range(2, isqrt(p) + 1):
            if p % i == 0:
                break
        else:
            return p
        p += 1


def get_sequence(n, start, end):
    """return the sequence of numbers that has `n` distinct primes,
    from the `start`-th smallest to the `end`-th smallest inclusive"""

    # collect first n primes
    primes = [2]
    for _ in range(n - 1):
        primes.append(next_prime(primes[-1]))

    # Each entry in hpq is `(num, signature)`, where `num` is the value of the number,
    # `signature` is a tuple of `(index, exponent)`, which represents `primes[index] ** exponent`
    # so that `num` is
    #   `primes[signature[0][0] ** signature[0][1]
    #  * primes[signature[1][0] ** signature[1][1]
    #  * ...
    #  * primes[signature[-1][0] ** signature[-1][1]`
    hpq = [(prod(primes[:n]), tuple((i, 1) for i in range(n)))]
    heapq.heapify(hpq)  # This does nothing, in fact.
    
    output = []
    cc = 0  # number of numbers we have got
    last_num = None
    while cc < end:
        num, signature = heapq.heappop(hpq)
        
        if num == last_num:
            continue
        cc += 1
        if cc >= start:
            output.append(num)
            last_num = num

        if signature[-1][0] == len(primes) - 1:
            primes.append(next_prime(primes[-1]))
        for i in range(n):
            # increase the exponent of `i`-th prime factor by 1
            heapq.heappush(hpq,
                           (num * primes[signature[i][0]],
                            tuple((signature[j][0], signature[j][1] + 1) if j == i else signature[j]
                                  for j in range(n))))
            # replace the `i`-th part by the next prime factor
            
            if signature[i][1] == 1 and (i < n -1 and signature[i + 1][0] > signature[i][0] + 1 or i == n - 1):
                heapq.heappush(hpq,
                               (num // pow(primes[signature[i][0]], signature[i][1]) * primes[signature[i][0] + 1],
                                tuple((signature[j][0] + 1, 1) if j == i else signature[j]
                                      for j in range(n))))
    return output


# print the first 12 smallest numbers that has 4 distinct prime factors
print(get_sequence(4, 1, 12))
# [210, 330, 390, 420, 462, 510, 546, 570, 630, 660, 690, 714]

Instead of returning a sequence, we can implement an iterator that given $n$, will yield the next smallest number that has $n$ distinct prime factors.

Moreover, we can use an auxiliary set to keep track of the numbers in the priority queue so that we will push a number into the queue only when it is not in the priority queue yet.

Furthermore, we can also keep a cap on the sum of all exponents. For $13$ distinct prime factors, the cap can be something like $2\log_2n+4$

Then we will be able to compute the first millions of numbers that have 13 distinct numbers without going out of memory.