Hamming numbers
From Rosetta Code
Hamming numbers is a programming task. Visitors like you are encouraged to solve it according to the task description, using any language they may happen to know.
Hamming numbers are numbers of the form
-
.
Generate the sequence of Hamming numbers, in increasing order. In particular:
- Show the first twenty Hamming numbers.
- Show the 1691st Hamming number (the last one below 231).
- Show the one millionth Hamming number (if the language – or a convenient library – supports arbitrary-precision integers).
References
- wp:Hamming_numbers
- wp:Smooth_number
- Hamming problem from Dr. Dobb's CodeTalk
Contents |
[edit] ALGOL 68
Works with: Algol 68 Genie 1.19.0
Hamming numbers are generated in a trivial iterative way as in the Python version below. This program keeps the series needed to generate the numbers as short as possible using flexible rows; on the downside, it spends considerable time on garbage collection.
PR precision=100 PR
MODE SERIES = FLEX [1 : 0] UNT, # Initially, no elements #
UNT = LONG LONG INT; # A 100-digit unsigned integer #
PROC hamming number = (INT n) UNT: # The n-th Hamming number #
CASE n
IN 1, 2, 3, 4, 5, 6, 8, 9, 10, 12 # First 10 in a table #
OUT # Additional operators #
OP MIN = (INT i, j) INT: (i < j | i | j), MIN = (UNT i, j) UNT: (i < j | i | j);
PRIO MIN = 9;
OP LAST = (SERIES h) UNT: h[UPB h]; # Last element of a series #
OP +:= = (REF SERIES s, UNT elem) VOID:
# Extend a series by one element, only keep the elements you need #
(INT lwb = (i MIN j) MIN k, upb = UPB s;
REF SERIES new s = HEAP FLEX [lwb : upb + 1] UNT;
(new s[lwb : upb] := s[lwb : upb], new s[upb + 1] := elem);
s := new s
);
# Determine the n-th hamming number iteratively #
SERIES h := 1, # Series, initially one element #
UNT m2 := 2, m3 := 3, m5 := 5, # Multipliers #
INT i := 1, j := 1, k := 1; # Counters #
TO n - 1
DO h +:= (m2 MIN m3) MIN m5;
(LAST h = m2 | m2 := 2 * h[i +:= 1]);
(LAST h = m3 | m3 := 3 * h[j +:= 1]);
(LAST h = m5 | m5 := 5 * h[k +:= 1])
OD;
LAST h
ESAC;
FOR k TO 20
DO print ((whole (hamming number (k), 0), blank))
OD;
print ((newline, whole (hamming number (1 691), 0)));
print ((newline, whole (hamming number (1 000 000), 0)))
Sample output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] C
C does not natively support long numbers. The following will print the first 20 and the 1690th hamming numbers#include <stdio.h>
typedef unsigned vlong;
typedef void (*Callback)(int,vlong);
void vextend( vlong *list, int *ix, vlong val)
{
vlong *val_l = list + *ix-1;
vlong *val_h = list + *ix;
while (*val_l > val) --val_l;
if (*val_l == val) return;
val_l = list + *ix-1;
if (*ix<8000-1)
*ix+=1;
else
if (val >= *val_h) return;
while ( *val_l > val)
*val_h-- = *val_l--;
*val_h = val;
}
unsigned MAX_UL = 0XFFFFFFFF;
void hamming( int n, Callback cbak)
{
vlong explist[8000];
int ix1, ix2;
int k2=0,k3=0,k5=0;
vlong v;
explist[0] = 1;
ix1 = 0; ix2 = 1;
while ((ix1 < ix2) && (ix1 < n) ) {
v = explist[ix1];
(*cbak)(ix1,v);
if ( v < MAX_UL/(unsigned)2)
vextend(explist, &ix2, 2*v);
if ( v < MAX_UL/(unsigned)3)
vextend(explist, &ix2, 3*v);
if ( v < MAX_UL/(unsigned)5)
vextend(explist, &ix2, 5*v);
ix1++;
}
}
void hprint(int ix, vlong val)
{
if ((ix < 20) || (ix=1689) || (ix==1847))
printf("h[%4d]= %u\n", 1+ix,val);
}
int main(int argc, char *argv[])
{
printf("%u\n", MAX_UL);
hamming(5002, &hprint);
return 0;
}
[edit] Clojure
This version implements Dijkstra's merge solution, so is closely related to the Haskell version.
(defn smerge [xs ys]
(lazy-seq
(let [x (first xs),
y (first ys),
[z xs* ys*]
(cond
(< x y) [x (rest xs) ys]
(> x y) [y xs (rest ys)]
:else [x (rest xs) (rest ys)])]
(cons z (smerge xs* ys*)))))
(defn smerge3 [xs ys zs]
(smerge xs (smerge ys zs)))
(defn map*n [n ks] (map #(* n %) ks))
(def hamming
(lazy-seq
(cons 1 (smerge3 (map*n 2 hamming) (map*n 3 hamming) (map*n 5 hamming)))))
Note that this version uses a lot of space and time after calculating a few hundred thousand elements of the sequence. This is no doubt due to its "holding on to the head": it maintains the entire generated sequence in memory.
[edit] Factor
Translation of: Scala
USING: accessors deques dlists fry kernel make math math.order
;
IN: rosetta.hamming
TUPLE: hamming-iterator 2s 3s 5s ;
: <hamming-iterator> ( -- hamming-iterator )
hamming-iterator new
1 1dlist >>2s
1 1dlist >>3s
1 1dlist >>5s ;
: enqueue ( n hamming-iterator -- )
[ [ 2 * ] [ 2s>> ] bi* push-back ]
[ [ 3 * ] [ 3s>> ] bi* push-back ]
[ [ 5 * ] [ 5s>> ] bi* push-back ] 2tri ;
: next ( hamming-iterator -- n )
dup [ 2s>> ] [ 3s>> ] [ 5s>> ] tri
3dup [ peek-front ] tri@ min min
[
'[
dup peek-front _ =
[ pop-front* ] [ drop ] if
] tri@
] [ swap enqueue ] [ ] tri ;
: next-n ( hamming-iterator n -- seq )
swap '[ _ [ _ next , ] times ] { } make ;
: nth-from-now ( hamming-iterator n -- m )
1 - over '[ _ next drop ] times next ;
<hamming-iterator> 20 next-n . <hamming-iterator> 1691 nth-from-now . <hamming-iterator> 1000000 nth-from-now .
Translation of: Haskell
Lazy lists aren quite slow in Factor, but still.
USING: combinators fry kernel lists lists.lazy locals math ;
IN: rosetta.hamming-lazy
:: sort-merge ( xs ys -- result )
xs car :> x
ys car :> y
{
{ [ x y < ] [ [ x ] [ xs cdr ys sort-merge ] lazy-cons ] }
{ [ x y > ] [ [ y ] [ ys cdr xs sort-merge ] lazy-cons ] }
[ [ x ] [ xs cdr ys cdr sort-merge ] lazy-cons ]
} cond ;
:: hamming ( -- hamming )
f :> h!
[ 1 ] [
h 2 3 5 [ '[ _ * ] lazy-map ] tri-curry@ tri
sort-merge sort-merge
] lazy-cons h! h ;
20 hamming ltake list>array . 1690 hamming lnth . 999999 hamming lnth .
[edit] Haskell
hamming = 1 : map (2*) hamming `merge` map (3*) hamming `merge` map (5*) hamming
where merge (x:xs) (y:ys)
| x < y = x : xs `merge` (y:ys)
| x > y = y : (x:xs) `merge` ys
| otherwise = x : xs `merge` ys
main = do
print $ take 20 hamming
print $ hamming !! 1690
print $ hamming !! 1000000
[edit] J
Solution:
A concise tacit expression using a (right) fold:
hamming=: {. (/:~@~.@], 2 3 5 * {)/ @ (1x,~i.@-)
Example usage:
hamming 20
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36
{: hamming 1691
2125764000
For the millionth (and billionth (1e9)) Hamming number see the hn verb from Hamming Number essay on the J wiki.
Explanation:
I'll explain this J-sentence by dividing it in three parts from left to right omitting the leftmost {.:
- sort and remove duplicates
/:~@~.@]
- produce (the next) 3 elements by selection and multiplication:
2 3 5 * {
or
2 3 5 * LHA { RHA
- the RH part forms an array of descending indices and the initial Hamming number 1
(1x,~i.@-)
e.g. if we want the first 5 Hamming numbers, it produces the array:
4 3 2 1 0 1
[edit] JavaScript
Works with: JavaScript version 1.7 Works with: Firefox version 2
Translation of: Ruby
This does not calculate the 1,000,000th Hamming number.
Note the use of for (x in obj) to iterate over the properties of an object, versus for each (y in obj) to iterate over the values of the properties of an object.
function hamming() {
var queues = {2: [], 3: [], 5: []};
var base;
var next_ham = 1;
while (true) {
yield next_ham;
for (base in queues) {queues[base].push(next_ham * base)}
next_ham = [ queue[0] for each (queue in queues) ].reduce(function(min, val) {
return Math.min(min,val)
});
for (base in queues) {if (queues[base][0] == next_ham) queues[base].shift()}
}
}
var ham = hamming();
var first20=[], i=1;
for (; i <= 20; i++)
first20.push(ham.next());
print(first20.join(', '));
print('...');
for (; i <= 1690; i++)
ham.next();
print(i + " => " + ham.next());
output
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36 ... 1691 => 2125764000
[edit] Logo
to init.ham
; queues
make "twos [1]
make "threes [1]
make "fives [1]
end
to next.ham
localmake "ham first :twos
if less? first :threes :ham [make "ham first :threes]
if less? first :fives :ham [make "ham first :fives]
if equal? :ham first :twos [ignore dequeue "twos]
if equal? :ham first :threes [ignore dequeue "threes]
if equal? :ham first :fives [ignore dequeue "fives]
queue "twos :ham * 2
queue "threes :ham * 3
queue "fives :ham * 5
output :ham
end
init.ham
repeat 20 [print next.ham]
repeat 1690-20 [ignore next.ham]
print next.ham
[edit] Lua
function hiter()
hammings = {1}
prev, vals = {1, 1, 1}
index = 1
local function nextv()
local n, v = 1, hammings[prev[1]]*2
if hammings[prev[2]]*3 < v then n, v = 2, hammings[prev[2]]*3 end
if hammings[prev[3]]*5 < v then n, v = 3, hammings[prev[3]]*5 end
prev[n] = prev[n] + 1
if hammings[index] == v then return nextv() end
index = index + 1
hammings[index] = v
return v
end
return nextv
end
j = hiter()
for i = 1, 20 do
print(j())
end
n, l = 0, 0
while n < 2^31 do n, l = j(), n end
print(l)
[edit] MUMPS
Hamming(n) New count,ok,next,number,which
For which=2,3,5 Set number=1
For count=1:1:n Do
. Set ok=0 Set:count<21 ok=1 Set:count=1691 ok=1 Set:count=n ok=1
. Write:ok !,$Justify(count,5),": ",number
. For which=2,3,5 Set next(number*which)=which
. Set number=$Order(next(""))
. Kill next(number)
. Quit
Quit
Do Hamming(2000)
1: 1
2: 2
3: 3
4: 4
5: 5
6: 6
7: 8
8: 9
9: 10
10: 12
11: 15
12: 16
13: 18
14: 20
15: 24
16: 25
17: 27
18: 30
19: 32
20: 36
1691: 2125764000
2000: 8062156800
[edit] Oz
Translation of: Haskell
declare
fun lazy {HammingFun}
1|{FoldL1 [{MultHamming 2} {MultHamming 3} {MultHamming 5}] LMerge}
end
Hamming = {HammingFun}
fun {MultHamming N}
{LMap Hamming fun {$ X} N*X end}
end
fun lazy {LMap Xs F}
case Xs
of nil then nil
[] X|Xr then {F X}|{LMap Xr F}
end
end
fun lazy {LMerge Xs=X|Xr Ys=Y|Yr}
if X < Y then X|{LMerge Xr Ys}
elseif X > Y then Y|{LMerge Xs Yr}
else X|{LMerge Xr Yr}
end
end
fun {FoldL1 X|Xr F}
{FoldL Xr F X}
end
in
{ForAll {List.take Hamming 20} System.showInfo}
{System.showInfo {Nth Hamming 1690}}
{System.showInfo {Nth Hamming 1000000}}
[edit] Python
from itertools import islice
def hamming2():
'''\
This version is based on a snippet from:
http://dobbscodetalk.com/index.php?option=com_content&task=view&id=913&Itemid=85
When expressed in some imaginary pseudo-C with automatic
unlimited storage allocation and BIGNUM arithmetics, it can be
expressed as:
hamming = h where
array h;
n=0; h[0]=1; i=0; j=0; k=0;
x2=2*h[ i ]; x3=3*h[j]; x5=5*h[k];
repeat:
h[++n] = min(x2,x3,x5);
if (x2==h[n]) { x2=2*h[++i]; }
if (x3==h[n]) { x3=3*h[++j]; }
if (x5==h[n]) { x5=5*h[++k]; }
'''
h = 1
_h=[h] # memoized
multipliers = (2, 3, 5)
multindeces = [0 for i in multipliers] # index into _h for multipliers
multvalues = [x * _h[i] for x,i in zip(multipliers, multindeces)]
yield h
while True:
h = min(multvalues)
_h.append(h)
for (n,(v,x,i)) in enumerate(zip(multvalues, multipliers, multindeces)):
if v == h:
i += 1
multindeces[n] = i
multvalues[n] = x * _h[i]
# cap the memoization
mini = min(multindeces)
if mini >= 1000:
del _h[:mini]
multindeces = [i - mini for i in multindeces]
#
yield h
Sample output:
>>> list(islice(hamming2(), 20)) [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] >>> list(islice(hamming2(), 1689, 1690)) [2123366400] >>> list(islice(hamming2(), 999999, 1000000)) [519312780448388736089589843750000000000000000000000000000000000000000000000000000000]
[edit] Alternate version using "Cyclic Iterators"
The original author is Raymond Hettinger and the code was first published here under the MIT license.
This implementation is the fastest of the Python implementations and is quite memory efficient.
from itertools import tee, chain, groupby
from heapq import merge
def raymonds_hamming():
# Generate "5-smooth" numbers, also called "Hamming numbers"
# or "Regular numbers". See: http://en.wikipedia.org/wiki/Regular_number
# Finds solutions to 2**i * 3**j * 5**k for some integers i, j, and k.
def no_repeats(it):
for k,g in groupby(it):
yield k
def deferred_output():
for i in output:
yield i
result, p2, p3, p5 = tee(deferred_output(), 4)
m2 = (2*x for x in p2) # multiples of 2
m3 = (3*x for x in p3) # multiples of 3
m5 = (5*x for x in p5) # multiples of 5
merged = merge(m2, m3, m5)
combined = chain([1], merged) # prepend a starting point
output = no_repeats(combined) # eliminate duplicates
return result
Results are the same as before.
[edit] R
Recursively find the Hamming numbers below 231. Works for larger limits, however to find the 1000000th value, one needs guess the correct limit
hamming=function(hamms,limit)
{
tmp=hamms
for(h in c(2,3,5))
{
tmp=c(tmp,h*hamms)
}
tmp=unique(tmp[tmp<=limit])
if(length(tmp)>length(hamms))
{
hamms=hamming(tmp,limit)
}
hamms
}
sort(hamming(1,limit=2^31)[-1])
[edit] Ruby
Translation of: Scala
This method Works with: Ruby version 1.8.6
require 'generator'
# the Hamming number generator
hamming = Generator.new do |generator|
next_ham = 1
queues = { 2 => [], 3 => [], 5 => [] }
loop do
generator.yield next_ham
[2,3,5].each {|m| queues[m] << (next_ham * m)}
next_ham = [2,3,5].collect {|m| queues[m][0]}.min
[2,3,5].each {|m| queues[m].shift if queues[m][0] == next_ham}
end
end
This method, Works with: Ruby version 1.9.1, does not require a library module.
hamming = Enumerator.new do |yielder|
next_ham = 1
queues = { 2 => [], 3 => [], 5 => [] }
loop do
yielder << next_ham # or: yielder.yield(next_ham)
[2,3,5].each {|m| queues[m]<< (next_ham * m)}
next_ham = [2,3,5].collect {|m| queues[m][0]}.min
[2,3,5].each {|m| queues[m].shift if queues[m][0]== next_ham}
end
end
And the "main" part of the task
start = Time.now
idx = 1
hamming.each do |ham|
case idx
when (1..20), 1691
p [idx, ham]
when 1_000_000
p [idx, ham]
break
end
idx += 1
end
puts "elapsed: #{Time.now - start} seconds"
output:
[1, 1] [2, 2] [3, 3] [4, 4] [5, 5] [6, 6] [7, 8] [8, 9] [9, 10] [10, 12] [11, 15] [12, 16] [13, 18] [14, 20] [15, 24] [16, 25] [17, 27] [18, 30] [19, 32] [20, 36] [1691, 2125764000] [1000000, 519312780448388736089589843750000000000000000000000000000000000000000000000000000000] elapsed: 143.96875 seconds
[edit] Scala
class Hamming extends Iterator[BigInt] {
import scala.collection.mutable.Queue
val qs = Seq.fill(3)(new Queue[BigInt])
def enqueue(n: BigInt) = qs zip Seq(2, 3, 5) foreach { case (q, m) => q enqueue n * m }
def next = {
val n = qs map (_.head) min;
qs foreach { q => if (q.head == n) q.dequeue }
enqueue(n)
n
}
def hasNext = true
qs foreach (_ enqueue 1)
}
However, the usage of closures adds a significant amount of time. The code below, though a bit uglier because of the repetitions, is twice as fast:
class Hamming extends Iterator[BigInt] {
import scala.collection.mutable.Queue
val q2 = new Queue[BigInt]
val q3 = new Queue[BigInt]
val q5 = new Queue[BigInt]
def enqueue(n: BigInt) = {
q2 enqueue n * 2
q3 enqueue n * 3
q5 enqueue n * 5
}
def next = {
val n = q2.head min q3.head min q5.head
if (q2.head == n) q2.dequeue
if (q3.head == n) q3.dequeue
if (q5.head == n) q5.dequeue
enqueue(n)
n
}
def hasNext = true
List(q2, q3, q5) foreach (_ enqueue 1)
}
Usage:
scala> new Hamming take 20 toList res87: List[BigInt] = List(1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36) scala> new Hamming drop 1690 next res88: BigInt = 2125764000 scala> new Hamming drop 999999 next res89: BigInt = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Scheme
(define-syntax lons
(syntax-rules ()
((_ lar ldr) (delay (cons lar (delay ldr))))))
(define (lar lons)
(car (force lons)))
(define (ldr lons)
(force (cdr (force lons))))
(define (lap proc . llists)
(lons (apply proc (map lar llists)) (apply lap proc (map ldr llists))))
(define (take n llist)
(if (zero? n)
(list)
(cons (lar llist) (take (- n 1) (ldr llist)))))
(define (llist-ref n llist)
(if (= n 1)
(lar llist)
(llist-ref (- n 1) (ldr llist))))
(define (merge llist-1 . llists)
(define (merge-2 llist-1 llist-2)
(cond ((null? llist-1) llist-2)
((null? llist-2) llist-1)
((< (lar llist-1) (lar llist-2))
(lons (lar llist-1) (merge-2 (ldr llist-1) llist-2)))
((> (lar llist-1) (lar llist-2))
(lons (lar llist-2) (merge-2 llist-1 (ldr llist-2))))
(else (lons (lar llist-1) (merge-2 (ldr llist-1) (ldr llist-2))))))
(if (null? llists)
llist-1
(apply merge (cons (merge-2 llist-1 (car llists)) (cdr llists)))))
(define hamming
(lons 1
(merge (lap (lambda (x) (* x 2)) hamming)
(lap (lambda (x) (* x 3)) hamming)
(lap (lambda (x) (* x 5)) hamming))))
(display (take 20 hamming))
(newline)
(display (llist-ref 1691 hamming))
(newline)
(display (llist-ref 1000000 hamming))
(newline)
Output:
(1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36)
2125764000
out of memory
[edit] Smalltalk
Works with: GNU Smalltalk
This is a straightforward implementation of the pseudocode snippet found in the Python section. Smalltalk supports arbitrary-precision integers, but the implementation is too slow to try it with 1 million.
Object subclass: Hammer [
Hammer class >> hammingNumbers: howMany [
|h i j k x2 x3 x5|
h := OrderedCollection new.
i := 0. j := 0. k := 0.
h add: 1.
x2 := 2. x3 := 2. x5 := 5.
[ ( h size) < howMany ] whileTrue: [
|m|
m := { x2. x3. x5 } sort first.
(( h indexOf: m ) = 0) ifTrue: [ h add: m ].
( x2 = (h last) ) ifTrue: [ i := i + 1. x2 := 2 * (h at: i) ].
( x3 = (h last) ) ifTrue: [ j := j + 1. x3 := 3 * (h at: j) ].
( x5 = (h last) ) ifTrue: [ k := k + 1. x5 := 5 * (h at: k) ].
].
^ h sort
]
].
(Hammer hammingNumbers: 20) displayNl.
(Hammer hammingNumbers: 1690) last displayNl.
[edit] Tcl
This uses coroutines to simplify the description of what's going on.
Works with: Tcl version 8.6
package require Tcl 8.6
# Simple helper: Tcl-style list "map"
proc map {varName list script} {
set l {}
upvar 1 $varName v
foreach v $list {lappend l [uplevel 1 $script]}
return $l
}
# The core of a coroutine to compute the product of a hamming sequence.
#
# Tricky bit: we don't automatically advance to the next value, and instead
# wait to be told that the value has been consumed (i.e., is the result of
# the [yield] operation).
proc ham {key multiplier} {
global hammingCache
set i 0
yield [info coroutine]
# Cannot use [foreach]; that would take a snapshot of the list in
# the hammingCache variable, so missing updates.
while 1 {
set n [expr {[lindex $hammingCache($key) $i] * $multiplier}]
# If the number selected was ours, we advance to compute the next
if {[yield $n] == $n} {
incr i
}
}
}
# This coroutine computes the hamming sequence given a list of multipliers.
# It uses the [ham] helper from above to generate indivdual multiplied
# sequences. The key into the cache is the list of multipliers.
#
# Note that it is advisable for the values to be all co-prime wrt each other.
proc hammingCore args {
global hammingCache
set hammingCache($args) 1
set hammers [map x $args {coroutine ham$x,$args ham $args $x}]
yield
while 1 {
set n [lindex $hammingCache($args) [incr i]-1]
lappend hammingCache($args) \
[tcl::mathfunc::min {*}[map h $hammers {$h $n}]]
yield $n
}
}
# Assemble the pieces so as to compute the classic hamming sequence.
coroutine hamming hammingCore 2 3 5
# Print the first 20 values of the sequence
for {set i 1} {$i <= 20} {incr i} {
puts [format "hamming\[%d\] = %d" $i [hamming]]
}
for {} {$i <= 1690} {incr i} {set h [hamming]}
puts "hamming{1690} = $h"
for {} {$i <= 1000000} {incr i} {set h [hamming]}
puts "hamming{1000000} = $h"
Produces this output:
hamming{1} = 1
hamming{2} = 2
hamming{3} = 3
hamming{4} = 4
hamming{5} = 5
hamming{6} = 6
hamming{7} = 8
hamming{8} = 9
hamming{9} = 10
hamming{10} = 12
hamming{11} = 15
hamming{12} = 16
hamming{13} = 18
hamming{14} = 20
hamming{15} = 24
hamming{16} = 25
hamming{17} = 27
hamming{18} = 30
hamming{19} = 32
hamming{20} = 36
hamming{1690} = 2123366400
hamming{1000000} = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Ursala
Smooth is defined as a second order function taking a list of primes and returning a function that takes a natural number n to the n-th smooth number with respect to them. An elegant but inefficient formulation based on the J solution is the following.
#import std
#import nat
smooth"p" "n" = ~&z take/"n" nleq-< (rep(length "n") ^Ts/~& product*K0/"p") <1>
This test program
main = smooth<2,3,5>* nrange(1,20)
yields this list of the first 20 Hamming numbers.
<1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36>
Although all calculations are performed using unlimited precision, the version above is impractical for large numbers. A more hardcore approach is the following.
#import std
#import nat
smooth"p" "n" =
~&H\"p" *-<1>; @NiXS ~&/(1,1); ~&ll~="n"->lr -+
^\~&rlPrrn2rrm2Zlrrmz3EZYrrm2lNCTrrm2QAX*rhlPNhrnmtPA2XtCD ~&lrPrhl2E?/~&l ^|/successor@l ~&hl,
^|/~& nleq-<&l+ * ^\~&r ~&l|| product@rnmhPX+-
#cast %nL
main = smooth<2,3,5>* nrange(1,20)--<1691,1000000>
The program produces this output. The great majority of time is spent calculating the millionth Hamming number.
< 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 2125764000, 519312780448388736089589843750000000000000000000000000000000000000000000000000000000>
