From: hrubin@mentor.cc.purdue.edu (Herman Rubin)
Subject: Re: Goals, Cost, and Flexibility (was Re: "Training" of programmers)
Organization: Purdue University Statistics Department
Date: Thu, 1 Oct 1992 02:05:58 GMT

In article <1992Sep29.232852.22456@CS.ORST.EDU> crowl@jade.CS.ORST.EDU (Lawrence Crowl) writes:
>In article <Bv1G5B.AGn@mentor.cc.purdue.edu>
>hrubin@pop.stat.purdue.edu (Herman Rubin) writes:
>>In article <1992Sep22.212534.14185@CS.ORST.EDU>
>>crowl@jade.CS.ORST.EDU (Lawrence Crowl) writes:
>>>As near as I can tell, your buffer is only present to provide a cache of
>>>values and is not important to the correct operation of the algorithms.  As
>>>such, it is a separable issue and I will separate it.

>>Your assessment is correct.  While it is not important to the correct 
>>operation of the algorithms, it is extremely important to the running time,
>>especially on parallel processors.

>Independence of computations is more important to the execution time on
>parallel processors than is a fast buffer.  A shared buffer, no matter how
>fast, will eventually limit performance.

I said nothing about a shared buffer.  However, there are places where a
fast shared buffer can improve performance.  I hope that each processor
has enough memory for its own buffer.  If not, performance is likely to
be greatly reduced.

One of the essential problems with acceptance-rejection and related 
procedures is that vectorization, and especially SIMD operation, done
crudely, even with the best algorithms, can be deadly to speed.  Consider
the simple algorithm for generating random variables with density 6x(1-x);
it succeeds with probability 2/3 and fails with probability 1/3.  Now if
we are doing things right, all processors and trials are independent.
Suppose we use the crude scalar approach.  At each trial on each processor,
one would have to check whether the uniform buffer is not empty, and also
carry out the test.  This must be repeated until all processors attain a
successful value.  Now if there are 16,000 processors, this is likely to
take about 10 stages.  Doing the appropriate modifications can reduce the
time to about that of 2 stages per element produced, and if there is 
adequate spare buffer space, to less than 1.5.  

It is not what the values of the random variables are which is important,
but their distribution.  Different correct algorithms produce different
results, and a good choice of the algorithm definitely depends on the
architecture.

>>>What's more important is that your presentations of your algorithms are
>>>ill-defined.  There is a reason that programming languages are precise, the
>>>machines have no capacity to deal with imprecision.  You presentation assumes
>>>far too much about how the rest of us will interpret your examples, in
>>>particular, you assume that we'll "set everything up right and call stuff in
>>>the right way".  Such imprecision in your presentation leads us to believe you
>>>don't know what you're talking about.

The programming languages are far from precise.  CERTAIN situations are
handled moderately well, but unless exact integer arithmetic is used, all
such languages as Fortran and Pascal and C have ambiguities.

I believe that any mathematician, knowing what is involved in random number
generation, will understand my terminology.  A programmer, used to programming
terminology, may very well be able to copy the procedure, but will not have
the slightest idea how to modify it from one machine to another.  

>>This is one of the problems with the current HLLs, the inability to use
>>"mathematical" type notation.  I admit that I have not used these aspects
>>of C enough to have the syntax right.

>I argue just the opposite.  One of the _advantages_ of current HLLs is the
>inability to use "mathematical" type notation.  Mathematical notation is
>ambiguous and highly dependent on context.  Current HLLs don't have this
>problem.  Furthermore, we don't want to introduce ambiguity and bizzare
>context into HLLs.

Current HLLs have other ambiguities, like what is m*n, where m and n are
integers?  Sometimes different versions in the same series will give 
different answers.  

			[Much deleted.]

>>     extern fixed(0,1) uniform( ) ;
>>     extern int geometric( );

>>     fixed(0,1) nonuniform( ) {  
>>         fixed(0,1) X = uniform( ) ;
>>	M = geometric(); N = geometric();
>>	K = (M+1)/2;

>>Now I can write out this next step in C with a little work, but it will
>>be more easily understood if I just write it in words:

>>	Change the N+K-th bit of X to be the complement of the K-th bit;
>>	return X;
>>	}

>Is that the N+K-th bit counting from least significant or from most
>significant?  Your assertion that it would be more easily understood in words
>is incorrect because you've left a major ambiguity.  The C code probably
>wouldn't have had the ambiguity.

Counting is from the most significant.  Incidentally, the C code in the case
that all bits are used, here illustrated for 32 bit words, is subject to the
ambiguity, although it can be eliminated by two additional operations.

	B = 0x80000000 >> (K-1);
	H = (B&X) >> N;
	B >>= N;
	H = B &~ H;
	X = (X &~ B) | H;

The shifts must be logical shifts.  How many additional instructions do you
need in portable C to do this?

>>The use of a buffer cache seems to be the only way to avoid a call at each
>>trial.  

>Inline expansion of procedures is another way.  It has the advantage of
>maintaining the abstraction boundary, which is particularly important in
>large-scale programs.

There is no way that this can work here, where there is no certain knowledge
of how many items will be used or produced.  Randomness is randomness, and
attempts to get around them, which I have indicated, involve judicious use
of buffers, and expansions cannot work at all.  

>>The Crays have obvious problems.  So do SIMD machines.

>Argh!  Don't you listen?  What are those problems?  What are the costs?  What
>are the alternatives?  Are the problems due to a misapplication of the tools?
>(E.g. my washing machine makes a really lousy hammer.)

Look at the description of the problems.  We have only probabilities of how
many, and which, steps are successful, and which are failures, and there can
even be different types of failure.

>>If one has an acceptance-rejection procedure which may use up an exponential
>>random variable, and an acceptance-rejection procedure to generate the 
>>exponential random variable, one is doing a lot of pointer moving by not
>>using buffers compared to that by using buffers.

>Can you back this up with total pointer motions in comparative algorithms?
>Many people are often suprised at where the time goes in their programs.

The ones I have, using fancy buffering, greatly reduce pointer motions.
Scalar versions may be able to reduce them only a little more, but at the
cost of continual tests.  But with vectoring and especially with SIMD, the
cheapest way to handle things is, in many cases, to put the addresses of
the failures in an array, and then come back to handle those.  This does
avoid pointer motions, except for the filling stages.

However, there are places where judicious use of pointer motions can be
a great advantage.  Consider the computation of a function, where the
particular procedure to be used depends on the argument.  I suggest you
see about doing this on a machine with vectorization or parallelism without
juggling the vectors before and after.  Random number generation is more
of the same.

>In reguards to extracting the next set bit in a bit stream, you will _always_
>have to keep track of the appropriate bit and make limit checks.  The only
>question is whether the application is faster do such checks in software or
>doing them in hardware.  The RISC effort has shown that it is often faster to
>do it in software.

If you are doing the simple deterministic stuff, much of the time.  It there
is massive nondeterminacy, the best one can do in software is to use big
buffers, so that checks can be avoided much of the time.  Buffers again!

>>I have considered costs, and would not consider what you have suggested here.
>>The spacings between the 1's in a random bit stream have the geometric 
>>distribution required, and on the average only 2 bits are looked at to
>>get one geometric random variable.  That this can be computationally
>>costly is a matter of the lack of extremely simple hardware.  It is 
>>not comparable to sorting.

>You made a grossly oversimplifying assumption again.  Looking at only two bits
>in a random bit stream may be cheap, but it begs the cost of producing the
>random bit stream.  You aren't doing reasonable science or engineering if you
>measure the cost as the incremental overhead from some other subroutine.  I
>already showed you that doing this with sort turns the maximum problem into a
>simple index.

The random bit stream is produced as I have suggested above, using a procedure
in which all bits can be considered random.  On a machine with 32-bit words,
the cost of generating the bits stream runs 1/16 of a random word per geometric
random variable.  I would always produce the initial random numbers in such a
manner unless the machine prohibited it.  I am not currently aware of any which
prohibit it.

>>>I do not believe that generating a string of uniform random numbers and then
>>>processing them into nonuniform random numbers is more efficient than
>>>generating nonuniform numbers directly.  You're going to have to work to
>>>convince me otherwise.  Until you've convinced me that the two step process is
>>>faster, I'm likely to believe you're applying micro-optimizations to the wrong
>>>algorithm rather than finding the right algorithm.  The micro-optimizations
>>>you use for the wrong algorithm may or may not be suitable for the right
>>>algorithm.

>>The entire literature on generating nonuniform random numbers without starting
>>with uniform ones is minuscule, while there are numerous ways to generate 
>>them from uniform ones.  Most of the methods currently in use are acceptance-
>>rejection methods, like the first two algorithms above.  Most of the good
>>acceptance-rejection algorithms also use exponential input, and a few use
>>normal input.  The best ones I know for exponential use both uniform and
>>exponential buffers in a not too difficult manner.

>This sounds like a problem with the random number literature, not with the
>programming language or computer architecture literature.  Again, I don't
>believe that you'll make your best speed improvements without looking at the
>whole problem.

>In particular, I now suspect that the fastest nonuniform random number
>generator will not use the bit-twiddling that you advocate.  In particular,
>testing the high (or low) bit of the words coming out of your
>U[n-460]+U[n-607] uniform random number generator will probably be faster than
>keeping track of the distance between bits in each word.

Testing, on a vector or parallel machine, is very definitely not the way to
go.  On most vector machines (not the Cray 1), it is not that difficult to
get the spacings between the 1's with NO testing other than store on a 1 in
a control vector, and the count-testing for breaking things up into units
which fit in vector registers.  On parallel machines, the cost of any test
is usually at least the cost of completely doing both branches.

>It's time for you to put up some work.  Write up a C (or Fortran, or some
>other high-level language) subroutine, using whatever number representation
>you deem appropriate, that generates nonuniform random numbers.  I must be
>able to compile what you give me and run it, so make sure I have a compiler
>for whatever language you use.  Do not assume the existance of some other
>"helper" routine, though you may write one if you wish.  Once you've done
>that, then we can _start_ to discuss performance enhancement.  If you won't do
>that, then you aren't doing science and you aren't being helpful.

I might be able to get a graduate student to do this on a versatile semi-
portable "package" which I already have, but not programmed.  It is wall-
to-wall buffering, however.  I say semi-portable, because for very good
reasons, nothing which is not really complicated is fully portable; 
changing the size of the buffers changes the results.
-- 
Herman Rubin, Dept. of Statistics, Purdue Univ., West Lafayette IN47907-1399
Phone: (317)494-6054
hrubin@snap.stat.purdue.edu (Internet, bitnet)  
{purdue,pur-ee}!snap.stat!hrubin(UUCP)