From: hrubin@pop.stat.purdue.edu (Herman Rubin)
Subject: Re: Hardware Support for Numeric Algorithms
Summary: Detailed derivation of an algorithm
Keywords: n
Organization: Purdue University Statistics Department
Date: Sun, 8 Nov 1992 19:40:40 GMT


I received this by email.  

Herman: 
Please consider posting this to the net, as I feel it does help to 
understand the problem. Thanks. Bob

 Here is a detailed description of what the algorithm is trying to do,
 and how it is trying to do it.  The object is to produce two arrays
 of random variables, an array n of integers, and an array m of masks,
 such that the elements of n are distributed as the integer parts of
 exponential random variables with mean 2, and such that if one takes
 a uniform fixed-point random variable independent of a mask, and zeroes
 out the places where the mask has a 1, the resulting fixed point random
 variable is distributed as the fractional part of an exponential random
 variable with mean 2.  This information is enough to specify the distribution,
 but there are infinitely many ways to do this.
 
 Now the algorithm.  We start with the known result that, if we want 
 exponentials with any mean m, we can do this by taking Poisson random
 variables with mean 1/m until a non-zero result j occurs.  The integer
 part is the number of 0's until this happens, and the fractional part
 is distributed as the minimum of j uniform random variables.  This last
 distinction is important, as we can get this much more cheaply than
 getting j uniform random variables.  The technique used here is to
 do things bitwise; if p_jk is the usual probability that of j random
 bits k are 0 (or 1), with probability p_j0 + p_jj we do nothing about
 the bit being examined, but for other k's we set the bit in the mask
 and replace j by k.  Once k gets down to 1, this necessarily stops.
 
 We use the explicit form of p_jk for j < 6, and for simplicity I will
 use q_j for p_j0 + p_jj = 2*p_j0.  For some cases, it is more convenient
 to use r_jk = p_jk/(1-q_j).
 
 	q_2 = 1/2	r_21 = 1/2
 	q_3 = 1/4	r_31 = r_32 = 1/2
 	q_4 = 1/8	p_31 = p_33 = 1/4	p_32 = 3/8
 	q_5 = 1/16	r_31 = r_34 = 1/6	r_32 = r_33 = 1/3
 
 Now the details.  We start out with the detailed expansion of exp(.5)/2,
 which is given by a_j = .5^(j+1)/j!.  However, for 6 or more, we will
 first multiply a_6 by 1.5.  Bit patterns will be allocated to the 
 outcomes so that the probability of returning j is a_j, and the remaining
 time a rejection occurs.  We find that easy by-products are, in some cases,
 random variables g4 which have a geometric distribution with p = 3/4, which
 happens to be the distribution of the number of places from the current 
 place which at which the bit mask will be set for j=3, and g16 for j=5.
 A geometric (1/2) random variable gets the location of a random bit.
 
 Let us consider the expansion, in hex and in binary
 
 0	.8			.1000 ....
 1	.4			.0100 ....
 2	.1			.0001 0000 ....
 3	.02aaaa....		      0010 1010 1010 1010 1010 ...
 4	.005555....			   0101 0101 0101 0101 ...
 5	.0008888 ...				1000 1000 1000 ...
 6*1.5	.0001111 ...				0001 0001 0001 ...
 
 remainder			.0010 1100 1111 0110 0110 0110 ...
 so we allocate the bits:
 				 01R2 RR3x RRRR 5RRy 5RRy 5RRy ...
 g16                                            1111 2222 3333 ...
 
 where R means reject, x means 3 or 4, y means 6 or more or reject.
 Now look at 3 and x.  In the x case, we get a geometric to find the
 bit in the binary representation.
 
 	10 1111 1111 1111 1111 ...
 	3  3434 3434 3434 3434 ...
 	1  2132 4354 6576 8798 ...
 
 In the y situation the probability is 1.5*a_6.  We now get a geometric
 random variable g, use it to both get a 2/3 - 1/3 division and also a g4.
 The 2/3 gives us the right value for a_6.  The 1/3 is half of that, or
 7*a_7.  We use a 2/7 acceptance to gat 2*a_7.  We now loop, starting with
 2*a_j
 
 	with probability 1/2, accept j		probability a_j of acceptance
 	with probability j/j+1, reject		probability 2*a_(j+1) of 
 						still going
  	augment j and proceed.
 
 I could similarly go into details about the procedures for forming or extending
 the mask when the active j is 5, 4, 3, or 2.  In some cases, byproducts of 
 previous stages are used.
-- 
Herman Rubin, Dept. of Statistics, Purdue Univ., West Lafayette IN47907-1399
Phone: (317)494-6054
hrubin@pop.stat.purdue.edu (Internet, bitnet)  
{purdue,pur-ee}!pop.stat!hrubin(UUCP)