No, the algorthm Knuth mentions deals a uniform deviate, looks at the top 6 bits, and for most of the possibilities returns the answer right away; for a few it does some more computation to get the answer.
Henry Rich > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of > Viktor Cerovski > Sent: Thursday, August 14, 2008 5:01 PM > To: [email protected] > Subject: RE: [Jprogramming] generating standard normal variates > > > > Henry Rich wrote: > > > > I'm not following this thread closely, but Algorithm M in section > > 3.4.1 of TAOCP is a very fast way of generating > > normal deviates. It can be modified to support generating > > them in batches. > > > > Knuth says (p. 113 of my 1969 edition) "a general-purpose > > subroutine based on Algorithm M will be a valuable part of > > any subroutine library". > > > > Henry Rich > > > Presumably it is about the code that avoids sin and cos of > the Box-Mueller code by picking a random point within the unit > circle and then doing a transformation involving only ln > and square root. > > It can be done like this: > > F=:(1-~2*(?,?)@0:)^:(1<:+&*:/)^:_ > Nrnd=:({.,)(*([:%:@([EMAIL PROTECTED]:@^. % ])+&*:/))@F"0@($1:)@>[EMAIL > PROTECTED]: > > F picks the point, Nrnd does the rest. There is doubling of > computation of +&*:/ for every found point. > > The following mathematically equivalent routine done via > complex numbers is a bit faster: > > ReIm=:9 11&o. > G=:(1j1-~2*(?j.?)@0:)^:(1<:|)^:_ > Nrnd1=:({.,)(((%:@[EMAIL PROTECTED] % ])@| * ReIm)"[EMAIL > PROTECTED]"0@($1:))@>[EMAIL PROTECTED]: > > but both are much slower than other tacit routines from this thread. > > -- > View this message in context: > http://www.nabble.com/generating-standard-normal-variates-tp18 > 898408s24193p18989673.html > Sent from the J Programming mailing list archive at Nabble.com. > > ---------------------------------------------------------------------- > For information about J forums see > http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
