11 apr 2009

Hamming numbers

They are also called regular numbers and are described here Regular numbers. See also A051037.

Inspired by Groeneveld who translated the Haskell solution Hamming's problem (its links to Dijkstra don't work, they should be Dijkstra) in

next=: {. (>:@[ , ] /:~@~.@, 2 3 5 * {) }.
hamming =: {. [: }. next^:(]`(0 1x"_))

I tried to generate these numbers as efficient as possible and came up with

HN1=: 3 : 0
h0=. (^.2 3 5) <.@%~ h2=. 3 %: y * 6 * */h1=. ^.2 3 5
h3=. (#~h2 >: h1 +/@:*"1 ]) >,{<@i."0 h0
y{./:~2 3 5x */@:^"1 h3
)

which appeared to be the tetrahedron method in Regular numbers.

An other algorithm appeard to be very lean in J, be it somewhat slower than HN1, but overall much more efficient.

HN2=: 3 : 0
z=. 1
t=. 0,,:~2 3 5x
for_j. i.<:y do.
 z=. z,a=. <./{:t
 a=. a = {:t
 t=. ((1{t) * z {~ {.t) _1} t=. (a + {.t) 0} t
end.
z

Here are some results. See J-blog for dspl, rnkng and scores.

   dspl rnkng scores'HN1 2000';'HN2 2000'
+----+-------+------+------+
|rank|rlt_prf|rlt_et|rlt_sz|   NB. relative performance, execution time and -size repsectively.
+----+-------+------+------+
| 1  |15.40  |1.00  |112.99|
+----+-------+------+------+
| 0  | 1.00  |7.34  |  1.00|
+----+-------+------+------+

   (HN2 -:HN1)2000
1

   dspl rnkng scores'hamming 1000';'HN1 1000';'HN2 1000'
+----+-------+------+------+
|rank|rlt_prf|rlt_et|rlt_sz|
+----+-------+------+------+
| 2  |16.45  |37.54 | 1.88 |
+----+-------+------+------+
| 1  | 3.39  | 1.00 |14.50 |
+----+-------+------+------+
| 0  | 1.00  | 4.28 | 1.00 |
+----+-------+------+------+

   2-:/\(hamming 1000);(HN1 1000);HN2 1000
1 1

 HN0  was at the base of HN1 but was skipped as it was very slow. Some slight improvements changed the scene.

HN0=: 3 : 0
h0=. h1 <.@%~ 3 %: y * 6 * */h1=. ^. 2 3 5
y {. /:~2 3 5x */@:^"1 >,{<@i."0 h0
)

HN00=: 3 : 0
h1=. h0 <.@%~ 3 %: y * 6 * */h0=. 2^. 2 3 5
h4=. (i.y {. /:~) h3=. h0 +/@:*"1 h2=. >,{<@i."0 h1
h5=. 56 > h4 { h3
(x: 2 3 5 */@:^"1 h5 # h6) , 2 3 5x */@:^"1 (-.h5) # h6=. h4 { h2 
)

   ('HN0';'HN1';'HN2';'HN00') dspl rnkng scores ('HN0 2000');('HN1 2000');('HN2 2000');'HN00 2000'
+----+----+-------+------+------+
|verb|rank|rlt_prf|rlt_et|rlt_sz|
+----+----+-------+------+------+
|HN0 | 3  |40.85  | 5.69 |75.88 |
+----+----+-------+------+------+
|HN1 | 2  | 1.47  | 1.26 |12.33 |
+----+----+-------+------+------+
|HN2 | 0  | 1.00  |10.56 | 1.00 |
+----+----+-------+------+------+
|HN00| 1  | 1.14  | 1.00 |12.08 |
+----+----+-------+------+------+

   2-:/\(HN0;HN1;HN2;HN00) 2000
1 1 1

Thanks to Groeneveld who pointed out an error in previous version of HN00 .

The next improvement was remarkable since I had not expected that such a simple adaptation of HN2 would make such a difference in footprint.

HN21=: 3 : 0
z=. 1
't0 t1 t2'=. 0 , ,:~2 3 5x
for_j. i.<:y do.
 z=. z , a=. <./ t2
 a=. a = t2
 t0=. a + t0
 t2=. ((a { t1) * z {~ a { t0) (a=. a # 0 1 2)} t2
end.
z
)

   ('HN1';'HN00';'HN2';'HN21') dspl rnkng scores ('HN1 2000');('HN00 2000');('HN2 2000');'HN21 2000'
+----+----+-------+------+------+
|verb|rank|rlt_prf|rlt_et|rlt_sz|
+----+----+-------+------+------+
|HN1 | 3  |12.60  | 1.29 | 95.09|
+----+----+-------+------+------+
|HN00| 2  |10.98  | 1.00 |106.71|
+----+----+-------+------+------+
|HN2 | 1  | 8.33  |10.12 |  7.99|
+----+----+-------+------+------+
|HN21| 0  | 1.00  | 9.72 |  1.00|
+----+----+-------+------+------+

   (HN21-:HN2) 2000
1

   {: HN21 1e4
288325195312500000

See Hui for more elegant code with a slight improvement:

hn21b=: 3 : 0
i=. 0
m=. 2 3 5x
z=. 1x
for. i.<:y do.
 z=. z , a=. <./ m
 b=. a = m
 i=. i + b
 m=. (m*-.b) + (b*2 3 5) * i{z
end.
z
)

   dspl rnkng scores'HN21 2000';'hn21b 2000'
+----+-------+------+------+
|rank|rlt_prf|rlt_et|rlt_sz|
+----+-------+------+------+
| 1  |1.03   |1.03  |1.00  |
+----+-------+------+------+
| 0  |1.00   |1.00  |1.00  |
+----+-------+------+------+

See Essay/Hamming Numbers for a lightning fast verb nh .

   dspl rnkng scores'HN21 2000';'nh 2000'
+----+--------+-------+------+
|rank|rlt_prf |rlt_et |rlt_sz|
+----+--------+-------+------+
| 1  |13208.42|1363.03|9.69  |
+----+--------+-------+------+
| 0  |    1.00|   1.00|1.00  |
+----+--------+-------+------+

RE Boss/J-blog/HN (last edited 2010-01-08 08:35:43 by RE Boss)