home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / comp / arch / 10530 < prev    next >
Encoding:
Internet Message Format  |  1992-11-08  |  4.8 KB
  1. Xref: sparky comp.arch:10530 comp.lang.misc:3546
  2. Newsgroups: comp.arch,comp.lang.misc
  3. Path: sparky!uunet!ferkel.ucsb.edu!taco!gatech!darwin.sura.net!spool.mu.edu!news.cs.indiana.edu!noose.ecn.purdue.edu!mentor.cc.purdue.edu!pop.stat.purdue.edu!hrubin
  4. From: hrubin@pop.stat.purdue.edu (Herman Rubin)
  5. Subject: Re: Hardware Support for Numeric Algorithms
  6. Message-ID: <BxExBt.C32@mentor.cc.purdue.edu>
  7. Summary: Detailed derivation of an algorithm
  8. Keywords: n
  9. Sender: news@mentor.cc.purdue.edu (USENET News)
  10. Organization: Purdue University Statistics Department
  11. References: <1992Nov5.202412.7266@linus.mitre.org> <1992Nov6.075459.13348@coe.montana.edu> <BxAJMA.DKy@mentor.cc.purdue.edu>
  12. Date: Sun, 8 Nov 1992 19:40:40 GMT
  13. Lines: 91
  14.  
  15.  
  16. I received this by email.  
  17.  
  18. Herman: 
  19. Please consider posting this to the net, as I feel it does help to 
  20. understand the problem. Thanks. Bob
  21.  
  22.  Here is a detailed description of what the algorithm is trying to do,
  23.  and how it is trying to do it.  The object is to produce two arrays
  24.  of random variables, an array n of integers, and an array m of masks,
  25.  such that the elements of n are distributed as the integer parts of
  26.  exponential random variables with mean 2, and such that if one takes
  27.  a uniform fixed-point random variable independent of a mask, and zeroes
  28.  out the places where the mask has a 1, the resulting fixed point random
  29.  variable is distributed as the fractional part of an exponential random
  30.  variable with mean 2.  This information is enough to specify the distribution,
  31.  but there are infinitely many ways to do this.
  32.  
  33.  Now the algorithm.  We start with the known result that, if we want 
  34.  exponentials with any mean m, we can do this by taking Poisson random
  35.  variables with mean 1/m until a non-zero result j occurs.  The integer
  36.  part is the number of 0's until this happens, and the fractional part
  37.  is distributed as the minimum of j uniform random variables.  This last
  38.  distinction is important, as we can get this much more cheaply than
  39.  getting j uniform random variables.  The technique used here is to
  40.  do things bitwise; if p_jk is the usual probability that of j random
  41.  bits k are 0 (or 1), with probability p_j0 + p_jj we do nothing about
  42.  the bit being examined, but for other k's we set the bit in the mask
  43.  and replace j by k.  Once k gets down to 1, this necessarily stops.
  44.  
  45.  We use the explicit form of p_jk for j < 6, and for simplicity I will
  46.  use q_j for p_j0 + p_jj = 2*p_j0.  For some cases, it is more convenient
  47.  to use r_jk = p_jk/(1-q_j).
  48.  
  49.      q_2 = 1/2    r_21 = 1/2
  50.      q_3 = 1/4    r_31 = r_32 = 1/2
  51.      q_4 = 1/8    p_31 = p_33 = 1/4    p_32 = 3/8
  52.      q_5 = 1/16    r_31 = r_34 = 1/6    r_32 = r_33 = 1/3
  53.  
  54.  Now the details.  We start out with the detailed expansion of exp(.5)/2,
  55.  which is given by a_j = .5^(j+1)/j!.  However, for 6 or more, we will
  56.  first multiply a_6 by 1.5.  Bit patterns will be allocated to the 
  57.  outcomes so that the probability of returning j is a_j, and the remaining
  58.  time a rejection occurs.  We find that easy by-products are, in some cases,
  59.  random variables g4 which have a geometric distribution with p = 3/4, which
  60.  happens to be the distribution of the number of places from the current 
  61.  place which at which the bit mask will be set for j=3, and g16 for j=5.
  62.  A geometric (1/2) random variable gets the location of a random bit.
  63.  
  64.  Let us consider the expansion, in hex and in binary
  65.  
  66.  0    .8            .1000 ....
  67.  1    .4            .0100 ....
  68.  2    .1            .0001 0000 ....
  69.  3    .02aaaa....              0010 1010 1010 1010 1010 ...
  70.  4    .005555....               0101 0101 0101 0101 ...
  71.  5    .0008888 ...                1000 1000 1000 ...
  72.  6*1.5    .0001111 ...                0001 0001 0001 ...
  73.  
  74.  remainder            .0010 1100 1111 0110 0110 0110 ...
  75.  so we allocate the bits:
  76.                   01R2 RR3x RRRR 5RRy 5RRy 5RRy ...
  77.  g16                                            1111 2222 3333 ...
  78.  
  79.  where R means reject, x means 3 or 4, y means 6 or more or reject.
  80.  Now look at 3 and x.  In the x case, we get a geometric to find the
  81.  bit in the binary representation.
  82.  
  83.      10 1111 1111 1111 1111 ...
  84.      3  3434 3434 3434 3434 ...
  85.      1  2132 4354 6576 8798 ...
  86.  
  87.  In the y situation the probability is 1.5*a_6.  We now get a geometric
  88.  random variable g, use it to both get a 2/3 - 1/3 division and also a g4.
  89.  The 2/3 gives us the right value for a_6.  The 1/3 is half of that, or
  90.  7*a_7.  We use a 2/7 acceptance to gat 2*a_7.  We now loop, starting with
  91.  2*a_j
  92.  
  93.      with probability 1/2, accept j        probability a_j of acceptance
  94.      with probability j/j+1, reject        probability 2*a_(j+1) of 
  95.                          still going
  96.       augment j and proceed.
  97.  
  98.  I could similarly go into details about the procedures for forming or extending
  99.  the mask when the active j is 5, 4, 3, or 2.  In some cases, byproducts of 
  100.  previous stages are used.
  101. -- 
  102. Herman Rubin, Dept. of Statistics, Purdue Univ., West Lafayette IN47907-1399
  103. Phone: (317)494-6054
  104. hrubin@pop.stat.purdue.edu (Internet, bitnet)  
  105. {purdue,pur-ee}!pop.stat!hrubin(UUCP)
  106.