How to iterate the Hardy-Ramanujan integers quickly

Question

The Hardy-Ramanujan integers, A025487 - OEIS, are integers which when factorized, have their exponents for all the primes starting from 2, in decreasing (not strictly) order. The first few terms are:

$$ \begin{array}{lll} 1 & = & 1\\ 2^1 & = & 2\\ 2^2 & = & 4\\ 2^1 \times 3^1 & = & 6\\ 2^3 & = & 8\\ 2^2 \times 3^1 & = & 12\\ 2^4 & = & 16\\ 2^3 \times 3^1 & = & 24\\ 2^1 \times 3^1 \times 5^1 & = & 30\\ 2^5 & = & 32\\ 2^2 \times 3^2 & = & 36\\ 2^4 \times 3^1 & = & 48\\ 2^2 \times 3^1 \times 5^1 & = & 60\\ \vdots \end{array} $$

As you can see, the exponents don't really follow much of an order that I can see. I thought perhaps the exponents increased in some way, but we have $288 = 2^5 \times 3^2$ shortly followed by $480 = 2^5 \times 3^1 \times 5^1$.

Is there a method to iterate these integers quickly?

I can only see a simple way of listing all of them below some upper bound $u$. I find the maximum exponent for $2$ by $\lfloor \log_2(u) \rfloor$, find the maximum # of primes, and iterate all decreasing sequences of exponents $\leq \lfloor \log_2(u) \rfloor$.

Answer 1

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 < 3 x, y < 5 y, z < 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 $< p$. So we only need to introduce a new prime when this occurs.

Before this occurs, The exponent of $p$ is 0. Any prime $> 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 < m:
            # Add a new prime & 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.

Answer 2

For a prime p, let hm(p) be the sequence of Hardy-Ramanujan numbers with no prime factor greater than p.

hm(2) is the sequence 1, 2, 4, 8, 16, 32 etc.

hm(3) is created by merging hm(2) and 6 x hm(3).

hm(5) is created by merging hm(3) and 30 x hm(5)

and so on. The numbers up to 210-1 are in hm(5), the numbers up to 2310-1 are in hm(7), the numbers up to 30,030-1 are in hm(11) etc.