Essays/Hamming Number
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.