Hamming numbers
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
Hamming numbers are numbers of the form
-
.
Hamming numbers are also known as ugly numbers and also 5-smooth numbers (numbers whose prime divisors are less or equal to 5).
Generate the sequence of Hamming numbers, in increasing order. In particular:
- Show the first twenty Hamming numbers.
- Show the 1691st Hamming number (the last one below 231).
- Show the one millionth Hamming number (if the language – or a convenient library – supports arbitrary-precision integers).
References
- wp:Hamming_numbers
- wp:Smooth_number
- Hamming problem from Dr. Dobb's CodeTalk (dead link as of Sep 2011; parts of the thread here and here).
[edit] Ada
GNAT provides the datatypes Integer, Long_Integer and Long_Long_Integer.
Values for GNAT Pro 6.3.1, 64 bit Linux version:
- Integer covers the range -2**31 .. 2**31-1 (-2147483648 .. 2147483647).
- Long_Integer and Long_Long_Integer cover the range -2**63 .. 2**63-1 (-9223372036854775808 .. 9223372036854775807).
Using your own modular integer type (for example type My_Unsigned_Integer is mod 2**64;), you can expand the range to 0 .. 18446744073709551615, but this still is not enough for the millionth Hamming number.
For bigger numbers, you have to use an external library, for example Big_Number.
The code for calculating the Hamming numbers is kept generic, to easily expand the range by changing the concrete type.
with Ada.Text_IO;
procedure Hamming is
generic
type Int_Type is private;
Zero : Int_Type;
One : Int_Type;
Two : Int_Type;
Three : Int_Type;
Five : Int_Type;
with function "mod" (Left, Right : Int_Type) return Int_Type is <>;
with function "/" (Left, Right : Int_Type) return Int_Type is <>;
with function "+" (Left, Right : Int_Type) return Int_Type is <>;
function Get_Hamming (Position : Positive) return Int_Type;
function Get_Hamming (Position : Positive) return Int_Type is
function Is_Hamming (Number : Int_Type) return Boolean is
Temporary : Int_Type := Number;
begin
while Temporary mod Two = Zero loop
Temporary := Temporary / Two;
end loop;
while Temporary mod Three = Zero loop
Temporary := Temporary / Three;
end loop;
while Temporary mod Five = Zero loop
Temporary := Temporary / Five;
end loop;
return Temporary = One;
end Is_Hamming;
Result : Int_Type := One;
Previous : Positive := 1;
begin
while Previous /= Position loop
Result := Result + One;
if Is_Hamming (Result) then
Previous := Previous + 1;
end if;
end loop;
return Result;
end Get_Hamming;
-- up to 2**32 - 1
function Integer_Get_Hamming is new Get_Hamming
(Int_Type => Integer,
Zero => 0,
One => 1,
Two => 2,
Three => 3,
Five => 5);
-- up to 2**64 - 1
function Long_Long_Integer_Get_Hamming is new Get_Hamming
(Int_Type => Long_Long_Integer,
Zero => 0,
One => 1,
Two => 2,
Three => 3,
Five => 5);
begin
Ada.Text_IO.Put ("1) First 20 Hamming numbers: ");
for I in 1 .. 20 loop
Ada.Text_IO.Put (Integer'Image (Integer_Get_Hamming (I)));
end loop;
Ada.Text_IO.New_Line;
Ada.Text_IO.Put_Line ("2) 1_691st Hamming number: " &
Integer'Image (Integer_Get_Hamming (1_691)));
-- even Long_Long_Integer overflows here
Ada.Text_IO.Put_Line ("3) 1_000_000st Hamming number: " &
Long_Long_Integer'Image (Long_Long_Integer_Get_Hamming (1_000_000)));
end Hamming;
- Output:
1) First 20 Hamming numbers: 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2) 1_691 st Hamming number: 2125764000 Execution terminated by unhandled exception Exception name: CONSTRAINT_ERROR Message: hamming.adb:34 overflow check failed Call stack traceback locations: 0x403212 0x402fd7 0x402a87 0x7f8b99517584 0x4026d7
For using Big_Number, you just have to add this to the code (additional to with Big_Number; and with Ada.Strings.Unbounded; in context clause):
type My_Index is mod 2**8;
package My_Big_Numbers is new Big_Number (Index_type => My_Index, Nb_Item => 64);
function Int2Big is new My_Big_Numbers.Generic_Conversion.Int_Number2Big_Unsigned (Integer);
function Big_Get_Hamming is new Get_Hamming
(Int_Type => My_Big_Numbers.Big_Unsigned,
Zero => My_Big_Numbers.Big_Unsigned_Zero,
One => My_Big_Numbers.Big_Unsigned_One,
Two => My_Big_Numbers.Big_Unsigned_Two,
Three => Int2Big(3),
Five => Int2Big(5),
"mod" => My_Big_Numbers.Unsigned_Number."mod",
"+" => My_Big_Numbers.Unsigned_Number."+",
"/" => My_Big_Numbers.Unsigned_Number."/");
and then use it like this:
Ada.Text_IO.Put_Line ("3) 1_000_000st Hamming number: " &
Ada.Strings.Unbounded.To_String (My_Big_Numbers.String_Conversion.Big_Unsigned2UString (Big_Get_Hamming (1_000_000))));
[edit] AutoHotkey
SetBatchLines, -1
Msgbox % hamming(1,20)
Msgbox % hamming(1690)
return
hamming(first,last=0)
{
if (first < 1)
ans=ERROR
if (last = 0)
last := first
i:=0, j:=0, k:=0
num1 := ceil((last * 20)**(1/3))
num2 := ceil(num1 * ln(2)/ln(3))
num3 := ceil(num1 * ln(2)/ln(5))
loop
{
H := (2**i) * (3**j) * (5**k)
if (H > 0)
ans = %H%`n%ans%
i++
if (i > num1)
{
i=0
j++
if (j > num2)
{
j=0
k++
}
}
if (k > num3)
break
}
Sort ans, N
Loop, parse, ans, `n, `r
{
if (A_index > last)
break
if (A_index < first)
continue
Output = %Output%`n%A_LoopField%
}
return Output
}
[edit] ALGOL 68
Hamming numbers are generated in a trivial iterative way as in the Python version below. This program keeps the series needed to generate the numbers as short as possible using flexible rows; on the downside, it spends considerable time on garbage collection.
PR precision=100 PR
MODE SERIES = FLEX [1 : 0] UNT, # Initially, no elements #
UNT = LONG LONG INT; # A 100-digit unsigned integer #
PROC hamming number = (INT n) UNT: # The n-th Hamming number #
CASE n
IN 1, 2, 3, 4, 5, 6, 8, 9, 10, 12 # First 10 in a table #
OUT # Additional operators #
OP MIN = (INT i, j) INT: (i < j | i | j), MIN = (UNT i, j) UNT: (i < j | i | j);
PRIO MIN = 9;
OP LAST = (SERIES h) UNT: h[UPB h]; # Last element of a series #
OP +:= = (REF SERIES s, UNT elem) VOID:
# Extend a series by one element, only keep the elements you need #
(INT lwb = (i MIN j) MIN k, upb = UPB s;
REF SERIES new s = HEAP FLEX [lwb : upb + 1] UNT;
(new s[lwb : upb] := s[lwb : upb], new s[upb + 1] := elem);
s := new s
);
# Determine the n-th hamming number iteratively #
SERIES h := 1, # Series, initially one element #
UNT m2 := 2, m3 := 3, m5 := 5, # Multipliers #
INT i := 1, j := 1, k := 1; # Counters #
TO n - 1
DO h +:= (m2 MIN m3) MIN m5;
(LAST h = m2 | m2 := 2 * h[i +:= 1]);
(LAST h = m3 | m3 := 3 * h[j +:= 1]);
(LAST h = m5 | m5 := 5 * h[k +:= 1])
OD;
LAST h
ESAC;
FOR k TO 20
DO print ((whole (hamming number (k), 0), blank))
OD;
print ((newline, whole (hamming number (1 691), 0)));
print ((newline, whole (hamming number (1 000 000), 0)))
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] C
Using a min-heap to keep track of numbers. Does not handle big integers.
#include <stdio.h>
#include <stdlib.h>
typedef unsigned long long ham;
size_t alloc = 0, n = 1;
ham *q = 0;
void qpush(ham h)
{
int i, j;
if (alloc <= n) {
alloc = alloc ? alloc * 2 : 16;
q = realloc(q, sizeof(ham) * alloc);
}
for (i = n++; (j = i/2) && q[j] > h; q[i] = q[j], i = j);
q[i] = h;
}
ham qpop()
{
int i, j;
ham r, t;
/* outer loop for skipping duplicates */
for (r = q[1]; n > 1 && r == q[1]; q[i] = t) {
/* inner loop is the normal down heap routine */
for (i = 1, t = q[--n]; (j = i * 2) < n;) {
if (j + 1 < n && q[j] > q[j+1]) j++;
if (t <= q[j]) break;
q[i] = q[j], i = j;
}
}
return r;
}
int main()
{
int i;
ham h;
for (qpush(i = 1); i <= 1691; i++) {
/* takes smallest value, and queue its multiples */
h = qpop();
qpush(h * 2);
qpush(h * 3);
qpush(h * 5);
if (i <= 20 || i == 1691)
printf("%6d: %llu\n", i, h);
}
/* free(q); */
return 0;
}
[edit] Alternative
Standard algorithm. Numbers are stored as exponents of factors instead of big integers, while GMP is only used for display. It's much more efficient this way.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gmp.h>
/* number of factors. best be mutually prime -- duh. */
#define NK 3
#define MAX_HAM 10000000
int n_hams = 0, idx[NK] = {0}, fac[] = { 2, 3, 5, 7, 11};
/* k-smooth numbers are stored as their exponents of each factor;
v is the log of the number, for convenience. */
typedef struct {
unsigned short e[NK];
double v;
} ham_t, *ham;
ham_t *hams, values[NK] = {{{0}, 0}};
double inc[NK];
/* most of the time v can be just incremented, but eventually
* floating point precision will bite us, so better recalculate */
inline
void _setv(ham x) {
int i;
for (x->v = 0, i = 0; i < NK; i++)
x->v += inc[i] * x->e[i];
}
inline
int _eq(ham x, ham y) {
int i;
for (i = 0; i < NK; i++)
if (x->e[i] != y->e[i])
return 0;
return 1;
}
ham get_ham(int n)
{
int i, ni;
ham h;
n--;
while (n_hams < n) {
for (ni = 0, i = 1; i < NK; i++)
if (values[i].v < values[ni].v)
ni = i;
*(h = hams + ++n_hams) = values[ni];
_setv(h);
for (ni = 0; ni < NK; ni++)
if ( _eq(values + ni, h) ) {
values[ni] = hams[++idx[ni]];
values[ni].e[ni]++;
_setv(values + ni);
}
}
return hams + n;
}
void show_ham(ham h)
{
static mpz_t das_ham, tmp;
int i;
mpz_init_set_ui(das_ham, 1);
mpz_init_set_ui(tmp, 1);
for (i = 0; i < NK; i++) {
mpz_ui_pow_ui(tmp, fac[i], h->e[i]);
mpz_mul(das_ham, das_ham, tmp);
}
gmp_printf("%Zu\n", das_ham);
}
int main()
{
int i;
hams = malloc(sizeof(ham_t) * MAX_HAM);
for (i = 0; i < NK; i++) {
values[i].e[i] = 1;
values[i].v = inc[i] = log(fac[i]);
}
printf(" 1,691: "); show_ham(get_ham(1691));
printf(" 1,000,000: "); show_ham(get_ham(1e6));
printf("10,000,000: "); show_ham(get_ham(1e7));
return 0;
}
- Output:
1,691: 2125764000 1,000,000: 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 10,000,000: 16244105063830431823239 ..<a gadzillion digits>.. 000000000000000000000
[edit] C#
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
public class LazyList : IEnumerable<int>
{
int First;
Func<LazyList> Rest;
public LazyList(int first, Func<LazyList> rest)
{
First = first;
Rest = rest;
}
public LazyList Map(Func<int, int> function)
{
return new LazyList(function(First), () => Rest().Map(function));
}
public LazyList Merge(LazyList list)
{
if (this == null)
{
return list;
}
if (list == null)
{
return this;
}
if (First < list.First)
{
return new LazyList(First, () => Rest().Merge(list));
}
if (First > list.First)
{
return new LazyList(list.First, () => Merge(list.Rest()));
}
return new LazyList(First, () => Rest().Merge(list.Rest()));
}
public IEnumerator<int> GetEnumerator()
{
yield return First;
if (Rest() != null)
{
foreach (var first in Rest())
{
yield return first;
}
}
}
IEnumerator IEnumerable.GetEnumerator()
{
return (IEnumerator)GetEnumerator();
}
}
class Program
{
static void Main()
{
var hamming = default(LazyList);
hamming = new LazyList(1, () => hamming.Map(x => x * 2).Merge(hamming.Map(x => x * 3).Merge(hamming.Map(x => x * 5))));
foreach (var number in hamming.Take(20))
{
Console.WriteLine(number);
}
}
}
[edit] Clojure
This version implements Dijkstra's merge solution, so is closely related to the Haskell version.
(defn smerge [xs ys]
(lazy-seq
(let [x (first xs),
y (first ys),
[z xs* ys*]
(cond
(< x y) [x (rest xs) ys]
(> x y) [y xs (rest ys)]
:else [x (rest xs) (rest ys)])]
(cons z (smerge xs* ys*)))))
(defn smerge3 [xs ys zs]
(smerge xs (smerge ys zs)))
(defn map*n [n ks] (map #(* n %) ks))
(def hamming
(lazy-seq
(cons 1 (smerge3 (map*n 2 hamming) (map*n 3 hamming) (map*n 5 hamming)))))
Note that this version uses a lot of space and time after calculating a few hundred thousand elements of the sequence. This is no doubt due to its "holding on to the head": it maintains the entire generated sequence in memory.
[edit] CoffeeScript
# Generate hamming numbers in order. Hamming numbers have the
# property that they don't evenly divide any prime numbers outside
# a given set, such as [2, 3, 5].
generate_hamming_sequence = (primes, max_n) ->
# We use a lazy algorithm, only ever keeping N candidates
# in play, one for each of our seed primes. Let's say
# primes is [2,3,5]. Our virtual streams are these:
#
# hammings: 1,2,3,4,5,6,8,10,12,15,16,18,20,...
# hammings*2: 2,4,6,9.10,12,16,20,24,30,32,36,40...
# hammings*3: 3,6,9,12,15,18,24,30,36,45,...
# hammings*5: 5,10,15,20,25,30,40,50,...
#
# After encountering 40 for the last time, our candidates
# will be
# 50 = 2 * 25
# 45 = 3 * 15
# 50 = 5 * 10
# Then, after 45
# 50 = 2 * 25
# 48 = 3 * 16 <= new
# 50 = 5 * 10
hamming_numbers = [1]
candidates = ([p, p, 1] for p in primes)
last_number = 1
while hamming_numbers.length < max_n
# Get the next candidate Hamming Number tuple.
i = min_idx(candidates)
candidate = candidates[i]
[n, p, seq_idx] = candidate
# Add to sequence unless it's a duplicate.
if n > last_number
hamming_numbers.push n
last_number = n
# Replace the candidate with its successor (based on
# p = 2, 3, or 5).
#
# This is the heart of the algorithm. Let's say, over the
# primes [2,3,5], we encounter the hamming number 32 based on it being
# 2 * 16, where 16 is the 12th number in the sequence.
# We'll be passed in [32, 2, 12] as candidate, and
# hamming_numbers will be [1,2,3,4,5,6,8,9,10,12,16,18,...]
# by now. The next candidate we need to enqueue is
# [36, 2, 13], where the numbers mean this:
#
# 36 - next multiple of 2 of a Hamming number
# 2 - prime number
# 13 - 1-based index of 18 in the sequence
#
# When we encounter [36, 2, 13], we will then enqueue
# [40, 2, 14], based on 20 being the 14th hamming number.
q = hamming_numbers[seq_idx]
candidates[i] = [p*q, p, seq_idx+1]
hamming_numbers
min_idx = (arr) ->
# Don't waste your time reading this--it just returns
# the index of the smallest tuple in an array, respecting that
# the tuples may contain integers. (CS compiles to JS, which is
# kind of stupid about sorting. There are libraries to work around
# the limitation, but I wanted this code to be standalone.)
less_than = (tup1, tup2) ->
i = 0
while i < tup2.length
return true if tup1[i] <= tup2[i]
return false if tup1[i] > tup2[i]
i += 1
min_i = 0
for i in [1...arr.length]
if less_than arr[i], arr[min_i]
min_i = i
return min_i
primes = [2, 3, 5]
numbers = generate_hamming_sequence(primes, 10000)
console.log numbers[1690]
console.log numbers[9999]
[edit] Common Lisp
Maintaining three queues, popping the smallest value every time.
(defun next-hamm (factors seqs)
(let ((x (apply #'min (map 'list #'first seqs))))
(loop for s in seqs
for f in factors
for i from 0
with add = t do
(if (= x (first s)) (pop s))
;; prevent a value from being added to multiple lists
(when add
(setf (elt seqs i) (nconc s (list (* x f))))
(if (zerop (mod x f)) (setf add nil)))
finally (return x))))
(loop with factors = '(2 3 5)
with seqs = (loop for i in factors collect '(1))
for n from 1 to 1000001 do
(let ((x (next-hamm factors seqs)))
(if (or (< n 21)
(= n 1691)
(= n 1000000)) (format t "~d: ~d~%" n x))))
A much faster method:
(defun hamming (n)
(let ((fac '(2 3 5))
(idx (make-array 3 :initial-element 0))
(h (make-array (1+ n)
:initial-element 1
:element-type 'integer)))
(loop for i from 1 to n
with e with x = '(1 1 1) do
(setf e (setf (aref h i) (apply #'min x))
x (loop for y in x
for f in fac
for j from 0
collect (if (= e y) (* f (aref h (incf (aref idx j)))) y))))
(aref h n)))
(loop for i from 1 to 20 do
(format t "~2d: ~d~%" i (hamming i)))
(loop for i in '(1691 1000000) do
(format t "~d: ~d~%" i (hamming i)))
- Output:
1: 1 2: 2 3: 3 4: 4 5: 5 6: 6 7: 8 8: 9 9: 10 10: 12 11: 15 12: 16 13: 18 14: 20 15: 24 16: 25 17: 27 18: 30 19: 32 20: 36 1691: 2125764000 1000000: 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] D
This version keeps all numbers in memory and runs in a couple of seconds, computing all the Hamming numbers up to the needed one. Performs constant number of operations per Hamming number produced.
import std.stdio, std.bigint, std.algorithm, std.range;
auto hamming(int n) {
BigInt two = 2, three = 3, five = 5;
auto h = new BigInt[n];
h[0] = 1;
BigInt x2 = 2, x3 = 3, x5 = 5;
int i, j, k;
foreach (ref el; h[1 .. $]) {
el = min(x2, x3, x5);
if (x2 == el) x2 = two * h[++i];
if (x3 == el) x3 = three * h[++j];
if (x5 == el) x5 = five * h[++k];
}
return h[$ - 1];
}
void main() {
writeln(map!hamming(iota(1, 21)));
writeln(hamming(1691));
writeln(hamming(1_000_000));
}
- Output:
[1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Alternative version 1
This keeps O(n2 / 3) numbers in memory (5.4 seconds run time) but overcomputes a sequence by a factor of about Θ(n2 / 3), calculating extra multiples past that as well. Incurs an extra O(logn) factor of operations per each number produced (reinserting its multiples into a tree). Doesn't stop when the target number is reached, instead continuing until it is no longer needed:
import std.stdio,std.bigint,std.container,std.algorithm,std.range;
BigInt hamming(int n)
in {
assert(n > 0);
} body {
auto frontier = redBlackTree(BigInt(2), BigInt(3), BigInt(5));
auto lowest = BigInt(1);
foreach (_; 1 .. n) {
lowest = frontier.front();
frontier.removeFront();
frontier.insert(lowest * 2);
frontier.insert(lowest * 3);
frontier.insert(lowest * 5);
}
return lowest;
}
void main() {
writeln("First 20 Hamming numbers: ", map!hamming(iota(1, 21)));
writeln("hamming(1691) = ", hamming(1691));
writeln("hamming(1_000_000) = ", hamming(1_000_000));
}
- Output:
First 20 Hamming numbers: [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] hamming(1691) = 2125764000 hamming(1_000_000) = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Alternative version 2
Does exactly what the first version does, creating an array and filling it with Hamming numbers, keeping the three back pointers into the sequence for next multiples calculations, except that it represents the numbers as their coefficients triples and their logarithm values (for comparisons), thus saving on BigInt calculations.
import std.stdio: writefln;
import std.bigint: BigInt, toDecimalString;
import std.numeric: gcd;
import std.algorithm: copy, map;
import std.math; // log, ^^
// number of factors
enum NK = 3;
enum MAX_HAM = 10_000_000;
static assert(gcd(NK, MAX_HAM) == 1);
enum int[NK] fac = [2, 3, 5];
/// k-smooth numbers (stored as their exponents of each factor).
struct Hamming {
double v; // log of the number, for convenience.
ushort[NK] e; // exponents of each factor.
// Compile-time constant, map!log(fac)
// log can't be used in CTFE yet
public static __gshared /*const*/ double[fac.length] inc;
/*nothrow*/ static this() {
// not nothrow because map isn't nothrow
copy(map!log(fac[]), inc[]);
}
bool opEquals(const ref Hamming y) const pure nothrow {
//return this.e == y.e; // too much slow
foreach (size_t i; 0 .. this.e.length)
if (this.e[i] != y.e[i])
return false;
return true;
}
void update() /*pure*/ nothrow {
// can't be pure because inc is not enum/const
//this.v = dotProduct(inc, this.e); // too much slow
this.v = 0.0;
foreach (size_t i; 0 .. this.e.length)
this.v += inc[i] * this.e[i];
}
string toString() const {
BigInt result = 1;
foreach (size_t i, f; fac)
result *= BigInt(f) ^^ this.e[i];
return toDecimalString(result);
}
}
// global variables
__gshared Hamming[] hams;
__gshared Hamming[NK] values;
nothrow static this() {
// Slower than malloc if you don't use all the MAX_HAM items
hams = new Hamming[MAX_HAM];
foreach (i, ref v; values) {
v.e[i] = 1;
v.v = Hamming.inc[i];
}
}
ref Hamming getHam(in size_t n) nothrow
in {
assert(n <= MAX_HAM);
} body {
// most of the time v can be just incremented, but eventually
// floating point precision will bite us, so better recalculate
__gshared static size_t[NK] idx;
__gshared static int n_hams;
for (; n_hams < n; n_hams++) {
{
// find the index of the minimum v
size_t ni = 0;
foreach (size_t i; 1 .. NK)
if (values[i].v < values[ni].v)
ni = i;
hams[n_hams] = values[ni];
hams[n_hams].update();
}
foreach (size_t i; 0 .. NK)
if (values[i] == hams[n_hams]) {
values[i] = hams[idx[i]];
idx[i]++;
values[i].e[i]++;
values[i].update();
}
}
return hams[n - 2];
}
void main() {
foreach (n; [1691, 10 ^^ 6, MAX_HAM])
writefln("%8d: %s", n, getHam(n));
}
The output is similar to the second C version. Runtime is about 0.2 seconds if MAX_HAM = 1_000_000 (as the task requires), and 1.7 seconds if MAX_HAM = 10_000_000.
[edit] Factor
USING: accessors deques dlists fry kernel make math math.order
;
IN: rosetta.hamming
TUPLE: hamming-iterator 2s 3s 5s ;
: <hamming-iterator> ( -- hamming-iterator )
hamming-iterator new
1 1dlist >>2s
1 1dlist >>3s
1 1dlist >>5s ;
: enqueue ( n hamming-iterator -- )
[ [ 2 * ] [ 2s>> ] bi* push-back ]
[ [ 3 * ] [ 3s>> ] bi* push-back ]
[ [ 5 * ] [ 5s>> ] bi* push-back ] 2tri ;
: next ( hamming-iterator -- n )
dup [ 2s>> ] [ 3s>> ] [ 5s>> ] tri
3dup [ peek-front ] tri@ min min
[
'[
dup peek-front _ =
[ pop-front* ] [ drop ] if
] tri@
] [ swap enqueue ] [ ] tri ;
: next-n ( hamming-iterator n -- seq )
swap '[ _ [ _ next , ] times ] { } make ;
: nth-from-now ( hamming-iterator n -- m )
1 - over '[ _ next drop ] times next ;
<hamming-iterator> 20 next-n . <hamming-iterator> 1691 nth-from-now . <hamming-iterator> 1000000 nth-from-now .
Lazy lists aren quite slow in Factor, but still.
USING: combinators fry kernel lists lists.lazy locals math ;
IN: rosetta.hamming-lazy
:: sort-merge ( xs ys -- result )
xs car :> x
ys car :> y
{
{ [ x y < ] [ [ x ] [ xs cdr ys sort-merge ] lazy-cons ] }
{ [ x y > ] [ [ y ] [ ys cdr xs sort-merge ] lazy-cons ] }
[ [ x ] [ xs cdr ys cdr sort-merge ] lazy-cons ]
} cond ;
:: hamming ( -- hamming )
f :> h!
[ 1 ] [
h 2 3 5 [ '[ _ * ] lazy-map ] tri-curry@ tri
sort-merge sort-merge
] lazy-cons h! h ;
20 hamming ltake list>array . 1690 hamming lnth . 999999 hamming lnth .
[edit] Forth
This version uses a compact representation of Hamming numbers: each 64-bit cell represents a number 2^l*3^m*5^n, where l, n, and m are bitfields in the cell (20 bits for now). It also uses a fixed-point logarithm to compare the Hamming numbers and prints them in factored form. This code has been tested up to the 10^9th Hamming number.
\ manipulating and computing with Hamming numbers:
: extract2 ( h -- l )
40 rshift ;
: extract3 ( h -- m )
20 rshift $fffff and ;
: extract5 ( h -- n )
$fffff and ;
' + alias h* ( h1 h2 -- h )
: h. { h -- }
." 2^" h extract2 0 .r
." *3^" h extract3 0 .r
." *5^" h extract5 . ;
\ the following numbers have been produced with bc -l as follows
1 62 lshift constant ldscale2
7309349404307464679 constant ldscale3 \ 2^62*l(3)/l(2) (rounded up)
10708003330985790206 constant ldscale5 \ 2^62*l(5)/l(2) (rounded down)
: hld { h -- ud }
\ ud is a scaled fixed-point representation of the logarithm dualis of h
h extract2 ldscale2 um*
h extract3 ldscale3 um* d+
h extract5 ldscale5 um* d+ ;
: h<= ( h1 h2 -- f )
2dup = if
2drop true exit
then
hld rot hld assert( 2over 2over d<> )
du>= ;
: hmin ( h1 h2 -- h )
2dup h<= if
drop
else
nip
then ;
\ actual algorithm
0 value seq
variable seqlast 0 seqlast !
: lastseq ( -- u )
\ last stored number in the sequence
seq seqlast @ th @ ;
: genseq ( h1 "name" -- )
\ h1 is the factor for the sequence
create , 0 , \ factor and index of element used for last return
does> ( -- u2 )
\ u2 is the next number resulting from multiplying h1 with numbers
\ in the sequence that is larger than the last number in the
\ sequence
dup @ lastseq { h1 l } cell+ dup @ begin ( index-addr index )
seq over th @ h1 h* dup l h<= while
drop 1+ repeat
>r swap ! r> ;
$10000000000 genseq s2
$00000100000 genseq s3
$00000000001 genseq s5
: nextseq ( -- )
s2 s3 hmin s5 hmin , 1 seqlast +! ;
: nthseq ( u1 -- h )
\ the u1 th element in the sequence
dup seqlast @ u+do
nextseq
loop
1- 0 max cells seq + @ ;
: .nseq ( u1 -- )
dup seqlast @ u+do
nextseq
loop
0 u+do
seq i th @ h.
loop ;
here to seq
0 , \ that's 1
20 .nseq
cr 1691 nthseq h.
cr 1000000 nthseq h.
- Output:
2^0*3^0*5^0 2^1*3^0*5^0 2^0*3^1*5^0 2^2*3^0*5^0 2^0*3^0*5^1 2^1*3^1*5^0 2^3*3^0*5^0 2^0*3^2*5^0 2^1*3^0*5^1 2^2*3^1*5^0 2^0*3^1*5^1 2^4*3^0*5^0 2^1*3^2*5^0 2^2*3^0*5^1 2^3*3^1*5^0 2^0*3^0*5^2 2^0*3^3*5^0 2^1*3^1*5^1 2^5*3^0*5^0 2^2*3^2*5^0 2^5*3^12*5^3 2^55*3^47*5^64
[edit] Fortran
Using big_integer_module from here[1]
program Hamming_Test
use big_integer_module
implicit none
call Hamming(1,20)
write(*,*)
call Hamming(1691)
write(*,*)
call Hamming(1000000)
contains
subroutine Hamming(first, last)
integer, intent(in) :: first
integer, intent(in), optional :: last
integer :: i, n, i2, i3, i5, lim
type(big_integer), allocatable :: hnums(:)
if(present(last)) then
lim = last
else
lim = first
end if
if(first < 1 .or. lim > 2500000 ) then
write(*,*) "Invalid input"
return
end if
allocate(hnums(lim))
i2 = 1 ; i3 = 1 ; i5 = 1
hnums(1) = 1
n = 1
do while(n < lim)
n = n + 1
hnums(n) = mini(2*hnums(i2), 3*hnums(i3), 5*hnums(i5))
if(2*hnums(i2) == hnums(n)) i2 = i2 + 1
if(3*hnums(i3) == hnums(n)) i3 = i3 + 1
if(5*hnums(i5) == hnums(n)) i5 = i5 + 1
end do
if(present(last)) then
do i = first, last
call print_big(hnums(i))
write(*, "(a)", advance="no") " "
end do
else
call print_big(hnums(first))
end if
deallocate(hnums)
end subroutine
function mini(a, b, c)
type(big_integer) :: mini
type(big_integer), intent(in) :: a, b, c
if(a < b ) then
if(a < c) then
mini = a
else
mini = c
end if
else if(b < c) then
mini = b
else
mini = c
end if
end function mini
end program
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Go
[edit] Concise version using dynamic-programming
package main
import (
"fmt"
"math/big"
)
func min(a, b *big.Int) *big.Int {
if a.Cmp(b) < 0 {
return a
}
return b
}
func hamming(n int) []*big.Int {
h := make([]*big.Int, n)
h[0] = big.NewInt(1)
two, three, five := big.NewInt(2), big.NewInt(3), big.NewInt(5)
next2, next3, next5 := big.NewInt(2), big.NewInt(3), big.NewInt(5)
i, j, k := 0, 0, 0
for m := 1; m < len(h); m++ {
h[m] = new(big.Int).Set(min(next2, min(next3, next5)))
if h[m].Cmp(next2) == 0 { i++; next2.Mul( two, h[i]) }
if h[m].Cmp(next3) == 0 { j++; next3.Mul(three, h[j]) }
if h[m].Cmp(next5) == 0 { k++; next5.Mul( five, h[k]) }
}
return h
}
func main() {
h := hamming(1e6)
fmt.Println(h[:20])
fmt.Println(h[1691-1])
fmt.Println(h[len(h)-1])
}
- Output:
[1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36] 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Longer version using dynamic-programming and logarithms
More than 10 times faster.
package main
import (
"flag"
"fmt"
"math"
"math/big"
"os"
)
var ordinal int // ordinal of last sequence element to compute
var sequenceMode bool // print the whole sequence or just one element?
var lg3, lg5 float64 // precomputed base-2 logarithms for 3 and 5
var table [][3]int16 // table for dynamic-programming stored results
var front [3]cursor // state of the three multiplied sequences
type cursor struct {
f int // index (0, 1, 2) corresponding to factor (2, 3, 5)
i int // index into table for the entry being multiplied
lg float64 // base-2 logarithm of the multiple (for ordering)
}
func (c *cursor) value() [3]int16 {
x := table[c.i]
x[c.f]++ // multiply by incrementing the exponent
return x
}
func (c *cursor) advance() {
c.i++
// skip entries that would produce duplicates
for (c.f < 2 && table[c.i][2] > 0) || (c.f < 1 && table[c.i][1] > 0) {
c.i++
}
x := c.value()
c.lg = float64(x[0]) + lg3*float64(x[1]) + lg5*float64(x[2])
}
func step() {
table = append(table, front[0].value())
front[0].advance()
// re-establish sorted order
if front[0].lg > front[1].lg {
front[0], front[1] = front[1], front[0]
if front[1].lg > front[2].lg {
front[1], front[2] = front[2], front[1]
}
}
}
func show(elem [3]int16) {
z := big.NewInt(1)
for i, base := range []int64{2, 3, 5} {
b := big.NewInt(base)
x := big.NewInt(int64(elem[i]))
z.Mul(z, b.Exp(b, x, nil))
}
fmt.Println(z)
}
func fail(msg string) {
fmt.Fprintf(os.Stderr, "%s: %s\n", os.Args[0], msg)
os.Exit(1)
}
func parse() {
flag.Parse()
if flag.NArg() != 1 {
fail("need one argument")
}
_, err := fmt.Sscan(flag.Arg(0), &ordinal)
if err != nil || ordinal <= 0 {
fail("argument must be a positive integer")
}
}
func init() {
flag.BoolVar(&sequenceMode, "s", false, "sequence mode")
lg3 = math.Log2(3)
lg5 = math.Log2(5)
front = [3]cursor{
{0, 0, 1}, // 2
{1, 0, lg3}, // 3
{2, 0, lg5}, // 5
}
}
func main() {
parse()
table = make([][3]int16, 1, ordinal)
for i, n := 1, ordinal; i < n; i++ {
if sequenceMode {
show(table[i-1])
}
step()
}
show(table[ordinal-1])
}
- Output:
$ ./hamming -s 20 | xargs 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 $ time ./hamming 1000000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 real 0m0.110s user 0m0.090s sys 0m0.020s $ uname -a Linux lance 3.0-ARCH #1 SMP PREEMPT Sat Aug 6 16:18:35 CEST 2011 x86_64 Intel(R) Core(TM)2 Duo CPU P8400 @ 2.26GHz GenuineIntel GNU/Linux
[edit] Haskell
[edit] The classic version
hamming = 1 : map (2*) hamming `union` map (3*) hamming `union` map (5*) hamming
union a@(x:xs) b@(y:ys) = case compare x y of
LT -> x : union xs b
EQ -> x : union xs ys
GT -> y : union a ys
main = do
print $ take 20 hamming
print (hamming !! (1691-1), hamming !! (1692-1))
print $ hamming !! (1000000-1)
-- Output:
-- [1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36]
-- (2125764000,2147483648)
-- 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Runs in about a second on Ideone.com. The nested unions' effect is to produce the minimal value at each step, without duplicates. As Haskell evaluation model is on-demand, the three map expressions are in effect iterators, maintaining hidden pointers back into the shared named storage with which they were each created (a name is a pointer/handle in Haskell; to name is to point at, to refer to, to take a handle on).
The amount of operations is constant for each number produced, so the time complexity should be O(n). Empirically, it is slightly above that and worsening, suggestive of extra cost of bignum arithmetics. Using triples representation with logarithm values for comparisons amends that problem, but runs ~ 1.5x slower for 1,000,000.
This is what that DDJ blog "pseudo-C" code was transcribing, mentioned at the Python entry that started this task. D, Go, PARI/GP, Prolog all implement the same idea of back-pointers into shared storage. A Haskell run-time system can actually free up the storage automatically at the start of the shared list and only keep the needed portion of it, from the (5*) back-pointer, behind the scenes – which is about O(n2 / 3) in length – as long as there's no re-use evident in the code.
[edit] Explicit multiples reinserting
This is a common approach which explicitly maintains an internal buffer of O(n2 / 3) elements, removing the numbers from its front and reinserting their 2- 3- and 5-multiples in order. It overproduces the sequence, stopping when the n-th number is no longer needed instead of when it's first found. Also overworks by maintaining this buffer in total order where just heap would be sufficient. Worse, this particular version uses a sequential list for its buffer. That means O(n2 / 3) operations for each number, instead of O(1) of the above version. Translation of Java (which does use priority queue though, so should have O (n logn) operations overall). Uses union from the "classic" version above:
hamm n = drop n $ iterate (\(_,(a:t))-> (a,union t [2*a,3*a,5*a])) (0,[1])
- Output:
*Main> map fst $ take 20 $ hamm 1
[1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36]
*Main> map fst $ take 2 $ hamm 1691
[2125764000,2147483648]
*Main> mapM_ print $ take 10 $ hamm 1
(1,[2,3,5])
(2,[3,4,5,6,10])
(3,[4,5,6,9,10,15])
(4,[5,6,8,9,10,12,15,20])
(5,[6,8,9,10,12,15,20,25])
(6,[8,9,10,12,15,18,20,25,30])
(8,[9,10,12,15,16,18,20,24,25,30,40])
(9,[10,12,15,16,18,20,24,25,27,30,40,45])
(10,[12,15,16,18,20,24,25,27,30,40,45,50])
(12,[15,16,18,20,24,25,27,30,36,40,45,50,60])
*Main> map (length.snd.head.hamm) [2000,4000,8000,16000]
[402,638,1007,1596]
Runs too slowly to reach 1,000,000, with empirical complexity above O (n 1.7 ) and worsening. Last two lines show the internal buffer's length for several sample n s.
[edit] Direct calculation through triples enumeration
It is also possible to more or less directly calculate the n-th Hamming number by enumerating (and counting) all the (i,j,k) triples below its estimated value – with ordering according to their exponents, i*ln2 + j*ln3 + k*ln5 – while storing only the "band" of topmost triples close enough to the target value (more in the original post on DDJ).
The total count of the enumerated triples is then the band's topmost value's position in the Hamming sequence, 1-based. The nth number is then found by a simple lookup in the sorted band, if it was wide enough. This produces the 1,000,000-th value in a few hundredths of a second on Ideone.com, running at about O(n0.7) empirical time and space complexity:
-- directly find n-th Hamming number, in ~ O(n^{2/3}) time
-- by Will Ness, based on "top band" idea by Louis Klauder, from DDJ discussion
-- http://drdobbs.com/blogs/architecture-and-design/228700538
{-# OPTIONS -O2 -XBangPatterns #-}
import Data.List (sortBy)
import Data.Function (on)
main = do let (r,t) = nthHam 1000000
sequence_ [print t, print $ trival t]
lb3 = logBase 2 3; lb5 = logBase 2 5
logval (i,j,k) = fromIntegral i + fromIntegral j*lb3 + fromIntegral k*lb5
trival (i,j,k) = 2^i * 3^j * 5^k
estval n = (6*lb3*lb5* fromIntegral n)**(1/3) -- estimated logval, base 2
rngval n
| n > 500000 = (2.4496 , 0.0076 ) -- empirical estimation
| n > 50000 = (2.4424 , 0.0146 ) -- correction, base 2
| n > 500 = (2.3948 , 0.0723 ) -- (dist,width)
| n > 1 = (2.2506 , 0.2887 ) -- around (log $ sqrt 30),
| otherwise = (2.2506 , 0.5771 ) -- says WP
nthHam n -- n: 1-based 1,2,3...
| w >= 1 = error $ "Breach of contract: (w < 1): " ++ show w
| m < 0 = error $ "Not enough triples generated: " ++ show (c,n)
| m >= nb = error $ "Generated band is too narrow: " ++ show (m,nb)
| True = res
where
(d,w) = rngval n -- correction dist, width
hi = estval n - d -- hi > logval > hi-w
(m,nb) = ( fromInteger $ c - n, length b ) -- m 0-based from top, |band|
(s,res) = ( sortBy (flip compare `on` fst) b, s!!m ) -- sorted decreasing, result
(c,b) = -- f 0 == prod (sum,concat) . unzip ; prod(f,g)(x,y) = (f x,g y)
f 0 -- total count, the band
[ ( i+1, -- total triples w/ this (j,k)
[ (r,(i,j,k)) | frac < w ] ) -- store it, if inside band
| k <- [ 0 .. floor ( hi /lb5) ], let p = fromIntegral k*lb5,
j <- [ 0 .. floor ((hi-p)/lb3) ], let q = fromIntegral j*lb3 + p,
let (i,frac) = pr (hi-q) ; r = hi - frac ] -- r = i + q
where pr = properFraction
f !c [] = (c,[])
f !c ((c1,b1):r) = let (cr,br) = f (c+c1) r
in case b1 of { [v] -> (cr,v:br)
; _ -> (cr, br) }
- Output:
-- time: 0.01s memory: 3688 kB (55,47,64) 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Icon and Unicon
This solution uses Unicon's object oriented extensions. An Icon only version has not been provided.
Lazy evaluation is used to improve performance.
# Lazily generate the three Hamming numbers that can be derived directly
# from a known Hamming number h
class Triplet : Class (cv, ce)
method nextVal()
suspend cv := @ce
end
initially (baseNum)
cv := 2*baseNum
ce := create (3|5)*baseNum
end
# Generate Hamming numbers, in order. Default is first 30
# But an optional argument can be used to generate more (or less)
# e.g. hamming 5000 generates the first 5000.
procedure main(args)
limit := integer(args[1]) | 30
every write("\t", generateHamming() \ limit)
end
# Do the work. Start with known Hamming number 1 and maintain
# a set of triplet Hamming numbers as they get derived from that
# one. Most of the code here is to figure out which Hamming
# number is next in sequence (while removing duplicates)
procedure generateHamming()
triplers := set()
insert(triplers, Triplet(1))
suspend 1
repeat {
# Pick a Hamming triplet that *may* have the next smallest number
t1 := !triplers # any will do to start
every t1 ~=== (t2 := !triplers) do {
if t1.cv > t2.cv then {
# oops we were wrong, switch assumption
t1 := t2
}
else if t1.cv = t2.cv then {
# t2's value is a duplicate, so
# advance triplet t2, if none left in t2, remove it
t2.nextVal() | delete(triplers, t2)
}
}
# Ok, t1 has the next Hamming number, grab it
suspend t1.cv
insert(triplers, Triplet(t1.cv))
# Advance triplet t1, if none left in t1, remove it
t1.nextVal() | delete(triplers, t1)
}
end
[edit] J
Solution:
A concise tacit expression using a (right) fold:
hamming=: {. (/:~@~.@], 2 3 5 * {)/ @ (1x,~i.@-)
Example usage:
hamming 20
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36
{: hamming 1691
2125764000
For the millionth (and billionth (1e9)) Hamming number see the hn verb from Hamming Number essay on the J wiki.
Explanation:
I'll explain this J-sentence by dividing it in three parts from left to right omitting the leftmost {.:
- sort and remove duplicates
/:~@~.@]
- produce (the next) 3 elements by selection and multiplication:
2 3 5 * {
or
2 3 5 * LHA { RHA
- the RH part forms an array of descending indices and the initial Hamming number 1
(1x,~i.@-)
e.g. if we want the first 5 Hamming numbers, it produces the array:
4 3 2 1 0 1
in other words, we compute a sequence which begins with the desired hamming sequence and then take the first n elements (which will be our desired hamming sequence)
({. (/:~@~.@], 2 3 5 * {)/ @ (1x,~i.@-)) 7
1 2 3 4 5 6 8</lang j>
This starts using a descending sequence with 1 appended:
<lang j> (1x,~i.@-) 7
6 5 4 3 2 1 0 1
and then the fold expression is inserted between these list elements and the result computed:
6(/:~@~.@], 2 3 5 * {) 5(/:~@~.@], 2 3 5 * {) 4(/:~@~.@], 2 3 5 * {) 3(/:~@~.@], 2 3 5 * {) 2(/:~@~.@], 2 3 5 * {) 1(/:~@~.@], 2 3 5 * {) 0(/:~@~.@], 2 3 5 * {) 1
1 2 3 4 5 6 8 9 10 12 15 18 20 25 30 16 24 40
(Note: A train of verbs in J is evaluated by supplying arguments to the every other verb (counting from the right) and the combining these results with the remaining verbs. Also: { has been implemented so that an index of 0 will select the only item from an array with no dimensions.)
[edit] Java
Has a common shortcoming of overproducing the sequence by about O(n2 / 3) entries, until the n-th number is no longer needed, instead of stopping as soon as it is reached. See Haskell for an illustration.
Inserting the top number's three multiples deep into the priority queue as it does, incurs extra cost for each number produced. To not worsen the expected algorithm complexity, the priority queue should have (amortized) O(1) implementation for both insertion and deletion operations but it looks like it's O(logn) in Java.
import java.math.BigInteger;
import java.util.PriorityQueue;
final class Hamming {
private static BigInteger THREE = BigInteger.valueOf(3);
private static BigInteger FIVE = BigInteger.valueOf(5);
private static void updateFrontier(BigInteger x,
PriorityQueue<BigInteger> pq) {
pq.offer(x.shiftLeft(1));
pq.offer(x.multiply(THREE));
pq.offer(x.multiply(FIVE));
}
public static BigInteger hamming(int n) {
if (n <= 0)
throw new IllegalArgumentException("Invalid parameter");
PriorityQueue<BigInteger> frontier = new PriorityQueue<BigInteger>();
updateFrontier(BigInteger.ONE, frontier);
BigInteger lowest = BigInteger.ONE;
for (int i = 1; i < n; i++) {
lowest = frontier.poll();
while (frontier.peek().equals(lowest))
frontier.poll();
updateFrontier(lowest, frontier);
}
return lowest;
}
public static void main(String[] args) {
System.out.print("Hamming(1 .. 20) =");
for (int i = 1; i < 21; i++)
System.out.print(" " + hamming(i));
System.out.println("\nHamming(1691) = " + hamming(1691));
System.out.println("Hamming(1000000) = " + hamming(1000000));
}
}
- Output:
Hamming(1 .. 20) = 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 Hamming(1691) = 2125764000 Hamming(1000000) = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] JavaScript
This does not calculate the 1,000,000th Hamming number.
Note the use of for (x in obj) to iterate over the properties of an object, versus for each (y in obj) to iterate over the values of the properties of an object.
function hamming() {
var queues = {2: [], 3: [], 5: []};
var base;
var next_ham = 1;
while (true) {
yield next_ham;
for (base in queues) {queues[base].push(next_ham * base)}
next_ham = [ queue[0] for each (queue in queues) ].reduce(function(min, val) {
return Math.min(min,val)
});
for (base in queues) {if (queues[base][0] == next_ham) queues[base].shift()}
}
}
var ham = hamming();
var first20=[], i=1;
for (; i <= 20; i++)
first20.push(ham.next());
print(first20.join(', '));
print('...');
for (; i <= 1690; i++)
ham.next();
print(i + " => " + ham.next());
- Output:
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36 ... 1691 => 2125764000
[edit] Liberty BASIC
LB has unlimited precision integers.
dim h( 1000000)
for i =1 to 20
print hamming( i); " ";
next i
print "H( 1691)", hamming( 1691)
print "H( 1000000)", hamming( 1000000)
end
function hamming( limit)
h( 0) =1
x2 =2: x3 =3: x5 =5
i =0: j =0: k =0
for n =1 to limit
h( n) = min( x2, min( x3, x5))
if x2 = h( n) then i = i +1: x2 =2 *h( i)
if x3 = h( n) then j = j +1: x3 =3 *h( j)
if x5 = h( n) then k = k +1: x5 =5 *h( k)
next n
hamming =h( limit -1)
end function
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 H( 1691) 2125764000 H( 1000000) 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Logo
to init.ham
; queues
make "twos [1]
make "threes [1]
make "fives [1]
end
to next.ham
localmake "ham first :twos
if less? first :threes :ham [make "ham first :threes]
if less? first :fives :ham [make "ham first :fives]
if equal? :ham first :twos [ignore dequeue "twos]
if equal? :ham first :threes [ignore dequeue "threes]
if equal? :ham first :fives [ignore dequeue "fives]
queue "twos :ham * 2
queue "threes :ham * 3
queue "fives :ham * 5
output :ham
end
init.ham
repeat 20 [print next.ham]
repeat 1690-20 [ignore next.ham]
print next.ham
[edit] Lua
function hiter()
hammings = {1}
prev, vals = {1, 1, 1}
index = 1
local function nextv()
local n, v = 1, hammings[prev[1]]*2
if hammings[prev[2]]*3 < v then n, v = 2, hammings[prev[2]]*3 end
if hammings[prev[3]]*5 < v then n, v = 3, hammings[prev[3]]*5 end
prev[n] = prev[n] + 1
if hammings[index] == v then return nextv() end
index = index + 1
hammings[index] = v
return v
end
return nextv
end
j = hiter()
for i = 1, 20 do
print(j())
end
n, l = 0, 0
while n < 2^31 do n, l = j(), n end
print(l)
[edit] MUMPS
Hamming(n) New count,ok,next,number,which
For which=2,3,5 Set number=1
For count=1:1:n Do
. Set ok=0 Set:count<21 ok=1 Set:count=1691 ok=1 Set:count=n ok=1
. Write:ok !,$Justify(count,5),": ",number
. For which=2,3,5 Set next(number*which)=which
. Set number=$Order(next(""))
. Kill next(number)
. Quit
Quit
Do Hamming(2000)
1: 1
2: 2
3: 3
4: 4
5: 5
6: 6
7: 8
8: 9
9: 10
10: 12
11: 15
12: 16
13: 18
14: 20
15: 24
16: 25
17: 27
18: 30
19: 32
20: 36
1691: 2125764000
2000: 8062156800
[edit] Oz
declare
fun lazy {HammingFun}
1|{FoldL1 [{MultHamming 2} {MultHamming 3} {MultHamming 5}] LMerge}
end
Hamming = {HammingFun}
fun {MultHamming N}
{LMap Hamming fun {$ X} N*X end}
end
fun lazy {LMap Xs F}
case Xs
of nil then nil
[] X|Xr then {F X}|{LMap Xr F}
end
end
fun lazy {LMerge Xs=X|Xr Ys=Y|Yr}
if X < Y then X|{LMerge Xr Ys}
elseif X > Y then Y|{LMerge Xs Yr}
else X|{LMerge Xr Yr}
end
end
fun {FoldL1 X|Xr F}
{FoldL Xr F X}
end
in
{ForAll {List.take Hamming 20} System.showInfo}
{System.showInfo {Nth Hamming 1690}}
{System.showInfo {Nth Hamming 1000000}}
[edit] PARI/GP
This is a basic implementation; finding the millionth term requires 3 seconds and 54 MB. Much better algorithms exist.
Hupto(n)={
my(v=vector(n),x2=2,x3=3,x5=5,i=1,j=1,k=1,t);
v[1]=1;
for(m=2,n,
v[m]=t=min(x2,min(x3,x5));
if(x2 == t, x2 = v[i++] << 1);
if(x3 == t, x3 = 3 * v[j++]);
if(x5 == t, x5 = 5 * v[k++]);
);
v
};
H(n)=Hupto(n)[n];
Hupto(20)
H(1691)
H(10^6)
- Output:
%1 = [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] %2 = 2125764000 %3 = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Perl
use List::Util 'min';
sub ham_gen {
my @s = ([1], [1], [1]);
my @m = (2, 3, 5);
return sub {
# use bigint;
my $n = min($s[0][0], $s[1][0], $s[2][0]);
for (0 .. 2) {
shift @{$s[$_]} if $s[$_][0] == $n;
push @{$s[$_]}, $n * $m[$_]
}
return $n
}
}
my ($h, $i) = ham_gen;
++$i, print $h->(), " " until $i > 20;
print "...\n";
++$i, $h->() until $i == 1690;
print ++$i, "-th: ", $h->(), "\n";
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 40 ... 1691-th: 2125764000 <perl's bigint is so horribly slow that I didn't have the patience for the 1000000th output> <some big number presumably>
[edit] Perl 6
The limit scaling is not required, but it cuts down on a bunch of unnecessary calculation.
my $limit = 32;
sub powers_of ($radix) { 1, [\*] $radix xx * }
my @hammings =
( powers_of(2)[^ $limit ] X*
( powers_of(3)[^($limit * 2/3)] X*
powers_of(5)[^($limit * 1/2)]
)
).sort;
say ~@hammings[^20];
say @hammings[1690]; # zero indexed
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000
[edit] PicoLisp
(de hamming (N)
(let (L (1) H)
(do N
(for (X L X (cadr X)) # Find smallest result
(setq H (car X)) )
(idx 'L H NIL) # Remove it
(for I (2 3 5) # Generate next results
(idx 'L (* I H) T) ) )
H ) )
(println (make (for N 20 (link (hamming N)))))
(println (hamming 1691))
(println (hamming 1000000))
- Output:
(1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36) 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Prolog
[edit] Generator idiom
%% collect N elements produced by a generator in a row
take( 0, Next, Z-Z, Next).
take( N, Next, [A|B]-Z, NZ):- N>0, !, next(Next,A,Next1),
N1 is N-1,
take(N1,Next1,B-Z,NZ).
%% a generator provides specific {next} implementation
next( hamm( A2,B,C3,D,E5,F,[H|G] ), H, hamm(X,U,Y,V,Z,W,G) ):-
H is min(A2, min(C3,E5)),
( A2 =:= H -> B=[N2|U],X is N2*2 ; (X,U)=(A2,B) ),
( C3 =:= H -> D=[N3|V],Y is N3*3 ; (Y,V)=(C3,D) ),
( E5 =:= H -> F=[N5|W],Z is N5*5 ; (Z,W)=(E5,F) ).
mkHamm( hamm(1,X,1,X,1,X,X) ). % Hamming numbers generator init state
main(N) :-
mkHamm(G),take(20,G,A-[],_), write(A), nl,
take(1691-1,G,_,G2),take(2,G2,B-[],_), write(B), nl,
take( N -1,G,_,G3),take(2,G3,[C1|_]-_,_), write(C1), nl.
Run on Ideone.com, which uses Prolog (swi) (swipl 5.6.64), it produces:
?- main(1000000). [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] [2125764000, 2147483648] ERROR: */2: Arithmetic: evaluation error: `int_overflow'
[edit] Laziness flavor
Works with SWI-Prolog. Laziness is simulate with freeze/2 and ground/2.
Took inspiration from this code : http://chr.informatik.uni-ulm.de/~webchr (click on hamming.pl: Solves Hamming Problem).
hamming(N) :-
% to stop cleanly
nb_setval(go, 1),
% display list
( N = 20 -> watch_20(20, L); watch(1,N,L)),
% go
L=[1|L235],
multlist(L,2,L2),
multlist(L,3,L3),
multlist(L,5,L5),
merge_(L2,L3,L23),
merge_(L5,L23,L235).
%% multlist(L,N,LN)
%% multiply each element of list L with N, resulting in list LN
%% here only do multiplication for 1st element, then use multlist recursively
multlist([X|L],N,XLN) :-
% the trick to stop
nb_getval(go, 1) ->
% laziness flavor
when(ground(X),
( XN is X*N,
XLN=[XN|LN],
multlist(L,N,LN)));
true.
merge_([X|In1],[Y|In2],XYOut) :-
% the trick to stop
nb_getval(go, 1) ->
% laziness flavor
( X < Y -> XYOut = [X|Out], In11 = In1, In12 = [Y|In2]
; X = Y -> XYOut = [X|Out], In11 = In1, In12 = In2
; XYOut = [Y|Out], In11 = [X | In1], In12 = In2),
freeze(In11,freeze(In12, merge_(In11,In12,Out)));
true.
%% display nth element
watch(Max, Max, [X|_]) :-
% laziness flavor
when(ground(X),
(format('~w~n', [X]),
% the trick to stop
nb_linkval(go, 0))).
watch(N, Max, [_X|L]):-
N1 is N + 1,
watch(N1, Max, L).
%% display nth element
watch_20(1, [X|_]) :-
% laziness flavor
when(ground(X),
(format('~w~n', [X]),
% the trick to stop
nb_linkval(go, 0))).
watch_20(N, [X|L]):-
% laziness flavor
when(ground(X),
(format('~w ', [X]),
N1 is N - 1,
watch_20(N1, L))).
- Output:
?- hamming(20). 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 true . ?- hamming(1691). 2125764000 true . ?- hamming(1000000). 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 true .
[edit] Python
[edit] Version based on example from Dr. Dobb's CodeTalk
from itertools import islice
def hamming2():
'''\
This version is based on a snippet from:
http://dobbscodetalk.com/index.php?option=com_content&task=view&id=913&Itemid=85
When expressed in some imaginary pseudo-C with automatic
unlimited storage allocation and BIGNUM arithmetics, it can be
expressed as:
hamming = h where
array h;
n=0; h[0]=1; i=0; j=0; k=0;
x2=2*h[ i ]; x3=3*h[j]; x5=5*h[k];
repeat:
h[++n] = min(x2,x3,x5);
if (x2==h[n]) { x2=2*h[++i]; }
if (x3==h[n]) { x3=3*h[++j]; }
if (x5==h[n]) { x5=5*h[++k]; }
'''
h = 1
_h=[h] # memoized
multipliers = (2, 3, 5)
multindeces = [0 for i in multipliers] # index into _h for multipliers
multvalues = [x * _h[i] for x,i in zip(multipliers, multindeces)]
yield h
while True:
h = min(multvalues)
_h.append(h)
for (n,(v,x,i)) in enumerate(zip(multvalues, multipliers, multindeces)):
if v == h:
i += 1
multindeces[n] = i
multvalues[n] = x * _h[i]
# cap the memoization
mini = min(multindeces)
if mini >= 1000:
del _h[:mini]
multindeces = [i - mini for i in multindeces]
#
yield h
- Output:
>>> list(islice(hamming2(), 20)) [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] >>> list(islice(hamming2(), 1689, 1690)) [2123366400] >>> list(islice(hamming2(), 999999, 1000000)) [519312780448388736089589843750000000000000000000000000000000000000000000000000000000]
[edit] Another implementation of same approach
This is the fastest of the Python implementations, it uses a lot of memory. Does not try to limit memory usage.
import psyco
def hamming(limit):
h = [1] * limit
x2, x3, x5 = 2, 3, 5
i = j = k = 0
for n in xrange(1, limit):
h[n] = min(x2, x3, x5)
if x2 == h[n]:
i += 1
x2 = 2 * h[i]
if x3 == h[n]:
j += 1
x3 = 3 * h[j]
if x5 == h[n]:
k += 1
x5 = 5 * h[k]
return h[-1]
psyco.bind(hamming)
print [hamming(i) for i in xrange(1, 21)]
print hamming(1691)
print hamming(1000000)
[edit] "Cyclical Iterators"
The original author is Raymond Hettinger and the code was first published here under the MIT license. Uses iterators dubbed "cyclical" in a sense that they are referring back (explicitly, with p2, p3, p5 iterators) to the previously produced values, same as the above versions (through indecies into shared storage) and the classic Haskell version (implicitly timed by lazy evaluation).
Memory is efficiently maintained automatically by the tee function for each of the three generator expressions, i.e. only that much is maintained as needed to produce the next value (although it looks like the storage is not shared so three copies are maintained implicitly there).
from itertools import tee, chain, groupby, islice
from heapq import merge
def raymonds_hamming():
# Generate "5-smooth" numbers, also called "Hamming numbers"
# or "Regular numbers". See: http://en.wikipedia.org/wiki/Regular_number
# Finds solutions to 2**i * 3**j * 5**k for some integers i, j, and k.
def deferred_output():
for i in output:
yield i
result, p2, p3, p5 = tee(deferred_output(), 4)
m2 = (2*x for x in p2) # multiples of 2
m3 = (3*x for x in p3) # multiples of 3
m5 = (5*x for x in p5) # multiples of 5
merged = merge(m2, m3, m5)
combined = chain([1], merged) # prepend a starting point
output = (k for k,g in groupby(combined)) # eliminate duplicates
return result
print list(islice(raymonds_hamming(), 20))
print islice(raymonds_hamming(), 1689, 1690).next()
print islice(raymonds_hamming(), 999999, 1000000).next()
Results are the same as before. Another formulation along the same lines, but greatly simplified, can be found here.
[edit] Qi
(define smerge
[X|Xs] [Y|Ys] -> [X | (freeze (smerge (thaw Xs) [Y|Ys]))] where (< X Y)
[X|Xs] [Y|Ys] -> [Y | (freeze (smerge [X|Xs] (thaw Ys)))] where (> X Y)
[X|Xs] [_|Ys] -> [X | (freeze (smerge (thaw Xs) (thaw Ys)))])
(define smerge3
Xs Ys Zs -> (smerge Xs (smerge Ys Zs)))
(define smap
F [S|Ss] -> [(F S)|(freeze (smap F (thaw Ss)))])
(set hamming [1 | (freeze (smerge3 (smap (* 2) (value hamming))
(smap (* 3) (value hamming))
(smap (* 5) (value hamming))))])
(define stake
_ 0 -> []
[S|Ss] N -> [S|(stake (thaw Ss) (1- N))])
(stake (value hamming) 20)
- Output:
[1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36]
[edit] R
Recursively find the Hamming numbers below 231. Works for larger limits, however to find the 1000000th value, one needs guess the correct limit
hamming=function(hamms,limit) {
tmp=hamms
for(h in c(2,3,5)) {
tmp=c(tmp,h*hamms)
}
tmp=unique(tmp[tmp<=limit])
if(length(tmp)>length(hamms)) {
hamms=hamming(tmp,limit)
}
hamms
}
sort(hamming(1,limit=2^31)[-1])
[edit] REXX
numeric digits 100
call hamming 1,20
call hamming 1691
call hamming 1000000
exit
hamming: procedure; parse arg x,y
if y=='' then y=x
p2=1
p3=1
p5=1
a.=0
a.1=1
do n=2 for y-1
a.n=min(2*a.p2, 3*a.p3, 5*a.p5)
if 2*a.p2==a.n then p2=p2+1
if 3*a.p3==a.n then p3=p3+1
if 5*a.p5==a.n then p5=p5+1
end
say
do j=x to y
say 'hamming('j")=" a.j
end
say
say right('length of last hamming number='length(a.y),70)
say
return
- Output:
hamming(1)= 1
hamming(2)= 2
hamming(3)= 3
hamming(4)= 4
hamming(5)= 5
hamming(6)= 6
hamming(7)= 8
hamming(8)= 9
hamming(9)= 10
hamming(10)= 12
hamming(11)= 15
hamming(12)= 16
hamming(13)= 18
hamming(14)= 20
hamming(15)= 24
hamming(16)= 25
hamming(17)= 27
hamming(18)= 30
hamming(19)= 32
hamming(20)= 36
length of last hamming number=2
hamming(1691)= 2125764000
length of last hamming number=10
hamming(1000000)= 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
length of last hamming number=84
[edit] Ruby
require 'generator'
# the Hamming number generator
hamming = Generator.new do |generator|
next_ham = 1
queues = { 2 => [], 3 => [], 5 => [] }
loop do
generator.yield next_ham
[2,3,5].each {|m| queues[m] << (next_ham * m)}
next_ham = [2,3,5].collect {|m| queues[m][0]}.min
[2,3,5].each {|m| queues[m].shift if queues[m][0] == next_ham}
end
end
This method does not require a library module.
hamming = Enumerator.new do |yielder|
next_ham = 1
queues = { 2 => [], 3 => [], 5 => [] }
loop do
yielder << next_ham # or: yielder.yield(next_ham)
[2,3,5].each {|m| queues[m]<< (next_ham * m)}
next_ham = [2,3,5].collect {|m| queues[m][0]}.min
[2,3,5].each {|m| queues[m].shift if queues[m][0]== next_ham}
end
end
And the "main" part of the task
start = Time.now
idx = 1
hamming.each do |ham|
case idx
when (1..20), 1691
p [idx, ham]
when 1_000_000
p [idx, ham]
break
end
idx += 1
end
puts "elapsed: #{Time.now - start} seconds"
- Output:
[1, 1] [2, 2] [3, 3] [4, 4] [5, 5] [6, 6] [7, 8] [8, 9] [9, 10] [10, 12] [11, 15] [12, 16] [13, 18] [14, 20] [15, 24] [16, 25] [17, 27] [18, 30] [19, 32] [20, 36] [1691, 2125764000] [1000000, 519312780448388736089589843750000000000000000000000000000000000000000000000000000000] elapsed: 143.96875 seconds
[edit] Scala
class Hamming extends Iterator[BigInt] {
import scala.collection.mutable.Queue
val qs = Seq.fill(3)(new Queue[BigInt])
def enqueue(n: BigInt) = qs zip Seq(2, 3, 5) foreach { case (q, m) => q enqueue n * m }
def next = {
val n = qs map (_.head) min;
qs foreach { q => if (q.head == n) q.dequeue }
enqueue(n)
n
}
def hasNext = true
qs foreach (_ enqueue 1)
}
However, the usage of closures adds a significant amount of time. The code below, though a bit uglier because of the repetitions, is twice as fast:
class Hamming extends Iterator[BigInt] {
import scala.collection.mutable.Queue
val q2 = new Queue[BigInt]
val q3 = new Queue[BigInt]
val q5 = new Queue[BigInt]
def enqueue(n: BigInt) = {
q2 enqueue n * 2
q3 enqueue n * 3
q5 enqueue n * 5
}
def next = {
val n = q2.head min q3.head min q5.head
if (q2.head == n) q2.dequeue
if (q3.head == n) q3.dequeue
if (q5.head == n) q5.dequeue
enqueue(n)
n
}
def hasNext = true
List(q2, q3, q5) foreach (_ enqueue 1)
}
Usage:
scala> new Hamming take 20 toList res87: List[BigInt] = List(1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36) scala> new Hamming drop 1690 next res88: BigInt = 2125764000 scala> new Hamming drop 999999 next res89: BigInt = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
There's also a fairly mechanical translation from Haskell using purely functional lazy streams
val hamming : Stream[BigInt] = {
def merge(inx : Stream[BigInt], iny : Stream[BigInt]) : Stream[BigInt] = {
if (inx.head < iny.head) inx.head #:: merge(inx.tail, iny) else
if (iny.head < inx.head) iny.head #:: merge(inx, iny.tail) else
merge(inx, iny.tail)
}
1 #:: merge(hamming map (_ * 2), merge(hamming map (_ * 3), hamming map (_ * 5)))
}
Use of "force" ensures that the stream is computed before being printed, otherwise it would just be left suspended and you'd see "Stream(1, ?)"
scala> (hamming take 20).force res0: scala.collection.immutable.Stream[BigInt] = Stream(1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36)
To get the nth code find the n-1th element because indexes are 0 based
scala> hamming(1690) res1: BigInt = 2125764000
To calculate the 1000000th code I had to increase the JVM heap from the default
scala> hamming(999999) res2: BigInt = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Scheme
(define-syntax lons
(syntax-rules ()
((_ lar ldr) (delay (cons lar (delay ldr))))))
(define (lar lons)
(car (force lons)))
(define (ldr lons)
(force (cdr (force lons))))
(define (lap proc . llists)
(lons (apply proc (map lar llists)) (apply lap proc (map ldr llists))))
(define (take n llist)
(if (zero? n)
(list)
(cons (lar llist) (take (- n 1) (ldr llist)))))
(define (llist-ref n llist)
(if (= n 1)
(lar llist)
(llist-ref (- n 1) (ldr llist))))
(define (merge llist-1 . llists)
(define (merge-2 llist-1 llist-2)
(cond ((null? llist-1) llist-2)
((null? llist-2) llist-1)
((< (lar llist-1) (lar llist-2))
(lons (lar llist-1) (merge-2 (ldr llist-1) llist-2)))
((> (lar llist-1) (lar llist-2))
(lons (lar llist-2) (merge-2 llist-1 (ldr llist-2))))
(else (lons (lar llist-1) (merge-2 (ldr llist-1) (ldr llist-2))))))
(if (null? llists)
llist-1
(apply merge (cons (merge-2 llist-1 (car llists)) (cdr llists)))))
(define hamming
(lons 1
(merge (lap (lambda (x) (* x 2)) hamming)
(lap (lambda (x) (* x 3)) hamming)
(lap (lambda (x) (* x 5)) hamming))))
(display (take 20 hamming))
(newline)
(display (llist-ref 1691 hamming))
(newline)
(display (llist-ref 1000000 hamming))
(newline)
- Output:
(1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36) 2125764000 out of memory
[edit] Seed7
$ include "seed7_05.s7i";
include "bigint.s7i";
const func bigInteger: min (in bigInteger: a, in bigInteger: b, in bigInteger: c) is func
result
var bigInteger: min is 0_;
begin
if a < b then
min := a;
else
min := b;
end if;
if c < min then
min := c;
end if;
end func;
const func bigInteger: hamming (in integer: n) is func
result
var bigInteger: hammingNum is 1_;
local
var array bigInteger: hammingNums is 0 times 0_;
var integer: index is 0;
var bigInteger: x2 is 2_;
var bigInteger: x3 is 3_;
var bigInteger: x5 is 5_;
var integer: i is 1;
var integer: j is 1;
var integer: k is 1;
begin
hammingNums := n times 1_;
for index range 2 to n do
hammingNum := min(x2, x3, x5);
hammingNums[index] := hammingNum;
if x2 = hammingNum then
incr(i);
x2 := 2_ * hammingNums[i];
end if;
if x3 = hammingNum then
incr(j);
x3 := 3_ * hammingNums[j];
end if;
if x5 = hammingNum then
incr(k);
x5 := 5_ * hammingNums[k];
end if;
end for;
end func;
const proc: main is func
local
var integer: n is 0;
begin
for n range 1 to 20 do
write(hamming(n) <& " ");
end for;
writeln;
writeln(hamming(1691));
writeln(hamming(1000000));
end func;
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
[edit] Smalltalk
This is a straightforward implementation of the pseudocode snippet found in the Python section. Smalltalk supports arbitrary-precision integers, but the implementation is too slow to try it with 1 million.
Object subclass: Hammer [
Hammer class >> hammingNumbers: howMany [
|h i j k x2 x3 x5|
h := OrderedCollection new.
i := 0. j := 0. k := 0.
h add: 1.
x2 := 2. x3 := 2. x5 := 5.
[ ( h size) < howMany ] whileTrue: [
|m|
m := { x2. x3. x5 } sort first.
(( h indexOf: m ) = 0) ifTrue: [ h add: m ].
( x2 = (h last) ) ifTrue: [ i := i + 1. x2 := 2 * (h at: i) ].
( x3 = (h last) ) ifTrue: [ j := j + 1. x3 := 3 * (h at: j) ].
( x5 = (h last) ) ifTrue: [ k := k + 1. x5 := 5 * (h at: k) ].
].
^ h sort
]
].
(Hammer hammingNumbers: 20) displayNl.
(Hammer hammingNumbers: 1690) last displayNl.
[edit] Tcl
This uses coroutines to simplify the description of what's going on.
package require Tcl 8.6
# Simple helper: Tcl-style list "map"
proc map {varName list script} {
set l {}
upvar 1 $varName v
foreach v $list {lappend l [uplevel 1 $script]}
return $l
}
# The core of a coroutine to compute the product of a hamming sequence.
#
# Tricky bit: we don't automatically advance to the next value, and instead
# wait to be told that the value has been consumed (i.e., is the result of
# the [yield] operation).
proc ham {key multiplier} {
global hammingCache
set i 0
yield [info coroutine]
# Cannot use [foreach]; that would take a snapshot of the list in
# the hammingCache variable, so missing updates.
while 1 {
set n [expr {[lindex $hammingCache($key) $i] * $multiplier}]
# If the number selected was ours, we advance to compute the next
if {[yield $n] == $n} {
incr i
}
}
}
# This coroutine computes the hamming sequence given a list of multipliers.
# It uses the [ham] helper from above to generate indivdual multiplied
# sequences. The key into the cache is the list of multipliers.
#
# Note that it is advisable for the values to be all co-prime wrt each other.
proc hammingCore args {
global hammingCache
set hammingCache($args) 1
set hammers [map x $args {coroutine ham$x,$args ham $args $x}]
yield
while 1 {
set n [lindex $hammingCache($args) [incr i]-1]
lappend hammingCache($args) \
[tcl::mathfunc::min {*}[map h $hammers {$h $n}]]
yield $n
}
}
# Assemble the pieces so as to compute the classic hamming sequence.
coroutine hamming hammingCore 2 3 5
# Print the first 20 values of the sequence
for {set i 1} {$i <= 20} {incr i} {
puts [format "hamming\[%d\] = %d" $i [hamming]]
}
for {} {$i <= 1690} {incr i} {set h [hamming]}
puts "hamming{1690} = $h"
for {} {$i <= 1000000} {incr i} {set h [hamming]}
puts "hamming{1000000} = $h"
- Output:
hamming{1} = 1
hamming{2} = 2
hamming{3} = 3
hamming{4} = 4
hamming{5} = 5
hamming{6} = 6
hamming{7} = 8
hamming{8} = 9
hamming{9} = 10
hamming{10} = 12
hamming{11} = 15
hamming{12} = 16
hamming{13} = 18
hamming{14} = 20
hamming{15} = 24
hamming{16} = 25
hamming{17} = 27
hamming{18} = 30
hamming{19} = 32
hamming{20} = 36
hamming{1690} = 2123366400
hamming{1000000} = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
A faster version can be built that also works on Tcl 8.5 (or earlier, if only small hamming numbers are being computed):
variable hamming 1 hi2 0 hi3 0 hi5 0
proc hamming {n} {
global hamming hi2 hi3 hi5
set h2 [expr {[lindex $hamming $hi2]*2}]
set h3 [expr {[lindex $hamming $hi3]*3}]
set h5 [expr {[lindex $hamming $hi5]*5}]
while {[llength $hamming] < $n} {
lappend hamming [set h [expr {
$h2<$h3
? $h2<$h5 ? $h2 : $h5
: $h3<$h5 ? $h3 : $h5
}]]
if {$h==$h2} {
set h2 [expr {[lindex $hamming [incr hi2]]*2}]
}
if {$h==$h3} {
set h3 [expr {[lindex $hamming [incr hi3]]*3}]
}
if {$h==$h5} {
set h5 [expr {[lindex $hamming [incr hi5]]*5}]
}
}
return [lindex $hamming [expr {$n - 1}]]
}
# Print the first 20 values of the sequence
for {set i 1} {$i <= 20} {incr i} {
puts [format "hamming\[%d\] = %d" $i [hamming $i]]
}
puts "hamming{1690} = [hamming 1690]"
puts "hamming{1691} = [hamming 1691]"
puts "hamming{1692} = [hamming 1692]"
puts "hamming{1693} = [hamming 1693]"
puts "hamming{1000000} = [hamming 1000000]"
[edit] Ursala
Smooth is defined as a second order function taking a list of primes and returning a function that takes a natural number n to the n-th smooth number with respect to them. An elegant but inefficient formulation based on the J solution is the following.
#import std
#import nat
smooth"p" "n" = ~&z take/"n" nleq-< (rep(length "n") ^Ts/~& product*K0/"p") <1>
This test program
main = smooth<2,3,5>* nrange(1,20)
yields this list of the first 20 Hamming numbers.
<1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36>
Although all calculations are performed using unlimited precision, the version above is impractical for large numbers. A more hardcore approach is the following.
#import std
#import nat
smooth"p" "n" =
~&H\"p" *-<1>; @NiXS ~&/(1,1); ~&ll~="n"->lr -+
^\~&rlPrrn2rrm2Zlrrmz3EZYrrm2lNCTrrm2QAX*rhlPNhrnmtPA2XtCD ~&lrPrhl2E?/~&l ^|/successor@l ~&hl,
^|/~& nleq-<&l+ * ^\~&r ~&l|| product@rnmhPX+-
#cast %nL
main = smooth<2,3,5>* nrange(1,20)--<1691,1000000>
- Output:
The great majority of time is spent calculating the millionth Hamming number.
< 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 2125764000, 519312780448388736089589843750000000000000000000000000000000000000000000000000000000>
Categories:
- Programming Tasks
- Solutions by Programming Task
- Ada
- AutoHotkey
- ALGOL 68
- C
- C sharp
- Clojure
- CoffeeScript
- Common Lisp
- D
- Factor
- Forth
- Fortran
- Go
- Haskell
- Icon
- Unicon
- J
- Java
- JavaScript
- Liberty BASIC
- Logo
- Lua
- MUMPS
- Oz
- PARI/GP
- Perl
- Perl 6
- PicoLisp
- Prolog
- Python
- Qi
- Qi examples needing attention
- Examples needing attention
- R
- REXX
- Ruby
- Scala
- Scheme
- Seed7
- Smalltalk
- Tcl
- Ursala
