How to iterate the Hardy-Ramanujan integers quickly
I'm very happy that the site linked me to How can I generate first n elements of the sequence 3^i * 5^j * 7^k?. It was a critical stepping stone, which enabled me to solve my own question after lots of thinking.

First I want to elaborate on that answer a bit, as it took me a while to understand and code it.
We want to iterate $3^i 5^j 7^k$.
The next element must be one of $3 x, 5 y, 7 z$, where $x, y, z$ is a previous number in the sequence. This is because $x &lt; 3 x, y &lt; 5 y, z &lt; 7 z$, and $3 x, 5 y, 7 z$ satisfy the constraints.
For $x$, we start with the first element in the sequence. We increment it's position whenever $3 x$ is the smallest out of $3 x, 5 y, 7 z$. To see why, we've already included $3 x$ in the sequence, for all $x$s in the sequence so far. So the only possible $3 x$ that can be inserted in the sequence, is if $x$ is the new element we just inserted.
Similarly for $y$ and $z$.
The following code iterates this sequence:
def main():
    x = 1
    y = 1
    z = 1
    S = []
    x_iter = iter(S)
    y_iter = iter(S)
    z_iter = iter(S)
    for _ in range(20):
        m = min(3 * x, 5 * y, 7 * z)
        S.append(m)
        if m == 3 * x:
            x = next(x_iter)
        if m == 5 * y:
            y = next(y_iter)
        if m == 7 * z:
            z = next(z_iter)
    print(S)


The Hardy-Ramanujan Integers can be defined as the integers $2^{e_1} 3^{e_2} 5^{e_3} \cdots$, s.t. $e_1 \geqslant e_2 \geqslant e_3 \geqslant \cdots \geqslant 0$.
It seems that these two problems are related, and indeed they are the same, if we re-write the Hardy-Ramanujan Integers by removing the decreasing exponents constraint, as $2^{e_1'} (2^{e_2'} 3^{e_2'}) (2^{e_3'} 3^{e_3'} 5^{e_3'}) \cdots$.
Now the only issue is that compared to the previous problem, our list of bases is infinite. But note that a new prime $p$ can only be included in the sequence, if it's smallest form, $2^1 3^1 \cdots p^1$, is less than the next sequence element, produced with primes $&lt; p$. So we only need to introduce a new prime when this occurs.
Before this occurs, The exponent of $p$ is 0. Any prime $&gt; p$ will give a sequence element bigger than $2^1 3^1 \cdots p^1$, so does not yet need to be considered.
This gives the following code:
import math

from sympy import nextprime

def main():
    S = [1]
    primes = [2]
    next_prime = nextprime(primes[0])
    # The smallest Hardy-Ramanujan integer that includes `next_prime`
    next_prime_product = primes[0] * next_prime
    candidates = [1]
    candidate_S_indexes = [0]
    for _ in range(20):
        m_options = [
            math.prod(primes[:i + 1]) * candidate
            for i, candidate in enumerate(candidates)
        ]
        m = min(m_options)
        if next_prime_product &lt; m:
            # Add a new prime &amp; candidate
            m = next_prime_product
            primes.append(next_prime)
            next_prime = nextprime(next_prime)
            next_prime_product *= next_prime
            candidates.append(m)
            candidate_S_indexes.append(len(S))
        S.append(m)
        for i, m_option in enumerate(m_options):
            if m_option == m:
                candidates[i] = S[candidate_S_indexes[i] + 1]
                candidate_S_indexes[i] += 1
    print(S)


P.S. I saw some answer in Haskell in German before at:

A good algorithm for A094348?
A025487Haskell - seqcomp

but could not understand the German at all, nor the Google-translated version, nor the Haskell code. But I'm quite satisfied with my algorithm here. It feels pretty optimal.