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.

About these ads

Pages: 1 2

31 Responses to “Hamming Numbers”

  1. Graham said

    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 hammingNumbers
    

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

  2. Graham said

    Looks like I had some redundant parentheses in my definition of merge, according to hlint. Apologies.

  3. Jussi Piitulainen said

    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:

    > (aq 10)
    #(1 2 3 4 5 6 8 9 10 12)
    

    (Cut-and-paste from my editor window to be safer against silly typos. I called my translation of Dijkstra’s code aq.)

  4. programmingpraxis said

    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.

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

  6. 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)))
    
  7. Philipp Siegmantel said


    hemming = [1] ++ concatMap (\ x -> [x*2, x*3, x*5]) hemming
    first1k = sort . take 1000 $ hemming

    My Haskell variant

  8. Jussi Piitulainen said

    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.

  9. Graham said

    @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?

  10. programmingpraxis said

    nub

  11. Lautaro Pecile said
    
    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**c
    
    
  12. Continuing 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
    
  13. Adolfo said
    #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))))
    
  14. Adolfo said

    In case you’re curious, my solution is in Lazy Racket :-)

  15. 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 hammingNumbers
    
  16. fa said

    HammingNumbers.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 
    25 
    

    hamm.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 
    
  17. Graham said

    @programmingpraxis: I’ve tried throwing nub into the solution, but it either (1) comes after take 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 :-(

  18. Mike said

    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]

  19. Mike said

    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]
    
    
  20. Axio said

    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 hn
    
  21. nbm said

    Inelegant but fast to write in Python:

    >>> from collections import defaultdict
    ... numbers = defaultdict(lambda: False)
    ... numbers[1] = False
    ... while len(numbers) < 1000:
    ...     tbd = [key for key in numbers.keys() if not numbers[key]]
    ...     tbd.sort()
    ...     me = tbd.pop(0)
    ...     numbers[me] = True
    ...     numbers[me * 2] = numbers[me * 2]
    ...     numbers[me * 3] = numbers[me * 3]
    ...     numbers[me * 5] = numbers[me * 5]
    ... hamming = numbers.keys()
    ... hamming.sort()
    
  22. CyberSpace17 said

    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

  23. programmingpraxis said

    The suggested solution makes n trips through a single loop, at each iteration calculating a minimum of three integers, performing three if tests comparing two integers at each iteration, and making a single addition and a single multiplication.

    Your code has 5 for loops, two of them with if tests inside them, and 1 while loop containing an if which has a nested if inside it. And the inner loop of your code performs a modulo operation (a very expensive operation) and a compare for divisorSize × sequenceIndex iterations.

    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

  24. [...] 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 [...]

  25. [...] This question is copyright of programmingpraxis.com. Click here to see the original [...]

  26. CyberSpace17 said

    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?”.

  27. cosmin said

    Another solution based upon Dijkstra’s paper:

    def hammingSeq(N):
    	h = [1]
    	i2, i3, i5 = 0, 0, 0
    	for i in xrange(1, N):
    		x = min(2*h[i2], 3*h[i3], 5*h[i5])
    		h.append(x)
    		if 2*h[i2] <= x: i2 += 1
    		if 3*h[i3] <= x: i3 += 1
    		if 5*h[i5] <= x: i5 += 1
    	return h
    
    print hammingSeq(1000)
    
  28. Kaos Engr said

    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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 609 other followers

%d bloggers like this: