Hamming Numbers Variant
 |
||
|
Hamming Numbers Variant |
|
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:
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.
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
