Hamming numbers: Difference between revisions
m (→{{header|REXX}}: changed two comments.) |
m (→{{header|REXX}}: changed the version names to be more descriptive.) |
||
Line 4,856: | Line 4,856: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
=== |
===idiomatic=== |
||
This REXX program was a direct translation from my old REXX subroutine to compute '''UGLY''' numbers, |
This REXX program was a direct translation from my old REXX subroutine to compute '''UGLY''' numbers, |
||
<br>it computes just enough Hamming numbers (one Hamming number after the current number). |
<br>it computes just enough Hamming numbers (one Hamming number after the current number). |
||
Line 4,911: | Line 4,911: | ||
</pre> |
</pre> |
||
===unrolled |
===unrolled=== |
||
This REXX version is roughly twice as fast as the 1<sup>st</sup> REXX version. |
This REXX version is roughly twice as fast as the 1<sup>st</sup> REXX version. |
||
<lang rexx>/*REXX program computes Hamming numbers: 1 ──► 20, # 1691, the one millionth.*/ |
<lang rexx>/*REXX program computes Hamming numbers: 1 ──► 20, # 1691, the one millionth.*/ |
Revision as of 20:48, 29 October 2015
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 ).
- Show the one millionth Hamming number (if the language – or a convenient library – supports arbitrary-precision integers).
References
- Hamming numbers
- Smooth number
- Hamming problem from Dr. Dobb's CodeTalk (dead link as of Sep 2011; parts of the thread here and here).
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. <lang Ada>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;</lang>
- 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): <lang Ada> 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."/");</lang>
and then use it like this: <lang Ada> 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))));</lang>
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. <lang algol68>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)))</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
ATS
<lang ATS> // // How to compile: // patscc -DATS_MEMALLOC_LIBC -o hamming hamming.dats //
- include
"share/atspre_staload.hats"
fun min3 (
A: arrayref(int, 3)
) : natLt(3) = i where {
var x: int = A[0] var i: natLt(3) = 0 val () = if A[1] < x then (x := A[1]; i := 1) val () = if A[2] < x then (x := A[2]; i := 2)
} (* end of [min3] *)
fun hamming {n:pos} (
n: int(n)
) : int = let // var A = @[int](2, 3, 5) val A = $UNSAFE.cast{arrayref(int, 3)}(addr@A) var I = @[int](1, 1, 1) val I = $UNSAFE.cast{arrayref(int, 3)}(addr@I) val H = arrayref_make_elt<int> (i2sz(succ(n)), 0) val () = H[0] := 1 // fun loop{k:pos}
(k: int(k)) : void =
( // if k < n then let
val i = min3(A) val k = ( if A[i] > H[k-1] then (H[k] := A[i]; k+1) else k ) : intBtwe(k, k+1) val ii = I[i] val () = I[i] := ii+1 val ii = $UNSAFE.cast{natLte(n)}(ii) val () = if i = 0 then A[i] := 2*H[ii] val () = if i = 1 then A[i] := 3*H[ii] val () = if i = 2 then A[i] := 5*H[ii]
in
loop(k)
end // end of [then] else () // end of [else] // ) (* end of [loop] *) // in
loop (1); H[n-1]
end (* end of [hamming] *)
implement main0 () = { val () = loop(1) where { fun loop {n:pos} (
n: int(n)
) : void = if n <= 20 then let
val () = println! ("hamming(",n,") = ", hamming(n))
in
loop(n+1)
end // end of [then] // end of [if] } (* end of [val] *) val n = 1691 val () = println! ("hamming(",n,") = ", hamming(n)) // } (* end of [main0] *) </lang>
- 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(1691) = 2125764000
AutoHotkey
<lang 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 }</lang>
AWK
<lang AWK>
- syntax: GAWK -f HAMMING_NUMBERS.AWK
BEGIN {
for (i=1; i<=20; i++) { printf("%d ",hamming(i)) } printf("\n1691: %d\n",hamming(1691)) exit(0)
} function hamming(limit, h,i,j,k,n,x2,x3,x5) {
h[0] = 1 x2 = 2 x3 = 3 x5 = 5 for (n=1; n<=limit; n++) { h[n] = min(x2,min(x3,x5)) if (h[n] == x2) { x2 = 2 * h[++i] } if (h[n] == x3) { x3 = 3 * h[++j] } if (h[n] == x5) { x5 = 5 * h[++k] } } return(h[limit-1])
} function min(x,y) {
return((x < y) ? x : y)
} </lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 1691: 2125764000
BBC BASIC
<lang bbcbasic> @% = &1010
FOR h% = 1 TO 20 PRINT "H("; h% ") = "; FNhamming(h%) NEXT PRINT "H(1691) = "; FNhamming(1691) END DEF FNhamming(l%) LOCAL i%, j%, k%, n%, m, x2, x3, x5, h%() DIM h%(l%) : h%(0) = 1 x2 = 2 : x3 = 3 : x5 = 5 FOR n% = 1 TO l%-1 m = x2 IF m > x3 m = x3 IF m > x5 m = x5 h%(n%) = m IF m = x2 i% += 1 : x2 = 2 * h%(i%) IF m = x3 j% += 1 : x3 = 3 * h%(j%) IF m = x5 k% += 1 : x5 = 5 * h%(k%) NEXT = h%(l%-1)</lang>
- Output:
H(1) = 1 H(2) = 2 H(3) = 3 H(4) = 4 H(5) = 5 H(6) = 6 H(7) = 8 H(8) = 9 H(9) = 10 H(10) = 12 H(11) = 15 H(12) = 16 H(13) = 18 H(14) = 20 H(15) = 24 H(16) = 25 H(17) = 27 H(18) = 30 H(19) = 32 H(20) = 36 H(1691) = 2125764000
Bracmat
<lang bracmat>( ( hamming
= x2 x3 x5 n i j k min . tbl$(h,!arg) { This creates an array. Arrays are always global in Bracmat. } & 1:?(0$h) & 2:?x2 & 3:?x3 & 5:?x5 & 0:?n:?i:?j:?k & whl ' ( !n+1:<!arg:?n & !x2:?min & (!x3:<!min:?min|) & (!x5:<!min:?min|) & !min:?(!n$h) { !n is index into array h } & ( !x2:!min & 2*!((1+!i:?i)$h):?x2 | ) & ( !x3:!min & 3*!((1+!j:?j)$h):?x3 | ) & ( !x5:!min & 5*!((1+!k:?k)$h):?x5 | ) ) & !((!arg+-1)$h) (tbl$(h,0)&) { We delete the array by setting its size to 0 } )
& 0:?I & whl'(!I+1:~>20:?I&put$(hamming$!I " ")) & out$ & out$(hamming$1691) & out$(hamming$1000000) );</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
C
Using a min-heap to keep track of numbers. Does not handle big integers. <lang c>#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; }</lang>
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. <lang c>#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 (1 << 24)
- define MAX_POW 1024
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 { int e[NK]; double v; } ham_t, *ham;
ham_t *hams, values[NK] = {{{0}, 0}}; double inc[NK][MAX_POW];
/* 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 a, ham b) { int i; for (i = 0; i < NK && a->e[i] == b->e[i]; i++);
return i == NK; }
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];
for (ni = 0; ni < NK; ni++) { if (! _eq(values + ni, h)) continue; 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, j; hams = malloc(sizeof(ham_t) * MAX_HAM);
for (i = 0; i < NK; i++) { values[i].e[i] = 1; inc[i][1] = log(fac[i]); _setv(values + i);
for (j = 2; j < MAX_POW; j++) inc[i][j] = j * inc[i][1]; }
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; }</lang>
- Output:
1,691: 2125764000 1,000,000: 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 10,000,000: 16244105063830431823239 ..<a gadzillion digits>.. 000000000000000000000
C#
<lang csharp>using System; using System.Numerics; using System.Linq;
namespace Hamming {
class MainClass {
public static BigInteger Hamming(int n) { BigInteger two = 2, three = 3, five = 5; var h = new BigInteger[n]; h[0] = 1; BigInteger x2 = 2, x3 = 3, x5 = 5; int i = 0, j = 0, k = 0; for (int index = 1; index < n; index++) { h[index] = BigInteger.Min(x2, BigInteger.Min(x3, x5)); if (h[index] == x2) x2 = two * h[++i]; if (h[index] == x3) x3 = three * h[++j]; if (h[index] == x5) x5 = five * h[++k]; } return h[n - 1]; }
public static void Main(string[] args) { Console.WriteLine(string.Join(" ", Enumerable.Range(1, 20).ToList().Select(x => Hamming(x)))); Console.WriteLine(Hamming(1691)); Console.WriteLine(Hamming(1000000)); } }
}</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Generic version for any set of numbers
The algorithm is similar to the one above. <lang csharp>using System; using System.Numerics; using System.Linq;
namespace Hamming {
class MainClass {
public static BigInteger[] Hamming(int n, int[] a) { var primes = a.Select(x => (BigInteger)x).ToArray(); var values = a.Select(x => (BigInteger)x).ToArray(); var indexes = new int[a.Length]; var results = new BigInteger[n]; results[0] = 1; for (int iter = 1; iter < n; iter++) { results[iter] = values[0]; for (int p = 1; p < primes.Length; p++) if (results[iter] > values[p]) results[iter] = values[p]; for (int p = 0; p < primes.Length; p++) if (results[iter] == values[p]) values[p] = primes[p] * results[++indexes[p]]; } return results; } public static void Main(string[] args) { foreach (int[] primes in new int[][] { new int[] {2,3,5}, new int[] {2,3,5,7} }) { Console.WriteLine("{0}-Smooth:", primes.Last()); Console.WriteLine(string.Join(" ", Hamming(20, primes))); Console.WriteLine(Hamming(1691, primes).Last()); Console.WriteLine(Hamming(1000000, primes).Last()); Console.WriteLine(); } } }
}</lang>
- Output:
5-Smooth: 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 7-Smooth: 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 3317760 4157409948433216829957008507500000000
Fast version
Like some of the other implementations on this page, this version represents each number as a list of exponents which would be applied to each prime number. So the number 60 would be represented as int[3] { 2, 1, 1 } which is interpreted as 2^2 * 3^1 * 5^1.
As often happens, optimizing for speed caused a marked increase in code size and complexity. Clearly the versions I wrote above are easier to read & understand. They were also much quicker to write. But the generic version above runs in 3+ seconds for the 1000000th 5-smooth number whereas this version does it in 0.35 seconds, 8-10 times faster.
I've tried to comment it as best I could, without bloating the code too much.
--Mike Lorenz
<lang csharp>using System; using System.Linq; using System.Numerics;
namespace HammingFast {
class MainClass {
private static int[] _primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };
public static BigInteger Big(int[] exponents) { BigInteger val = 1; for (int i = 0; i < exponents.Length; i++) for (int e = 0; e < exponents[i]; e++) val = val * _primes[i]; return val; }
public static int[] Hamming(int n, int nprimes) { var hammings = new int[n, nprimes]; // array of hamming #s we generate var hammlogs = new double[n]; // log values for above var primelogs = new double[nprimes]; // pre-calculated prime log values var indexes = new int[nprimes]; // intermediate hamming values as indexes into hammings var listheads = new int[nprimes, nprimes]; // intermediate hamming list heads var listlogs = new double[nprimes]; // log values of list heads for (int p = 0; p < nprimes; p++) { listheads[p, p] = 1; // init list heads to prime values primelogs[p] = Math.Log(_primes[p]); // pre-calc prime log values listlogs[p] = Math.Log(_primes[p]); // init list head log values } for (int iter = 1; iter < n; iter++) { int min = 0; // find index of min item in list heads for (int p = 1; p < nprimes; p++) if (listlogs[p] < listlogs[min]) min = p; hammlogs[iter] = listlogs[min]; // that's the next hamming number for (int i = 0; i < nprimes; i++) hammings[iter, i] = listheads[min, i]; for (int p = 0; p < nprimes; p++) { // update each list head if it matches new value bool equal = true; // test each exponent to see if number matches for (int i = 0; i < nprimes; i++) { if (hammings[iter, i] != listheads[p, i]) { equal = false; break; } } if (equal) { // if it matches... int x = ++indexes[p]; // set index to next hamming number for (int i = 0; i < nprimes; i++) // copy each hamming exponent listheads[p, i] = hammings[x, i]; listheads[p, p] += 1; // increment exponent = mult by prime listlogs[p] = hammlogs[x] + primelogs[p]; // add log(prime) to log(value) = mult by prime } } }
var result = new int[nprimes]; for (int i = 0; i < nprimes; i++) result[i] = hammings[n - 1, i]; return result; }
public static void Main(string[] args) { foreach (int np in new int[] { 3, 4, 5 }) { Console.WriteLine("{0}-Smooth:", _primes[np - 1]); Console.WriteLine(string.Join(" ", Enumerable.Range(1, 20).Select(x => Big(Hamming(x, np))))); Console.WriteLine(Big(Hamming(1691, np))); Console.WriteLine(Big(Hamming(1000000, np))); Console.WriteLine(); } } }
}</lang>
- Output:
5-Smooth: 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 7-Smooth: 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 3317760 4157409948433216829957008507500000000 11-Smooth: 1 2 3 4 5 6 7 8 9 10 11 12 14 15 16 18 20 21 22 24 296352 561912530929780078125000
C# Enumerator Version
I wanted to fix the enumerator (old) version, as it wasn't working. It became a bit of an obsession... after a few iterations I came up with the following, which is the fastest C# version on my computer - your mileage may vary. It combines the speed of the Log method; Log(2)+Log(3)=Log(2*3) to help determine which is the next one to use. Then I have added some logic (using the series property) to ensure that exponent sets are never duplicated - which speeds the calculations up a bit.... Adding this trick to the Fast Version will probably result in the fastest version, but I'll leave that to someone else to implement. Finally it's all enumerated through a crazy one-way-linked-list-type-structure that only exists as long as the enumerator and is left up to the garbage collector to remove the bits no longer needed... I hope it's commented enough... follow it if you dare!
<lang csharp>using System; using System.Collections.Generic; using System.Linq; using System.Numerics;
namespace HammingTest {
class HammingNode { public double log; public int[] exponents; public HammingNode next; public int series; } class HammingListEnumerator : IEnumerable<BigInteger> { private int[] primes; private double[] primelogs; private HammingNode next; private HammingNode[] values; private HammingNode[] indexes;
public HammingListEnumerator(IEnumerable<int> seeds) { // Ensure our seeds are properly ordered, and generate their log values primes = seeds.OrderBy(x => x).ToArray(); primelogs = primes.Select(x => Math.Log10(x)).ToArray(); // Start at 1 (log(1)=0, exponents are all 0, series = none) next = new HammingNode { log = 0, exponents = new int[primes.Length], series = primes.Length }; // Set all exponent sequences to the start, and calculate the first value for each exponent indexes = new HammingNode[primes.Length]; values = new HammingNode[primes.Length]; for(int i = 0; i < primes.Length; ++i) { indexes[i] = next; values[i] = AddExponent(next, i); } }
// Make a copy of a node, and increment the specified exponent value private HammingNode AddExponent(HammingNode node, int i) { HammingNode ret = new HammingNode { log = node.log + primelogs[i], exponents = (int[])node.exponents.Clone(), series = i }; ++ret.exponents[i]; return ret; }
private void GetNext() { // Find which exponent value is the lowest int min = 0; for(int i = 1; i < values.Length; ++i) if(values[i].log < values[min].log) min = i; // Add it to the end of the 'list', and move to it next.next = values[min]; next = values[min];
// Find the next node in an allowed sequence (skip those that would be duplicates) HammingNode val = indexes[min].next; while(val.series < min) val = val.next;
// Keep the current index, and calculate the next value in the series for that exponent indexes[min] = val; values[min] = AddExponent(val, min); }
// Skip values without having to calculate the BigInteger value from the exponents public HammingListEnumerator Skip(int count) { for(int i = count; i > 0; --i) GetNext();
return this; }
// Calculate the BigInteger value from the exponents internal BigInteger ValueOf(HammingNode n) { BigInteger val = 1; for(int i = 0; i < n.exponents.Length; ++i) for(int e = 0; e < n.exponents[i]; e++) val = val * primes[i]; return val; }
public IEnumerator<BigInteger> GetEnumerator() { while(true) { yield return ValueOf(next); GetNext(); } }
System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator() { return this.GetEnumerator(); } }
class Program { static void Main(string[] args) { foreach(int[] primes in new int[][] { new int[] { 2, 3, 5 }, new int[] { 2, 3, 5, 7 }, new int[] { 2, 3, 5, 7, 9}}) { HammingListEnumerator hammings = new HammingListEnumerator(primes); System.Diagnostics.Debug.WriteLine("{0}-Smooth:", primes.Last()); System.Diagnostics.Debug.WriteLine(String.Join(" ", hammings.Take(20).ToArray())); System.Diagnostics.Debug.WriteLine(hammings.Skip(1691 - 20).First()); System.Diagnostics.Debug.WriteLine(hammings.Skip(1000000 - 1691).First()); System.Diagnostics.Debug.WriteLine(""); } } }
} </lang>
- Output:
5-Smooth: 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 7-Smooth: 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 3317760 4157409948433216829957008507500000000 11-Smooth: 1 2 3 4 5 6 7 8 9 10 11 12 14 15 16 18 20 21 22 24 296352 561912530929780078125000
C++
C++11 For Each Generator
<lang cpp>
- include <iostream>
- include <vector>
// Hamming like sequences Generator // // Nigel Galloway. August 13th., 2012 // class Ham { private: std::vector<unsigned int> _H, _hp, _hv, _x; public: bool operator!=(const Ham& other) const {return true;} Ham begin() const {return *this;}
Ham end() const {return *this;}
unsigned int operator*() const {return _x.back();} Ham(const std::vector<unsigned int> &pfs):_H(pfs),_hp(pfs.size(),0),_hv({pfs}),_x({1}){} const Ham& operator++() { for (int i=0; i<_H.size(); i++) for (;_hv[i]<=_x.back();_hv[i]=_x[++_hp[i]]*_H[i]); _x.push_back(_hv[0]); for (int i=1; i<_H.size(); i++) if (_hv[i]<_x.back()) _x.back()=_hv[i]; return *this; } }; </lang>
5-Smooth
<lang cpp> int main() {
int count = 1; for (unsigned int i : Ham({2,3,5})) { if (count <= 62) std::cout << i << ' '; if (count++ == 1691) { std::cout << "\nThe one thousand six hundred and ninety first Hamming Number is " << i << std::endl; break; } } return 0;
} </lang> Produces:
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 108 120 125 128 135 144 150 160 162 180 192 200 216 225 240 243 250 256 270 288 300 320 324 360 375 384 400 405 The one thousand six hundred and ninety first Hamming Number is 2125764000
7-Smooth
<lang cpp> int main() {
int count = 1; for (unsigned int i : Ham({2,3,5,7})) { std::cout << i << ' '; if (count++ == 64) break; } std::cout << std::endl; return 0;
} </lang> Produces:
1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 28 30 32 35 36 40 42 45 48 49 50 54 56 60 63 64 70 72 75 80 81 84 90 96 98 100 105 108 112 120 125 126 128 135 140 144 147 150 160 162 168 175 180 189
Clojure
This version implements Dijkstra's merge solution, so is closely related to the Haskell version. <lang clojure>(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*)))))
(def hamming
(lazy-seq (->> (map #(*' 5 %) hamming) (smerge (map #(*' 3 %) hamming)) (smerge (map #(*' 2 %) hamming)) (cons 1))))</lang>
Note that the above 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.
In order to fix the problems with the above program as to memory use and extra time expended, the following code is implemented as a function so that it does not retain the pointers to the streams used so that they can be garbage collected from the beginning as they are consumed; it also avoids duplicate number generation by using intermediate streams for each of the multiples and building each on the results of the last; finally, it orders the stream use from least dense to most so that the intermediate streams retained are as short as possible, with the "s5" stream only from one fifth to a third of the current value, the "s35" stream only between a third and a half of the current output value, and the s235 stream only between a half and the current output - as the sequence is not very dense with increasing range, mot many values need be retained: <lang clojure> (defn hamming
"Computes the unbounded sequence of Hamming 235 numbers." [] (letfn [(merge [xs ys] (let [xv (first xs), yv (first ys)] (if (< xv yv) (cons xv (lazy-seq (merge (next xs) ys))) (cons yv (lazy-seq (merge xs (next ys))))))), (smult [m s] ;; equiv to map (* m) s -- faster (cons (*' m (first s)) (lazy-seq (smult m (next s)))))] (do (def s5 (cons 5 (lazy-seq (smult 5 s5)))) (def s35 (cons 3 (lazy-seq (merge s5 (smult 3 s35))))) (def s235 (cons 2 (lazy-seq (merge s35 (smult 2 s235))))) (cons 1 (lazy-seq s235)))))
</lang>
Much of the time expended for larger ranges (say 10 million or more) is due to the time doing extended precision arithmetic, with also a significant percentage spent in garbage collection.
After compiling code in REPL:
- Output:
(take 20 (hamming)) (1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36) (->> (hamming) (drop 1690) (first) (time)) "Elapsed time: 1.105582 msecs" 2125764000 (->> (hamming) (drop 999999) (first) (time)) "Elapsed time: 447.561128 msecs" 519312780448388736089589843750000000000000000000000000000000000000000000000000000000N
CoffeeScript
<lang 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]</lang>
Common Lisp
Maintaining three queues, popping the smallest value every time. <lang lisp>(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))))</lang> A much faster method: <lang lisp>(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)))</lang>
- 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
D
Basic Version
This version keeps all numbers in memory, computing all the Hamming numbers up to the needed one. Performs constant number of operations per Hamming number produced. <lang d>import std.stdio, std.bigint, std.algorithm, std.range, core.memory;
auto hamming(in uint n) pure nothrow /*@safe*/ {
immutable BigInt two = 2, three = 3, five = 5; auto h = new BigInt[n]; h[0] = 1; BigInt x2 = 2, x3 = 3, x5 = 5; size_t i, j, k;
foreach (ref el; h.dropOne) { el = min(x2, x3, x5); if (el == x2) x2 = two * h[++i]; if (el == x3) x3 = three * h[++j]; if (el == x5) x5 = five * h[++k]; } return h.back;
}
void main() {
GC.disable; iota(1, 21).map!hamming.writeln; 1_691.hamming.writeln; 1_000_000.hamming.writeln;
}</lang>
- Output:
[1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Runtime is about 1.6 seconds with LDC2.
Alternative Version 1
This keeps numbers in memory, but over-computes a sequence by a factor of about , calculating extra multiples past that as well. Incurs an extra 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:
<lang d>import std.stdio, std.bigint, std.container, std.algorithm, std.range,
core.memory;
BigInt hamming(in int n) in {
assert(n > 0);
} body {
auto frontier = redBlackTree(2.BigInt, 3.BigInt, 5.BigInt); auto lowest = 1.BigInt; foreach (immutable _; 1 .. n) { lowest = frontier.front; frontier.removeFront; frontier.insert(lowest * 2); frontier.insert(lowest * 3); frontier.insert(lowest * 5); } return lowest;
}
void main() {
GC.disable; writeln("First 20 Hamming numbers: ", iota(1, 21).map!hamming); writeln("hamming(1691) = ", 1691.hamming); writeln("hamming(1_000_000) = ", 1_000_000.hamming);
}</lang>
- 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
About 3.2 seconds run time with LDC2.
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.
<lang d>import std.stdio: writefln; import std.bigint: BigInt; import std.conv: text; import std.numeric: gcd; import std.algorithm: copy, map; import std.array: array; import core.stdc.stdlib: calloc; 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] factors = [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.
public static __gshared immutable double[factors.length] inc = factors[].map!log.array;
bool opEquals(in ref Hamming y) const pure nothrow @nogc { //return this.e == y.e; // Too much slow. foreach (immutable i; 0 .. this.e.length) if (this.e[i] != y.e[i]) return false; return true; }
void update() pure nothrow @nogc { //this.v = dotProduct(inc, this.e); // Too much slow. this.v = 0.0; foreach (immutable i; 0 .. this.e.length) this.v += inc[i] * this.e[i]; }
string toString() const { BigInt result = 1; foreach (immutable i, immutable f; factors) result *= f.BigInt ^^ this.e[i]; return result.text; }
}
// Global variables. __gshared Hamming[] hams; __gshared Hamming[NK] values;
nothrow @nogc static this() {
// Slower than calloc if you don't use all the MAX_HAM items. //hams = new Hamming[MAX_HAM];
auto ptr = cast(Hamming*)calloc(MAX_HAM, Hamming.sizeof); static const err = new Error("Not enough memory."); if (!ptr) throw err; hams = ptr[0 .. MAX_HAM];
foreach (immutable i, ref v; values) { v.e[i] = 1; v.v = Hamming.inc[i]; }
}
ref Hamming getHam(in size_t n) nothrow @nogc
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 (immutable i; 1 .. NK) if (values[i].v < values[ni].v) ni = i;
hams[n_hams] = values[ni]; hams[n_hams].update; }
foreach (immutable 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 (immutable n; [1691, 10 ^^ 6, MAX_HAM]) writefln("%8d: %s", n, n.getHam);
}</lang> The output is similar to the second C version. Runtime is about 0.11 seconds if MAX_HAM = 1_000_000 (as the task requires), and 0.90 seconds if MAX_HAM = 10_000_000.
Alternative Version 3
This version is similar to the precedent, but frees unused values. It's a little slower than the precedent version, but it uses much less RAM, so it allows to compute the result for larger n. <lang d>import std.stdio: writefln; import std.bigint: BigInt; import std.conv: text; import std.algorithm: map; import std.array: array; import core.stdc.stdlib: malloc, calloc, free; import std.math: log; // ^^
// Number of factors. enum NK = 3;
__gshared immutable int[NK] primes = [2, 3, 5]; __gshared immutable double[NK] lnPrimes = primes[].map!log.array;
/// K-smooth numbers (stored as their exponents of each factor).
struct Hamming {
double ln; // Log of the number. ushort[NK] e; // Exponents of each factor. Hamming* next; size_t n;
// Recompute the logarithm from the exponents. void recalculate() pure nothrow @safe @nogc { this.ln = 0.0; foreach (immutable i, immutable ei; this.e) this.ln += lnPrimes[i] * ei; }
string toString() const { BigInt result = 1; foreach (immutable i, immutable f; primes) result *= f.BigInt ^^ this.e[i]; return result.text; }
}
Hamming getHam(in size_t n) nothrow @nogc in {
assert(n && n != size_t.max);
} body {
static struct Candidate { typeof(Hamming.ln) ln; typeof(Hamming.e) e;
void increment(in size_t n) pure nothrow @safe @nogc { e[n] += 1; ln += lnPrimes[n]; }
bool opEquals(T)(in ref T y) const pure nothrow @safe @nogc { // return this.e == y.e; // Slow. return !((this.e[0] ^ y.e[0]) | (this.e[1] ^ y.e[1]) | (this.e[2] ^ y.e[2])); }
int opCmp(T)(in ref T y) const pure nothrow @safe @nogc { return (ln > y.ln) ? 1 : (ln < y.ln ? -1 : 0); } }
static struct HammingIterator { // Not a Range. Candidate cand; Hamming* base; size_t primeIdx;
this(in size_t i, Hamming* b) pure nothrow @safe @nogc { primeIdx = i; base = b; cand.e = base.e; cand.ln = base.ln; cand.increment(primeIdx); }
void next() pure nothrow @safe @nogc { base = base.next; cand.e = base.e; cand.ln = base.ln; cand.increment(primeIdx); } }
HammingIterator[NK] its; Hamming* head = cast(Hamming*)calloc(Hamming.sizeof, 1); Hamming* freeList, cur = head; Candidate next;
foreach (immutable i, ref it; its) it = HammingIterator(i, cur);
for (size_t i = cur.n = 1; i < n; ) { auto leastReferenced = size_t.max; next.ln = double.max; foreach (ref it; its) { if (it.cand == *cur) it.next; if (it.base.n < leastReferenced) leastReferenced = it.base.n; if (it.cand < next) next = it.cand; }
// Collect unferenced numbers. while (head.n < leastReferenced) { auto tmp = head; head = head.next; tmp.next = freeList; freeList = tmp; }
if (!freeList) { cur.next = cast(Hamming*)malloc(Hamming.sizeof); } else { cur.next = freeList; freeList = freeList.next; }
cur = cur.next; version (fastmath) { cur.ln = next.ln; cur.e = next.e; } else { cur.e = next.e; cur.recalculate; // Prevent FP error accumulation. }
cur.n = i++; cur.next = null; }
auto result = *cur; version (leak) {} else { while (head) { auto tmp = head; head = head.next; tmp.free; }
while (freeList) { auto tmp = freeList; freeList = freeList.next; tmp.free; } }
return result;
}
void main() {
foreach (immutable n; [1691, 10 ^^ 6, 10_000_000]) writefln("%8d: %s", n, n.getHam);
}</lang> The output is the same as the second alternative version.
DCL
<lang DCL>$ limit = p1 $ $ n = 0 $ h_'n = 1 $ x2 = 2 $ x3 = 3 $ x5 = 5 $ i = 0 $ j = 0 $ k = 0 $ $ n = 1 $ loop: $ x = x2 $ if x3 .lt. x then $ x = x3 $ if x5 .lt. x then $ x = x5 $ h_'n = x $ if x2 .eq. h_'n $ then $ i = i + 1 $ x2 = 2 * h_'i $ endif $ if x3 .eq. h_'n $ then $ j = j + 1 $ x3 = 3 * h_'j $ endif $ if x5 .eq. h_'n $ then $ k = k + 1 $ x5 = 5 * h_'k $ endif $ n = n + 1 $ if n .le. limit then $ goto loop $ $ i = 0 $ loop2: $ write sys$output h_'i $ i = i + 1 $ if i .lt. 20 then $ goto loop2 $ $ n = limit - 1 $ write sys$output h_'n</lang>
- Output:
Here's the output; $ @hamming 1691 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000
Eiffel
<lang Eiffel> note description : "Initial part, in order, of the sequence of Hamming numbers" math : "[ Hamming numbers, also known as regular numbers and 5-smooth numbers, are natural integers that have 2, 3 and 5 as their only prime factors. ]" computer_arithmetic : "[ This version avoids integer overflow and stops at the last representable number in the sequence. ]" output : "[
Per requirements of the RosettaCode example, execution will produce items of indexes 1 to 20 and 1691. The algorithm (procedure `hamming') is more general and will produce the first `n' Hamming numbers for any `n'. ]"
source : "This problem was posed in Edsger W. Dijkstra, A Discipline of Programming, Prentice Hall, 1978" date : "8 August 2012" authors : "Bertrand Meyer", "Emmanuel Stapf" revision : "1.0" libraries : "Relies on SORTED_TWO_WAY_LIST from EiffelBase" implementation : "[ Using SORTED_TWO_WAY_LIST provides an elegant illustration of how to implement a lazy scheme in Eiffel through the use of object-oriented data structures. ]" warning : "[ The formatting (<lang>) specifications for Eiffel in RosettaCode are slightly obsolete: `note' and other newer keywords not supported, red color for manifest strings. This should be fixed soon. ]"
class APPLICATION
create make
feature {NONE} -- Initialization
make -- Print first 20 Hamming numbers, in order, and the 1691-st one. local Hammings: like hamming -- List of Hamming numbers, up to 1691-st one. do Hammings := hamming (1691) across 1 |..| 20 as i loop io.put_natural (Hammings.i_th (i.item)); io.put_string (" ") end io.put_new_line; io.put_natural (Hammings.i_th (1691)); io.put_new_line end
feature -- Basic operations
hamming (n: INTEGER): ARRAYED_LIST [NATURAL] -- First `n' elements (in order) of the Hamming sequence, -- or as many of them as will not produce overflow. local sl: SORTED_TWO_WAY_LIST [NATURAL] overflow: BOOLEAN first, next: NATURAL do create Result.make (n); create sl.make sl.extend (1); sl.start across 1 |..| n as i invariant -- "The numbers output so far are the first `i' - 1 Hamming numbers, in order". -- "Result.first is the `i'-th Hamming number." until sl.is_empty loop first := sl.first; sl.start Result.extend (first); sl.remove across << 2, 3, 5 >> as multiplier loop next := multiplier.item * first overflow := overflow or next <= first if not overflow and then not sl.has (next) then sl.extend (next) end end end end end </lang>
- Output:
1 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000
Elixir
<lang elixir>defmodule Hamming do
def generater do queues = [{2, queue}, {3, queue}, {5, queue}] Stream.unfold({1, queues}, fn {n, q} -> next(n, q) end) end defp next(n, queues) do queues = Enum.map(queues, fn {m, queue} -> {m, push(queue, m*n)} end) min = Enum.map(queues, fn {_, queue} -> top(queue) end) |> Enum.min queues = Enum.map(queues, fn {m, queue} -> {m, (if min==top(queue), do: erase_top(queue), else: queue)} end) {n, {min, queues}} end defp queue, do: {[], []} defp push({input, output}, term), do: {[term | input], output} defp top({input, []}), do: List.last(input) defp top({_, [h|_]}), do: h defp erase_top({input, []}), do: erase_top({[], Enum.reverse(input)}) defp erase_top({input, [_|t]}), do: {input, t}
end
IO.puts "first twenty Hamming numbers:" IO.inspect Hamming.generater |> Enum.take(20) IO.puts "1691st Hamming number:" IO.puts Hamming.generater |> Enum.take(1691) |> List.last IO.puts "one millionth Hamming number:" IO.puts Hamming.generater |> Enum.take(1_000_000) |> List.last</lang>
- Output:
first twenty Hamming numbers: [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] 1691st Hamming number: 2125764000 one millionth Hamming number: 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
ERRE
For bigger numbers, you have to use an external program, like MULPREC.R <lang ERRE>PROGRAM HAMMING
!$DOUBLE
DIM H[2000]
PROCEDURE HAMMING(L%->RES)
LOCAL I%,J%,K%,N%,M,X2,X3,X5 H[0]=1 X2=2 X3=3 X5=5 FOR N%=1 TO L%-1 DO M=X2 IF M>X3 THEN M=X3 END IF IF M>X5 THEN M=X5 END IF H[N%]=M IF M=X2 THEN I%+=1 X2=2*H[I%] END IF IF M=X3 THEN J%+=1 X3=3*H[J%] END IF IF M=X5 THEN K%+=1 X5=5*H[K%] END IF END FOR RES=H[L%-1]
END PROCEDURE
BEGIN
FOR H%=1 TO 20 DO HAMMING(H%->RES) PRINT("H(";H%;")=";RES) END FOR HAMMING(1691->RES) PRINT("H(1691)=";RES)
END PROGRAM</lang>
- Output:
H( 1 )= 1H( 2 )= 2 H( 3 )= 3 H( 4 )= 4 H( 5 )= 5 H( 6 )= 6 H( 7 )= 8 H( 8 )= 9 H( 9 )= 10 H( 10 )= 12 H( 11 )= 15 H( 12 )= 16 H( 13 )= 18 H( 14 )= 20 H( 15 )= 24 H( 16 )= 25 H( 17 )= 27 H( 18 )= 30 H( 19 )= 32 H( 20 )= 36 H(1691)= 2125764000
F#
This version implements Dijkstra's merge solution, so is closely related to the Haskell classic version. <lang fsharp> type LazyList<'a> = Cons of 'a * Lazy<LazyList<'a>>
let rec hamming =
let rec (-|-) (Cons(x, nxf) as xs) (Cons(y, nyf) as ys) = if x < y then Cons(x, lazy(nxf.Force() -|- ys)) elif x > y then Cons(y, lazy(xs -|- nyf.Force())) else Cons(x, lazy(nxf.Force() -|- nyf.Force())) let rec inf_map f (Cons(x, nxf)) = Cons(f x, lazy(inf_map f (nxf.Force()))) Cons(1I, lazy(let x = inf_map ((*) 2I) hamming let y = inf_map ((*) 3I) hamming let z = inf_map ((*) 5I) hamming x -|- y -|- z))
// testing... [<EntryPoint>] let main args =
let rec iterLazyListFor f n (Cons(v, rf)) = if n > 0 then f v; iterLazyListFor f (n - 1) (rf.Force()) let rec nthLazyList n ((Cons(v, rf)) as ll) = if n <= 1 then v else nthLazyList (n - 1) (rf.Force()) printf "( "; iterLazyListFor (printf "%A ") 20 (hamming); printfn ")" printfn "%A" (hamming |> nthLazyList 1691) printfn "%A" (hamming |> nthLazyList 1000000) 0
</lang>
The above code memory residency is quite high as it holds the entire lazy sequence in memory due to the reference preventing garbage collection as the sequence is consumed,
The following code reduces that high memory residency by making the routine a function and using internal local stream references for the intermediate streams so that they can be collected as the stream is consumed as long as no reference is held to the main results stream (which is not in the sample test functions); it also avoids duplication of factors by successively building up streams and further reduces memory use by ordering of the streams so that the least dense are determined first: <lang fsharp> let rec hamming() =
let rec (-|-) (Cons(x, nxf) as xs) (Cons(y, nyf) as ys) = if x < y then Cons(x, lazy(nxf.Force() -|- ys)) else Cons(y, lazy(xs -|- nyf.Force())) let rec smult m (Cons(x, nxf)) = // like map but faster Cons(m * x, lazy(smult m (nxf.Force()))) let rec s5 = Cons(5I, lazy(smult 5I s5)) let rec s35 = Cons(3I, lazy(s5 -|- smult 3I s35)) let rec s235 = Cons(2I, lazy(s35 -|- smult 2I s235)) Cons(1I, lazy s235)
</lang>
The above code can by used just by substituting it for the "hamming" binding and substituting "hamming()" for "hamming" in the main testing function calls (three places).
Both codes output the same results as follows but the second is over three times faster:
- Output:
( 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 ) 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Factor
<lang 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 ;</lang>
<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. <lang factor>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 ;</lang>
20 hamming ltake list>array . 1690 hamming lnth . 999999 hamming lnth .
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. <lang Forth>\ 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.</lang>
- 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
A smaller, less capable solution is presented here. It solves two out of three requirements and is ANS-Forth compliant. <lang Forth>2000 cells constant /hamming create hamming /hamming allot
( n1 n2 n3 n4 n5 n6 n7 -- n3 n4 n5 n6 n1 n2 n8)
- min? >r dup r> min >r 2rot r> ;
- hit? ( n1 n2 n3 n4 n5 n6 n7 n8 -- n3 n4 n9 n10 n1 n2 n7)
>r 2dup = \ compare number with found minimum if -rot drop 1+ hamming over cells + @ r@ * rot then r> drop >r 2rot r>
- \ if so, increment and rotate
- hamming# ( n1 -- n2)
1 hamming ! >r \ set first cell and initialize parms 0 5 over 3 over 2 r@ 1 ?do \ determine minimum and set cell dup min? min? min? dup hamming i cells + ! 2 hit? 5 hit? 3 hit? drop loop \ find if minimum equals value 2drop 2drop 2drop hamming r> 1- cells + @
- \ clean up stack and fetch hamming number
- test
cr 21 1 ?do i . i hamming# . cr loop 1691 hamming# . cr
- </lang>
Fortran
Using big_integer_module from here [1] <lang fortran>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</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
FunL
<lang funl>native scala.collection.mutable.Queue
val hamming =
q2 = Queue() q3 = Queue() q5 = Queue() def enqueue( n ) = q2.enqueue( n*2 ) q3.enqueue( n*3 ) q5.enqueue( n*5 )
def stream = val n = min( min(q2.head(), q3.head()), q5.head() ) if q2.head() == n then q2.dequeue() if q3.head() == n then q3.dequeue() if q5.head() == n then q5.dequeue() enqueue( n ) n # stream() for q <- [q2, q3, q5] do q.enqueue( 1 ) stream()</lang>
<lang funl>val hamming = 1 # merge( map((2*), hamming), merge(map((3*), hamming), map((5*), hamming)) )
def
merge( inx@x:_, iny@y:_ ) | x < y = x # merge( inx.tail(), iny ) | x > y = y # merge( inx, iny.tail() ) | otherwise = merge( inx, iny.tail() )
println( hamming.take(20) ) println( hamming(1690) ) println( hamming(2000) )</lang>
- Output:
[1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36] 2125764000 8100000000
Go
Concise version using dynamic-programming
<lang go>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])
}</lang>
- Output:
[1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36] 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Longer version using dynamic-programming and logarithms
More than 10 times faster. <lang go>package main
import ( "flag" "fmt" "log" "math" "math/big" "os" )
var ( // print the whole sequence or just one element?
seqMode = flag.Bool("s", false, "sequence mode") // precomputed base-2 logarithms for 3 and 5 lg3, lg5 float64 = math.Log2(3), math.Log2(5)
// state of the three multiplied sequences front = [3]cursor{ {0, 0, 1}, // 2 {1, 0, lg3}, // 3 {2, 0, lg5}, // 5 }
// table for dynamic-programming stored results table [][3]int16 )
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) val() [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.val() c.lg = float64(x[0]) + lg3*float64(x[1]) + lg5*float64(x[2]) }
func step() { table = append(table, front[0].val()) 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 main() { log.SetPrefix(os.Args[0] + ": ") log.SetOutput(os.Stderr) flag.Parse() if flag.NArg() != 1 { log.Fatalln("need one positive integer argument") } var ordinal int // ordinal of last sequence element to compute _, err := fmt.Sscan(flag.Arg(0), &ordinal) if err != nil || ordinal <= 0 { log.Fatalln("argument must be a positive integer") } table = make([][3]int16, 1, ordinal) for i, n := 1, ordinal; i < n; i++ { if *seqMode { show(table[i-1]) } step() } show(table[ordinal-1]) }</lang>
- 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
Haskell
The classic version
<lang haskell>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</lang>
Runs in about a second on Ideone.com.
The nested union
s' effect is to produce the minimal value at each step,
with duplicates removed. 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 . 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 ( curiously, it is in almost word-for-word correspondence with Edsger Dijkstra's code from his book A Discipline of Programming, p. 132 ). 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, – which is about in length – behind the scenes, as long as there's no re-use evident in the code.
Avoiding generation of duplicates
The classic version can be sped up quite a bit (above twice, and with better orders of growth, even closer to 1.0) by avoiding generation of duplicate values: <lang haskell>hamming = 1:foldl u [] [5,3,2] where u s n = ar where ar = merge s (n:map (n*) ar) merge [] b = b merge a@(x:xs) b@(y:ys) | x < y = x:merge xs b | otherwise = y:merge a ys
main = do print $ take 20 hamming print $ hamming !! 1690 print $ hamming !! (1000000-1)</lang>
Explicit multiples reinserting
This is a common approach which explicitly maintains an internal buffer of 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 operations for each number, instead of of the above version (and thus overall). Translation of Java (which does use priority queue though, so should have O (n logn) operations overall). Uses union
from the "classic" version above:
<lang Haskell>hamm n = drop n $ iterate (\(_,(a:t))-> (a,union t [2*a,3*a,5*a])) (0,[1])</lang>
- Output:
<lang Haskell>*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]</lang> Runs too slowly to reach 1,000,000, with empirical orders of growth above ~ (n 1.7 ) and worsening. Last two lines show the internal buffer's length for several sample n s.
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 empirical time and space complexity: <lang haskell>-- 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 = let (r,t) = nthHam 1000000 in print t >> print (trival t)
lg3 = logBase 2 3; lg5 = logBase 2 5 logval (i,j,k) = fromIntegral i + fromIntegral j*lg3 + fromIntegral k*lg5 trival (i,j,k) = 2^i * 3^j * 5^k estval n = (6*lg3*lg5* 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 :: Int -> (Double, (Int, Int, Int)) 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) | otherwise = res where (d,w) = rngval n -- correction dist, width hi = estval n - d -- hi > logval > hi-w (m,nb) = ( fromIntegral $ 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 -- 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 /lg5) ], let p = fromIntegral k*lg5, j <- [ 0 .. floor ((hi-p)/lg3) ], let q = fromIntegral j*lg3 + p, let (i,frac) = pr (hi-q) ; r = hi-frac ] -- r = i + q -- f 0 z == (sum $ map fst z, concat $ map snd z) where pr = properFraction f !c [] = (c,[]) -- code as a loop f !c ((c1,b1):r) = let (cr,br) = f (c+c1) r -- to prevent space leak in case b1 of { [v] -> (cr,v:br) ; _ -> (cr, br) }</lang>
- Output:
-- time: 0.01s memory: 3688 kB (55,47,64) 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
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. <lang Unicon># 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</lang>
J
Solution:
A concise tacit expression using a (right) fold:
<lang j>hamming=: {. (/:~@~.@] , 2 3 5 * {)/@(1x ,~ i.@-)</lang>
Example usage:
<lang j> 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</lang>
For the millionth (and billionth (1e9)) Hamming number see the nh
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
<lang j> /:~@~.@]</lang>
- produce (the next) 3 elements by selection and multiplication:
<lang j> 2 3 5 * {</lang> note that LHA (2 3 5 * {) RHA is equivalent to <lang j> 2 3 5 * LHA { RHA</lang>
- the RH part forms an array of descending indices and the initial Hamming number 1
<lang j> (1x ,~ i.@-)</lang> 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)
<lang j> ({. (/:~@~.@] , 2 3 5 * {)/@(1x ,~ i.@-)) 7
1 2 3 4 5 6 8</lang>
This starts using a descending sequence with 1 appended:
<lang j> (1x ,~ i.@-) 7
6 5 4 3 2 1 0 1</lang>
and then the fold expression is inserted between these list elements and the result computed:
<lang j> 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</lang>
(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.)
Java
Has a common shortcoming of overproducing the sequence by about 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) implementation for both insertion and deletion operations but it looks like it's in Java. <lang java5>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)); }
}</lang>
- 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
Another possibility is to realize that Hamming numbers can be represented and stored as triples of nonnegative integers which are the powers of 2, 3 and 5, and do all operations needed by the algorithms directly on these triples without converting to , which saves tremendous memory and time. Also, the search frontier through this three-dimensional grid can be generated in an order that never reaches the same state twice, so we don't need to keep track which states have already been encountered, saving even more memory. The objects of encode Hamming numbers in three fields , and . Multiplying by 2, 3 and 5 can now be done just by incrementing that field. The order comparison of triples needed by the priority queue is implemented with simple logarithm formulas of multiplication and addition, resorting to conversion to only in the rare cases that the floating point arithmetic produces too close results.
<lang java5> import java.math.BigInteger; import java.util.*;
// ilkka.kokkarinen@gmail.com
public class HammingTriple implements Comparable<HammingTriple> {
// Precompute a couple of constants that we need all the time private static final BigInteger two = BigInteger.valueOf(2); private static final BigInteger three = BigInteger.valueOf(3); private static final BigInteger five = BigInteger.valueOf(5); private static final double logOf2 = Math.log(2); private static final double logOf3 = Math.log(3); private static final double logOf5 = Math.log(5); // The powers of this triple private int a, b, c; public HammingTriple(int a, int b, int c) { this.a = a; this.b = b; this.c = c; } public String toString() { return "[" + a + ", " + b + ", " + c + "]"; } public BigInteger getValue() { return two.pow(a).multiply(three.pow(b)).multiply(five.pow(c)); } public boolean equals(Object other) { if(other instanceof HammingTriple) { HammingTriple h = (HammingTriple) other; return this.a == h.a && this.b == h.b && this.c == h.c; } else { return false; } } // Return 0 if this == other, +1 if this > other, and -1 if this < other public int compareTo(HammingTriple other) { // equality if(this.a == other.a && this.b == other.b && this.c == other.c) { return 0; } // this dominates if(this.a >= other.a && this.b >= other.b && this.c >= other.c) { return +1; } // other dominates if(this.a <= other.a && this.b <= other.b && this.c <= other.c) { return -1; } // take the logarithms for comparison double log1 = this.a * logOf2 + this.b * logOf3 + this.c * logOf5; double log2 = other.a * logOf2 + other.b * logOf3 + other.c * logOf5; // are these different enough to be reliable? if(Math.abs(log1 - log2) > 0.0000001) { return (log1 < log2) ? -1: +1; } // oh well, looks like we have to do this the hard way return this.getValue().compareTo(other.getValue()); // (getting this far should be pretty rare, though) } public static BigInteger computeHamming(int n, boolean verbose) { if(verbose) { System.out.println("Hamming number #" + n); } long startTime = System.currentTimeMillis(); // The elements of the search frontier PriorityQueue<HammingTriple> frontierQ = new PriorityQueue<HammingTriple>(); int maxFrontierSize = 1; // Initialize the frontier frontierQ.offer(new HammingTriple(0, 0, 0)); // 1 while(true) { if(frontierQ.size() > maxFrontierSize) { maxFrontierSize = frontierQ.size(); } // Pop out the next Hamming number from the frontier HammingTriple curr = frontierQ.poll(); if(--n == 0) { if(verbose) { System.out.println("Time: " + (System.currentTimeMillis() - startTime) + " ms"); System.out.println("Frontier max size: " + maxFrontierSize); System.out.println("As powers: " + curr.toString()); System.out.println("As value: " + curr.getValue()); } return curr.getValue(); } // Current times five, if at origin in (a,b) plane if(curr.a == 0 && curr.b == 0) { frontierQ.offer(new HammingTriple(curr.a, curr.b, curr.c + 1)); } // Current times three, if at line a == 0 if(curr.a == 0) { frontierQ.offer(new HammingTriple(curr.a, curr.b + 1, curr.c)); } // Current times two, unconditionally curr.a++; frontierQ.offer(curr); // reuse the current HammingTriple object } }
} </lang>
Hamming number #1000000 Time: 650 ms Frontier max size: 10777 As powers: [55, 47, 64] As value: 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 Hamming number #1000000000 Time: 1763306 ms Frontier max size: 1070167 As powers: [1334, 335, 404] As value: 62160757555652448616308163328720720039470565190896527065916324096423370220027531418244175407 772567327803701726166152919355404186200255249167295000868314547113136940786355040041603128729517887 0364794838245609107270160079056207179759030665476588225699039176388785014115448224991592743918456282 8227449023750262318234797192076792208033475638322151983772515798004125909334741121595323950448656375 1044570269974247729669174417794061727369755885568000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000
Alternative
This uses memoized streams - similar in principle to the classic lazy-evaluated sequence, but with a java flavor. Hope you like this recipe!
<lang java> import java.math.BigInteger;
public class Hamming {
public static void main(String args[]) { Stream hamming = makeHamming(); System.out.print("H[1..20] "); for (int i=0; i<20; i++) { System.out.print(hamming.value()); System.out.print(" "); hamming = hamming.advance(); } System.out.println();
System.out.print("H[1691] "); hamming = makeHamming(); for (int i=1; i<1691; i++) { hamming = hamming.advance(); } System.out.println(hamming.value());
hamming = makeHamming(); System.out.print("H[10^6] "); for (int i=1; i<1000000; i++) { hamming = hamming.advance(); } System.out.println(hamming.value()); }
public interface Stream { BigInteger value(); Stream advance(); }
public static class MultStream implements Stream { MultStream(int mult) { m_mult = BigInteger.valueOf(mult); } MultStream setBase(Stream s) { m_base = s; return this; } public BigInteger value() { return m_mult.multiply(m_base.value()); } public Stream advance() { return setBase(m_base.advance()); }
private final BigInteger m_mult; private Stream m_base; }
private final static class RegularStream implements Stream { RegularStream(Stream[] streams, BigInteger val) { m_streams = streams; m_val = val; } public BigInteger value() { return m_val; }
public Stream advance() { // memoized value for the next stream instance. if (m_advance != null) { return m_advance; }
int minidx = 0 ; BigInteger next = nextStreamValue(0); for (int i=1; i<m_streams.length; i++) { BigInteger v = nextStreamValue(i); if (v.compareTo(next) < 0) { next = v; minidx = i; } } RegularStream ret = new RegularStream(m_streams, next); // memoize the value! m_advance = ret; m_streams[minidx].advance(); return ret; } private BigInteger nextStreamValue(int streamidx) { // skip past duplicates in the streams we're merging. BigInteger ret = m_streams[streamidx].value(); while (ret.equals(m_val)) { m_streams[streamidx] = m_streams[streamidx].advance(); ret = m_streams[streamidx].value(); } return ret; } private final Stream[] m_streams; private final BigInteger m_val; private RegularStream m_advance = null; }
private final static Stream makeHamming() { MultStream nums[] = new MultStream[] { new MultStream(2), new MultStream(3), new MultStream(5) }; Stream ret = new RegularStream(nums, BigInteger.ONE); for (int i=0; i<nums.length; i++) { nums[i].setBase(ret); } return ret; }
} </lang>
$ java Hamming H[1..20] 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 H[1691] 2125764000 H[10^6] 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 $
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.
<lang javascript>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());</lang>
- Output:
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36 ... 1691 => 2125764000
Fast & complete version
A translation of my fast C# version. I was curious to see how much slower JavaScript is. The result: it runs about 5x times slower than C#, though YMMV. You can try it yourself here: http://jsfiddle.net/N7AFN/
--Mike Lorenz
<lang javascript><html> <head></head> <body>
</body> <script src="http://code.jquery.com/jquery-latest.min.js"></script> <script src="http://peterolson.github.com/BigInteger.js/BigInteger.min.js"></script> <script type="text/javascript">
var _primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
function log(text) { $('#main').append(text + "\n"); }
function big(exponents) { var i, e, val = bigInt.one; for (i = 0; i < exponents.length; i++) for (e = 0; e < exponents[i]; e++) val = val.times(_primes[i]); return val.toString(); }
function hamming(n, nprimes) { var i, iter, p, q, min, equal, x;
var hammings = new Array(n); // array of hamming #s we generate hammings[0] = new Array(nprimes); for (p = 0; p < nprimes; p++) { hammings[0][p] = 0; }
var hammlogs = new Array(n); // log values for above hammlogs[0] = 0;
var primelogs = new Array(nprimes); // pre-calculated prime log values var listlogs = new Array(nprimes); // log values of list heads for (p = 0; p < nprimes; p++) { primelogs[p] = listlogs[p] = Math.log(_primes[p]); }
var indexes = new Array(nprimes); // intermediate hamming values as indexes into hammings for (p = 0; p < nprimes; p++) { indexes[p] = 0; }
var listheads = new Array(nprimes); // intermediate hamming list heads for (p = 0; p < nprimes; p++) { listheads[p] = new Array(nprimes); for (q = 0; q < nprimes; q++) { listheads[p][q] = 0; } listheads[p][p] = 1; }
for (iter = 1; iter < n; iter++) { min = 0; for (p = 1; p < nprimes; p++) if (listlogs[p] < listlogs[min]) min = p; hammlogs[iter] = listlogs[min]; // that's the next hamming number hammings[iter] = listheads[min].slice(); for (p = 0; p < nprimes; p++) { // update each list head if it matches new value equal = true; // test each exponent to see if number matches for (i = 0; i < nprimes; i++) { if (hammings[iter][i] != listheads[p][i]) { equal = false; break; } } if (equal) { // if it matches... x = ++indexes[p]; // set index to next hamming number listheads[p] = hammings[x].slice(); // copy hamming number listheads[p][p] += 1; // increment exponent = mult by prime listlogs[p] = hammlogs[x] + primelogs[p]; // add log(prime) to log(value) = mult by prime } } }
return hammings[n - 1]; }
$(document).ready(function() { var i, nprimes; var t = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,1691,1000000];
for (nprimes = 3; nprimes <= 4; nprimes++) { var start = new Date();
log('
' + _primes[nprimes - 1] + '-Smooth:' + '
'); log('
'); for (i = 0; i < t.length; i++) log('' + '');var end = new Date();log('' + ''); log('
' + t[i] + ':' + ' | ' + big(hamming(t[i], nprimes)) + ' |
' + 'Elapsed time:' + ' | ' + (end-start)/1000 + ' seconds' + ' |
');
} });
</script> </html></lang>
- Output:
5-Smooth: 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 time: 1.73 seconds 7-Smooth: 1: 1 2: 2 3: 3 4: 4 5: 5 6: 6 7: 7 8: 8 9: 9 10: 10 11: 12 12: 14 13: 15 14: 16 15: 18 16: 20 17: 21 18: 24 19: 25 20: 27 1691: 3317760 1000000: 4157409948433216829957008507500000000 Elapsed time: 1.989 seconds
jq
We take the primary challenge here to be to write a Hamming number generator that can generate a given number of Hamming numbers, or the n-th Hamming number, without storing previously generated numbers.
To motivate a more complex version, in Part 1 of this section hamming(n) is defined as a generator of Hamming numbers, as numbers. This function uses an efficient algorithm and can run indefinitely, but it has one disadvantage: currently, jq converts large integers to floating point approximations, and thus precision is lost. For example, it reports the millionth Hamming number as 1.926511252902403e+44.
In Part 2, the algorithm in the first part is modified to use the [p,q,r] representation of Hamming numbers, where p, q, and r are the relevant exponents respectively of 2, 3, and 5.
The task description focuses on finding the n-th element of an infinite sequence and so it should be mentioned that using jq versions greater than 1.4, it would be possible to simply the generator so that is always unbounded, and then harness it with new builtins such as "limit" and "nth".
Hamming number generator
<lang jq># Return the index in the input array of the min_by(f) value def index_min_by(f):
. as $in | if length == 0 then null else .[0] as $first | reduce range(0; length) as $i ([0, $first, ($first|f)]; # state: [ix; min; f|min] ($in[$i]|f) as $v | if $v < .[2] then [ $i, $in[$i], $v ] else . end) | .[0] end;
- Emit n Hamming numbers if n>0; the nth if n<0
def hamming(n):
# input: [twos, threes, fives] of which at least one is assumed to be non-empty # output: the index of the array holding the min of the firsts def next: map( .[0] ) | index_min_by(.);
# input: [value, [twos, threes, fives] ....] # ix is the index in [twos, threes, fives] of the array to be popped # output: [popped, updated_arrays ...] def pop(ix): .[1] as $triple | setpath([0]; $triple[ix][0]) | setpath([1,ix]; $triple[ix][1:]);
# input: [x, [twos, threes, fives], count] # push value*2 to twos, value*3 to threes, value*5 to fives and increment count def push(v): [.[0], [.[1][0] + [2*v], .[1][1] + [3*v], .[1][2] + [5*v]], .[2] + 1];
# _hamming is the workhorse # input: [previous, [twos, threes, fives], count] def _hamming: .[0] as $previous | if (n > 0 and .[2] == n) or (n<0 and .[2] == -n) then $previous else (.[1]|next) as $ix # $ix cannot be null | pop($ix) | .[0] as $next | (if $next == $previous then empty elif n>=0 then $previous else empty end), (if $next == $previous then . else push($next) end | _hamming) end; [1, [[2],[3],[5]], 1] | _hamming;
. as $n | hamming($n)</lang> Examples: <lang jq># First twenty: hamming(20)
- See elsewhere for output
- 1691st Hamming number:
hamming(-1691)
- => 2125764000
- Millionth:
hamming(-1000000)
- => 1.926511252902403e+44
</lang>
Hamming numbers as triples
In this section, Hamming numbers are represented as triples, [p,q,r], where p, q and r are the relevant powers of 2, 3, and 5 respectively. We therefore begin with some functions for managing Hamming numbers represented in this manner: <lang jq># The log (base e) of a Hamming triple: def ln_hamming:
if length != 3 then error("ln_hamming: \(.)") else . end | (.[0] * (2|log)) + (.[1] * (3|log)) + (.[2] * (5|log));
- The numeric value of a Hamming triple:
def hamming_tof: ln_hamming | exp;
def hamming_toi:
def pow(n): . as $in | reduce range(0;n) as $i (1; . * $in); . as $in | (2|pow($in[0])) * (3|pow($in[1])) * (5|pow($in[2]));
- Return the index in the input array of the min_by(f) value
def index_min_by(f):
. as $in | if length == 0 then null else .[0] as $first | reduce range(0; length) as $i ([0, $first, ($first|f)]; # state: [ix; min; f|min] ($in[$i]|f) as $v | if $v < .[2] then [ $i, $in[$i], $v ] else . end) | .[0] end;
- Emit n Hamming numbers (as triples) if n>0; the nth if n<0; otherwise indefinitely.
def hamming(n):
# n must be 2, 3 or 5 def hamming_times(n): n as $n | if $n==2 then .[0] += 1 elif $n==3 then .[1] += 1 else .[2] += 1 end;
# input: [twos, threes, fives] of which at least one is assumed to be non-empty # output: the index of the array holding the min of the firsts def next: map( .[0] ) | index_min_by( ln_hamming );
# input: [value, [twos, threes, fives] ....] # ix is the index in [twos, threes, fives] of the array to be popped # output: [popped, updated_arrays ...] def pop(ix): .[1] as $triple | setpath([0]; $triple[ix][0]) | setpath([1,ix]; $triple[ix][1:]);
# input: [x, [twos, threes, fives], count] # push value*2 to twos, value*3 to threes, value*5 to fives and increment count def push(v): [.[0], [.[1][0] + [v|hamming_times(2)], .[1][1] + [v|hamming_times(3)], .[1][2] + [v|hamming_times(5)]], .[2] + 1];
# _hamming is the workhorse # input: [previous, [twos, threes, fives], count] def _hamming: .[0] as $previous | if (n > 0 and .[2] == n) or (n<0 and .[2] == -n) then $previous else (.[1]|next) as $ix # $ix cannot be null | pop($ix) | .[0] as $next | (if $next == $previous then empty elif n>=0 then $previous else empty end), (if $next == $previous then . else push($next) end | _hamming) end; [[0,0,0], [ 1,0,0 ,0,1,0, 0,0,1 ], 1] | _hamming;
</lang> Examples <lang jq># The first twenty Hamming numbers as integers: hamming(-20) | hamming_toi
- => (see elsewhere)
- 1691st as a Hamming triple:
hamming(-1691)
- => [5,12,3]
- The millionth:
hamming(-1000000)
- => [55,47,64]</lang>
Julia
The n parameter was chosen by trial and error. You have to pick an n large enough that the powers of 2, 3 and 5 will all be greater than n at the 1691st Hamming number.
<lang julia>n = 40
powers_2 = 2.^[0:n-1] powers_3 = 3.^[0:n-1] powers_5 = 5.^[0:n-1]
matrix = powers_2 * powers_3' powers_23 = sort(reshape(matrix,length(matrix),1),1)
matrix = powers_23 * powers_5' powers_235 = sort(reshape(matrix,length(matrix),1),1)
- Remove the integer overflow values.
powers_235 = powers_235[powers_235 .> 0]
println(powers_235[1:20]) println(powers_235[1691])</lang>
Kotlin
<lang kotlin>import java.math.BigInteger import java.util.PriorityQueue
val Three = BigInteger.valueOf(3) val Five = BigInteger.valueOf(5)
fun updateFrontier(x : BigInteger, pq : PriorityQueue<BigInteger>) {
pq add(x shiftLeft(1)) pq add(x multiply(Three)) pq add(x multiply(Five))
}
fun hamming(n : Int) : BigInteger {
val frontier = PriorityQueue<BigInteger>() updateFrontier(BigInteger.ONE, frontier) var lowest = BigInteger.ONE for (i in 1 .. n-1) { lowest = frontier.poll() ?: lowest while (frontier.peek() equals(lowest)) frontier.poll() updateFrontier(lowest, frontier) } return lowest
}
fun main(args : Array<String>) {
System.out print("Hamming(1 .. 20) =") for (i in 1 .. 20) System.out print(" ${hamming(i)}") System.out println("\nHamming(1691) = ${hamming(1691)}") System.out println("Hamming(1000000) = ${hamming(1000000)}")
}</lang>
- 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
Overloaded method:
<lang scala>import java.math.BigInteger import java.util.PriorityQueue
val One = BigInteger.ONE val Three = BigInteger.valueOf(3) val Five = BigInteger.valueOf(5)
fun PriorityQueue<BigInteger>.update(x: BigInteger) {
add(x shiftLeft(1)) add(x multiply(Three)) add(x multiply(Five))
}
fun hamming(n: Int): BigInteger {
val frontier = PriorityQueue<BigInteger>() frontier.update(One) var lowest = One repeat(n - 1) { lowest = frontier.poll() ?: lowest while (frontier.peek() == lowest) frontier.poll() frontier.update(lowest) } return lowest
}
fun hamming(i : Iterable<Int>) : Iterable<BigInteger> = i.map { hamming(it) }
fun main(args: Array<String>) {
val r = 1..20 println("Hamming($r) = " + hamming(r)) arrayOf(1691, 1000000).forEach { println("Hamming(${it}) = " + hamming(it)) }
}</lang>
Recursive method:
<lang scala>import java.math.BigInteger import java.util.PriorityQueue
val One = BigInteger.ONE val Three = BigInteger.valueOf(3) val Five = BigInteger.valueOf(5)
fun PriorityQueue<BigInteger>.update(x: BigInteger) {
add(x shiftLeft 1) add(x multiply Three) add(x multiply Five)
}
fun hamming(a: Any?): Any = when (a) {
is Number -> { val pq = PriorityQueue<BigInteger>() pq update One var lowest = One repeat(a.toInt() - 1) { lowest = pq.poll() ?: lowest while (pq.peek() == lowest) pq.poll() pq update lowest } lowest } is Iterable<*> -> a.map { hamming(it) } else -> throw IllegalArgumentException("cannot parse argument")
}
fun main(args: Array<String>) {
arrayOf(1..20, 1691, 1000000).forEach { println("Hamming($it) = " + hamming(it)) }
}</lang>
- 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
Liberty BASIC
LB has unlimited precision integers. <lang lb> dim h( 1000000)
for i =1 to 20
print hamming( i); " ";
next i
print 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</lang>
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
Logo
<lang 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</lang>
Lua
<lang 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)</lang>
Mathematica
<lang mathematica>HammingList[N_] := Module[{A, B, C}, {A, B, C} = (N^(1/3))*{2.8054745679851933, 1.7700573778298891, 1.2082521307023026} - {1, 1, 1};
Take[ Sort@Flatten@Table[ 2^x * 3^y * 5^z ,
{x, 0, A}, {y, 0, (-B/A)*x + B}, {z, 0, C - (C/A)*x - (C/B)*y}], N]];</lang>
HammingList[20] -> {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36} HammingList[1691] // Last -> 2125764000 HammingList[1000000] // Last ->519312780448388736089589843750000000000000000000000000000000000000000000000000000000
MATLAB / Octave
The n parameter was chosen by trial and error. You have to pick an n large enough that the powers of 2, 3 and 5 will all be greater than n at the 1691st Hamming number.
<lang Matlab>n = 40;
powers_2 = 2.^[0:n-1]; powers_3 = 3.^[0:n-1]; powers_5 = 5.^[0:n-1];
matrix = powers_2' * powers_3; powers_23 = sort(reshape(matrix,n*n,1));
matrix = powers_23 * powers_5;
powers_235 = sort(reshape(matrix,n*n*n,1));
% % Remove the integer overflow values. % powers_235 = powers_235(powers_235 > 0);
disp(powers_235(1:20)) disp(powers_235(1691))</lang>
MUMPS
<lang 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</lang>
Nim
<lang nim>import bigints, math
proc hamming(limit): BigInt =
var h = newSeq[BigInt](limit) x2 = initBigInt(2) x3 = initBigInt(3) x5 = initBigInt(5) i, j, k = 0 for i in 0..h.high: h[i] = initBigInt(1)
for n in 1 .. < limit: h[n] = min(x2, x3, x5) if x2 == h[n]: inc i x2 = h[i] * 2 if x3 == h[n]: inc j x3 = h[j] * 3 if x5 == h[n]: inc k x5 = h[k] * 5
result = h[h.high]
for i in 1 .. 20:
echo hamming(i)
echo hamming(1691) echo hamming(1_000_000)</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Oz
Lazy Version
<lang 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}}</lang>
Strict Version
The strict version uses iterators and a priority queue. Note that it can calculate other variations of the hamming numbers too. By changing K, it will calculate the p(K)-smooth numbers. (E.g. K = 3, it will use the first three primes 2,3 and 5, thus resulting in the 5-smooth numbers, see [2]) <lang oz> functor import Application System define
class Multiplier attr lst factor current
meth init(Factor Lst) lst := Lst factor := Factor {self next} end meth next local A AS in A|AS = @lst current := A*@factor lst := AS end end meth peek(?X) X = @current end
meth dump {System.showInfo "DUMP"} {System.showInfo "Factor: "#@factor} {System.showInfo "current: "#@current} end end
% a priority queue of multipliers. The one which currently holds the smallest value is put on front class PriorityQueue attr mults current % for duplicate detection
meth init(Mults) mults := Mults current := 0 end
meth insert(Mult) local fun {Insert M Lst} local Av Mv in case Lst of nil then M|Lst [] A|AS then {A peek(Av)} {M peek(Mv)} if Av < Mv then A|{Insert M AS} else M|A|AS end end end end in mults := {Insert Mult @mults} end end
meth next(Tail NextTail) local M Ms X Curr in M|Ms = @mults {M peek(X)} % gets value of lowest iterator Curr = @current if Curr == X then skip else Tail = X|NextTail % if we found a new value: append end {M next} mults := Ms {self insert(M)} if Curr == X then {self next(Tail NextTail)} else current := X end end end
end
local
% Sieve of erasthothenes, adapted from http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Oz fun {Sieve N} S = {Array.new 2 N true} M = {Float.toInt {Sqrt {Int.toFloat N}}} in for I in 2..M do if S.I then for J in I*I..N;I do S.J := false end end end S end
fun {Primes N} S = {Sieve N} in for I in 2..N collect:C do if S.I then {C I} end end end
% help method to extract args
proc {GetNK ArgList N K}
case ArgList of
A|B|_ then
N={StringToInt A}
K={StringToInt B}
end
end
proc {Generate N PriorQ Tail}
local
NewTail
in
if N == 0 then
Tail = nil
else
{PriorQ next(Tail NewTail)}
{Generate (N-1) PriorQ NewTail}
end
end
end
K = 3 PrimeFactors Lst Tail in ArgList = {Application.getArgs plain} Lst = 1|Tail PrimeFactors = {List.take {Primes K*K} K} Mults = {List.map PrimeFactors fun {$ A} {New Multiplier init(A Lst) } end} PriorQ = {New PriorityQueue init(Mults)} {Generate 20 PriorQ Tail} {ForAll Lst System.showInfo} {Application.exit 0} end end </lang> Strict version made by pietervdvn; do what you want with the code.
PARI/GP
This is a basic implementation; finding the millionth term requires 3 seconds and 54 MB. Much better algorithms exist. <lang parigp>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)</lang>
- 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
Pascal
Simple brute force til 2^32-1.I was astonished by the speed.The inner loop is taken 2^32 -1 times.DIV by constant is optimized to Mul and shift. Using FPC_64 3.1.1, i4330 3.5 Ghz <lang pascal> program HammNumb; {$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} { type
NativeUInt = longWord;
} var
pot : array[0..2] of NativeUInt;
function NextHammNumb(n:NativeUInt):NativeUInt; var
q,p,nr : NativeUInt;
begin
repeat nr := n+1; n := nr;
p := 0; while NOT(ODD(nr)) do begin inc(p); nr := nr div 2; end; Pot[0]:= p;
p := 0; q := nr div 3; while q*3=nr do Begin inc(P); nr := q; q := nr div 3; end; Pot[1] := p;
p := 0; q := nr div 5; while q*5=nr do Begin inc(P); nr := q; q := nr div 5; end; Pot[2] := p;
until nr = 1; result:= n;
end;
procedure Check; var
i,n: NativeUint;
begin
n := 1; for i := 1 to 20 do begin n := NextHammNumb(n); write(n,' '); end; writeln; writeln; n := 1; for i := 1 to 1690 do n := NextHammNumb(n); writeln('No ',i:4,' | ',n,' = 2^',Pot[0],' 3^',Pot[1],' 5^',Pot[2]);
end;
Begin
Check;
End.</lang> Output
2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 40 No 1690 | 2125764000 = 2^5 3^12 5^3 real 0m17.328s user 0m17.310s
a fast alternative
The Pascal code above is by far slower.Easily to use for smooth-3 .. smooth-37.
Big(O) is nearly linear to sub-linear . 1E7-> 0.028s => x10 =>1e8 ->0.273s => x1000 => 100'200'300'400 ~ 1e11 35.907s // estimated 270 s! This depends extreme on sorting speed.
http://rosettacode.org/wiki/Hamming_numbers#Direct_calculation_through_triples_enumeration is head to head, but still faster for very big numbers >1e8 (10^8: 4 MB 0.27 sec)
100'200'300'400 calculates in 8.33 s
For fpc 3.1.1_64 linux on 3.5 Ghz i4330, depends on 64-Bit by a factor of 4 slower on 32-Bit
/* For 12 primes "smooth-37" 1e8 it takes 02.807 s */
I collect only the factors between p^n and p^(n+1), in a recursive way in different lists
5 is a list consisting only 5^? = 1 factor
3 is a sorted list 3^?..3^?+1 and inserted values of 5
2 is a sorted list 2^?..2^?+1 and inserted values of list 3
Changing sizeOf(tElem) to 32 {maxPrimFakCnt = 3+8} instead of 16 ( x2) {maxPrimFakCnt = 3} results in increasing the runtime by x4 ( 2^2 )
<lang pascal>program hammNumb; {$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ASMCSE,CSE,PEEPHOLE} {$ALIGN 16}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} uses
sysutils;
const
maxPrimFakCnt = 3;//3 or 3+8 if tNumber= double, else -1 for extended to keep data aligned minElemCnt = 10;
type
tPrimList = array of NativeUint; tnumber = double; tpNumber= ^tnumber; tElem = record n : tnumber;//ln(prime[0]^Pots[0]*... Pots: array[0..maxPrimFakCnt] of word; end; tpElem = ^tElem; tElems = array of tElem; tElemArr = array [0..0] of tElem; tpElemArr = ^tElemArr;
tpFaktorRec = ^tFaktorRec; tFaktorRec = record frElems : tElems; frInsElems: tElems; frAktIdx : NativeUint; frMaxIdx : NativeUint; frPotNo : NativeUint; frActPot : NativeUint; frNextFr : tpFaktorRec; frActNumb: tElem; frLnPrime: tnumber; end; tArrFR = array of tFaktorRec;
var
Pl : tPrimList; ActIndex : NativeUint; ArrInsert : tElems;
procedure PlInit(n: integer); const
cPl : array[0..11] of byte=(2,3,5,7,11,13,17,19,23,29,31,37);
var
i : integer;
Begin
IF n>High(cPl)+1 then n := High(cPl) else IF n < 0 then n := 1; setlength(Pl,n); dec(n); For i := 0 to n do Pl[i] := cPl[i];
end;
procedure AusgabeElem(pElem: tElem); var
i : integer;
Begin
with pElem do Begin IF n < 23 then write(round(exp(n)):16) else write('ln ',n:13:7); For i := 0 to maxPrimFakCnt do write(' ',PL[i]:2,'^',Pots[i]); end; writeln
end;
//LoE == List of Elements function LoEGetNextNumber(pFR :tpFaktorRec):tElem;forward;
procedure LoECreate(const Pl: tPrimList;var FA:tArrFR); var
i : integer;
Begin
setlength(ArrInsert,100); setlength(FA,Length(PL)); For i := 0 to High(PL) do with FA[i] do Begin //automatic zeroing IF i < High(PL) then Begin setlength(frElems,minElemCnt); setlength(frInsElems,minElemCnt); frNextFr := @FA[i+1] end else Begin setlength(frElems,2); setlength(frInsElems,0); frNextFr := NIL; end; frPotNo := i; frLnPrime:= ln(PL[i]); frMaxIdx := 0; frAktIdx := 0; frActPot := 1; With frElems[0] do Begin n := frLnPrime; Pots[i]:= 1; end; frActNumb := frElems[0]; end;
end;
procedure LoEFree(var FA:tArrFR);
var
i : integer;
Begin
For i := High(FA) downto Low(FA) do setlength(FA[i].frElems,0); setLength(FA,0);
end;
function LoEGetActElem(pFr:tpFaktorRec):tElem; Begin
with pFr^ do result := frElems[frAktIdx];
end;
function LoEGetActLstNumber(pFr:tpFaktorRec):tpNumber; Begin
with pFr^ do result := @frElems[frAktIdx].n;
end;
procedure LoEIncInsArr(var a:tElems); Begin
setlength(a,Length(a)*8 div 5);
end;
procedure LoEIncreaseElems(pFr:tpFaktorRec;minCnt:NativeUint); var
newLen: NativeUint;
Begin
with pFR^ do begin newLen := Length(frElems); minCnt := minCnt+frMaxIdx; repeat newLen := newLen*8 div 5 +1; until newLen > minCnt; setlength(frElems,newLen); end;
end;
procedure LoEInsertNext(pFr:tpFaktorRec;Limit:tnumber); var
pNum : tpNumber; pElems : tpElemArr; cnt,i,u : NativeInt;
begin
with pFr^ do Begin //collect numbers of heigher primes cnt := 0; pNum := LoEGetActLstNumber(frNextFr); while Limit > pNum^ do Begin frInsElems[cnt] := LoEGetNextNumber(frNextFr);
// writeln( 'Ins ',frInsElems[cnt].n:10:8,' < ',pNum^:10:8);
inc(cnt); IF cnt > High(frInsElems) then LoEIncInsArr(frInsElems); pNum := LoEGetActLstNumber(frNextFr); end;
if cnt = 0 then EXIT;
i := frMaxIdx; u := frMaxIdx+cnt+1;
IF u > High(frElems) then LoEIncreaseElems(pFr,cnt);
IF frPotNo = 0 then inc(ActIndex,u); //Merge pElems := @frElems[0]; dec(cnt); dec(u); frMaxIdx:= u; repeat
// writeln(i:10,cnt:10,u:10); writeln( pElems^[i].n:10:8,' < ',frInsElems[cnt].n:10:8);
IF pElems^[i].n < frInsElems[cnt].n then Begin pElems^[u] := frInsElems[cnt]; dec(cnt); end else Begin pElems^[u] := pElems^[i]; dec(i); end; dec(u); until (i<0) or (cnt<0); IF i < 0 then For u := cnt downto 0 do pElems^[u] := frInsElems[u];
end;
end;
procedure LoEAppendNext(pFr:tpFaktorRec;Limit:tnumber); var
pNum : tpNumber; pElems : tpElemArr; i : NativeInt;
begin
with pFr^ do Begin i := frMaxIdx+1; pElems := @frElems[0]; pNum := LoEGetActLstNumber(frNextFr); while Limit > pNum^ do Begin IF i > High(frElems) then Begin LoEIncreaseElems(pFr,10); pElems := @frElems[0]; end; pElems^[i] := LoEGetNextNumber(frNextFr); inc(i); pNum := LoEGetActLstNumber(frNextFr); end; inc(ActIndex,i); frMaxIdx:= i-1; end;
end;
procedure LoENextList(pFr:tpFaktorRec); var
pElems : tpElemArr; j : NativeUint;
begin
with pFR^ do Begin //increase Elements by factor pElems := @frElems[0]; for j := frMaxIdx Downto 0 do with pElems^[j] do Begin n := n+frLnPrime; inc(Pots[frPotNo]); end; //x^j -> x^(j+1) j := frActPot+1; with frActNumb do begin n:= j*frLnPrime; Pots[frPotNo]:= j; end; frActPot := j; //if something follows IF frNextFr <> NIL then LoEInsertNext(pFR,frActNumb.n); frAktIdx := 0; end;
end;
function LoEGetNextNumber(pFR :tpFaktorRec):tElem; Begin
with pFr^ do Begin result := frElems[frAktIdx]; inc(frAktIdx); IF frMaxIdx < frAktIdx then LoENextList(pFr); end;
end;
procedure LoEGetNumber(pFR :tpFaktorRec;no:NativeUint); Begin
dec(no); while ActIndex < no do LoENextList(pFR); with pFr^ do frAktIdx := (no-(ActIndex-frMaxIdx)-1);
end;
var
T1,T0: tDateTime; FA: tArrFR; i : integer;
Begin
PlInit(3);// 3 -> 2,3,5 LoECreate(Pl,FA); i := 1; i := 1; T0 := time;
For i := 1 to 20 do AusgabeElem(LoEGetNextNumber(@FA[0]));
LoEGetNumber(@FA[0],1691); AusgabeElem(LoEGetNextNumber(@FA[0]));
LoEGetNumber(@FA[0],1000*1000); AusgabeElem(LoEGetNextNumber(@FA[0])); LoEGetNumber(@FA[0],100*1000*1000); T1 := time; AusgabeElem(LoEGetNextNumber(@FA[0])); Writeln('Timed 100*1000*1000 in ',FormatDateTime('HH:NN:SS.ZZZ',T1-T0));
Writeln('Actual Index ',ActIndex ); AusgabeElem(LoEGetNextNumber(@FA[0])); For i := 0 to High(FA) do writeln(pL[i]:2, ' elemcount ',FA[i].frMaxIdx+1:7,' out of',length(FA[i].frElems):7); LoEFree(FA);
End.</lang> Output
2,3,4,5....,36,40 ... shortened 2125764000 2^5 3^12 5^3 0^0 ln 192.7618989 2^55 3^47 5^64 0^0 ln 900.9063136 2^2 3^454 5^249 0^0 Timed 100*1000*1000 in 00:00:00.276 Actual Index 100061507 ln 900.9063159 2^142 3^80 5^444 0^0 2 elemcount 230300 out of 348159 3 elemcount 561 out of 772 5 elemcount 1 out of 2 real 0m0.278s user 0m0.273s sys 0m0.003s 2125764000 2^5 3^12 5^3 0^0 ln 192.7618989 2^55 3^47 5^64 0^0 ln 417.2530468 2^80 3^92 5^162 0^0 00:00:00.028 Actual Index 10201068944--> hamming Nr: 100200300400 see http://ideone.com/q3fma ln 4215.6152353 2^942 3^2276 5^660 0^0 2 elemcount 5028911 out of 5841156 3 elemcount 2620 out of 3165 5 elemcount 1 out of 2 real 0m35.963s user 0m35.907s sys 0m0.023s ... change zu use 12 primes [2..37] ( 32 bit ) -> 2.2x runtime over using 3 primes Begin PlInit(12) ln 40.8834947 2^14 3^0 5^6 7^4 11^2 13^1 17^0 19^1 23^0 29^0 31^1 37^0 Actual Index 100269652 Timed 100000000 in 00:00:02.807 2 elemcount 14322779 out of 14953361 3 elemcount 3387290 out of 3650722 5 elemcount 891236 out of 891289 7 elemcount 289599 out of 348159 11 elemcount 92240 out of 135999 13 elemcount 28272 out of 33202 17 elemcount 9394 out of 12969 19 elemcount 2639 out of 3165 23 elemcount 676 out of 772 29 elemcount 119 out of 188 31 elemcount 15 out of 17 37 elemcount 1 out of 2
Perl
<lang perl>use List::Util 'min';
- If you want the large output, uncomment either the one line
- marked (1) or the two lines marked (2)
- use Math::GMP qw/:constant/; # (1) uncomment this to use Math::GMP
- use Math::GMPz; # (2) uncomment this plus later line for Math::GMPz
sub ham_gen { my @s = ([1], [1], [1]); my @m = (2, 3, 5); #@m = map { Math::GMPz->new($_) } @m; # (2) uncomment for Math::GMPz
return sub { 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";
- You will need to pick one of the bigint choices
- ++$i, $h->() until $i == 999999;
- print ++$i, "-th: ", $h->(), "\n";
</lang>
- 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 1000000-th: 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
The core module bigint (Math::BigInt) is very slow, even with the GMP backend. There are some common alternatives. Math::GMP is handy and takes about 15 seconds. Math::GMPz takes slightly more work but finishes in about 5 seconds.
Perl 6
The limit scaling is not required, but it cuts down on a bunch of unnecessary calculation. <lang perl6>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</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000
PicoLisp
<lang 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))</lang>
- Output:
(1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36) 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
PL/I
<lang PL/I>(subscriptrange): Hamming: procedure options (main); /* 14 November 2013 */
declare (H(3000), t) fixed (15); declare (i, j, k, m, n) fixed binary; declare swaps bit (1);
on underflow ;
m = 0; n = 12; do k = 0 to n; do j = 0 to n; do i = 0 to n; m = m + 1; H(m) = 2**i * 3**j * 5**k; end; end; end; /* sort */ swaps = '1'b; do while (swaps); /* Cocktail-shaker sort is adequate, because values are largely sorted */ swaps = '0'b; do i = 1 to m-1, i-1 to 1 by -1; if H(i) > H(i+1) then /* swap */ do; t = H(i); H(i) = H(i+1); H(i+1) = t; swaps = '1'b; end; end; end; do i = 1 to m; put skip data (H(i)); end; put skip data (H(1653));
end Hamming;</lang> Results:
H(1)= 1; H(2)= 2; H(3)= 3; H(4)= 4; H(5)= 5; H(6)= 6; H(7)= 8; H(8)= 9; H(9)= 10; H(10)= 12; H(11)= 15; H(12)= 16; H(13)= 18; H(14)= 20; H(15)= 24; H(16)= 25; H(17)= 27; H(18)= 30; H(19)= 32; H(20)= 36;
Prolog
Generator idiom
<lang Prolog>%% 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.</lang>
SWI Prolog 6.2.6 produces (in about 7 ideone seconds):
?- time( main(1000000) ). [1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36] [2125764000,2147483648] 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 % 10,017,142 inferences
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).
<lang Prolog>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))).</lang>
- 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 .
Python
Version based on example from Dr. Dobb's CodeTalk
<lang python>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</lang>
- 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(), 1690, 1691)) [2125764000] >>> list(islice(hamming2(), 999999, 1000000)) [519312780448388736089589843750000000000000000000000000000000000000000000000000000000]
Another implementation of same approach
This version uses a lot of memory, it doesn't try to limit memory usage. <lang python>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)</lang>
"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).
<lang python>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()</lang> Results are the same as before.
Non-sharing recursive generator
Another formulation along the same lines, but greatly simplified, found here. Lacks data sharing, i.e. calls self recursively thus creating a separate copy of the data stream fed to the tee() call, again and again, instead of using its own output. This gravely impacts the efficiency. Not to be used.
<lang python>from heapq import merge from itertools import tee
def hamming_numbers():
last = 1 yield last
a,b,c = tee(hamming_numbers(), 3)
for n in merge((2*i for i in a), (3*i for i in b), (5*i for i in c)): if n != last: yield n last = n</lang>
Cyclic generator method #2.
Cyclic generator method #2. Considerably faster due to early elimination (before merge) of duplicates. Currently the faster Python version. Direct copy of Haskell code. <lang python>from itertools import islice, chain, tee
def merge(r, s):
# This is faster than heapq.merge. rr = r.next() ss = s.next() while True: if rr < ss: yield rr rr = r.next() else: yield ss ss = s.next()
def p(n):
def gen(): x = n while True: yield x x *= n return gen()
def pp(n, s):
def gen(): for x in (merge(s, chain([n], (n * y for y in fb)))): yield x r, fb = tee(gen()) return r
def hamming(a, b = None):
if not b: b = a + 1 seq = (chain([1], pp(5, pp(3, p(2))))) return list(islice(seq, a - 1, b - 1))
print hamming(1, 21) print hamming(1691)[0] print hamming(1000000)[0]</lang>
Qi
<lang 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)</lang>
- Output:
[1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36]
R
Recursively find the Hamming numbers below . Shown are results for tasks 1 and 2. Arbitrary precision integers are not supported natively. <lang R>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
} h <- sort(hamming(1,limit=2^31-1)) print(h[1:20]) print(h[length(h)])</lang>
- Output:
[1] 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 [1] 2125764000
Racket
<lang racket>
- lang racket
(require racket/stream) (define first stream-first) (define rest stream-rest)
(define (merge s1 s2)
(define x1 (first s1)) (define x2 (first s2)) (cond [(= x1 x2) (merge s1 (rest s2))] [(< x1 x2) (stream-cons x1 (merge (rest s1) s2))] [else (stream-cons x2 (merge s1 (rest s2)))]))
(define (mult k) (λ(x) (* x k)))
(define hamming
(stream-cons 1 (merge (stream-map (mult 2) hamming) (merge (stream-map (mult 3) hamming) (stream-map (mult 5) hamming)))))
(for/list ([i 20] [x hamming]) x) (stream-ref hamming 1690) (stream-ref hamming 999999) </lang>
- Output:
'(1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36) 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Raven
<lang raven>define hamming use $limit
[ ] as $h 1 $h 0 set 2 as $x2 3 as $x3 5 as $x5 0 as $i 0 as $j 0 as $k 1 $limit 1 + 1 range each as $n $x3 $x5 min $x2 min $h $n set $h $n get $x2 = if $i 1 + as $i $h $i get 2 * as $x2 $h $n get $x3 = if $j 1 + as $j $h $j get 3 * as $x3 $h $n get $x5 = if $k 1 + as $k $h $k get 5 * as $x5 $h $limit 1 - get
1 21 1 range each as $lim
$lim hamming print " " print
"\n" print
"Hamming(1691) is: " print 1691 hamming print "\n" print
- Raven can't handle > 2^31 using integers
- "Hamming(1000000) is: " print 1000000 hamming print "\n" print</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 Hamming(1691) is: 2125764000
REXX
idiomatic
This REXX program was a direct translation from my old REXX subroutine to compute UGLY numbers,
it computes just enough Hamming numbers (one Hamming number after the current number).
<lang rexx>/*REXX program computes Hamming numbers: 1 ──► 20, # 1691, the one millionth.*/
numeric digits 100 /*ensure enough decimal digits. */
call hamming 1, 20 /*show the 1st ──► twentieth Hamming #s*/
call hamming 1691 /*show the 1,691st Hamming number. */
call hamming 1000000 /*show the 1 millionth Hamming number.*/
exit /*stick a fork in it, we're all done. */
/*────────────────────────────────────────────────────────────────────────────*/
hamming: procedure; parse arg x,y; if y== then y=x; w=length(y)
#2=1; #3=1; #5=1; @.=0; @.1=1 do n=2 for y-1 @.n = min(2*@.#2, 3*@.#3, 5*@.#5) /*pick the minimum of 3 Hamming numbers*/ if 2*@.#2 == @.n then #2 = #2+1 /*number already defined? Use next #. */ if 3*@.#3 == @.n then #3 = #3+1 /* " " " " " " */ if 5*@.#5 == @.n then #5 = #5+1 /* " " " " " " */ end /*n*/ /* [↑] assign next 3 Hamming numbers. */ do j=x to y /*W is used to align the (output) index*/ say 'Hamming('right(j,w)") =" @.j /*display 'em, Dano.*/ end /*j*/
say right( 'length of last Hamming number =' length(@.y), 70); say return</lang> output when using the default input(s):
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
unrolled
This REXX version is roughly twice as fast as the 1st REXX version. <lang rexx>/*REXX program computes Hamming numbers: 1 ──► 20, # 1691, the one millionth.*/ numeric digits 100 /*ensure enough decimal digits. */ call hamming 1, 20 /*show the 1st ──► twentieth Hamming #s*/ call hamming 1691 /*show the 1,691st Hamming number. */ call hamming 1000000 /*show the 1 millionth Hamming number.*/ exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ hamming: procedure; parse arg x,y; if y== then y=x; w=length(y)
#2=1; #3=1; #5=1; @.=0; @.1=1 do n=2 for y-1 _2 = @.#2 + @.#2 /*this is faster than: 2 * @.#2 */ _3 = 3 * @.#3 _5 = 5 * @.#5 m =_2 /*assume a minimum of the 3 Hamming #s.*/ if _3 < m then m =_3 /*is this number less than the minimum?*/ if _5 < m then m =_5 /* " " " " " " " */ @.n = m /*now, assign the next Hamming number. */ if _2 == m then #2 = #2 + 1 /*number already defined? Use next #. */ if _3 == m then #3 = #3 + 1 /* " " " " " " */ if _5 == m then #5 = #5 + 1 /* " " " " " " */ end /*n*/ /* [↑] assign next 3 Hamming numbers. */ do j=x to y /*W is used to align the (output) index*/ say 'Hamming('right(j,w)") =" @.j /*display 'em, Dano.*/ end /*j*/
say right( 'length of last Hamming number =' length(@.y), 70); say
return</lang>
output is identical to the 1st REXX version.
Ruby
<lang ruby>hamming = Enumerator.new do |yielder|
next_ham = 1 queues = [[ 2, []], [3, []], [5, []] ] loop do yielder << next_ham # or: yielder.yield(next_ham) queues.each {|m,queue| queue << next_ham * m} next_ham = queues.collect{|m,queue| queue.first}.min queues.each {|m,queue| queue.shift if queue.first==next_ham} end
end</lang> And the "main" part of the task <lang ruby>start = Time.now
hamming.each.with_index(1) do |ham, idx|
case idx when (1..20), 1691 puts "#{idx} => #{ham}" when 1_000_000 puts "#{idx} => #{ham}" break end
end
puts "elapsed: #{Time.now - start} seconds"</lang>
- 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: 6.522811 seconds
Run BASIC
<lang runbasic> dim h(1000000) for i =1 to 20
print hamming(i);" ";
next i
print print "Hamming List First(1691) =";chr$(9);hamming(1691) print "Hamming List Last(1000000) =";chr$(9);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</lang>
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 Hamming List First(1691) = 2125764000 Hamming List Last(1000000) = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Scala
<lang 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)
}</lang> 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: <lang scala>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)
}</lang> 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
<lang scala>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)))
}</lang> 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
Translation of Haskell code avoiding duplicates
One can fix the problems of the memory use of the above code resulting from the entire stream being held in memory due to the use a "val hamming: Stream[BigInt]" by using a function "def hamming(): Stream[BigInt]" and making temporary local variables for intermediate streams so that the beginnings of the streams are garbage collected as the output stream is consumed; one can also implement the other Haskell algorithm to avoid factor duplication by building each stream on successive streams, again with memory conserved by building the least dense first: <lang scala>
def hamming(): Stream[BigInt] = { def merge(a: Stream[BigInt], b: Stream[BigInt]): Stream[BigInt] = { val av = a.head; val bv = b.head if (av < bv) av #:: merge(a.tail, b) else bv #:: merge(a, b.tail) } def smult(m:BigInt, s: Stream[BigInt]): Stream[BigInt] = (m * s.head) #:: smult(m, s.tail) // equiv to map (m *) s - faster lazy val s5: Stream[BigInt] = 5 #:: smult(5, s5) lazy val s35: Stream[BigInt] = 3 #:: merge(s5, smult(3, s35)) lazy val s235: Stream[BigInt] = 2 #:: merge(s35, smult(2, s235)) 1 #:: s235 }
</lang> Usage:
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) scala> hamming() drop 1690 head res1: BigInt = 2125764000 scala> hamming() drop 999999 head res2: BigInt = 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
It only takes under a half second to find the millionth number in the sequence in the last output.
Scheme
<lang 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)</lang>
- 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
Avoiding Generation of Duplicates, including reduced memory use
Although the algorithm above is true to the classic Dijkstra version and although the algorithm does require a form of lazy list/stream processing in order to utilize memoization and avoid repeated recalculations/comparisons, the stream implementation can be simplified, and the modified algorithm as per the Haskell code avoids duplicate generations of factors. As well, the following code implements the algorithm as a procedure/function so that it restarts the calculation from the beginning on every new call and so that internal stream variables are not top level so that the garbage collector can collect the beginning of all intermediate and final streams when they are no longer referenced; in this way total memory used (after interspersed garbage collections) is almost zero for a sequence of the first million numbers: <lang scheme>(define (hamming)
(define (merge a b) (let ((x (car a)) (y (car b))) (if (< x y) (cons x (delay (merge (force (cdr a)) b))) (cons y (delay (merge a (force (cdr b)))))))) (define (smult m s) (cons (* m (car s)) (delay (smult m (force (cdr s)))))) ;; equiv to map (* m) s (define s5 (cons 5 (delay (smult 5 s5)))) (define s35 (cons 3 (delay (merge s5 (smult 3 s35))))) (define s235 (cons 2 (delay (merge s35 (smult 2 s235))))) (cons 1 (delay s235)))
- test...
(define (stream-take->list n strm)
(if (= n 0) (list) (cons (car strm) (stream-take->list (- n 1) (force (cdr strm))))))
(define (stream-ref strm nth)
(do ((nxt strm (force (cdr nxt))) (cnt 0 (+ cnt 1))) ((>= cnt nth) (car nxt))))
(display (stream-take->list 20 (hamming))) (newline) (display (stream-ref (hamming) (- 1691 1))) (newline) (display (stream-ref (hamming) (- 1000000 1))) (newline)</lang>
- Output:
{1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36} 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
The "stream-ref" procedure is zero based as is the Scheme standard for array indices, thus the subtraction of one from the desired nth number in the sequence.
Seed7
<lang 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;</lang>
- Output:
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 2125764000 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
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. <lang smalltalk>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.</lang>
<lang smalltalk> limit := 10 raisedToInteger: 84. tape := Set new.
hammingProcess := [:newHamming| (newHamming <= limit) ifTrue: [| index | index := tape scanFor: newHamming. (tape array at: index) ifNil: [tape atNewIndex: index put: newHamming asSetElement. hammingProcess value: newHamming * 2. hammingProcess value: newHamming * 3. hammingProcess value: newHamming * 5]]].
hammingProcess value: 1.
sc := tape asSortedCollection. sc first: 20. "a SortedCollection(1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36)" sc at: 1691. "2125764000" sc at: 1000000. "519312780448388736089589843750000000000000000000000000000000000000000000000000000000" </lang>
Tcl
This uses coroutines to simplify the description of what's going on.
<lang tcl>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"</lang>
- 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): <lang tcl>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]"</lang>
uBasic/4tH
uBasic's single array does not have the required size to calculate the 1691st number, let alone the millionth. <lang>For H = 1 To 20
Print "H("; H; ") = "; Func (_FnHamming(H))
Next
End
_FnHamming Param (1)
@(0) = 1
X = 2 : Y = 3 : Z = 5 I = 0 : J = 0 : K = 0
For N = 1 To a@ - 1 M = X If M > Y Then M = Y If M > Z Then M = Z @(N) = M
If M = X Then I = I + 1 : X = 2 * @(I) If M = Y Then J = J + 1 : Y = 3 * @(J) If M = Z Then K = K + 1 : Z = 5 * @(K) Next
Return (@(a@-1))</lang>
- Output:
H(1) = 1 H(2) = 2 H(3) = 3 H(4) = 4 H(5) = 5 H(6) = 6 H(7) = 8 H(8) = 9 H(9) = 10 H(10) = 12 H(11) = 15 H(12) = 16 H(13) = 18 H(14) = 20 H(15) = 24 H(16) = 25 H(17) = 27 H(18) = 30 H(19) = 32 H(20) = 36 0 OK, 0:379
UNIX Shell
Large numbers are not supported. <lang bash>typeset -a hamming=(1) function nextHamming {
typeset -Sa q2 q3 q5 integer h=${hamming[${#hamming[@]}-1]} q2+=( $(( h*2 )) ) q3+=( $(( h*3 )) ) q5+=( $(( h*5 )) ) h=$( min3 ${q2[0]} ${q3[0]} ${q5[0]} ) (( ${q2[0]} == h )) && ashift q2 >/dev/null (( ${q3[0]} == h )) && ashift q3 >/dev/null (( ${q5[0]} == h )) && ashift q5 >/dev/null hamming+=($h)
}
function ashift {
nameref ary=$1 print -- "${ary[0]}" ary=( "${ary[@]:1}" )
}
function min3 {
if (( $1 < $2 )); then (( $1 < $3 )) && print -- $1 || print -- $3 else (( $2 < $3 )) && print -- $2 || print -- $3 fi
}
for ((i=1; i<=20; i++)); do
nextHamming printf "%d\t%d\n" $i ${hamming[i-1]}
done for ((; i<=1690; i++)); do nextHamming; done nextHamming printf "%d\t%d\n" $i ${hamming[i-1]} print "elapsed: $SECONDS"</lang>
- 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 1690 2125764000 elapsed: 0.568
Ursala
Smooth is defined as a second order function taking a list of primes and returning a function that takes a natural number to the -th smooth number with respect to them. An elegant but inefficient formulation based on the J solution is the following. <lang Ursala>#import std
- import nat
smooth"p" "n" = ~&z take/"n" nleq-< (rep(length "n") ^Ts/~& product*K0/"p") <1></lang> This test program <lang Ursala>main = smooth<2,3,5>* nrange(1,20)</lang> 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. <lang Ursala>#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></lang>
- 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>
VBScript
<lang vb> For h = 1 To 20 WScript.StdOut.Write "H(" & h & ") = " & Hamming(h) WScript.StdOut.WriteLine Next WScript.StdOut.Write "H(" & 1691 & ") = " & Hamming(1691) WScript.StdOut.WriteLine
Function Hamming(l) Dim h() : Redim h(l) : h(0) = 1 i = 0 : j = 0 : k = 0 x2 = 2 : x3 = 3 : x5 = 5 For n = 1 To l-1 m = x2 If m > x3 Then m = x3 End If If m > x5 Then m = x5 End If h(n) = m If m = x2 Then i = i + 1 : x2 = 2 * h(i) End If If m = x3 Then j = j + 1 : x3 = 3 * h(j) End If If m = x5 Then k = k + 1 : x5 = 5 * h(k) End If Next Hamming = h(l-1) End Function </lang>
- Output:
H(1) = 1 H(2) = 2 H(3) = 3 H(4) = 4 H(5) = 5 H(6) = 6 H(7) = 8 H(8) = 9 H(9) = 10 H(10) = 12 H(11) = 15 H(12) = 16 H(13) = 18 H(14) = 20 H(15) = 24 H(16) = 25 H(17) = 27 H(18) = 30 H(19) = 32 H(20) = 36 H(1691) = 2125764000
zkl
<lang zkl>var BN=Import("zklBigNum"); // only needed for large N fcn hamming(N){
h:=List.createLong(N+1); (0).pump(N+1,h.write,Void); // fill list with stuff h[0]=1;
- if 1 // regular (64 bit) ints
x2:=2; x3:=3; x5:=5; i:=j:=k:=0;
- else // big ints
x2:=BN(2); x3:=BN(3); x5:=BN(5); i:=j:=k:=0;
- endif
foreach n in ([1..N]){ z:=(x2<x3) and x2 or x3; z=(z<x5) and z or x5; h[n]=z; if (h[n] == x2) { x2 = h[i+=1]*2 } if (h[n] == x3) { x3 = h[j+=1]*3 } if (h[n] == x5) { x5 = h[k+=1]*5 } } return(h[N-1])
} [1..20].apply(hamming).println(); hamming(1691).println();</lang>
- Output:
L(1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36) 2125764000
While the other algorithms save [lots of] space, run time still sucks when n > 100,000 so memory usage might as well too. Change the #if 0 to 1 and
- Output:
hamming(0d1_000_000).println(); 519312780448388736089589843750000000000000000000000000000000000000000000000000000000
Direct calculation through triples enumeration
OK, I was wrong, calculating the nth Hamming number can be fast and efficient.
as direct a translation as I can, except using a nested for loop instead of list comprehension (which makes it easier to keep the count).
<lang zkl>
- -- 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
var BN=Import("zklBigNum"); var lg3 = (3.0).log()/(2.0).log(), lg5 = (5.0).log()/(2.0).log(); fcn logval(i,j,k){ lg5*k + lg3*j + i } fcn trival(i,j,k){ BN(2).pow(i) * BN(3).pow(j) * BN(5).pow(k) } fcn estval(n){ (6.0*lg3*lg5*n).pow(1.0/3) } #-- estimated logval, base 2 fcn rngval(n){
if(n > 500000) return(2.4496 , 0.0076); #-- empirical estimation if(n > 50000) return(2.4424 , 0.0146); #-- correction, base 2 if(n > 500) return(2.3948 , 0.0723); #-- (dist,width) if(n > 1) return(2.2506 , 0.2887); #-- around (log $ sqrt 30),
return(2.2506 , 0.5771); #-- says WP }
fcn nthHam(n){ // -> (Double, (Int, Int, Int)) #-- n: 1-based: 1,2,3...
d,w := rngval(n); #-- correction dist, width hi := estval(n.toFloat()) - d; #-- hi > logval > hi-w c,b := band(hi,w); #-- total count, the band s := b.sort(fcn(a,b){ a[0]>b[0] }); #-- sorted decreasing, result m := c - n; #-- m 0-based from top nb := b.len(); #-- |band| res := s[m]; #-- result
if(w >= 1) throw(Exception.Generic("Breach of contract: (w < 1): " + w)); if(m < 0) throw(Exception.Generic("Not enough triples generated: " +c+n)); if(m >= nb)throw(Exception.Generic("Generated band is too narrow: " +m+nb)); return(res);
}
fcn band(hi,w){ //--> #-- total count, the band
b := Sink(List); cnt := 0; foreach k in ([0 .. (hi/lg5).floor()]){ p := lg5*k; foreach j in ([0 .. ((hi-p)/lg3).floor()]){ q := lg3*j + p; i,frac := (hi-q).modf(); r := hi-frac; #-- r = i + q
cnt+=(i+1); #-- total count if(frac<w) b.write(T(r,T(i,j,k))); #-- store it, if inside band
} } return(cnt,b.close());
}</lang> <lang>fcn printHam(n){
r,t:=nthHam(n); i,j,k:=t; h:=trival(i,j,k); println("Hamming(%,d)-->2^%d * 3^%d * 5^%d-->\n%s".fmt(n,i,j,k,h));
}
printHam(1691); //(5,12,3), 10 digits printHam(0d1_000_000); //(55,47,64), 84 digits printHam(0d10_000_000); //(80,92,162), 182 digits, 80 zeros at end printHam(0d1_000_000_000); //(1334,335,404), 845 digits</lang>
- Output:
Hamming(1,691)-->2^5 * 3^12 * 5^3--> 2125764000 Hamming(1,000,000)-->2^55 * 3^47 * 5^64--> 519312780448388736089589843750000000000000000000000000000000000000000000000000000000 Hamming(10,000,000)-->2^80 * 3^92 * 5^162--> 162441050638304318232392153117595750351085388205966408633356724833252116013682098127901554107666015625 <80 zeros> Hamming(1,000,000,000)-->2^1334 * 3^335 * 5^404--> 621607575556524486163081633287207200394705651908965270659163240.......
- Programming Tasks
- Solutions by Programming Task
- Ada
- ALGOL 68
- ATS
- AutoHotkey
- AWK
- BBC BASIC
- Bracmat
- C
- C sharp
- C++
- Clojure
- CoffeeScript
- Common Lisp
- D
- DCL
- Eiffel
- Elixir
- ERRE
- F Sharp
- Factor
- Forth
- Fortran
- FunL
- Go
- Haskell
- Icon
- Unicon
- J
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- Liberty BASIC
- Logo
- Lua
- Mathematica
- MATLAB
- Octave
- MUMPS
- Nim
- Bigints
- Oz
- PARI/GP
- Pascal
- Perl
- Perl 6
- PicoLisp
- PL/I
- Prolog
- Python
- Qi
- Qi examples needing attention
- Examples needing attention
- R
- Racket
- Raven
- REXX
- Ruby
- Run BASIC
- Scala
- Scheme
- Seed7
- Smalltalk
- Tcl
- UBasic/4tH
- UNIX Shell
- Ursala
- VBScript
- Zkl