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