-- numbers kept as {log,{pow2,pow3,pow5}},
--  value is ~= exp(log), == (2^pow2)*(3^pow3)*(5^pow5)
enum LOG, POWS
enum POW2, POW3, POW5

function lnmin(sequence a, sequence b)
    return iff(a[LOG]<b[LOG]?a:b)
end function

constant ln1 = log(1), ln2 = log(2), ln3 = log(3), ln5 = log(5)

function hamming(integer N)
sequence h = repeat(0,N)
sequence x2 = {ln2,{1,0,0}},
         x3 = {ln3,{0,1,0}},
         x5 = {ln5,{0,0,1}}
integer i = 1, j = 1, k = 1
    h[1] = {ln1,{0,0,0}}
    for n=2 to N do
        h[n] = lnmin(x2,lnmin(x3,x5))
        sequence p = h[n][POWS]
        if p=x2[POWS] then i += 1 x2 = h[i] x2[LOG] += ln2 x2[POWS][POW2] += 1 end if
        if p=x3[POWS] then j += 1 x3 = h[j] x3[LOG] += ln3 x3[POWS][POW3] += 1 end if
        if p=x5[POWS] then k += 1 x5 = h[k] x5[LOG] += ln5 x5[POWS][POW5] += 1 end if
    end for
    return h[N]
end function

function hint(sequence hm)
-- (obviously not accurate above 53 bits on a 32-bit system, or 64 bits on a 64 bit system)
    sequence p = hm[POWS]
    return power(2,p[POW2])*power(3,p[POW3])*power(5,p[POW5])
end function

sequence s = {}
for i=1 to 20 do
    s = append(s,hint(hamming(i)))
end for
?s
?hint(hamming(1691))
?hint(hamming(1000000))
printf(1," %d (approx)\n",hint(hamming(1000000)))

include builtins\mpfr.e

function mpz_hint(sequence hm)
-- (as accurate as you like)
    integer {p2,p3,p5} = hm[POWS]
    mpz {tmp2,tmp3,tmp5} = mpz_inits(3)
    mpz_ui_pow_ui(tmp2,2,p2)
    mpz_ui_pow_ui(tmp3,3,p3)
    mpz_ui_pow_ui(tmp5,5,p5)
    mpz_mul(tmp3,tmp3,tmp5)
    mpz_mul(tmp2,tmp2,tmp3)
    return mpz_get_str(tmp2)
end function

?mpz_hint(hamming(1000000))