QualityCentral

Public Report

Report From: Delphi-BCB/RTL/Delphi/Math Unit [ Add a report in this area ]

Report #: 6619

Status: Open

Expand the Random function

Project: Delphi

Build #: 4.453

Version: 7.0

Submitted By: John Herbster

Report Type: Suggestion / Enhancement Request

Date Reported: 12/12/2003 8:23:11 AM

Severity: Commonly encountered problem

Last Updated: 5/26/2008 10:47:38 AM

Platform: All platforms

Internal Tracking #: 225315

Resolution: None (Resolution Comments)

Resolved in Build: : None

Duplicate of: None

Voting and Rating

Overall Rating: (4 Total Ratings)

4.75 out of 5

Total Votes: 3

Description

Machines are getting so fast, I wonder if it is time to slow down the Random function by making it use 64-bit seeds and results. Maybe call it Random64( or 63 or 62 if we need to waste a couple of bits)?

Programmers are starting to be concerned about the limited randomness that can be gotten from a 32-bit seed.

Please review the recommendations and code provided as comments of Sven Pran in this QC report's comments.

Also see these reports:
  #4368 -- Making the random seed thread-local
  #4394 -- Calling Randomize at program start.
  #4493 -- Providing an overloaded call with the seed as a parameter.

** The following are suggestions by John Stockton, Surrey, UK posted with permission **

The Random and Randomize functions should not be changed, for compatibility reasons.

A new Unit should instead be provided, containing thoughtfully- named and -written routines including versions of Random, Deal, and Shuffle (and, if there is not one already, a Sorting Unit).  The task could be subcontracted to skilled volunteer labour, such as is found here.

The documentation for the existing routines should be changed to say that they are OK for most normal uses but that users should look at the new stuff before deciding.

New randomisers should be capable of initialising their seeds to any value within range, at least if called a small but arbitrary time after the start of a program which was started at an arbitrary time of day on a Windows booted an arbitrary time ago.

If a reasonably efficient reverse algorithm is known, it should be provided; it is necessarily an alternative, and reversibility might be of use.

The new routines should not use a specific RandSeed; instead, each should have a first parameter   var seed : <type> .

For completeness, the unit should contain the old algorithms but with the new parameter.

Steps to Reproduce:

Please look in the Comments tab for Discussion and Implementations (by Sven Pran, Dr. John Stockton, and myself (Herbster))

Workarounds

None

Attachment

None

Comments

John Herbster at 12/12/2003 8:29:13 AM -

=== Googling finds these previous discussions of PRNGs ===

From about a year ago in objectpascal forum:
http://tinyurl.com/cf2y short for
http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&threadm=3C62D7FA.10906%40streamsec.se&rnum=5&prev=/groups%3Fhl%3Den%26lr%3D%26ie%3DUTF-8%26oe%3DUTF-8%26q%3D%2B32-bit%2B%2BPRNG%2BOR%2BRNG%2BOR%2Brandom%2B%252263%2Bbit%2B%2522%2B-cypher%2B-penknife%2B-encryption%26sa%3DN%26tab%3Dwg

Good discussion of 64-bit PRNGs:
http://tinyurl.com/cf3p short for
http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&threadm=%25_P99.180496%24SS.7379596%40bin3.nnrp.aus1.giganews.com&rnum=1&prev=/groups%3Fq%3DPRNG%2BOR%2BRNG%2BOR%2Brandom%2B%252264%2Bbit%2B%2522%2B-cypher%2B-penknife%2B-encryption%26hl%3Den%26lr%3D%26ie%3DUTF-8%26oe%3DUTF-8%26selm%3D%2525_P99.180496%2524SS.7379596%2540bin3.nnrp.aus1.giganews.com%26rnum%3D1

From the above thread in a post by George Marsaglia (geo@stat.fsu.edu) on Aug. 24, 2002:
/*
A C version of a very very good 64-bit RNG is given below.
You should be able to adapt it to your particular needs.

It is based on the complimentary-multiple-with-carry sequence
         x(n)=a*x(n-4)+carry mod 2^64-1,
which works as follows:
Assume a certain multiplier 'a' and a base 'b'.
Given a current x value and a current carry 'c',
form:               t=a*x+c
Then the new carry is     c=floor(t/b)
and the new x value is    x = b-1-(t mod b).

Ordinarily, for 32-bit mwc or cmwc sequences, the value t=a*x+c can be formed in 64 bits, then the new c is the top and the new x the bottom 32 bits (with a little fiddling when b=2^32-1 and cmwc rather than mwc.) ...
*/

A simplier faster method by Marsaglia on Aug 25, 2002:
http://tinyurl.com/cf4f short for
http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&threadm=Jc9a9.203434%242p2.8301471%40bin4.nnrp.aus1.giganews.com&rnum=4&prev=/groups%3Fq%3DPRNG%2BOR%2BRNG%2BOR%2Brandom%2B%252264%2Bbit%2B%2522%2B-cypher%2B-penknife%2B-encryption%26hl%3Den%26lr%3D%26ie%3DUTF-8%26oe%3DUTF-8%26selm%3DJc9a9.203434%25242p2.8301471%2540bin4.nnrp.aus1.giganews.com%26rnum%3D4

From above thread Marsaglia writes:
"In response to a posting of Aug 24, in which I suggested a very good method for producing 64-bit random integers, I was asked if I had any simpler and faster methods. The following is simpler, faster
and seems to provide 64-bit integers that will satisfy all but extreme tests of independence among successive values.

It is an extension to 64-bits of the widely used 32-bit Super Duper RNG that I developed in the 60's, after showing that (congruential)  "random numbers fall mainly in the planes".

The 64-bit version of SuperDuper, say supdup64( ) combines (by addition mod 2^64)  the sequence
  x(n)=6906969069*x(n-1)+1234567 mod 2^64
with a 3-shift shift register generator, (in C:)
      y^=(y<<13);  y^=(y>>17);  y^(y<<43);
The period of supdup64( ) is 2^128-2^64, and it is very fast, about 30 nanosecs per number, even on machines that have to carry out 64-bit integer arithmetic with paired 32-bit registers, (as most currently do).

Other than a=6906969069 you may choose any multiplier that is 3 or 5 mod 8, and for shifts other than 13 left,17 right,43 left,  you may choose any of some 2200 triples that will provide y's with period 2^64-1.
The full set is in a forthcoming paper; you may send for a preprint.

Thus a selection of one from some 2200*2^62 versions of supdup64( ) could provide an excellent 64-bit RNG for potential security purposes, and you could double that number by choosing exclusive-or for combining  x and y.   (Randomly choosing between  return (x+y);  and return ( x^y);  might provide a difficult-to-crack output.)

Another thread by Marsaglia on Jan. 8, 2002:
http://tinyurl.com/cf5d short for
http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&threadm=20020110.1244.47432snz%40genesis.demon.co.uk&rnum=18&prev=/groups%3Fq%3DPRNG%2BOR%2BRNG%2BOR%2Brandom%2B%252264%2Bbit%2B%2522%2B-cypher%2B-penknife%2B-encryption%26hl%3Den%26lr%3D%26ie%3DUTF-8%26oe%3DUTF-8%26start%3D10%26sa%3DN

Here is a very long thread on the sugject started around Jan. 1992:
http://tinyurl.com/cf5r short for
http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&threadm=%25XSe6.3214%24lS6.80400%40news1.oke.nextra.no&rnum=34&prev=/groups%3Fq%3DPRNG%2BOR%2BRNG%2BOR%2Brandom%2B%252264%2Bbit%2B%2522%2B-cypher%2B-penknife%2B-encryption%26hl%3Den%26lr%3D%26ie%3DUTF-8%26oe%3DUTF-8%26start%3D30%26sa%3DN

And another reference on PRNGs from Mark Vaughan:
http://www-rab.larc.nasa.gov/nmp/nmpIndex.htm#RandomNumbers

Sven Pran at 12/13/2003 11:27:49 PM -

The following is my own implentation of a 64 bit PRNG which is given by Knuth as an exceptional good generator. It uses the standard Lehmer algorithm:  Rand64Seed := ((Rand64Seed * a) + c) mod (M)

with M = 2^64, c = 1 and a = $5851F42D4C957F2D;

Global definition:

type

TRand64 = packed record
              case integer of
              1: (RSeed : int64);

              2: (Rft   : TFILETIME);

              3: (RILo  : longWord;
                  RIHi  : longWord);

              4: (RWrk  : array[0..1] of longWord);
             end;

TExtended = packed record
              case integer of
              1: (Rext : extended);

              2: (Rmant : int64;
                  Rexp  : word);
             end;

var

  R64          : TRand64;

procedure NextR64;  {internal procedure}
const
MHI = $5851F42D;
MLO = $4C957F2D;
var
XLO, XHI, YHI, YLO, ZLO, ZHI: longword;
begin
YLO := MLO;
YHI := MHI;
XLO := R64.RiLo;
XHI := R64.RiHi;
asm
  mov EAX,XLO
  mul YLO
  add EAX,1
  adc EDX,0
  mov ZLO,EAX
  mov ZHI,EDX
  mov EAX,XHI
  mul YLO
  add EAX,ZHI
  mov ZHI,EAX
  mov EAX,XLO
  mul YHI
  add EAX,ZHI
  mov ZHI,EAX
end;
  R64.RiLo := ZLO;
  R64.RiHi := ZHI;
end;

function  RandFrac : extended;
var
R : Textended;
begin
NextR64;
if R64.Rseed = 0
  then
   begin
    Result := 0.0;
   end
  else
   begin
    R.Rmant := R64.RSeed;
    R.Rexp := $3FFE;
    while R.Rmant > 0 do
     begin
      R.Rmant := R.Rmant shl 1;
      dec(R.Rexp);
     end;
    Result := R.Rext;
   end;
end;

function  RandWord (const Limit: longWord) : longWord;
var
XLO, XHI, LIM, VAL: longword;
begin
NextR64;
XLO := R64.RiLo;
XHI := R64.RiHi;
LIM := Limit;
asm
  mov EAX,XLO
  mul LIM
  mov VAL,EDX
  mov EAX,XHI
  mul LIM
  add EAX,VAL
  adc EDX,0
  mov VAL,EDX
end;
Result := VAL
end;

Sven Pran at 12/13/2003 7:50:33 AM -

I have never seen the fundamental restrictions on pseudo random generators been published, not even by David E Knuth in his book "The Art of Computer Programming" Vol 2, second edition which was published as late as in 1980 but they do exist (assuming the Lehmer algorithm: X := (a*X + C) mod M ):

Let M be the number of "atomic values" inside the generator (the "cycle" length) and let these numbers be numbered from 0 (zero) to M-1
Let a be the multiplier used in the generator
Let N be the number of "desired alternatives" among which to select the output from the generator (e.g. used as integer parameter in the calls to the "random" function). The "atomic values" in the first segment is then numbered from 0 to approximately M/N and the "atomic values" within the last segment is numbered from approximately (N - 1) * M / N to (M-1)

Each of the N "desired alternatives" is represented by M/N "atomic values" forming a "segment" of contiguous such "atomic values". (The M "atomic values" are grouped into N "desired alternatives" numbered from 0 to N-1)

1: In the call to Random(N) N should never be greater than (a)

Proof: For an "atomic value" in the last segment to be reachable from an "atomic value" in the first segment we must have:

a * M / N > (N - 1) * M / N

This condition resolves to the condition N - 1 < a QED

(Note that the introduction of different values for c in the algorithm does not affect the reality of this proof. If for instance a great c is used to compensate for a too small a so that you will still reach the last segment you will instead be unable to reach the first segment. This fact is most easily appreciated when considering the M "atomic values" on a ring where value 0 immediately follows value M-1. The introduction of different vcalues for c will then only shift all "next" "atomic values" the same amount around the ring).

2: In the call to Random(N) N should never be greater than (SQRT(M)

Proof: As the algorithm results in a definite transition from one "atomic element" to the next we must have at least N "atomic values" in each of the N "segments" or there will always be some segments that cannot be reached from each of the N segments. This gives the condition:

M/N >= N which resolves to N <= SQRT(M) QED

3: To guarantee a statistical error less than p% the number of calls to Random(N) for a single application run must be less than p% of M

Proof: Within any application run each call to the random function "uses" one of the "atomic values" which then cannot be used again until the entire cycle has been used. That means that for each of the N "desired alternatives" returned by the generator the probability to have this same value (between 0 and N-1) is reduced with the rate 1/M.

For subsequent calls to Random there will therefore be a smaller probability for the alternatives already returned than for the alternatives not yet returned and this error is maximally proportional to the number of calls to Random.

In order to guarantee that this statistical error shall not exceed p% we must therefore not have had more than p% of M previous calls to the generator within the same application run. (It is true that this "worst case" only occurs when all calls up to a point have returned the same value among the N alternatives, but there is still a certain however small probability for this to happen).

Sven Pran

Sven Pran at 12/13/2003 11:15:29 PM -

Adding an improved Randomize procedure (applicable for both 32 bits and 64 bits RandSeeds):

For my own use I have implemented the following routine which I consider a significant improvement to the Randomize procedure. The main feature is that it ensures the more significdant bits in RandSeed to become most random.

I use the following little routine to "reverse" the bits in a 32 bits variable:

function  Reversed(source: longint): longint;
asm
mov EDX,EAX
mov ECX,32
@@loop:
RCL EDX,1
RCR EAX,1
LOOP @@LOOP
end;

I have defined the following record type (for 64 biots RandSeed):

TRand64 = packed record
              case integer of
              1: (RSeed : int64);

              2: (Rft   : TFILETIME);

              3: (RILo  : longWord;
                  RIHi  : longWord);

              4: (RWrk  : array[0..1] of longWord);
             end;

I have defined the following global variables:

  fw           : LongWord;
  fr           : TFILETIME;
  fr64         : int64  absolute fr;
  R64          : TRand64;

and I use the following routine to initialize R64 (64 bits randseed):

R64.RIHi := GetTickCount;  {First to initialize only 32 bits in case the 64 bits routine below fails}
repeat
  fW := GetTickCount;
until fW <> R64.RIHi;   {That is until GetTickCount has changed}
fW := fW xor R64.RIHi;
while (fW and 1) = 0 do  {shift out low-order bits that did not change}
  begin
   R64.RIHi := R64.RIHi shr 1;
   fW := fW shr 1;
  end;
R64.RIHi := reversed(R64.RIHi);  {Now we have the "best" possible 32 bit RandSeed R64}
{Here we contine with an attempt to initialize all 64 bits if possible}
GetSystemTimeAsFileTime(fr);  {See if we can get even more bits in R64 randomized}
if fr64 <> 0
  then  {yes we can}
   begin
    repeat
     GetSystemTimeAsFileTime(R64.Rft);
    until R64.RSeed <> fr64;  {that is until GetSystemTimeAsFileTime has changed}
    fr64 := fr64 xor R64.RSeed;
    while (fr64 and 1) = 0 do  {again shift out low-order bits that did not change}
     begin
      R64.RSeed := R64.RSeed shr 1;
      fr64 := fr64 shr 1;
     end;
    fw       := reversed(R64.RiLo);
    R64.RiLo := reversed(R64.RiHi);
    R64.RiHi := fw;   {And now we have the "best" possible 64 bit RandSeed R64}
   end;

John Herbster at 12/15/2003 10:30:34 AM -

The following are suggestions by John Stockton, Surrey, UK posted by John Herbster, with permission:

IMHO, the Random and Randomize functions should not be changed, for compatibility reasons.

A new Unit should instead be provided, containing thoughtfully- named and -written routines including versions of Random, Deal, and Shuffle (and, if there is not one already, a Sorting Unit). The task could be subcontracted to skilled volunteer labour, such as is found here.

The documentation for the existing routines should be changed to say that they are OK for most normal uses but that users should look at the new stuff before deciding.

New randomisers should be capable of initialising their seeds to any value within range, at least if called a small but arbitrary time after the start of a program which was started at an arbitrary time of day on a Windows booted an arbitrary time ago.

If a reasonably efficient reverse algorithm is known, it should be provided; it is necessarily an alternative, and reversibility might be of use.

The new routines should not use a specific RandSeed; instead, each should have a first parameter var seed : <type> .

For completeness, the unit should contain the old algorithms but with the new parameter.

Brian Cook at 12/15/2003 3:36:42 PM -

John,

When the report is opened, the comments are not transferred. All the good ideas you've written won't make it. To ensure the people at Borland have a chance to read your comments, please put them in Description or Steps.

- Brian

John Herbster at 5/26/2008 10:49:26 AM -

Brian,
I tried, but there is not enough room in the Steps area for all of Comments and Implementations infomation.
Rgds, JohnH

John Herbster at 6/10/2008 11:03:49 AM -

I just happened to notice that if the type of "integer" changes from LongInt to Int64, that the present Random(i) function may fail because its header is defined as "function Random [ ( Range: Integer) ];"
--JohnH

Quality Central

Developer Tools

Database Tools