To Foot"Random numbers should not be generated with
a method chosen at random".
Donald E Knuth,
The Art of Computer Programming Volume 2,
Seminumerical algorithms.
This page was written originally for Pascal; it now also considers Delphi. On Random in general, see also JavaScript Random, which, inter alia, includes "To Generate a Random Slide Show".
Note : I use Borland Pascal 7 and Delphi 2006, and have Delphi 3 and Turbo Pascal 5.
The Turbo/Borland Pascal and Delphi function Random generates pseudo-random sequences by updating a 32-bit variable RandSeed, which is normally initialised by Randomize from a system clock. The return value is then calculated from RandSeed :-
X := Random { 0.0 <= X < 1.0 } ;
J := Random(N) { 0 <= J < N } ;
Note that Random(N) has N possible results.
In Pascal, N is a 16-bit word; in 32-bit Delphi, it is a 32-bit integer.
For a given argument N, all possible results are nominally equi-probable. Other distributions can be derived.
The sequence of numbers obtained from Random(N) repeats, for any N, with a period of 232, RandSeed being a 32-bit quantity. Longer sequences could be obtained by appropriate, similar, non-complex code; but the constants involved would need expert selection.
Under some circumstances, Randomize may be partly systematic; and there is a system behind the updating of RandSeed which can occasionaly give detectable effects. True random sequences or initialisations can be generated by special hardware; recent Intel CPUs include a (slow) generator of such, intended for use as seeds.
In generating random numbers, use of Round will often be an error; Trunc is more likely to be right. But not always.
The implementation of Randomize and Random may vary between language versions. That of Random has three aspects : the updating of RandSeed, and the conversion of that to a float or to an integer.
A programmer can supersede the functions provided with better ones; therefore, it may be useful for code using these functions to be of higher quality than the native functions justify.
For both Pascal and Delphi, my program delprand.pas gives pure Pascal/Delphi apparently equivalent to but independent of the Borland library routines, which it also demonstrates. For serious work, such code should be tested over a complete cycle of 232 values of the seed RS.
In Pascal and Delphi, Randomize should normally be called once near the start of any program that uses Random or Random(). But, in testing, it may be well to initialise RandSeed with a fixed number, for reproducibility. Without Randomize, RandSeed probably starts at zero.
Randomize should NOT be called immediately before every call of Random.
In Pascal Version 7, Randomize initialises the longint RandSeed with a 32-bit value derived by permuting values (H, M, S, C) from the DOS time. Previously thought to be from permuting the date and the time; could this have been another version of Pascal? The time has 55 ms resolution; be sure to consider whether this is adequate. There are only $1800B0 ticks per day (and fewer per working day), so the initial RandSeed cannot take all possible values of its type, by a factor of $100000000/$1800B0 ~ $AAA = 2730.
In Delphi, RandSeed is initialised differently (dependent on version). D3 (& D4, D6 apparently) uses milliseconds since midnight GMT, which would imply that the initial RandSeed cannot take all possible values of its type by a factor of $100000000/86400000 = $31 = 49. However, at least in my D3, the number of milliseconds is always divisible by 10, and evidently rises at 18.2 Hz; it must be derived from the DOS timer at $40:$6C, as for Pascal 7 above. The factor is again ~2730.
John Herbster reports : D7 code looks like this :-
if QueryPerformanceCounter(Counter)
then RandSeed := Counter {overflow checking off here}
else RandSeed := GetTickCount; {rollover in 49.5 days}
Turbo Delphi 2006 has the same code. Note that the increment of GetTickCount may not always be unity; see my Delphi Programming, and Win32 Help. Counter increments at CPU frequency.
If RandSeed were to be initialised by the result of an RDTSC instruction [permuted] (Delphi Programming), all values could be obtained on a system that didn't get rebooted too often. RDTSC gives 64 bits; maybe XOR the halves, and BSWAP.
Here is a test Pascal program, to be run at a known time :-
begin Randomize ; writeln(RandSeed:10) end.
In a program scheduled to be run at a specific time each day, initialisation from time alone could lead to non-negligibly imperfect randomness of sequence start; to improve on this, one might XOR the initial RandSeed with GetTickCount, or with an absolute day count (* Trunc(now)? *), or, say, with ((Year-2000)*16+Month)*32+Day. A program auto-run at system startup must equally not generate RandSeed from GetTickCount alone.
Risks Digest 20.56 item #1 is a MUST for users of Randomize; on-line casinos using Randomize before every deal. Chortle.
Remember that Random(N) generates values in [0..N-1]; N is of type word (16-bit unsigned) in Pascal and non-negative integer (31-bit unsigned) in Delphi 3. I suspect that, with RandSeed being only 32-bit, calls such as Random(232×3/8) will not give near-equal probabilities.
Random generates a floating-point random number R, in the range 0.0≤R<1.0.
For 0.0<R<1.0, perhaps repeat R := Random until R>0 ; ?
To get a particular, repeatable, sequence of pseudo-random numbers, just assign a particular longint value to RandSeed instead of using Randomize (RandSeed does start with a fixed value; I believe always zero). The sequence may depend on which version of BP/TP/D is being used. This is particularly useful at the testing stage of a program.
I gather that the generator behind Random is satisfactory, but there are better ones. I've read that a Web search for '+Intel +"random numbers" +generator +true' finds interesting URLs.
For Random(N), the value of RandSeed (4 bytes) is multiplied by N, and the upper part of the result is the answer. I don't know that all versions of Pascal and Delphi do this, but it seems likely.
Since RandSeed takes 232 = 3×K+1 possible values, Random(3) cannot generate 0, 1, & 2 with exactly equal frequency; there is one more 0 in every 232 samples. Similar applies whenever the argument of Random is not a power of 2. The point will only rarely be significant.
The case of N=3 has a simple solution, I think. At design
time, determine the RandSeed corresponding to one of those
zeroes; and at run time, if that RandSeed appears, try again :-
repeat Selection := Random(3) until RandSeed <>
TryAgain ;
For Random(N) with N greater than three, there may be more numbers to be avoided.
Another way of generating randoms uses a Hardware Shift-register approach, and at best generates 2N-1 values. Note, 232-1 = (216+1)×(216-1) and so must be a multiple of three; in fact it is (216+1)×(28+1)×(24+1)×(22+1)×(21+1)×(21-1) and so has other factors. Non-maximal sequences have different factors IIRC.
RAH Prins wrote that, in TP6, X := Random can generate 1.0. This has now been seen by me in TP5; starting with RandSeed not initialised, at exactly half-way through the sequence, a call of Random with RandSeed=-1498392781 gives 1.0 and sets RandSeed to -2147483648 = $80000000. This does not occur in BP7 or Delphi 3.
Duncan Murdoch wrote, in TP6BUGS7.LST, "The Random(N) function doesn't produce a uniform distribution unless N is a power of 2, especially for large values of N. (Fixed in BP 7.)".
Use of mod is sometimes suggested as a way of reducing a uniformly-distributed random number R to lie within the range (0..N-1). Unless the range of R is a + (0..bN-1) for integer a, b, the new R will not be evenly distributed, for a fairly obvious reason.
When Random[(N)] is called, RandSeed is updated first.
I have read that :- Random implementation in Delphi (and Borland Pascal) is as follows: Pseudorandom longint sequence is implemented as a congruential generator with the formula X[n+1] = 134775813*X[n] + 1 (mod 2^32); my tests support this, for BP7 & D3. It can be shown that the sequence has full period (its length is 2^32); the multiplier in Hex is $08088405. D5 and D7 are reported as using the same expression. So does TD2006.
And :-The Lehmer algorithm is : X := ((X * A) + C) modulo M
Tests also confirm expectation that the result of Random(N) is as if obtained by dividing RandSeed by 2^32, then multiplying by N and truncating to a whole number. BP7 and D3 give the same result.
Tests in BP7 and D3 show that the result of Random is as if obtained by dividing RandSeed by 2^32, except that in BP7 the upper and lower halves of the range appear swapped. This is probably somehow due to BP7 having longint but not dword, or being 16-bit rather than 32-bit.
I have read that the library code, in Pascal, is :-
NextRand PROC NEAR
imul eax, RandSeed, 08088405h
inc eax
mov RandSeed, eax
ret
NextRand ENDP
I have read that the library code, in Delphi 7, is equivalent to :-
function _RandInt(Range : LongInt) : LongInt ; var NewSeed : LongInt ; begin NewSeed := RandSeed * $08088405 + 1 ; RandSeed := NewSeed ; Result := Int64(EDX) * Int64(Range) shr 32 ; end ;
The following is the JavaScript FAQ way to get a Random(N) from a Random :-
function Random(N) { return Math.floor(N*Math.random()) }
and presumably could be used by a Delphi programmer.
And I have heard that X[n+1] = 1664525*X[n] + 1013904223 is good.
What constant[s] would be appropriate for a 48-bit or 64-bit version?
2002-10-26: I read that Knuth reports the multiplier
6364136223846793005 (Hex 5851F42D4C957F2D) as excellent
(but not necessarily best) for 64-bit :-
try X[n+1] = 6364136223846793005*X[n] + 1 (mod 2^64).
Program longcalc.pas could be scripted to generate that sequence; or ones with larger multipliers.
From mathutys.pas, in and tested by mathunit.zip :-
// 2008-03-18 Turbo Delphi 2006, OK
uses Windows ;
{$OVERFLOWCHECKS ON} {$RANGECHECKS ON}
var RandSeed64 : UInt64 ;
const
Factor = 6364136223846793005 { Said to be OK by Knuth } ;
HalfTo64 = 1.0 / (65536.0 * 65536.0 * 65536.0 * 65536.0) ;
procedure Randomise64 ;
begin QueryPerformanceCounter(Int64(RandSeed64)) end {Randomise64} ;
function Random64 : extended ; overload ;
begin
{$OVERFLOWCHECKS OFF} {$RANGECHECKS OFF}
RandSeed64 := RandSeed64*Factor + 1 ;
{$RANGECHECKS ON} {$OVERFLOWCHECKS ON}
Result := RandSeed64*HalfTo64 ;
end {Random64} ;
function Random64(const UI64 : UInt64) : UInt64 ; overload ;
begin Result := Trunc((Random64)*UI64) end {Random64} ;
But see, at CodeGear, Report #6619 - Status: Open Expand the Random function.
In articles by Sven Pran I read that
X[n+1] := (a*X[n] + c) mod M
is reversed by
X[n-1] := b*(X[n] - c) mod M
where (a * b) mod M = 1; note the parenthesis positioning.
The value of c can be found by setting Randseed:=0 then calling Random and reading RandSeed; the value of a is then obvious from the next RandSeed. The value of b can be found by looping through the full range and considering the value previous to Randseed=0.
The part of the BP7 RandSeed sequence around 0 is :-
-19094774, 649090867, 0, 1, 134775814, -596792289.
Using the forward sequence 0, 1, 134775814, c = 1 and then a = 134775813.
Using the reverse sequence 1, 0, 649090867, c = 1 and then b = -649090867.
So actually b is either -649090867 or +3645876429 (which are equal, mod 232) in Borland Pascal 7 and Delphi 3; those satisfy (a * b) mod M = 1.
With the Pascal/Delphi 32-bit PRNG, one can expect that the Extended value of Random will take values which are multiples of 2-32, and that those will store exactly in a Double. They will not store exactly in a Single. To get a good Single distribution, the conversion must truncate. If the values are rounded by the conversion, the Single may be set to 1.0. But Random(N) need not be affected.
If a 64-bit PRNG is used, similar considerations apply to Doubles. See also Generator Properties in JavaScript Random.
Of course, such refinement will not matter in most applications; but with the possibility of leaving a program running on a 3 GHz PC for a few days, during which RandSeed will cycle many times, one should not be too ready to disregard the possibility of an end result being affected.
{ For truncation towards zero of a longer float to a shorter, one should
only need to cast the longer float to an array of integer or word, mask
off the superfluous bits until it has a value of the original type that
the shorter float can hold, then assign. Something like :- }
function XtndToDble(X : extended) : double ;
type A5 = array [0..4] of word ;
begin A5(X)[4] := A5(X)[4] and $F800 ;
Result := X end ;
function DbleToSngl(X : double) : single ;
type A5 = array [0..1] of longint ;
begin A5(X)[1] := A5(X)[1] and $E0000000 ;
Result := X end ;
function XtndToSngl(const X : extended) : single ;
begin Result := DbleToSngl(XtndToDble(X)) end ; // or mask 40 bits
// UNTESTED. Independent of the FPU Control Word.
If a random position, to the nearest centimetre, for a particle in a ten-centimetre box is needed, Round(10*Random) should do; it gives 1..9 each with 10% probability, and 0 & 10 each with 5% probability.
The native Random function is intended to produce a flat distribution. It will usually be important that any extension has the same property.
The range of numbers required may well be wider than the 0..65534 range of TP/BP Random, which, alas, cannot even fully cover a 16-bit type.
Note that the available randomness in Random is only 32 bits; therefore when larger numbers are generated using it alone, these can have at most 232 possible values.
A good algorithm for generating large-range random numbers will be able also to generate small-range ones without fundamental change; any imperfections in the latter should be more obvious, and are liable to occur also in the former.
In TP/BP, Random(N:word) returns a random value in 0..N-1. This is IMHO a design fault, since 0..N would be more congenial. Therefore it cannot generate a random word, as 65535 is not possible (Random(0) also generates only zero values).
Use, therefore, the high word of the longint RandSeed, after calling Random.
The following program needs testing, but CR may return a random of type comp :-
function CR : comp ; var A : array [0..7] of byte ; B : byte ; begin for B := 0 to 7 do A[B] := Random(256) ; CR := comp(A) end ; var K : word ; begin for K := 1 to 20 do Writeln(CR) end.
Note, however, that $8000000000000000 is perhaps a NaN in Pascal; the stated range of type comp is every integer value from -(263-1) to +(263-1) inclusive. If this value, -263, appears, try again.
A random number of any length can be generated by using the above idea to generate a random of range 0..2N-1, where N is as small as possible, and repeatedly trying again if the number found is too big. Worst case, the possible repeating should on average less than double the time taken.
To cover a type of range, say, 0..Pred(2^24) will IMHO need something like
Ran := (longint((Random(256) shl 8) + Random(256)) shl 8) + Random(256) ;
or for the full longint range, as in this test program :-
type X = array [0..3] of byte ;
var Ran : longint ; K : byte ;
const XMin = -500000000 ; XMax = 1000000000 ;
begin repeat
repeat for K := 0 to 3 do X(Ran)[K] := Random(256)
until (Ran>=Xmin) and (Ran<=Xmax) ;
Write(Ran:14) ; Readln ; until false end.
One approach should be to use something like Trunc(N*Random), since Random returns 0.0≤float<1.0; do not overestimate the quality of the result. Simple scaling up of a random integer will give missing values.
Another should be to build up by concatenation of bytes & bits, as Random(N) produces random bit strings for N a power of 2 in [2..32768]; (Random(16384) shl 12)+Random(4096) should give a random 26-bit number (range 0..16384×4096). For a range not a power of two, use a "repeat until in-range" loop. Untested.
In TP7/BP7, I believe that RandSeed, declared as longint, takes all 2^32 possible values as one successively calls Random. It is, therefore, a better start than Random is; but beware of detectable correlation between successive values (ISTM that RandSeed is alternately odd and even); just use the high word of it).
With its limited number of internal state bits, native Random is entirely unsuited to generating globally unique identifiers; and Randomize makes it worse.
Native Random is essentially limited by having only a 32-bit internal state. Much better generators exist, and can be found either via Links or by a Web search, including versions for Delphi. In fact, for any serious work, the native 32-bit Lehmer Random should not be used.
A validated large-random-number function would be a useful addition to a library.
A different form of generator can be expressed as
XK+1 := XK shl 1 + Ord( (Q0 in BS(X)) xor ... xor (Qn in BS(X)) )
where BS casts the Z-bit variable X to a set of 0..Z-1, and [Q0..Qn] is a strategically chosen constant set of type BS. Clearly this sequence has length less than 2Z - because if XK=0, then XK+1:=0 - but it can have length 2Z-1.
Note that if Z is such as 8, 16, 24, 32 then several factors of 2Z-1 are easily identified.
There may be a better way of computing that Ord() if n is not small; but good sequences can be obtained with n=1. To be sensible, Q0=Z-1.
There is an obvious relationship between consecutive bit-patterns in the variable. This will vanish if the variable is read only after Z clocks. In practice, Z and 2Z-1 will have no common factors, so the sequence length will not be reduced thereby.
This form can be effectively implemented in simple hardware, but is not favoured for software.
See also in my JavaScript Random.
Delphi implementations should exist for all good methods.
A thread "Delphi V7 PRNG?" in newsgroup a.c.l.b-d in early November 2008 includes a Delphi version of KISS32.
My definitions : A Deal generates the items that it delivers. A Shuffle works with an existing collection; it can be emulated by using a Deal to index the collection. A Draw is a partial Deal, perhaps indexing a collection. My definitions.
To shuffle, one can just stir a set of items; but this either does not introduce essentially full randomness, or takes longer than necessary.
A perfect deal or shuffle introduces complete randomness; no trace of the original order is left. For N items, there is one chance in N! of ending up with the original order.
Up to the date of writing this paragraph, I have given no thought whatsoever to the shuffling, dealing, or drawing of linked lists.
Shuffling rearranges a pre-existing sequence.
To shuffle a pack of cards, or any similar problem, the correct method is to select the elements randomly, one at a time, from those elements as yet unselected. So, in an efficient manner using swapping in a single array A[1..Max] [I'm told : Donald E Knuth, The Art of Computer Programming (Chapter 3.4.2), and he also gives the original publishers of this algorithm: L.E.Moses and R.V.Oakford (1963) and R Durstenfeld (1964)] [also Fisher-Yates] [Wikipedia] :-
for J := Max downto 2 do Swap(A[J], A[Succ(Random(J))]) ;
By selecting in reverse order, the indexing of the random element is eased. I'm assuming that A is [1..Max]; using downto 2 saves a little time. Depending on details (such as range-checking) it may be better to "inline" Swap, and (if Max is not big) to prevent swapping an element with itself.
for J := 52 downto 2 do begin K := Succ(Random(J)) ;
if J<>K then begin p := @A[J] ; q := @A[K] ;
t := p^ ; p^ := q^ ; q^ := t end end ;
From "Subject: Cryptography FAQ (08/10: Technical Miscellany) ... Date: 26 Jun 1997 08:48:11 GMT", available in news.answers : "8.13. How do I shuffle cards?" :-
Card shuffling is a special case of the permutation of an array of values, using a random or pseudo-random function. All possible output permutations of this process should be equally likely. To do this, you need a random function (modran(x)) which will produce a uniformly distributed random integer in the interval [0..x-1]. Given that function, you can shuffle with the following [C] code: (assuming ARRLTH is the length of array arr[] and swap() interchanges values at the two addresses given)
for ( n=ARRLTH-1; n>0 ; n-- ) swap(&arr[modran(n+1)], &arr[n]) ;
I think that is equivalent <g>.
Note that, if an element of A consists of more than a few bytes, in any substantial reordering operation (randomisation, sorting, ...) it will be better to manipulate a list or array of pointer-to-element rather than a list or array of element itself.
After Gerald Brandt: A common error is to do :-
for J := 1 to Max do Swap(A[J], A[Succ(Random(Max))]) ;
As this has Max random selections each with Max possibilities, there are Max^Max equi-probable possibilities in doing this. However, there Max! different permutations of A, and (for Max>2) Max! is not a factor of Max^Max. Therefore it is impossible that every permutation is generated by the same number of random sequences. This condition is necessary and sufficient for the uniform distribution of the permutations. The recommended code does give Max! sequences, and (being also correct) so is maximally efficient.
I have a program RANDCARD.PAS which demonstrates good and bad shuffling.
To rearrange with no element left in its original place, select the elements, one at a time, at random from those other elements as yet unselected (tested) :-
for J := Max downto 2 do Swap(A[J], A[Succ(Random(Pred(J)))]) ;
Dealing generates a new random sequence, without using an old one. Since the result depends only on the initial RandSeed the method can generate only 232 sequences. To get more possible sequences, shuffle an already-disordered array.
To deal a pack of numbered cards, not already in an array, is slightly quicker by the method given by Horst Kraemer at the end of TSFAQP#140 :-
A[1] := 1 ; for J := 2 to 52 do begin K := Random(J)+1 ; A[J] := A[K] ; A[K] := J end ; = for J := 1 to 52 do begin K := Random(J)+1 ; A[J] := A[K] ; A[K] := J end ;
One can understand that by seeing that every value of J is assigned a random position in 1..J, and the existing occupant of that random position is placed at the new end before J is emplaced; thus, given that the original set of 0 or 1 element was properly random, each increased one will also be so. It may be worth modifying the above to use pointers, as in Shuffle.
New 2004-09-25, untested :-
for J := 1 to 52 do begin K := Random(J)+1 ; A[J] := A[K] ; A[K] := Deck[J] end ;
should fill A with a randomised copy of array Deck.
That is a complete deal. For a partial deal (also see "Drawing" below), it can be better to choose at random, and reject if in the set already chosen (untested) :-
SetGot := [] ; for j := 1 to Few do begin repeat k := Random(52)+1 until not (k in SetGot) ; Include(SetGot, k) ; a[j] := k end ;
Other methods of obtaining a non-repeating sequence exist; I'm told that Knuth, volume 2 chapter 3 refers (and uses: for X := 1 to M do pseudo := (x*A+C) mod M;}.
It's probably best here to use a one dimensional array of numbers 0..51, deal these, and interpret N mod 4 as the suit and N div 4 as the number. Suit can be an enumerated type, (Club, Diamond, Heart, Spade); in that case, use Suit(N mod 4).
Drawing produces a set containing a few members of an existing, generally complete, set. There is no need to use a full Deal; there is no significance in the order of the result.
David Rifkind has presented (1999-05-17, news:c.l.p.d.m) :-
type RandomSet = set of 1..MaxN ;
{ Generate M non-repeating random numbers from 1 to N }
procedure GenerateRandomSet(M, N : integer ; var aSet : RandomSet) ;
var J : integer ;
begin
if M < 1 then aSet := []
else begin
GenerateRandomSet(M - 1, N - 1, aSet) ;
J := Random(N) + 1 { J is in [1..N] } ;
if J in aSet then Include(aSet, N) else Include(aSet, J) ;
end ;
end {GenerateRandomSet} ;
It seems to work, & DR is trustworthy.
NOTE that it is easy enough to handle more than 256 elements, by indexing into an array[0..x] of set of 0..255.
DRAFT Proof by induction :-
Sets with more than 256 elements can be treated, fairly efficiently, as arrays of sets.
2001-11-21 : Rufus V Smith (news:c.l.p.m) has a non-recursive version :-
aSet := [] ; for K := N-M+1 to N do begin J := Random(K)+1 ; if J in aSet then aSet := aSet + [K] else aSet := aSet + [J] ; end ;
which my program dealtest.pas shows to give the same sets in the same order.
In translating it into JavaScript - formally, JS does not have Sets, but an array with elements that may or may not be defined can be used for a Set - it appears that the long line can be replaced by
if J in aSet then J := K ; aSet := aSet + [J] ;
which seems clearer and possibly gives better code.
In Delphi, one can generate numbers at random and insert them into a TStringList Sorted dupIgnore until Count is adequate; untested. The sort will be alphabetical, but one can add a large number to each entry.
Provided that the Random Number Generator is perfect, the above recommended methods give all possible sequences with equal probability. Actually, any generator of the general type used has (at least) three faults :-
Thus only for a very limited number of items can all orders actually be obtained using a 32-bit generator seeded from the DOS time-of-day.
It seems highly likely that repeatedly shuffling the same deck (provided that the use of the deck does not sort it) will generate many more possible sequences; and probable that it can generate them all with virtually equal probability. The system now has 232 × 52! (about 3.5E77) internal states.
Always check a generator for plausibility, even if trivial; in TP/BP, remember that Random(N) takes values in [0..N-1], and remember to Randomize once.
AFAIR, I have put into SWAG an algorithm for generating random numbers with any specified distribution curve (the example is Gaussian) - it should be in MATH, with a number around 100 or so (i.e. search backwards from the end) - seen at MATH0104. Or see distribs.pas.
The idea is to draw the distribution within a box, select a point at random within the box, and accept it's x-value if the point lies under the curve, else repeat.
The distribution is, of course, truncated; but in the real world many Gaussians are, at a few sigma; anyway, it makes little difference.
From a working program :-
function RandomGaussianSD1 : DatT { Gauss dist. SD=1, cutoff @±4.0 } ;
var X : DatT ;
begin repeat X := (Random-0.5)*8.0 until Exp(-Sqr(X)*0.5)>Random ;
RandomGaussianSD1 := X end {RandomGaussianSD1} ;
{ "Sqr" is "Square"; replace Exp(-Sqr(X)*0.5) for other distributions }
procedure GaussTest ;
var j : byte ; s, t : DatT ; const k = 80 ;
begin t := 0.0 ;
for j := 1 to k do begin
s := RandomGaussianSD1 ; t := t+Sqr(s) ;
Write(s:9:3) ; if (j and 7)=0 then Writeln ;
end ;
Writeln('RMS: ', Sqrt(t/k):10:3, ' _1? '^G) ; Delay(500) ;
Writeln('Does that look Gaussian, mean 0, RMS 1 ?'^G) end {GaussTest} ;
Imagine a box bounded by x=-4, x=+4, y=0, y=1, with a unit Gaussian drawn within it. Repeatedly choose a point x, y within it at random, until the point lies under the Gaussian; the x-value is a random from the required distribution. (A direct method for Gaussians is below.)
The scheme may be adapted for selecting an item from a weighted collection; ... until Weight > Random(MaxWeight) ; for discrete weights. Untested pseudo-code :-
type TWeights = array [0..NW-1] of word ;
const Weights : TW = ( ... ) ;
function RandomSelect(const Wts : TWeights) : word ;
var X : word;
begin repeat J := Random(NW) until Wts[J]>Random*MaxWeight ;
RandomSelect := J end {RandomSelect} ;
The "box" method will give a distribution which is statistically correct, but is not guaranteed to be any more than near to the nominal distribution.
To get an exact distribution, with a specified number of each element, it is sufficient, and may be necessary or optimum, to assemble the required number of each in an array, and then to give one full, efficient shuffle as above.
If a distribution, continuous or discrete, has very disparate weights, the above method becomes slow. Just divide it into a set of distributions, where the sets have non-disparate weights; select a set at random, as above, using its total weight; then select from the chosen set, with its lower maximum weight.
For a Gaussian distribution, Jyrki Lahtonen <lahtonen@utu.fi> has recommended (I have condensed it a little) :-
function RandomGaussian : real { mean=0, SD=1 } ;
var Phi, Rho, X, Y : real ;
begin
Phi := 2.0*Pi*Random ;
Rho := sqrt(-2.0*Ln(1.0-Random)) { uses 1.0-random, as Ln(0.0) fails } ;
X := cos(Phi)*Rho ; Y := sin(Phi)*Rho ;
{ A polar to cartesian coordinate transformation is taking place -
observe that here both X and Y will have normal Gaussian distribution }
RandomGaussian := X { or Y } end ;
This is fundamentally a better method for a Gaussian (I do not doubt, but have not verified, the theory behind it).
MRM wrote, 2001-07-24 : This is a standard Box-Muller transformation (1958). ... Eric Weisstein has a nice proof.
Bob Schor <bschor@vms.cis.pitt.edu> has a mathematically equivalent routine, adapted to deliver X and Y alternately, while only calling one pair of Randoms for each X, Y pair. He says that the proof involves integral(exp(-sqr(X))) which is equal to sqrt(integral(exp(-(sqr(X)+sqr(Y))))) and leads to an integral over a unit circle, which is known to be Pi.
Delphi, at least in version 3, has function RandG, maybe :-
{ Marsaglia-Bray algorithm: }
function RandG(Mean, StdDev : extended) : extended ;
var U1, S2 : extended ;
begin
repeat
U1 := 2*Random - 1 ;
S2 := Sqr(U1) + Sqr(2*Random-1) ;
until S2 < 1 ;
Result := Sqrt(-2*Ln(S2)/S2) * U1 * StdDev + Mean ;
end ;
JH has posted :-
// Very quick and dirty :
x := +random-random+random-random+random-random
+random-random+random-random+random-random ;
Adapted from an unverified (by me) post by Gary Williams on 2002-05-01 in news:borland.public.delphi.objectpascal :-
function WeightedRandom(const Weights : array of Double) : Integer ;
var R : Double ;
begin
R := Random * Sum(Weights) ;
for Result := 0 to High(Weights) do
if (R < Weights[Result]) then BREAK
else R := R - Weights[Result] ;
end ;
It takes an array of weights, and chooses an index into it. Check that BREAK is always executed, even with rounding errors.