Hamming Numbers
August 30, 2011
The sequence of Hamming numbers 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, … (A051037) consists of all numbers of the form 2i·3j·5k where i, j and k are non-negative integers. Edsger Dijkstra introduced this sequence to computer science in his book A Discipline of Programming, and it has been a staple of beginning programming courses ever since. Dijkstra wrote a program based on three axioms:
Axiom 1: The value 1 is in the sequence.
Axiom 2: If x is in the sequence, so are 2 * x, 3 * x, and 5 * x.
Axiom 3: The sequence contains no other values other than those that belong to it on account of Axioms 1 and 2.
.
Your task is to write a program to compute the first thousand terms of the Hamming sequence. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
I wrote a Haskell version of Python’s test of using lazy evaluation to generate them; as it turns out, it’s the same as the SRFI-41 version.
merge :: (Ord a) => [a] -> [a] -> [a] merge [] ys = ys merge xs [] = xs merge (x:xs) (y:ys) | x < y = x : (merge xs (y:ys)) | x > y = y : (merge (x:xs) ys) | x == y = y : merge xs ys hammingNumbers :: [Integer] hammingNumbers = 1 : merge hN2 (merge hN3 hN5) where hN2 = map (*2) hammingNumbers hN3 = map (*3) hammingNumbers hN5 = map (*5) hammingNumbers main :: IO () main = print $ take 1000 hammingNumbersSince I’m new to Haskell, I probably didn’t find the most elegant solution (and suggestions are welcome), but lazy evaluation is kind of amazing.
Looks like I had some redundant parentheses in my definition of merge, according to hlint. Apologies.
Re the imperative version, I think it is nicer to omit the 0 that does not belong to the sequence and simply produce a 0-based sequence of n elements:
(Cut-and-paste from my editor window to be safer against silly typos. I called my translation of Dijkstra’s code aq.)
Jussi: I tend to add a useless element to the beginning of arrays when porting code from languages that base arrays at 1 instead of 0. It’s mostly harmless, and saves a lot of x+1 and x-1 statements mapping between the two formats. I agree it’s ugly, but only a little bit. Where it actually matters, I do it right.
Here’s a shortish but quite efficient lazy Clojure version:
(defn hammings [initial-set]
(let [v (first initial-set)
others (rest initial-set)]
(lazy-seq
(cons
v
(hammings (into (sorted-set (* v 2) (* v 3) (* v 5)) others))))))
(take 1000 (hammings #{1}))
Runs in about 0.07 ms for the first 1000 hamming numbers on my machine.
Here’s my clojure solution. More details at my blog.
;;Hamming Numbers (defn hamming-generator ([] (hamming-generator (sorted-set 1))) ([some-set] (let [n (first some-set) s (disj some-set n)] (cons n (lazy-seq (hamming-generator (conj s (* n 2) (* n 3) (* n 5)))))))) (def hamming-numbers (take 1000 (hamming-numbers)))hemming = [1] ++ concatMap (\ x -> [x*2, x*3, x*5]) hemming
first1k = sort . take 1000 $ hemming
My Haskell variant
Dijkstra advocated 0-based indexing, so I think his algorithm in the chapter may actually be intended to be 0-based. It is hard to tell for sure. However, his initialization of aq to (1 1) is funny either way when only one 1 is meant to remain in the result, and the assignments to ik, xk seem to need to be either parallel or swapped when I implement it. Below is a 0-based version with parallel assignments. For example, (set*! (x y) y x) swaps x and y. (Originally I had tail-recursion instead of assignments, and I did not even notice that the second assignments in Dijkstra depend on the first ones.)
(define-syntax set*! (syntax-rules () ((set*! (var ...) exp ...) (set*! "gen" (var ...) () (var ...) exp ...)) ((set*! "gen" (x y ...) (t ...) (var ...) exp ...) (set*! "gen" (y ...) (u t ...) (var ...) exp ...)) ((set*! "gen" () (tmp ...) (var ...) exp ...) (set*! "map" (tmp ...) (var ...) (exp ...) ())) ((set*! "map" (x y ...) (t u ...) (d e ...) (w ...)) (set*! "map" (y ...) (u ...) (e ...) ((x t d) w ...))) ((set*! "map" () () () ((tmp var exp) ...)) (let ((tmp exp) ...) (set! var tmp) ...)))) (define (hamming2 n) (let ((aq (make-vector n 1)) (i2 1) (i3 1) (i5 1) (x2 2) (x3 3) (x5 5)) (do ((h 1 (+ h 1))) ((= h n) aq) (let ((t (min x2 x3 x5))) (vector-set! aq h t) (if (= x2 t) (set*! (i2 x2) (+ i2 1) (* 2 (vector-ref aq i2)))) (if (= x3 t) (set*! (i3 x3) (+ i3 1) (* 3 (vector-ref aq i3)))) (if (= x5 t) (set*! (i5 x5) (+ i5 1) (* 5 (vector-ref aq i5))))))))Cheers.
@Philipp: I like the brevity of your Haskell answer, but it seems to include a lot of repeats; as a Haskell newbie, I couldn’t figure out a way to remove them easily. Any ideas?
nub
from itertools import product from bisect import insort def hamming_numbers(): result = [] for n in hamming_generator(10): insort(result, n) return result def hamming_generator(n=0): for a, b, c in product(xrange(n), repeat=3): yield 2**a * 3**b * 5**cContinuing to attempt to become more familiar with Haskell. Here is what I came up with.
import qualified Data.Set as Set hammingNumbers = generate $ Set.fromList [1] where generate curr = case Set.minView curr of Just (min, rest) -> min : generate (foldl (\s x -> Set.insert x s) rest (map (*min) [2,3,5])) Nothing -> error "should not get here" answer = take 1000 hammingNumbers#lang lazy (define naturals (cons 1 (map add1 naturals))) (define (merge . lsts) (let ((ordered (sort lsts (λ (x y) (< (car x) (car y)))))) (cons (caar ordered) (apply merge (cons (cdar ordered) (map (λ (x) (remove (caar ordered) x)) (cdr ordered))))))) (define hamming-numbers (cons 1 (merge (map (λ (x) (* x 5)) hamming-numbers) (map (λ (x) (* x 3)) hamming-numbers) (map (λ (x) (* x 2)) hamming-numbers))))In case you’re curious, my solution is in Lazy Racket :-)
Ended up playing with my answer some more while waiting for some other code to compile. I like how this one better.
import qualified Data.Set as Set hammingNumbers = generate $ Set.fromList [1] where generate curr = min : generate (foldl (\s x -> Set.insert (min * x) s) rest [2,3,5]) where (min, rest) = Set.deleteFindMin curr hammingAnswer = take 1000 hammingNumbersHammingNumbers.java
1 public class HammingNumbers 2 { 3 public static int[] hammSeq(int aLen) { 4 int seq[] = new int[aLen]; 5 int i, i2, i3, i5, x, x2, x3, x5; 6 7 seq[0] = 1; i = 0; x = 1; 8 i2 = i3 = i5 = -1; x2 = 2; x3 = 3; x5 = 5; 9 while (++i < aLen) { 10 seq[i] = x = (x2 <= x3 && x2 <= x5) ? x2 : (x3 <= x5) ? x3 : x5; 11 while (x2 <= x) x2 = 2 * seq[++i2]; 12 while (x3 <= x) x3 = 3 * seq[++i3]; 13 while (x5 <= x) x5 = 5 * seq[++i5]; 14 } 15 return seq; 16 } 17 18 public static void main(String args[]) { 19 for (int i : hammSeq(1000)) 20 System.out.print(" " + i); 21 System.out.println(); 22 } 23 } 24 25hamm.go
1 package main 2 import "fmt" 3 4 func hamm(slen int, sfunc []func(int)int) (seq []int) { 5 seq = make([]int, slen); seq[0] = 1 // the calculated values 6 isf := make([]int, len(sfunc)); for i := 0; i < len(sfunc); i++ { isf[i] = -1 } // the last index for each function 7 xsf := make([]int, len(sfunc)); for i := 0; i < len(sfunc); i++ { xsf[i] = sfunc[i](seq[0]) } // the next value of each function 8 xmin := func()int { x := xsf[0]; for i := 1; i < len(sfunc); i++ { if xsf[i] < x { x = xsf[i] } }; return x } 9 c := make(chan int) 10 for is := 1; is < slen; is++ { 11 seq[is] = xmin() 12 for i, sf := range sfunc { // find the next value for each function, paralellize 13 go func(sf func(int)int, isf *int, xsf *int, seq []int, is int) { 14 for *xsf <= seq[is] { *isf++; *xsf = sf(seq[*isf]) } 15 c <- 1 // next value found 16 } (sf, &isf[i], &xsf[i], seq, is) // pass the arguments 17 } 18 for i := 0; i < len(sfunc); i++ { <-c } // wait for each calculation to end 19 } 20 return 21 } 22 23 func main() { // print the first 1000 hamming numbers calculated from the given functions 24 sfunc := []func(int)int{ func(x int)int{ return 2*x }, func(x int)int{ return 3*x }, func(x int)int{ return 5*x } } 25 for _, x := range hamm(1000, sfunc) { fmt.Printf("%v ", x) }; fmt.Printf("\n") 26 } 27 28@programmingpraxis: I’ve tried throwing
nubinto the solution, but it either (1) comes aftertake 1000, in which case there are no longer 1000 entries, or (2) comes before it, and the program just hangs. I am not very adept in Haskell :-(Hamming number generator in Python 3.
Follows directly from the definition. Yield 1, then yield 2*, 3* and 5* a number in the list. It’s lazy too.
from heapq import merge
from itertools import tee
def hamming_numbers():
last = 1
yield last
a,b,c = tee(hamming_numbers(), 3)
for n in merge((2*i for i in a), (3*i for i in b), (5*i for i in c)):
if n != last:
yield n
last = n
x = list(islice(hamming_numbers(),1000))
print(x[:10], x[-5:])
# output -> [1, 2, 3, 4, 5, 6, 8, 9, 10, 12] [50331648, 50388480, 50625000, 51018336, 51200000]
This time with proper formating:
from heapq import merge from itertools import tee def hamming_numbers(): last = 1 yield last a,b,c = tee(hamming_numbers(), 3) for n in merge((2*i for i in a), (3*i for i in b), (5*i for i in c)): if n != last: yield n last = n #test x = list(islice(hamming_numbers(),1000)) print(x[:10],x[-5:]) #output -> [1, 2, 3, 4, 5, 6, 8, 9, 10, 12] [50331648, 50388480, 50625000, 51018336, 51200000]Nothing magical here…
merge (x:xs) (y:ys) = case compare x y of LT -> x:(merge xs (y:ys)) EQ -> x:(merge xs ys) GT -> y:(merge (x:xs) ys) hn = 1: merge (map (* 2) hn) (merge (map (* 3) hn) (map (* 5) hn)) run = take 1000 hnInelegant but fast to write in Python:
my C++ solution. I’m disappointed that it takes so long for my implementation to calculate the sequence. (especially since someone got it in less thana tenth of a second.)Can someone tell me where I can improve the code?
http://codepad.org/WhHDxJAZ
The suggested solution makes n trips through a single loop, at each iteration calculating a minimum of three integers, performing three
iftests comparing two integers at each iteration, and making a single addition and a single multiplication.Your code has 5
forloops, two of them withiftests inside them, and 1whileloop containing anifwhich has a nestedifinside it. And the inner loop of your code performs amodulooperation (a very expensive operation) and a compare fordivisorSize×sequenceIndexiterations.And you forgot the 1 that starts the sequence.
By the way, the suggested solution runs in a thousandth of a second, not a tenth of a second:
> (define x #f)> (time (set! x (hamming 1000)))
(time (set! x ...))
no collections
1 ms elapsed cpu time
1 ms elapsed real time
12048 bytes allocated
> (vector-ref x 1000)
51200000
To improve your code, you should follow the suggested solution.
I’m sorry if that’s a little bit harsh. I’ll be happy to look at your improved solution.
Phil
[…] compute the first thousand terms of the Hamming sequence. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments […]
[…] This question is copyright of programmingpraxis.com. Click here to see the original […]
Here is my C implementation!
http://codepad.org/ilubfEfd
Thanks Phil for a review of my code. I now see how inefficient it is and was working on a revision but I don’t think I can make my algorithm as efficient as Supriya Koduri’s, to whom I want to ask “what train of thought allowed him to create such an elegant solution?”.
ruby solution : http://codepad.org/b87tK8qQ
Another solution based upon Dijkstra’s paper:
cosmin’s solution in Go: http://play.golang.org/p/qD1Tqntcy5
The results of the cosmin’s solution in Go posted by glnbrwn 12Dec2012 in Go do not match.
Python version posted October 30, 2012 by cosmin results are:
[1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, …
Go when executed at the link given above results in:
[1 2 3 4 6 8 9 12 16 18 24 27 32 36 48 54 64
Go is missing 10, 15, 20, 25, 30 when executing it. Multiples of 5 are not getting added to the result.