58 lines
1.6 KiB
Text
58 lines
1.6 KiB
Text
-- 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))
|