Hamming Numbers Variant

wiki

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/). 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 -> 3*x, hamming));

This, however, doesn't work as is.

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