gm2
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Gm2] Problem with urandom


From: Andreas Fischlin
Subject: Re: [Gm2] Problem with urandom
Date: Thu, 22 Oct 2009 02:29:23 +0200
User-agent: Thunderbird 2.0.0.23 (Macintosh/20090812)

Dear Martin,

That is correct, that you get all the time the same result, since you do only access the variable without actually calling any random number generator. AFAIK your implementation of RandInt is wrong, since urandom is not a call to the random number generator but merely a variable of type SeqFile.ChanId. It does not matter that you call it urandom nor that you call your procedure RandInt. This only misleads you; if you call the varaible someInteger and and the procedure TheInteger then it becomes clearer what you actually do. I find it also very confusing that you call procedure WholeIO.ReadInt, i.e. as if you would try to read from a file those random numbers. That would only work if you would have first generated and written the wanted random numbers to that file and then you would read them back from there. But such an approach makes bascially no sense. All practical random number generation is done such that you generate random numbers only in the memory by having a single global variable storing a sequence of integers. To get a decent random number generator you need a 32 bit integer or even better 3 integers used by 3 multiplicative linear congruential random number generators. Such random number generators follow the formula:

z(k+1) = A*z(k) MOD M

Having a single one research tells that it has to be:

M = 2**31 - 1 = 2147483647 = MAX(LONGINT)

with A = 950'706'376 (Fischman & Moore III, 1982). Implementation in Modula-2 is tricky unless you have 64-bit integer arithmetic available, then it is straightforward.

Having three multiplicative linear congruential random number generators is easier and considered even superior to above solution. The interface to a basic random number generator then looks similar to this:

DEFINITION MODULE Randoms;
  PROCEDURE Seed(z0: LONGINT); (*default z0 = 1*)      
  PROCEDURE U(): REAL;  (* returns within (0,1] uniformly distributed variates *)
END Randoms.

DEFINITION MODULE Randoms;
 
  PROCEDURE SetSeeds(z0,z1,z2: INTEGER);
    (*defaults:  z0 = 1, z1 = 10000, z2 = 3000 *)

  PROCEDURE U(): REAL;
    (*returns within (0,1) uniformly distributed variates*)

    (*
      Based on a combination of three multiplicative linear
      congruential random number generators of the form   z(k+1) =
      A*z(k) MOD M   with a prime modulus and a primitive root
      multiplier (=> individual generator full length period). The
      multipliers A are: 171, 172, and 170; the modulus' M are:
      30269, 30307, and 30323.
    *)
END Randoms.

implementation of key routine:

  PROCEDURE U() : REAL;
    VAR
      temp : REAL;
  BEGIN

    (* 1st generator: *)
    x:= 171*(x MOD 177) -  2*(x DIV 177);
    IF x < 0 THEN x:= x+m1 END;

    (* 2nd generator: *)
    y:= 172*(y MOD 176) - 35*(y DIV 176);
    IF y < 0 THEN y:= y+m2 END;

    (* 3rd generator: *)
    z:= 170*(z MOD 178) - 63*(z DIV 178);
    IF z < 0 THEN z:= z+m3 END;

    (* combine to give function value: *)
    temp:= FLOAT(x)/m1R + FLOAT(y)/m2R + FLOAT(z)/m3R;
    RETURN temp-FLOAT(TRUNC(temp));

  END U;

to generate decent integer random numbers in the interval [min, max] use an algorithm similar to this:

  PROCEDURE Jp(min,max: INTEGER): INTEGER;
    CONST absTolerance = 1.0E-7; (* assuming U() returns a 32-bit single precision real *)
  BEGIN
    RETURN min + INTEGER( TRUNC( (FLOAT(max-min+1)-absTolerance)*U() ) )
  END Jp;

You might also find Press et al. (1986, 2002) useful.

It may be that ISO Modula-2 provides a useful random number generator. But in general be aware, unless you know exactly what algorithm any random number generator is using, any scientific applications require careful investigation. Only games may work fine with using predefined random number generators using an unknown generator. Randomization, useful typically in the context of games but not scientific applications, instead of seeding is an entire other issue, I don't have the time to touch upon.

Regards,
Andreas


Cited references:
--------------------
Fishman, G. & L.R. Moore III, 1982. A statistical
                evaluation of multiplicative congruential random
                number generators with modulus 2^31-1.  J.  Amer.
                Statist.  Assoc., 77: 129ff.

Wichmann, B.A. & Hill, I.D., 1982.  An efficient and
                  portable pseudo-random number generator.  Algorithm
                  AS 183. Applied Statistics, 31(2): 188-190.

Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 2002. Numerical recipes multi-language source code CD-ROMs with Windows, DOS, or Mac single screen license. v 2.10, Cambridge University Press, Cambridge, 2 CD ROM pp.

Press, H.W., Flannery, B.P., Teukolsky, S.A. & Vetterling, W.T., 1986. Numerical recipes: the art of scientific computing. Cambridge University Press: New York, 818 pp.


 

Martin Kalbfuß wrote:
Hi,

I created the following module to produce random integers. But when I
call RandInt I always get the same value. Maybe I do something wrong
about The stream. Opening the file with fopen in C an read a value from
it produces a random number. Do you see anything wrong here? I tried it
with stream file, too. 

Thanks

IMPLEMENTATION MODULE Dice;

IMPORT SeqFile,
       WholeIO,
       IOChan,
       EXCEPTIONS,
       SYSTEM;

VAR urandom : SeqFile.ChanId; 
    res     : SeqFile.OpenResults;
    source  : EXCEPTIONS.ExceptionSource;

PROCEDURE RandInt() : INTEGER;
VAR result : INTEGER;
BEGIN
	WholeIO.ReadInt(urandom, result);
	RETURN result;
END RandInt;

BEGIN
	EXCEPTIONS.AllocateSource(source);

	SeqFile.OpenRead(urandom, '/dev/urandom', SeqFile.raw, res);
	IF urandom = IOChan.InvalidChan() THEN
				EXCEPTIONS.RAISE(source, 0, 'Error: Could not access /dev/urandom');
	END;
FINALLY
	SeqFile.Close(urandom);
END Dice.

The calling module:
MODULE Kniffel;

IMPORT Dice,
       SWholeIO,
       STextIO;
BEGIN
	SWholeIO.WriteInt(Dice.RandInt(), 0);
	STextIO.WriteLn;	
	SWholeIO.WriteInt(Dice.RandInt(), 0);
	STextIO.WriteLn;
	SWholeIO.WriteInt(Dice.RandInt(), 0);
	STextIO.WriteLn;		
END Kniffel.

The results:

+134653772
+134653772
+134653772

Always the same number. All the time. No matter what I do.




_______________________________________________
gm2 mailing list
address@hidden
http://lists.nongnu.org/mailman/listinfo/gm2
  

--
________________________________________________________________________
ETH Zurich
Prof. Dr. Andreas Fischlin
Systems Ecology - Institute of Integrative Biology
CHN E 21.1
Universitaetstrasse 16
8092 Zurich
SWITZERLAND

address@hidden
www.sysecol.ethz.ch/staff/af

+41 44 633-6090 phone
+41 44 633-1136 fax

Make it as simple as possible, but distrust it!
________________________________________________________________________
 

Attachment: andreas_fischlin.vcf
Description: Vcard


reply via email to

[Prev in Thread] Current Thread [Next in Thread]