Going back to the original Haskell source in http://dobbscodetalk.com/index.php?option=com_content&task=view&id=913&Itemid=85 I see that the code can be fixed readily to handle all valid arguments, viz., when y<:5e4, use 1.56 and 0.2 instead of 1.693 and 0.01 as the magic constants. Thus:
NB. the y-th Hamming number (1-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*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 */2 3 5x ^ (y-~(#s)++/s) { :: 0: z\:z+/ .*^.2 3 5 ) Check (using hn21b from a few days ago): (nh >:i.3000) -: hn21b 3000 1 ----- Original Message ----- From: Roger Hui <[email protected]> Date: Thursday, January 7, 2010 13:56 Subject: Re: [Jprogramming] Rosetta code: hamming numbers To: Programming forum <[email protected]> > I have convinced myself that the quantity d in nh4 > is boolean, enabling further simplifications. Thus: > > 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 > ) > > > > ----- Original Message ----- > From: Roger Hui <[email protected]> > Date: Thursday, January 7, 2010 12:19 > Subject: Re: [Jprogramming] Rosetta code: hamming numbers > To: Programming forum <[email protected]> > > > nh4 obtained by manipulation of the J code for nthHam > > without much understanding of the underlying mathematics. > > > > 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 > > ) > > > > > > > > ----- Original Message ----- > > From: Roger Hui <[email protected]> > > Date: Thursday, January 7, 2010 11:38 > > Subject: Re: [Jprogramming] Rosetta code: hamming numbers > > To: Programming forum <[email protected]> > > > > > The }.z (and }.t) in your code are necessary because > > > z and t are initialized incorrectly, engendering extra > > > leading rows of 0. If instead you initialize > > > t=. i.0 5 and z=. i.0 4 then drops > > > can be removed. > > > > > > Anyway I have a non-looping and clearer version > > > that provides a factor of 10 improvement for y=1e6. > > > I will post that version in a little while. > > > > > > > > > > > > ----- Original Message ----- > > > From: Aai <[email protected]> > > > Date: Thursday, January 7, 2010 11:21 > > > Subject: Re: [Jprogramming] Rosetta code: hamming numbers > > > To: Programming forum <[email protected]> > > > > > > > Correction to make the last assert work correctly: > > > > > > > > nthHam=: 3 : 0 > > > > lns=. 'ln2 ln3 ln5'=. ^. tdv=.2 3 5 > > > > lo=. 0.01 -~ hi=. 1.693-~(*/6,y,lns)^%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 + 3{"1 }. t > > > > z=.,:0$0 > > > > for_r. (#~(2&{ <: 3&{)"1) }. t do. > > > > for_i. ([+i.@>:@-~)/2 3{r do. z=.z, i,(2{.r), > > > > ({:r)+i*ln2 end. > > > > end. > > > > assert. 0 <: c =. c - y > > > > assert. c < # }. z > > > > __ q: inv tdv ,: x: 3 {. c { (\: {:"1) }. z > > > > ) > > > > > > > > It shouldn't be: > > > > > > > > nthHam 31 > > > > 1 > > > > > > > > but: > > > > > > > > nthHam 31 > > > > |assertion failure: nthHam > > > > | c<#}.z > > > > > > > > > > > > Anyway, if you dare using it, use nthHam for N>:50000 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
