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