\ Structure and Interpretation of Computer Programs, Exercise 3.56
\ enumerate, in ascending order with no repetitions, all positive
\ integers with no prime factors other than 2, 3, or 5.
\ This version represents Hamming numbers with their logarithm dualis
\ (in fixed-point form scaled by 2^42; testing using another, slower
\ representation has shown that these truncated logarithms don't
\ collide up to 10^9.
\ manipulating and computing with Hamming numbers:
' + alias h* ( h1 h2 -- h )
1 42 lshift constant scale-factor
: h. { h -- }
2e h 0 d>f scale-factor 0 d>f f/ f** 15 0 10 f.rdp space ;
\ the following numbers have been produced with bc -l as follows
1 42 lshift constant ldscale2
6970738796527 constant ldscale3 \ 2^42*l(3)/l(2) (rounded down)
10211947756754 constant ldscale5 \ 2^42*l(5)/l(2) (rounded up)
' u<= alias h<= ( h1 h2 -- f )
' umin alias hmin ( h1 h2 -- h )
\ actual algorithm
0 value seq
variable seqlast 0 seqlast !
: lastseq ( -- u )
\ last stored number in the sequence
seq seqlast @ th @ ;
: genseq ( h1 "name" -- )
\ h1 is the factor for the sequence
create , 0 , \ factor and index of element used for last return
does> ( -- u2 )
\ u2 is the next number resulting from multiplying h1 with numbers
\ in the sequence that is larger than the last number in the
\ sequence
dup @ lastseq { h1 l } cell+ dup @ begin ( index-addr index )
seq over th @ h1 h* dup l h<= while
drop 1+ repeat
>r swap ! r> ;
scale-factor genseq s2
6970738796527 genseq s3
10211947756754 genseq s5
: nextseq ( -- )
s2 s3 hmin s5 hmin , 1 seqlast +! ;
: nthseq ( u1 -- h )
\ the u1 th element in the sequence
dup seqlast @ u+do
nextseq
loop
1- 0 max cells seq + @ ;
: .nseq ( u1 -- )
dup seqlast @ u+do
nextseq
loop
0 u+do
seq i th @ h.
loop ;
here to seq
0 , \ that's 1
\ usage: 1000 .nseq
\ 1000 nthseq .