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
