From: hrubin@pop.stat.purdue.edu (Herman Rubin)
Newsgroups: comp.arch,comp.lang.misc
Subject: Re: Hardware Support for Numeric Algorithms
Date: 16 Oct 92 13:52:34 GMT
Organization: Purdue University Statistics Department

In article <1992Oct12.201841.22506@CS.ORST.EDU> crowl@jade.CS.ORST.EDU (Lawrence Crowl) writes:
>I apologize for this disjointed message, but I had a work crunch and
>missed some articles.  Herman Rubin and I were discussing the hardware
>support needed to execute his numeric algorithms efficiently.  My last
>message asked him to provide a complete, working subroutine for one of
>his non-uniform random number generators.  In the last message I saw,
>he said he might look at getting one of his graduate students to look
>into it.  Have I missed anything since then?

>At any rate, I don't find that answer satisfactory.  I'm not asking for
>lots of code, only one page.  It doesn't have to be "efficient", it
>just has to be correct.  Once we have working code, we can find ways to
>make it efficient given current hardware, and estimate the cost and
>benefit of new hardware.  So, Dr. Rubin, are you prepared to provide
>the subroutine?

Here are two routines.  The first is a "bootstrap" type routine which
can be used to produce results with the same distribution as the second.
The second uses less bits, and can be improved considerably from the
efficiency standpoint.  Note the significant use of single random bits
and spacings between the 1's in a random bit stream in both.

In the first one, I have arrays of sizes to be determined at run time.
This is not as great a problem as it seems, and can be handled if 
there are reasonable bounds.   This algorithm is vectorizable.    

In the second one, I did not expand some code for certain tests.  This
is a rare part of the program, but must be there.  Optimizing this 
algorithm requires a much more complicated version of the first part.

/* This procedure essentially starts out with an array of increasing   *
*  indices and an equal sized array of masks and produces a usually    *
*  larger array of the same type.  The array sizes are not fixed, so   *
*  the lengths are also variables.  There is also a "startup" number.  *
*  In this, b will be an unsigned integer which is 1 just to the left  *
*  of the binary point for the masks (this can be easily worked around *
*  if this is not possible), G() gets a geometric (.5) random variable *
*  at each call, which can be generated as the spacings between 1's in *
*  a random bit stream, and BIT() gets a random bit.  All variables    *
*  are integers.  For simplicity, all arrays start with element #1.    */

narr(lold, nold, mold, rold, lnew, nnew, mnew, rnew) 
{h = nold[lold];
while(i <= h) ntemp[i] = rold + (position of the i-th 1 in a random bit stream);
while(i <= h) mtemp[i] = 0;
k = ntemp[h];
while(i <= lold){
	j = nold[i];
	if (BIT()==1) ntemp[j] = -1;	#	any impossible value
	else
		c = b >> G();
		if (c&mold[i] != 0) ntemp[j] = -1;
	else
		mtemp[j] = (mold[i] &~ (c-1)) | c;
	}
lnew = 0;
while(i <= k) {if(ntemp[i] < 0) break;	# test for the same impossible value
	lnew += 1;
	nnew[lnew] = ntemp[i];
	mnew[lnew] = mtemp[i];}
rnew = k - nnew[lnew];
return;
}



/* This procedure returns an integer n and a mask m, such that  *
*	n + (x &~ m) is exponential with mean 2, when x is a    *
*  fixed-point uniform (0,1) random variable.  		        *
*  B is an unsigned integer, which is 1 just to the left of the *
*  binary point; G() returns geometric (.5) random variables,   *
*  which can be generated as the spacings between 1's in a      *
*  random bit stream, BIT() is a random bit, and CT(i) returns  *
*  the number of 1's in i random bits.	Where events with other *
*  specified probabilites are needed, the method used will not  *
*  be described.						*
*  All variables in this program are integers.	All breaks are  *
*  to the statement "a=G()".					*/

n=0; m=0; b=B;

{	a = G();
	if(a==1){n++; break;}
	if(a==2)goto case1;
	if(a==3)break;
	if(a==4)goto case2;
	if(a==5)break;
	if(a==6)break;
	if(a==7){g4=1; goto case3;}
	if(a==8){c=G(); g4=c>>1;
		if(EVEN(c))goto case4;
		else {g4+=2; goto case3;};}
	if(a==9)break;
	if(a==10)break;
	if(a==11)break;
	if(a==12)break;
	g16 = (a-9)>>2; a &= 3;
	if(a==1) goto case5;
	if(a==2)break;
	if(a==3)break;
	i = 6;
	c = G();
	g4 = (c+1)>>1;
	if(ODD(c))goto caseh;
	with probability 5/7 break;
	for(i=7;;i++){
		if(BIT()==1)goto caseh;
		with probability i/(i+1)break;}
}
		
caseh:	b >> 1;
	j = CT(i);
	if(j==i || j==0) goto caseh;
	i=j;
	m |= b;
	if(i >= 6)goto caseh;
	goto case^i;		 	# ^ is here used for concatenation

case5:	b >> g16;
	m |= b;
	c = G();
	if(BIT()==1){if(ODD(c)){goto case2; goto case1;}
	g4 = (c+1)>>1;
	if(ODD(c))goto case3;
case4:	c = G();
	if(c > g4){m |= (b>>c | b>>g4); goto case1;}
	m |= b>>c;
	if(BIT()==1)goto case1;
	b >>= 1;
case3:  b >>= g4;
	m |= b;
	b <<= 1;
case2:  m |= (b>>G());
case1:  exit;
	
-- 
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)
Newsgroups: comp.arch