jmitchel@tjhsst.vak12ed.edu (John Mitchell) (03/08/91)
How does TP 5.5 generate random numbers? We are using this statement: r := random; which returns a real number in the range [0,1). Does TP use some kind of formula to generate these values? Or is there an internal list of random numbers that is accessed? Does "randomize" seed a formula or find a spot in the list to start from? My students are writing a program to approximate pi using the Monte Carlo method. Our problem is that we are experiencing a relatively wide range of answers for pi (+/- 0.2) which is not improved by running the program with more iterations. Thanks in advance for your help. == John Mitchell== P.S. Does TP generate random integers using the same method? P.P.S. Would anyone be able to send me a formula for generating good random numbers that might be substituted for the TP function call?
ajburke@casbah.acns.nwu.edu (Austin Burke) (03/08/91)
In article 1624 John Mitchell writes: > How does TP 5.5 generate random numbers? The numbers generated are not really random. Each number is derived from the preceding "random" number by a complex math equation, and then forced (by another formula) to fit into the given range. > We are using this statement: r := random; which returns a >real number in the range [0,1). Does TP use some kind of formula >to generate these values? Or is there an internal list of random >numbers that is accessed? Does "randomize" seed a formula or >find a spot in the list to start from? There is no internal list. The formula used is, approximately, the following: seed = (multiplier * seed) + (incrementor mod modulus) where "seed" is the previous number generated, "multiplier" is some very large constant number, as is "increment," and "modulus" is 65536. The result of this formula is then forced to fit into the range given by the user. The "randomize" command gives TP its first "seed." This first seed is derived (I think) from the current time, measured in hundreths of a second (or even more exact - I'm not sure). Through this process, seemingly random numbers are generated. It is conceivable, but not likely, that the exact same sequence of "random" numbers may be generated if the initial seed is the same. This is why the current time is used - it makes this possiblity all the more unlikely. > My students are writing a program to approximate pi using the >Monte Carlo method. Our problem is that we are experiencing a >relatively wide range of answers for pi (+/- 0.2) which is not >improved by running the program with more iterations. Sorry, I'm not familiar with the Monte Carlo method. >P.S. Does TP generate random integers using the same method? >P.P.S. Would anyone be able to send me a formula for generating >good random numbers that might be substituted for the TP function >call? See above. -- __ _____ ___ ___________________________ / | | |__) Austin J. Burke III /---| | | \ ajburke@casbah.acns.nwu.edu / | \__/ |__/ Northwestern University
reino@cs.eur.nl (Reino de Boer) (03/08/91)
In <1991Mar8.031013.14902@casbah.acns.nwu.edu> ajburke@casbah.acns.nwu.edu (Austin Burke) writes: >In article 1624 John Mitchell writes: >> How does TP 5.5 generate random numbers? > ...... > There is no internal list. The formula used is, approximately, the >following: > seed = (multiplier * seed) + (incrementor mod modulus) This probably is: seed = ((multiplier * seed) + incrementor) mod modulus where "modulus == computer word size" to simplify the calculation. >>P.P.S. Would anyone be able to send me a formula for generating >>good random numbers that might be substituted for the TP function >>call? Check out Volume 2 of Donald E. Knuth -- The Art of Computer Programming The first half of this book is entirely dedicated to pseudo random numbers and their generation. Following are some routines I recently used to check some properties of prime numbers (this is sun-pascal, although I think it's even level 0): function random1 : float; const a = 1664525.0; c = 1.0; m = 4294967296.0; begin seed1 := modulo( seed1 * a + c, m ); random1 := seed1 / m end; { random1 } function randint1( max : int ) : int; begin randint1 := floor( random1 * max ) end; { randint1 } function random2 : float; const a = 17059465.0; c = 1.0; m = 34359738368.0; begin seed2 := modulo( seed2 * a + c, m ); random2 := seed2 / m end; { random2 } function randint2( max : int ) : int; begin randint2 := floor( random2 * max ) end; { randint2 } function randnumber( digits : integer ) : int; var result : int; i : integer; begin result := floor( random( 9.0 ) * 9.0 + 1.0 ); for i := 2 to digits do result := 10.0 * result + floor( random( 10.0 ) * 10.0 ); randnumber := result end; { randnumber } If needed, I can look up the implementation of modulo() and floor() for you, although they are fairly trivial. Hope this helps -- Reino -- Reino R. A. de Boer "We want to build the right product right, right?" Erasmus University Rotterdam ( Informatica ) e-mail: reino@cs.eur.nl
CDCKAB%EMUVM1.BITNET@cunyvm.cuny.edu ( Karl Brendel) (03/08/91)
In article <1991Mar8.031013.14902@casbah.acns.nwu.edu>, ajburke@casbah.acns.nwu.edu (Austin Burke) wrote: >In article 1624 John Mitchell writes: > >> How does TP 5.5 generate random numbers? > > The numbers generated are not really random. Each number is >derived from the preceding "random" number by a complex math >equation, and then forced (by another formula) to fit into the given >range. > >> We are using this statement: r := random; which returns a >>real number in the range [0,1). Does TP use some kind of formula >>to generate these values? Or is there an internal list of random >>numbers that is accessed? Does "randomize" seed a formula or >>find a spot in the list to start from? > > There is no internal list. The formula used is, >approximately, the following: > > seed = (multiplier * seed) + (incrementor mod modulus) > > where "seed" is the previous number generated, "multiplier" >is some very large constant number, as is "increment," and "modulus" >is 65536. The result of this formula is then forced to fit into the >range given by the user. > The "randomize" command gives TP its first "seed." This >first seed is derived (I think) from the current time, measured in >hundreths of a second (or even more exact - I'm not sure). Specifics of the random number generator have been discussed at length in the list by Duncan Murdoch <dmurdoch@watstat.waterloo.edu>. If he doesn't respond directly to this question, email to him or a check of the archives would be appropriate. See particularly articles <1991Jan14.213834.26624@maytag.waterloo.edu> and <1991Jan15.163811.15107@maytag.waterloo.edu>. (Yes, Duncan, I keep a _lot_ of your stuff. :) ) > Through this process, seemingly random numbers are >generated. It is conceivable, but not likely, that the exact same >sequence of "random" numbers may be generated if the initial seed is >the same. This is why the current time is used - it makes this >possiblity all the more unlikely. See the documentation for Randomize. A seed value can be _assigned_. Reuse of a seed value _will_ result in 'the exact same sequence of "random" numbers" being generated. > >> My students are writing a program to approximate pi using the >>Monte Carlo method. Our problem is that we are experiencing a >>relatively wide range of answers for pi (+/- 0.2) which is not >>improved by running the program with more iterations. I wonder if the source of the problem is really in the random number generator itself. Would you like to throw out some code for us to look at? :) > > Sorry, I'm not familiar with the Monte Carlo method. > >>P.S. Does TP generate random integers using the same method? >>P.P.S. Would anyone be able to send me a formula for generating >>good random numbers that might be substituted for the TP function >>call? This routine is based on "Abramowitz & Stegun (1965)", implemented in code received from a statistician whose name I have misplaced. You should be able to modify it for real numbers as well. var seed : longint; function randomInteger(least, most : integer):integer; begin seed := (seed * 23) mod (100000001); randomInteger := least + round((most - least)*seed/100000000); end; {of randomInteger} +--------------------------------------------------------------------+ | Karl Brendel Centers for Disease Control | | Internet: CDCKAB@EMUVM1.BITNET Epidemiology Program Office | | Bitnet: CDCKAB@EMUVM1 Atlanta, GA, USA | | Home of Epi Info 5.0 | +--------------------------------------------------------------------+
protonen@daimi.aau.dk (Lars J|dal) (03/15/91)
jmitchel@tjhsst.vak12ed.edu (John Mitchell) writes: > My students are writing a program to approximate pi using the >Monte Carlo method. Our problem is that we are experiencing a >relatively wide range of answers for pi (+/- 0.2) which is not >improved by running the program with more iterations. The Monte Carlo method of finding pi is trying to "shoot" at a 1 times 1 square and see how many is within distance 1.0 from origo, right? (That should give pi/4) >P.P.S. Would anyone be able to send me a formula for generating >good random numbers that might be substituted for the TP function >call? Here follows a little program which should be better than the one you describes (I have used it to find pi and it's about +/- 0.03 from pi after a few thousend "shots"). Unfortunately I don't have TP myself, so this hasn't been run in pascal, so I there'll probably be a few missing ;'s and so. var seed : array [1..3] of integer; procedure initrnd; begin seed[1] := 1; (* 0 < arbitrary number < 30269 *) seed[2] := 10000; (* 0 < arbitrary number < 30307 *) seed[3] := 5000; (* 0 < arbitrary number < 30323 *) end; procedure rndrandomize; begin (* This should call the clock and get the minutes, seconds and *) (* 1/100 seconds (say). - Check for 0 as this only gives 0's *) (* (if e.g. sec = 0 you could just give it the value 117 (say) ). *) end; function rnd : real; var temp : real; begin (* This is the actual random generator *) seed[1] := 171*(seed[1] mod 177) - 2*(seed[1] div 177); if seed[1] < 0 then seed[1] := seed[1]+30269; seed[2] := 172*(seed[2] mod 176) - 35*(seed[2] div 176); if seed[2] < 0 then seed[2] := seed[2]+30307; seed[3] := 170*(seed[3] mod 178) - 63*(seed[3] div 178); if seed[3] < 0 then seed[3] := seed[3]+30323; temp := seed[1]/30269 + seed[2]/30307 + seed[3]/30323; temp := temp-int(temp); (* Get fractional part *) rnd := temp; end; Don't ask me why this works (although it looks good :-) ). It's from an issue of BYTE from some years ago (I can't remember which). Lars J|dal protonen@daimi.aau.dk