# How to iterate the Hardy-Ramanujan integers quickly

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$$.

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)