Essays/Hamming Number

From J Wiki
Jump to: navigation, search

A Hamming number (aka regular number aka 5-smooth number) is an integer of the form (2^i)*(3^j)*(5^k) for non-negative integers i, j, and k. The first 20 Hamming numbers are 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36.


nh

The verb nh finds the y-th Hamming number (0-origin). It is translated from a Haskell function presented here.

nh=: 3 : 0 " 0
 hi=. (1.56 1.693{~y>5e4) -~ 3%:6*(1+y)**/'ln2 ln3 ln5'=. ^.2 3 5
 k=. i.1+<.hi%ln5              NB. exponents for 5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ;i.&.>m                   NB. exponents for 3
 k=. m#k
 s=. <.p=. ln2%~hi-(j*ln3)+k*ln5
 i=. >.p-ln2%~0.2 0.01{~y>5e4  NB. exponents for 2
 z=. (i=s)#i,.j,.k
 c=. ((#s)++/s)-1+y
 assert. 0<:c
 */2 3 5x ^ c { :: 0: z\:z+/ .*^.2 3 5
)

For example:

   nh i.20
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36

   nh 1e6
519381797917090766274082018159448243742493816603938969600000000000000000000000000000

   _50 ]\ ": t=: nh 1e9
62160789282005611071334298264867143857669458682016
16271891228224144262383715571325506849547578511309
06044868429977012580044918560955998407925870165832
00725201893976680899723970584153571737144573016210
05622628546074961245527554715641510294092346419243
91465509121412949352977345288230138936337661611574
05935601838226002405329313436790825046261524824703
63136000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000
   _ q: t
761 572 489
   t = */ 2 3 5x ^ 761 572 489
1

Derivation

Here we present the intermediate steps from the original Haskell function to nh .

nthHam

The initial J version was posted to the J Programming Forum on 2010-01-07. The argument y is the 1-origin index of the desired Hamming number, and must be greater than 5e4.

nthHam=: 3 : 0
 'ln2 ln3 ln5'=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~(*/6,y,ln2, ln3, ln5)^%3
 t=.,:0$0
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
 c=. +/ 1+_2{"1 }. t
 z=.,:0$0
 for_r. t do.
   for_i. ([+i.@>:@-~)/2{.2}.r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
 c=. c - y
 assert. 0 <: c
 assert. c <: # z
 __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z
)

A slightly improved version was subsequently posted, which processes immaterial rows of the pronoun t more efficiently. (An extra line just before the for_r. loop.)

nthHam1=: 3 : 0
 'ln2 ln3 ln5'=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~(*/6,y,ln2, ln3, ln5)^%3
 t=.,:0$0
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
 c=. +/ 1+_2{"1 }. t
 z=.,:0$0
 t=. (#~(2&{<:3&{)"1) t
 for_r. t do.
   for_i. ([+i.@>:@-~)/2{.2}.r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
 c=. c - y
 assert. 0 <: c
 assert. c <: # z
 __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z
)

nthHam is a translation of a Haskell function from here.

-- find n-th (base 1) Hamming number directly
-- by Will Ness, based on "band" idea by Louis Klauder

ln2 = log 2; ln3 = log 3; ln5 = log 5
trival (i,j,k) = 2^i * 3^j * 5^k          -- use system's bignums

nthHam n
  | m < 0   = error $ "Not enough triples generated: " ++ show (c,n)
  | m >= hn = error $ "Generated band too narrow: " ++ show (m,hn)
  | True    = ( (trival (fst x), fst x), (m, hn) )
  where
  hi = (6 * ln2 * ln3 * ln5 * n)**(1/3) - 1.693 -- good for n>50,000
  lo = hi - 0.01                                -- or else -1.56 -0.2
  ts = [ (imax + 1, [ triple | i <- [imin..imax],
                               let triple = ((i,j,k), q + i*ln2) ])
         | k <- [ 0 .. floor ( hi   /ln5) ],  let p =     k*ln5,
           j <- [ 0 .. floor ((hi-p)/ln3) ],  let q = p + j*ln3,
           let imax =  floor ((hi-q)/ln2)
               imin = ceiling((lo-q)/ln2) ]
  c  = foldr ((+).fst) 0 ts               -- total count
  m  = c – n                              -- target index in the band
  hn = length h                           -- band's length
  h  = concatMap snd ts                   -- band of triples
  hs = sortBy (flip compare `on` snd) h   -- in *descending* order
  x  = hs !! m                            -- the n-th Hamming number

nh2

nh2 differs from nthHam1 only in using ln=.^.2 3 5 instead of ln2,ln3,ln5 .

nh2=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 t=.,:0$0
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
 c=. +/ 1+_2{"1 }. t
 z=.,:0$0
 t=. (#~(2&{<:3&{)"1) t
 for_r. t do.
   for_i. ([+i.@>:@-~)/2{.2}.r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
 c=. c - y
 assert. 0 <: c
 assert. c <: # z
 __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z
)

nht

We focus on the first loops in nh2 , the for_k. and for_j. nested loops which compute the pronoun t .

nh2a simply copies the lines from nh2 , together with the lines necessary to initialize the required pronouns. In addition, t is initialized to i.0 5 instead of ,:0$0 (a 1-by-0 matrix). The latter initialization leads to t incorrectly having a leading row of 0s.

nh2a=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 t=. i.0 5
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
)

   $ nh2a 5e4
1437 5

   10 {. nh2a 5e4
0 0 101 100       0
1 0 100  99 1.09861
2 0  98  97 2.19722
3 0  97  96 3.29584
4 0  95  94 4.39445
5 0  93  92 5.49306
6 0  92  91 6.59167
7 0  90  89 7.69029
8 0  89  88  8.7889
9 0  87  86 9.88751

It is seen from the key line t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p+j*ln3 that entries 2 3 4 of t can be computed outside of the nested loops, from j and k and from pronouns that remain constant throughout. Thus:

nh2b=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 t=. i.0 2
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k
   end.
 end.
 'j k'=. |: t
 p=. k*ln5
 q=. p+j*ln3
 t,.(>.ln2%~lo-q),.(<.ln2%~hi-q),.q
)

   (nh2a -: nh2b) 5e4
1

It is a short step to see that j is the raze of i.&.>f k for a simple function f of k ,  and i{k itself is repeated f i{k times to make up column 0 of t .  Both loops are rendered unnecessary.

nht=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. 2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 k=. i.1+<.hi%ln5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ; i.&.> m
 k=. m#k
 p=. k*ln5
 q=. p+j*ln3
 r=. >.ln2%~lo-q
 s=. <.ln2%~hi-q
 j,.k,.r,.s,.q
)

   (nh2a -: nht) 50
1

nhz

We turn our attention to the computation of z in nh2 .  The verb nhz1 simply copies the for_r. and for_i. loops from nh2 ,  renaming t to y to account for that t is now an argument to the verb. z is initialized to 0 4$0 because it is a 4-column matrix, rather than ,:0$0 (a 1-by-0 matrix, which leads to z incorrectly having a leading row of 0).

nhz1=: 3 : 0
 ln2=. ^.2
 y=. (#~(2&{<:3&{)"1) y
 z=. 0 4$0
 for_r. y do.
   for_i. ([+i.@>:@-~)/2 3{r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
)

   $ nhz1 nht 5e4
21 4
   nhz1 nht 5e4
28 46  0 69.9443
97  1  1 69.9433
13 54  1 69.9454
82  9  2 69.9445
67 17  3 69.9456
52 25  4 69.9467
11 45  8 69.9377
65  8 10 69.9378
50 16 11  69.939
35 24 12 69.9401
20 32 13 69.9412
 5 40 14 69.9424
59  3 16 69.9425
44 11 17 69.9437
29 19 18 69.9448
14 27 19 69.9459
27 10 25  69.937
12 18 26 69.9382
21  5 31 69.9417
 6 13 32 69.9429
15  0 37 69.9464

We notice that the last column of z is a simple function of i (column 0 of z) and the last column of y (aka t aka r), and does not need to be carried in z . Thus:

nhz2=: 3 : 0
 ln2=. ^.2
 y=. (#~(2&{<:3&{)"1) y
 z=. 0 3$0
 for_r. y do.
   for_i. ([+i.@>:@-~)/2 3{r do. z=.z, i,2{.r end.
 end.
)

   (}:"1@nhz1 -: nhz2) nht 5e4
1

Columns 1 2 of z are simply (the rows of) columns 2 3 of y repeated an appropriate number of times. The key column is column 0, the quantity i , readily seen to be the raze of  y2 +&.> i.&.> 1+y3-y2 , where y2 and y3 are columns 2 and 3 of y . Both loops are rendered unnecessary. Thus:

nhz=: 3 : 0
 y=. (>:/"1@(3 2&({"1)) # ]) y
 i=. ; (2{"1 y) +&.> i.&.> 1+-/"1 ]3 2{"1 y
 i,.2{."1 y
)

   (}:"1@nhz1 -: nhz) nht 5e4
1

Actually, nhz assumes that 1 is the "appropriate number of times" that (the rows of) columns 2 3 of y are repeated. This (for all we know) could be a mistake. (On the positive side, if ever 1 is not the appropriate number, the verb would signal a length error rather than return an erroneous result.) More about that later.

nh3

Putting it all together:

nh3=: 3 : 0
 t=. nht y
 c=. y -~ +/ 1+_2{"1 t
 z=. nhz t
 assert. 0 <: c
 assert. c < # z
 */2 3 5 ^ x: c{z\:+/"1 z*"1 ^.2 3 5
)

   nh2 5e4
2379126835648701693226200858624
   nh3 5e4
2379126835648701693226200858624

   (nh2 = nh3)"0 ] 5e4 1e5 1e6
1 1 1

From nh2 to nh3 , we also changed last line: first, __ q: inv is replaced with a more direct and clearer computation; second, the last column of the original z computed in the inner loop discussed in nhz is computed using non-looping vector operations.

nh4

The component verbs nht and nhz were good intermediate steps, but their use actually obscures some of the goings-on. Merging everything together, we see readily that i (column 0 of z), j (column 0 of t and column 1 of z), and k (column 1 of t and column 2 of z) are the exponents for 2, 3, and 5, respectively; the 3 columns of z are the exponents i , j , and k . The phrase z\:+/"1 z*"1 ^.2 3 5 on the last line is ordering the rows of z by +/"1 z*"1 ^.2 3 5 , itself the logs of the numbers represented by the exponents. (Ordering by logs produces the same sequence as ordering by the numbers themselves.) We note that the exponents, when separated from other quantities, are machine integers, making their handling simpler and more efficient.

nh4=: 3 : 0
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,'ln2 ln3 ln5'=. ^. 2 3 5
 k=. i.1+<.hi%ln5       NB. exponents for 5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ;i.&.>m            NB. exponents for 3
 k=. m#k
 r=. >.ln2%~lo-q=. (j*ln3)+k*ln5
 s=. <.ln2%~hi-q
 d=. 0>.1+s-r
 i=. (d#r)+;i.&.>d-.0   NB. exponents for 2
 z=. i,.d#j,.k
 c=. y -~ (#s)++/s
 assert. 0 <: c
 assert. c < #z
 */2 3 5x^c{z\:z+/ .*^.2 3 5
)

   (nh2 = nh4)"0 ] 5e4 1e5 1e6
1 1 1

nh5

We focus on the question of "appropriate number of times" that first arose in nhz . This is now the quantity d in nh4.

 r=. >.ln2%~lo-q=. (j*ln3)+k*ln5
 s=. <.ln2%~hi-q
 d=. 0>.1+s-r

Simple algebraic manipulation gives:

 p=. ln2%~hi-(j*ln3)+k*ln5
 s=. <.p
 r=. >.p-0.01
 d=. 0>.1+s-r

From which we conclude that either s=r or s=r-1 , whence d is boolean. A few further simplifications are enabled thereby:

nh5=: 3 : 0
 hi=. 1.693-~ 3%:*/6,y,'ln2 ln3 ln5'=. ^.2 3 5
 k=. i.1+<.hi%ln5       NB. exponents for 5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ;i.&.>m            NB. exponents for 3
 k=. m#k
 s=. <.p=. ln2%~hi-(j*ln3)+k*ln5
 i=. >.p-0.01%ln2       NB. exponents for 2
 z=. (i=s)#i,.j,.k
 c=. y -~ (#s)++/s
 assert. (0<:c)*0<#z
 */2 3 5x^c{z\:z+/ .*^.2 3 5
)

   (nh2 = nh5)"0 ] 5e4 1e5 1e6
1 1 1

nh

The restriction that the argument must be greater than 5e4 is irksome. Going back to the original Haskell function, we see that the restriction can be lifted if we replace the magic constants 1.693 and 0.01 by 1.56 and 0.2 when y is bounded by 5e4. Moreover, using 1+y instead of y accommodates the 0-origin indexing used in the rest of the language. The verb nh obtains as a result.

In the last line of nh , the adversary in { :: 0: is invoked only when the argument y=0 .

Collected Definitions

NB. the y-th Hamming number (0-origin indexing)
NB. http://dobbscodetalk.com/index.php?option=com_content&task=view&id=913&Itemid=85

nh=: 3 : 0 " 0
 hi=. (1.56 1.693{~y>5e4) -~ 3%:6*(1+y)**/'ln2 ln3 ln5'=. ^.2 3 5
 k=. i.1+<.hi%ln5              NB. exponents for 5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ;i.&.>m                   NB. exponents for 3
 k=. m#k
 s=. <.p=. ln2%~hi-(j*ln3)+k*ln5
 i=. >.p-ln2%~0.2 0.01{~y>5e4  NB. exponents for 2
 z=. (i=s)#i,.j,.k
 c=. ((#s)++/s)-1+y
 assert. 0<:c
 */2 3 5x ^ c { :: 0: z\:z+/ .*^.2 3 5
)

NB. the y-th Hamming number (1-origin indexing), y>5e4

nthHam=: 3 : 0
 'ln2 ln3 ln5'=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~(*/6,y,ln2, ln3, ln5)^%3
 t=.,:0$0
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
 c=. +/ 1+_2{"1 }. t
 z=.,:0$0
 for_r. t do.
   for_i. ([+i.@>:@-~)/2{.2}.r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
 c=. c - y
 assert. 0 <: c
 assert. c <: # z
 __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z
)

nthHam1=: 3 : 0
 'ln2 ln3 ln5'=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~(*/6,y,ln2, ln3, ln5)^%3
 t=.,:0$0
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
 c=. +/ 1+_2{"1 }. t
 z=.,:0$0
 t=. (#~(2&{<:3&{)"1) t
 for_r. t do.
   for_i. ([+i.@>:@-~)/2{.2}.r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
 c=. c - y
 assert. 0 <: c
 assert. c <: # z
 __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z
)

nh2=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 t=.,:0$0
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
 c=. +/ 1+_2{"1 }. t
 z=.,:0$0
 t=. (#~(2&{<:3&{)"1) t
 for_r. t do.
   for_i. ([+i.@>:@-~)/2{.2}.r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
 c=. c - y
 assert. 0 <: c
 assert. c <: # z
 __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z
)

nh2a=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 t=. i.0 5
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k,(>.ln2%~lo-q),(<.ln2%~hi-q), q=. p + j * ln3
   end.
 end.
)

nh2b=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. tdv=.2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 t=. i.0 2
 for_k. i. 1+ <. hi%ln5 do.
   for_j. i. 1+ <. ln3 %~ hi - p=. k*ln5 do.
     t=. t, j,k
   end.
 end.
 'j k'=. |: t
 p=. k*ln5
 q=. p+j*ln3
 t,.(>.ln2%~lo-q),.(<.ln2%~hi-q),.q
)

nht=: 3 : 0
 'ln2 ln3 ln5'=. ln=. ^. 2 3 5
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,ln
 k=. i.1+<.hi%ln5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ; i.&.> m
 k=. m#k
 p=. k*ln5
 q=. p+j*ln3
 r=. >.ln2%~lo-q
 s=. <.ln2%~hi-q
 j,.k,.r,.s,.q
)

nhz1=: 3 : 0
 ln2=. ^.2
 y=. (#~(2&{<:3&{)"1) y
 z=. 0 4$0
 for_r. y do.
   for_i. ([+i.@>:@-~)/2 3{r do. z=.z, i,(2{.r), ({:r)+i*ln2 end.
 end.
)

nhz2=: 3 : 0
 ln2=. ^.2
 y=. (#~(2&{<:3&{)"1) y
 z=. 0 3$0
 for_r. y do.
   for_i. ([+i.@>:@-~)/2 3{r do. z=.z, i,2{.r end.
 end.
)

nhz=: 3 : 0
 y=. (>:/"1@(3 2&({"1)) # ]) y
 i=. ; (2{"1 y) +&.> i.&.> 1+-/"1 ]3 2{"1 y
 i,.2{."1 y
)

nh3=: 3 : 0
 t=. nht y
 c=. y -~ +/ 1+_2{"1 t
 z=. nhz t
 assert. 0 <: c
 assert. c < # z
 */2 3 5^x: c{z\:+/"1 z*"1 ^.2 3 5
)

nh4=: 3 : 0
 lo=. 0.01 -~ hi=. 1.693-~ 3%:*/6,y,'ln2 ln3 ln5'=. ^. 2 3 5
 k=. i.1+<.hi%ln5       NB. exponents for 5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ;i.&.>m            NB. exponents for 3
 k=. m#k
 r=. >.ln2%~lo-q=. (j*ln3)+k*ln5
 s=. <.ln2%~hi-q
 d=. 0>.1+s-r
 i=. (d#r)+;i.&.>d-.0   NB. exponents for 2
 z=. i,.d#j,.k
 c=. y -~ (#s)++/s
 assert. 0 <: c
 assert. c < #z
 */2 3 5x^c{z\:z+/ .*^.2 3 5
)

nh5=: 3 : 0
 hi=. 1.693-~ 3%:*/6,y,'ln2 ln3 ln5'=. ^.2 3 5
 k=. i.1+<.hi%ln5       NB. exponents for 5
 m=. 1+<.ln3%~hi-k*ln5
 j=. ;i.&.>m            NB. exponents for 3
 k=. m#k
 s=. <.p=. ln2%~hi-(j*ln3)+k*ln5
 i=. >.p-0.01%ln2       NB. exponents for 2
 z=. (i=s)#i,.j,.k
 c=. y -~ (#s)++/s
 assert. (0<:c)*0<#z
 */2 3 5x^c{z\:z+/ .*^.2 3 5
)



See also



Contributed by Roger Hui.