Hamming Numbers Variant

lua-users home
wiki

Inspired by HammingNumbers, this is another program that computes the sequence of Hamming Numbers. You may want to read that page for more explanation.

The Hamming numbers are those numbers that are not divisible by any other prime than 2, 3, and 5. This sequence is A051037 (http://www.research.att.com/~njas/sequences/A051037) in Sloane. The first few Hamming numbers are:

1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100,

Our task is to generate these numbers in increasing order. This task is a well-known example on why lazy lists are useful. In fact, both this code an the one on HammingNumbers uses lazy lists as a solution. A third implementation can be found in Higher Order Perl by Mark Jason Dominus (http://hop.perl.plover.com/). I'd like to thank the authors of both implementations for the ideas given there.

So, let's see how the lazy list solution goes.

We define a lazy list as a promise to a pair of the first element of the list and the rest of the list. (Most other people define lazy lists a bit differently: as a pair of the first element and a promise to the rest of the list. I prefer to use the first approach, but there isn't too much difference here.)

Now, classically we can generate the Hamming sequence like this:

hamming := cons(1, merge(map(x -> 2*x, Hamming), merge(map(x -> 3*x, Hamming), map(x -> 5*x, Hamming))));

Here, cons(a, d) makes a list from the first element a and the rest of the list d; map(f, x) returns a list each of whose element is the result of applying the function f to the respective elements of x; merge(a, b) merges two sorted lists to another sorted list; x -> n*x is a unary function that multiples its input by n; and the operations are sufficently lazy.

This solution uses the fact that a number is Hamming number iff it's equal to 1 or it's 2 times a Hamming number or it's 3 times a Hamming number or it's 5 times a Hamming number.

The problem with this approach is this. Some numbers are generated by multiple times. For example, 6 is both 2 times a Hamming number and 3 times a Hamming number. Thus, to make this approach work, you have to define merge in such a way that it puts common elements only once to the resulting list.

However, there's another way to generate the Hamming numbers. For this, let's call a number Very Hamming iff its only prime divisors are 2 and 3. Then, we generate the powers of twos, then then Very Hamming numbers from these, and then the Hamming numbers from these.

Powers of two are easy: a number is a power of two iff it's 1 or it's 2 times a power of two. Now notice that a number is Very Hamming iff it's either a power of two or 3 times a Very Hamming number. Similarly, a number is Hamming iff it's either Very Hamming or it's 5 times a Hamming number. Notice how this generates every number in each sequence exactly once.

The corresponding formulae are these.

two_power := cons(1, map(x -> 2*x, two_power));

very_hamming := merge(two_power, map(x -> 3*x, very_hamming));

hamming := merge(very_hamming, map(x -> 5*x, hamming));

This, however, doesn't work as is. Here's why. The equations for the sequences are true, but they are not enough to generate the sequences because when you try to calculate them, you get an infinite loop. Why is that?

The first sequence, two_power is fine. But suppose you want to try to get the first element of the very_hamming list. For this, you have to calculate the first element of two_power and the first element of map(x -> 3*x, very_hamming), and take the lesser of the two. But, to calculate the first element of map(x -> 3*x, very_hamming), you have to calculate the first element of very_hamming. That's an infinite recursion.

We know that the first element of very_hamming is the first element of two_power, 1. However, that isn't clear from the definition, because, for that, you'd have to know that the first element of map(x -> 3*x, very_hamming) is greater than 1. This, however, you can't know before you know the first element of very_hamming.

To solve this, you have to explicitly tell the program the first few elements of very_hamming and also those of hamming. After those changes, the algorithm indeed works.

Now let's see the program.

The first part is the bignum library I've taken verbatim from HammingNumbers. The second defines promise operations, the third defines pairs. These two are each treated as abstract data types. The fourth part defines some functions for lazy lists, including map and merge. The last section defines the sequence of Hamming numbers and prints some of them.

-- hamming.lua hamming numbers example
-- see http://lua-users.org/wiki/HammingNumbers
-- this code is unoptimized

-- bignums -----------------------------------------------------------
-- bignum library

do
  -- very limited bignum stuff; just enough for the examples here.
  -- Please feel free to improve it.
  local base = 1e15
  local fmt = "%015.0f"
  local meta = {}
  function meta:__lt(other)
    if self.n ~= other.n then return self.n < other.n end
    for i = 1, self.n do
      if self[i] ~= other[i] then return self[i] < other[i] end
    end
  end
  function meta:__eq(other)
    if self.n == other.n then
      for i = 1, self.n do
        if self[i] ~= other[i] then return false end
      end
      return true
    end
  end
  function meta:__mul(k)
    -- If the base where a multiple of all possible multipliers, then
    -- we could figure out the length of the result directly from the
    -- first "digit". On the other hand, printing the numbers would be
    -- difficult. So we accept the occasional overflow.
    local offset = 0
    if self[1] * k >= base then offset = 1 end
    local carry = 0
    local result = {}
    local n = self.n
    for i = n, 1, -1 do
      local tmp = self[i] * k + carry
      local digit = tmp % base
      carry = (tmp - digit) / base
      result[offset + i] = digit
    end
    if carry ~= 0 then
      n = n + 1
      if offset == 0 then
        for i = n, 2, -1 do
          result[i] = result[i - 1]
        end
      end
      result[1] = carry
    end
    result.n = n
    return setmetatable(result, meta)
  end
  -- Added so that Fib would work; other must be a Bignum
  function meta:__add(other)
    local result = {}
    if self.n < other.n then self, other = other, self end
    local n, m = self.n, other.n
    local diff = n - m
    local carry = 0
    local result = {}
    for i = m, 1, -1 do
      local tmp = self[i + diff] + other[i] + carry
      if tmp < base then
        carry = 0
      else
        tmp = tmp - base
        carry = 1
      end
      result[i + diff] = tmp
    end
    for i = diff, 1, -1 do
      local tmp = self[i] + carry
      if tmp < base then
        carry = 0
      else
        tmp = tmp - base
        carry = 1
      end
      result[i] = tmp
    end
    if carry > 0 then
      n = n + 1
      for i = n, 2, -1 do
        result[i] = result[i - 1]
      end
      result[1] = carry
    end
    result.n = n
    return setmetatable(result, meta)
  end

  function meta:__tostring()
    local tmp = {}
    tmp[1] = ("%.0f"):format(self[1])
    for i = 2, self.n do
      tmp[i] = fmt:format(self[i])
    end
    return table.concat(tmp)
  end
  function Bignum(k)
    return setmetatable({k, n = 1}, meta)
  end
end


-- promises ----------------------------------------------------------
function delay(f)
	return {promise_tag = true; promise_func = f};
end;
function force(p)
	if not p.promise_tag then error("forcing a non-promise object of type "..type(p)) end;
	local f = p.promise_func;
	if f then
		local x = f();
		p.promise_val, p.promise_func = x, nil;
		return x;
	else
		return p.promise_val;
	end;
end;


-- pairs -------------------------------------------------------------
-- we need only infinite lists so we won't mess with nulls here
function cons(a, d)
	return {pair_car = a, pair_cdr = d};
end;
function car(c)
	return c.pair_car;
end;
function cdr(c)
	return c.pair_cdr;
end;


-- lazy lists --------------------------------------------------------
-- a lazy list is a promise to a pair whose cdr is a lazy list

-- access fields (thus forcing the list)
function lcar(l)
	return car(force(l));
end;
function lcdr(l)
	return cdr(force(l));
end;

-- map a lazy list
function lmap(f, l)
	return delay(function()
		return cons(f(lcar(l)), lmap(f, lcdr(l)));
	end);
end;

-- merge two nondecreasing lazy lists to a lazy list
function lmerge(a, b)
	return delay(function()
		local x, y = lcar(a), lcar(b);
		if x <= y then
			return cons(x, lmerge(lcdr(a), b));
		else
			return cons(y, lmerge(a, lcdr(b)));
		end;
	end);
end;


-- iterate a lazy list
function lforeach(f, l)
	f(lcar(l));
	return lforeach(f, lcdr(l));
end;

function lany(f, l)
	x = f(lcar(l));
	if x then
		return x;
	else
		return lany(f, lcdr(l));
	end;
end;


-- main ------------------------------------------------------------
--dofile "dump.lua";

-- sequence of natural numbers
local N; 
N = delay(function() 
	return cons(1, lmap(function(x) return x + 1; end, N)) 
end);

-- sequence of Hamming numbers
local timeser = function(x) return function(y) return y * x; end; end;
local H2; H2 = delay(function()
	return cons(Bignum(1), lmap(timeser(2), H2));
end);
local H3; H3 = delay(function()
	return cons(Bignum(1), 
		delay(function() 
			return cons(Bignum(2), 
				lcdr(lcdr(lmerge(lmap(timeser(3), H3), H2)))
			);
		end)
	);	
end);
local H5; H5 = delay(function()
	return cons(Bignum(1), 
		delay(function() 
			return cons(Bignum(2), 
				lcdr(lcdr(lmerge(lmap(timeser(5), H5), H3)))
			);
		end)
	);	
end);

local n, m = 1, 500005;
lany(function(a)
	if 0 == n % 50000 or n <= 200 then
		print(n, a);
	end;
	n = n + 1;
	return m <= n;
end, H5);

-- END

RecentChanges · preferences
edit · history
Last edited March 23, 2009 7:28 pm GMT (diff)