home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / sigm / vols200 / vol243 / near.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1994-07-13  |  3.3 KB  |  92 lines
  1. {AUTHOR: Mark Aldon Weiss    donated to public domain}
  2.  
  3. CONST
  4.  
  5. MaxNumPts  =  2000;
  6.  
  7. MaxNumRuns =    500;
  8.  
  9.  
  10.  
  11. VAR
  12.  
  13. X,Y,Distance: Array[1..MaxNumPts] of Real;
  14.  
  15. NearestNeighbor: Array[1..MaxNumPts] of Integer;
  16.  
  17. Probability: Array[1..MaxNumRuns] of Real;
  18.  
  19. i,j,NumPts: 1..MaxNumPts;  NumRuns,Run:1..MaxNumRuns;
  20.  
  21. NumCommuting: 1..MaxNumPts;
  22.  
  23. ch: char;
  24.  
  25. Smallest,StDev,Avg: Real;
  26.  
  27.  
  28.  
  29. BEGIN  { M A I N    P R O G R A M }
  30.  
  31. Writeln;
  32. Writeln('This program is a Monte Carlo computation of the following problem.');
  33. Writeln('Given random points in the plane, what is the probability that a');
  34. Writeln('point''s nearest neighbor has as IT''s nearest neighbor, the first');
  35. Writeln('point?'#7);
  36. Writeln;
  37. Writeln('TURN THE PRINTER ON before we continue.');
  38. Writeln;
  39. REPEAT
  40.  Write('Did you turn the printer on?   ');  Readln(ch);  Writeln;
  41.  IF NOT( ch IN ['y','Y'] ) THEN Writeln('Well go turn it on'#7#7#7)
  42. UNTIL ch IN ['y','Y'];
  43. Writeln(lst);
  44. Writeln(lst,#27'E RESULTS of NEAREST NEIGHBOR TEST on TURBO PASCAL RNG');
  45. Writeln(lst,#27'F');
  46. Write('Good, now how many points per run do you want (',MaxNumPts:5,' max)? ');
  47. Readln(NumPts);  Writeln;
  48. Write('How many runs do you want (',MaxNumRuns:3,' max)?    ');
  49. Readln(NumRuns);  Writeln;
  50. Writeln(lst,#27'G EACH RUN HAS ',NumPts:3,' POINTS'#27'H');  Writeln(lst);
  51. FOR  Run := 1  to  NumRuns  DO
  52.      Begin
  53.      for i := 1 to NumPts do  begin  X[i] := RANDOM;  Y[i] := RANDOM  end;
  54.      NumCommuting := 0;
  55.      FOR i := 1 to NumPts DO
  56.          Begin
  57.          FOR j := 1 to NumPts DO
  58.          Distance[j] := SQR( X[i] - X[j] ) + SQR( Y[i] - Y[j]);
  59.          { whether its the square of the distance or not does't matter }
  60.          { why have the computer take the sqaure root when it's not needed }
  61.          Distance[i] := 20;  {any number > 1 so dist. pt. to self not smallest}
  62.          Smallest := Distance[1];  NearestNeighbor[i] := 1;
  63.          FOR j := 2 to Numpts DO
  64.              IF Distance[j] < Smallest THEN
  65.                 begin  Smallest := Distance[j];  NearestNeighbor[i] := j  end
  66.          End;
  67.      FOR i := 1 to NumPts DO
  68.          IF  NearestNeighbor[ NearestNeighbor[i] ] = i  THEN
  69.              NumCommuting := NumCommuting + 1;
  70.      Probability[Run] := NumCommuting/NumPts;
  71.      Writeln;
  72.      Writeln(lst);
  73.      Writeln('Run',Run:4,' There were ',NumCommuting:4,' commuting points out of');
  74.      Writeln(Numpts:5,' the probability was therfore ',Probability[Run]:7:5);
  75.      Writeln(lst,'Run',Run:4,' probability = ',Probability[Run]:11:8);
  76.      Writeln(lst);
  77.      Writeln;
  78.      Avg := 0;  StDev := 0;
  79.      for i := 1 to Run do Avg := Avg + Probability[i];
  80.      Avg := Avg/Run;
  81.      for i := 1 to Run do StDev := StDev + SQR( Probability[i] - Avg );
  82.      if Run <> 1 then StDev := SQRT( StDev/(Run-1) ) else StDev := 0;
  83.      Writeln('The average so far is     ',Avg:7:5);
  84.      Writeln('The standard deviation is ',StDev:7:5);
  85.      Writeln(lst,'The average so far is     ',Avg:11:8);
  86.      Writeln(lst,'The standard deviation is ',StDev:11:8)
  87.      End;
  88. Writeln;  Writeln;
  89. Writeln('The nearest neighbor problem may be solved analytically.');
  90. Writeln('The result is  6pi / ( 8pi + 3 SQRT(3) ) = 0.621505')
  91.  
  92. END.   { M A I N    P R O G R A M }