Path: utzoo!news-server.csri.toronto.edu!cs.utexas.edu!wuarchive!zaphod.mps.ohio-state.edu!rpi!uupsi!sunic!dkuug!daimi!protonen
From: protonen@daimi.aau.dk (Lars J|dal)
Newsgroups: comp.lang.pascal
Subject: Re: Random numbers
Message-ID: <1991Mar13.131738.16813@daimi.aau.dk>
Date: 13 Mar 91 13:17:38 GMT
References: <1991Mar7.193727.1629@tjhsst.vak12ed.edu>
Sender: protonen@daimi.aau.dk (Lars J|dal)
Organization: DAIMI: Computer Science Department, Aarhus University, Denmark
Lines: 64

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