{AUTHOR: Mark Aldon Weiss    donated to public domain}

CONST

MaxNumPts  =  2000;

MaxNumRuns =    500;



VAR

X,Y,Distance: Array[1..MaxNumPts] of Real;

NearestNeighbor: Array[1..MaxNumPts] of Integer;

Probability: Array[1..MaxNumRuns] of Real;

i,j,NumPts: 1..MaxNumPts;  NumRuns,Run:1..MaxNumRuns;

NumCommuting: 1..MaxNumPts;

ch: char;

Smallest,StDev,Avg: Real;



BEGIN  { M A I N    P R O G R A M }

Writeln;
Writeln('This program is a Monte Carlo computation of the following problem.');
Writeln('Given random points in the plane, what is the probability that a');
Writeln('point''s nearest neighbor has as IT''s nearest neighbor, the first');
Writeln('point?'#7);
Writeln;
Writeln('TURN THE PRINTER ON before we continue.');
Writeln;
REPEAT
 Write('Did you turn the printer on?   ');  Readln(ch);  Writeln;
 IF NOT( ch IN ['y','Y'] ) THEN Writeln('Well go turn it on'#7#7#7)
UNTIL ch IN ['y','Y'];
Writeln(lst);
Writeln(lst,#27'E RESULTS of NEAREST NEIGHBOR TEST on TURBO PASCAL RNG');
Writeln(lst,#27'F');
Write('Good, now how many points per run do you want (',MaxNumPts:5,' max)? ');
Readln(NumPts);  Writeln;
Write('How many runs do you want (',MaxNumRuns:3,' max)?    ');
Readln(NumRuns);  Writeln;
Writeln(lst,#27'G EACH RUN HAS ',NumPts:3,' POINTS'#27'H');  Writeln(lst);
FOR  Run := 1  to  NumRuns  DO
     Begin
     for i := 1 to NumPts do  begin  X[i] := RANDOM;  Y[i] := RANDOM  end;
     NumCommuting := 0;
     FOR i := 1 to NumPts DO
         Begin
         FOR j := 1 to NumPts DO
         Distance[j] := SQR( X[i] - X[j] ) + SQR( Y[i] - Y[j]);
         { whether its the square of the distance or not does't matter }
         { why have the computer take the sqaure root when it's not needed }
         Distance[i] := 20;  {any number > 1 so dist. pt. to self not smallest}
         Smallest := Distance[1];  NearestNeighbor[i] := 1;
         FOR j := 2 to Numpts DO
             IF Distance[j] < Smallest THEN
                begin  Smallest := Distance[j];  NearestNeighbor[i] := j  end
         End;
     FOR i := 1 to NumPts DO
         IF  NearestNeighbor[ NearestNeighbor[i] ] = i  THEN
             NumCommuting := NumCommuting + 1;
     Probability[Run] := NumCommuting/NumPts;
     Writeln;
     Writeln(lst);
     Writeln('Run',Run:4,' There were ',NumCommuting:4,' commuting points out of');
     Writeln(Numpts:5,' the probability was therfore ',Probability[Run]:7:5);
     Writeln(lst,'Run',Run:4,' probability = ',Probability[Run]:11:8);
     Writeln(lst);
     Writeln;
     Avg := 0;  StDev := 0;
     for i := 1 to Run do Avg := Avg + Probability[i];
     Avg := Avg/Run;
     for i := 1 to Run do StDev := StDev + SQR( Probability[i] - Avg );
     if Run <> 1 then StDev := SQRT( StDev/(Run-1) ) else StDev := 0;
     Writeln('The average so far is     ',Avg:7:5);
     Writeln('The standard deviation is ',StDev:7:5);
     Writeln(lst,'The average so far is     ',Avg:11:8);
     Writeln(lst,'The standard deviation is ',StDev:11:8)
     End;
Writeln;  Writeln;
Writeln('The nearest neighbor problem may be solved analytically.');
Writeln('The result is  6pi / ( 8pi + 3 SQRT(3) ) = 0.621505')

END.   { M A I N    P R O G R A M }
.