program Percolation2D_v4;

////////////////////////////////////
// Author: Matjaz Gomilsek        //
//  Email: matjaz.gomilsek@ijs.si //
//   Date: 2020/03/31             //
////////////////////////////////////

{$mode objfpc}

uses
  Math, Windows, SysUtils;

type
  PCluster = ^TCluster;
  TCluster = record
    Size: LongInt;
    Complete: Boolean;
  end;
//  TArrayPCluster = array of PCluster; // Not actually used
//  PArrayPCluster = ^TArrayPCluster;   // Not actually used
  TArray2PCluster = array of array of PCluster;
  PArray2PCluster = ^TArray2PCluster;

  TLatticeEnum = (laIsolated, laChain, laChain2, laSquare, laSquare2, laSquareBond, laTriangular, laTriangular2, laHoneycomb, laHoneycomb2, laKagome, laKagome2);
  TLatticeStr = record
    Str: AnsiString;
    LattEnum: TLatticeEnum;
  end;
  TLongIntArray = array of LongInt; //
  TLattice = class
    private
      // Primitive
      _NumB: LongInt;
      _NumN: array of LongInt;
//      _MaxDRAtB: array of LongInt;
      _MaxDR: LongInt;
      // Summary over primitive
      _TotalNumN: LongInt;
    public
      TotalNeighZ: LongInt;                               // Total number of actual neighbours. Only for information, not used at all in the percolation calculations!
      DI, DJ, B2: array of array of LongInt;              // Coordinates and basis indices of neighbours, but only for those accessible through a layer-by-layer scan. The first index is the basis index, the second is the neighbour index at that basis index
      property NumB: LongInt read _NumB;                  // Number of basis indices/vectors
      property NumN: TLongIntArray read _NumN;            // Number of accessible neighbours at a basis index (basis indexing starts from 0!!!) //function NumN(AB: LongInt): LongInt; //property NumN: array of LongInt read _NumN;
      property TotalNumN: LongInt read _TotalNumN;        //function TotalNumN: LongInt; // Total number of accessible neighbours (summed over all basis indices)
//      property MaxDRAtB: array of LongInt read _MaxDRAtB; // Maximum number of layers distance to the furthest accessible neighbour at a given basis index
      property MaxDR: LongInt read _MaxDR;                // Maximum number of layers distance to the furthest accessible neighbour

      constructor Create(ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ: LongInt);
      destructor Destroy; override;
      function Insert(AB, ADI, ADJ, AB2: LongInt): LongInt; // Add a new neighbour at basis index AB
      procedure FinishInserting;                            // Resizes arrays to the actual number of neighbours and the highest B and highest referenced B2. Just for saving memory.
  end;

const
  // Lattice names
  LatticeNames: array [1..65] of TLatticeStr = (
    (Str: 'i';             LattEnum: laIsolated),
    (Str: 'iso';           LattEnum: laIsolated),
    (Str: 'isolated';      LattEnum: laIsolated),
    (Str: 'none';          LattEnum: laIsolated),
    (Str: '0D';            LattEnum: laIsolated),

    (Str: 'c';             LattEnum: laChain),
    (Str: 'ch';            LattEnum: laChain),
    (Str: 'chain';         LattEnum: laChain),
    (Str: '1D';            LattEnum: laChain),

    (Str: 'c2';            LattEnum: laChain2),
    (Str: 'ch2';           LattEnum: laChain2),
    (Str: 'chain2';        LattEnum: laChain2),
    (Str: '1D2';           LattEnum: laChain2),
    (Str: 'cNNN';          LattEnum: laChain2),
    (Str: 'chNNN';         LattEnum: laChain2),
    (Str: 'chainNNN';      LattEnum: laChain2),
    (Str: '1DNNN';         LattEnum: laChain2),

    (Str: 's';             LattEnum: laSquare),
    (Str: 'sq';            LattEnum: laSquare),
    (Str: 'square';        LattEnum: laSquare),

    (Str: 's2';            LattEnum: laSquare2),
    (Str: 'sq2';           LattEnum: laSquare2),
    (Str: 'square2';       LattEnum: laSquare2),
    (Str: 'sNNN';          LattEnum: laSquare2),
    (Str: 'sqNNN';         LattEnum: laSquare2),
    (Str: 'squareNNN';     LattEnum: laSquare2),

    (Str: 'sBond';         LattEnum: laSquareBond),
    (Str: 'sqBond';        LattEnum: laSquareBond),
    (Str: 'squareBond';    LattEnum: laSquareBond),
    (Str: 'sB';            LattEnum: laSquareBond),
    (Str: 'sqB';           LattEnum: laSquareBond),
    (Str: 'squareB';       LattEnum: laSquareBond),

    (Str: 't';             LattEnum: laTriangular),
    (Str: 'tri';           LattEnum: laTriangular),
    (Str: 'triangular';    LattEnum: laTriangular),

    (Str: 't2';            LattEnum: laTriangular2),
    (Str: 'tri2';          LattEnum: laTriangular2),
    (Str: 'triangular2';   LattEnum: laTriangular2),
    (Str: 'tNNN';          LattEnum: laTriangular2),
    (Str: 'triNNN';        LattEnum: laTriangular2),
    (Str: 'triangularNNN'; LattEnum: laTriangular2),

    (Str: 'h';             LattEnum: laHoneycomb),
    (Str: 'hc';            LattEnum: laHoneycomb),
    (Str: 'hon';           LattEnum: laHoneycomb),
    (Str: 'honeycomb';     LattEnum: laHoneycomb),

    (Str: 'h2';            LattEnum: laHoneycomb2),
    (Str: 'hc2';           LattEnum: laHoneycomb2),
    (Str: 'hon2';          LattEnum: laHoneycomb2),
    (Str: 'honeycomb2';    LattEnum: laHoneycomb2),
    (Str: 'hNNN';          LattEnum: laHoneycomb2),
    (Str: 'hcNNN';         LattEnum: laHoneycomb2),
    (Str: 'honNNN';        LattEnum: laHoneycomb2),
    (Str: 'honeycombNNN';  LattEnum: laHoneycomb2),

    (Str: 'k';             LattEnum: laKagome),
    (Str: 'kg';            LattEnum: laKagome),
    (Str: 'kag';           LattEnum: laKagome),
    (Str: 'kagome';        LattEnum: laKagome),

    (Str: 'k2';            LattEnum: laKagome2),
    (Str: 'kg2';           LattEnum: laKagome2),
    (Str: 'kag2';          LattEnum: laKagome2),
    (Str: 'kagome2';       LattEnum: laKagome2),
    (Str: 'kNNN';          LattEnum: laKagome2),
    (Str: 'kgNNN';         LattEnum: laKagome2),
    (Str: 'kagNNN';        LattEnum: laKagome2),
    (Str: 'kagomeNNN';     LattEnum: laKagome2)
  );
  LatticeNiceName: array [TLatticeEnum] of AnsiString = (
    'None',                // laIsolated
    '1D chain (NN)',       // laChain
    '1D chain (NN+NNN)',   // laChain2
    'Square (NN)',         // laSquare
    'Square (NN+NNN)',     // laSquare2
    'Square bonds (NN)',   // laSquareBond
    'Triangular (NN)',     // laTriangular
    'Triangular (NN+NNN)', // laTriangular2
    'Honeycomb (NN)',      // laHoneycomb
    'Honeycomb (NN+NNN)',  // laHoneycomb2
    'Kagome (NN)',         // laKagome
    'Kagome (NN+NNN)'      // laKagome2
  );

var
  MaxR: LongInt;
  P: Extended;
  PrintSparse: Boolean;
  LattEnum: TLatticeEnum;
  Latt: TLattice;
  Hist: array of LongInt; // Histogram of cluster sizes
  HistMaxS, HistTotalClusters, HistTotalOccSites, HistTotalAllSites: LongInt;
  Prob: array of Extended; // Probabilities that an occupied site is in a cluster of a given size (w_s in "_Lecture - Christensen_K_2002 - Percolation Theory - Percolation.pdf")
  ProbWO, ProbWE, ProbSO, ProbW1, ProbWSum, ProbMeanClusterSize: Extended;
  StartAll, StartFindSolution, StopFindSolution: UInt64;

function NewDoublingLength(OldLength, MinLength: LongInt): LongInt; inline;
begin
  if OldLength<=0
    then
      Result := MinLength
    else
      begin
        Result := OldLength;
        while Result<MinLength do
          Result := 2*Result+1;
      end;
end;

//
// --- Lattice ---
//

constructor TLattice.Create(ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ: LongInt);
var
  B, N: LongInt;
begin
  TotalNeighZ := ATotalNeighZ;

   // Initialize
  _MaxDR := 1; //_MaxDR := 0; // Even with no neighbour connections (as in the laIsolated lattice) we can still keep 1 extra layer (it won't hurt too much)
  _NumB := 1;  //_NumB := 0;  // Even with no neighbour connections (as in the laIsolated lattice) there is still at least 1 basis index
  _TotalNumN := 0;

  ExpectedNumB := Max(ExpectedNumB, 0); // Safety
  ExpectedNumB := Max(ExpectedNumB, 0); // Safety

  SetLength(_NumN, ExpectedNumB);
  SetLength(_NumN, ExpectedNumB);
  SetLength(DI, ExpectedNumB, ExpectedNumNPerB);
  SetLength(DJ, ExpectedNumB, ExpectedNumNPerB);
  SetLength(B2, ExpectedNumB, ExpectedNumNPerB);
  for B:=Low(_NumN) to High(_NumN) do
    begin
      _NumN[B] := 0;
      for N:=Low(DI[B]) to High(DI[B]) do
        begin
          DI[B, N] := 0;
          DJ[B, N] := 0;
          B2[B, N] := 0;
        end;
    end;
end;

destructor TLattice.Destroy;
begin
  Assert((Length(_NumN)=Length(DI)) and (Length(DI)=Length(DJ)) and (Length(DJ)=Length(B2)), 'Internal error (lengths of _MaxN, DI, DJ, and B2 are not synchronized)!'); //
  SetLength(_NumN, 0);
  SetLength(DI, 0);
  SetLength(DJ, 0);
  SetLength(B2, 0);

  Inherited;
end;

function TLattice.Insert(AB, ADI, ADJ, AB2: LongInt): LongInt;
var
  B, N, OldNumB, OldLengthB, OldLengthN: LongInt;
begin
  Assert(AB2>=0, 'Internal error (tried to insert negative basis index into lattice)!');

  // Increase the basis number if needed and resize arrays accordingly if needed
  OldNumB := _NumB;
  _NumB := Max(_NumB, Max(AB+1, AB2+1)); // Update _NumB to the max of itself, AB+1, and AB2+1 (so that all references point to existent basis indices AB2)
  OldLengthB := Length(_NumN);
  if _NumB>OldLengthB //if (_NumB>OldNumB) and (_NumB>OldLengthB) // Did we just increase the maximum of basis index AND we have to resize arrays as a consequence of that?
    then
      begin
        SetLength(_NumN, NewDoublingLength(OldLengthB, _NumB));
        SetLength(DI, Length(_NumN));
        SetLength(DJ, Length(_NumN));
        SetLength(B2, Length(_NumN));

        if OldNumB>0
          then
            OldLengthN := _NumN[OldNumB-1] // Use the last basis index number of neighbours for a guess for how many neighbours this basis index will likely have
          else
            OldLengthN := 1;               // Default (start with 1)

        // Zero all the new entries
        for B:=OldLengthB to High(_NumN) do
          begin
            _NumN[B] := 0; // No neighbours yet
            SetLength(DI[B], OldLengthN);
            SetLength(DJ[B], OldLengthN);
            for N:=Low(DI[B]) to High(DI[B]) do
              begin
                DI[B, N] := 0;
                DJ[B, N] := 0;
              end;
          end;
      end;

  // Increase the number of accessible neighbours at this basis index
  _NumN[AB] := _NumN[AB] + 1;
  // Increase the summary total number of accessible neighbours to match
  _TotalNumN := _TotalNumN + 1;

  // Resize DI[AB], DJ[AB], and B2[AB] arrays if they are too small
  OldLengthN := Length(DI[AB]);
  if _NumN[AB]>OldLengthN
    then
      begin
        SetLength(DI[AB], NewDoublingLength(OldLengthN, _NumN[AB]));
        SetLength(DJ[AB], Length(DI[AB]));
        SetLength(B2[AB], Length(DI[AB]));

        // Zero all the new entries
        for N:=OldLengthN to High(DI[AB]) do
          begin
            DI[AB, N] := 0;
            DJ[AB, N] := 0;
            B2[AB, N] := 0;
          end;
      end;

  // Actually insert the new neighbour coordinates
  DI[AB, _NumN[AB]-1] := ADI;
  DJ[AB, _NumN[AB]-1] := ADJ;
  B2[AB, _NumN[AB]-1] := AB2;

  // Update _MaxDR
  _MaxDR := Max(_MaxDR, Max(Abs(ADI), Abs(ADJ)));

  // Return the index of the newly added neighbour in the DI, DJ, and B2 arrays
  Result := (_NumN[AB]-1);
end;

procedure TLattice.FinishInserting;
var
  B: LongInt;
begin
  // Cut down the lengths of all arrays (if needed)
  if Length(_NumN)<>_NumB
    then
      begin
        SetLength(_NumN, _NumB);
        SetLength(DI, _NumB);
        SetLength(DJ, _NumB);
        SetLength(B2, _NumB);
        for B:=Low(_NumN) to High(_NumN) do
          if Length(DI)<>_NumN[B]
            then
              begin
                SetLength(DI, _NumN[B]);
                SetLength(DJ, _NumN[B]);
                SetLength(B2, _NumN[B]);
              end;
      end;
end;

function CreateLattice(ALattEnum: TLatticeEnum): TLattice;
begin
  Result := nil; // Safety

  case ALattEnum of
    laIsolated:    begin
                     Result := TLattice.Create(0, 0, 0); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     // Nothing to insert
                   end;
    laChain:       begin
                     Result := TLattice.Create(1, 1, 2); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
{
DI:    ( -1,NIC,NIC,NIC,NIC,NIC,NIC,NIC), // laChain
DJ:    (  0,NIC,NIC,NIC,NIC,NIC,NIC,NIC), // laChain
}
                         // --- NN (1) ---
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laChain2:      begin
                     Result := TLattice.Create(1, 2, 4); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
{
DI:    ( -1, -2,NIC,NIC,NIC,NIC,NIC,NIC), // laChain2
DJ:    (  0,  0,NIC,NIC,NIC,NIC,NIC,NIC), // laChain2
}
                         // --- NN (1) ---
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         // --- NNN (1) ---
                         Insert(0, -2,  0, 0); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laSquare:      begin
                     Result := TLattice.Create(1, 3, 4); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
{
DI:    (  0, -1,  0,NIC,NIC,NIC,NIC,NIC), // laSquare
DJ:    ( +1,  0, -1,NIC,NIC,NIC,NIC,NIC), // laSquare
}
                         // --- NN (3) ---
                         Insert(0,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 0); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laSquare2:     begin
                     Result := TLattice.Create(1, 6, 8); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
{
DI:    (  0, -1,  0, -1, -1, +1,NIC,NIC), // laSquare2
DJ:    ( +1,  0, -1, +1, -1, -1,NIC,NIC), // laSquare2
}
                         // --- NN (3) ---
                         Insert(0,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 0); // AB, ADI, ADJ, AB2
                         // --- NNN (3) ---
                         Insert(0, -1, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, +1, -1, 0); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laSquareBond:  begin
                     Result := TLattice.Create(2, 4, 6); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
//                         Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//                         Insert(0, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(0,  0, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 1); // AB, ADI, ADJ, AB2

                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(1,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(1,  0, +1, 1); // AB, ADI, ADJ, AB2
//                         Insert(1, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
//                         Insert(1, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(1,  0, -1, 1); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laTriangular:  begin
                     Result := TLattice.Create(1, 4, 6); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
{
DI:    (  0, -1, -1,  0,NIC,NIC,NIC,NIC), // laTriangular
DJ:    ( +1,  0, -1, -1,NIC,NIC,NIC,NIC), // laTriangular
}
                         // --- NN (4) ---
                         Insert(0,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 0); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laTriangular2: begin
                     Result := TLattice.Create(1, 8, 12); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
{
DI:    (  0, -1, -1,  0, +1, -1, -2, -1)  // laTriangular2
DJ:    ( +1,  0, -1, -1, -1, +1, -1, -2)  // laTriangular2
}
                         // --- NN (4) ---
                         Insert(0,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 0); // AB, ADI, ADJ, AB2
                         // --- NNN (4) ---
                         Insert(0, +1, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -2, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -2, 0); // AB, ADI, ADJ, AB2
                       end;
                   end;
    laHoneycomb:   begin
                     Result := TLattice.Create(2, 2, 3); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
                         { {=== OLD BUT CORRECT ===}
                         // --- NN (4) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(0, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(1,  0, +1, 0); // AB, ADI, ADJ, AB2
                         }

                         {} {=== NEW AND OPTIMIZED (newDI = oldDI + oldDJ; newDJ = oldDJ) ===}
                         // --- NN (3) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(0, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
//Insert(1, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
                         {}
                       end;
                   end;
    laHoneycomb2:  begin
                     //Result := TLattice.Create(2, 7, 9); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ // OLD BUT CORRECT
                     Result := TLattice.Create(2, 6, 9); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ // NEW AND OPTIMIZED (newDI = oldDI + oldDJ; newDJ = oldDJ)
                     with Result do
                       begin
                         { {=== OLD BUT CORRECT ===}
                         // --- NN (4) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(0, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(1,  0, +1, 0); // AB, ADI, ADJ, AB2
                         // --- NNN (10) ---
//Insert(0, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(0,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, +1, -1, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(1,  0, +1, 1); // AB, ADI, ADJ, AB2
                         Insert(1, -1, +1, 1); // AB, ADI, ADJ, AB2
                         Insert(1, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(1,  0, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(1, +1, -1, 1); // AB, ADI, ADJ, AB2
                         }

                         {} {=== NEW AND OPTIMIZED (newDI = oldDI + oldDJ; newDJ = oldDJ) ===}
                         // --- NN (3) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(0, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
//Insert(1, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
                         // --- NNN (8) ---
//Insert(0, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
//Insert(0, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
                         Insert(0,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
//Insert(1, +1, +1, 1); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
                         Insert(1,  0, +1, 1); // AB, ADI, ADJ, AB2
                         Insert(1, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(1, -1, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(1,  0, -1, 1); // AB, ADI, ADJ, AB2
                         {}
                       end;
                   end;
    laKagome:      begin
                     //Result := TLattice.Create(3, 3, 4); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ // OLD BUT CORRECT
                     Result := TLattice.Create(3, 4, 4); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ // NEW AND OPTIMIZED (newDI = -oldDI; newDJ = -oldDJ)
                     with Result do
                       begin
                         { {=== OLD BUT CORRECT ===}
                         // --- NN (7) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//Insert(0,  0,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//Insert(0, +1, +1, 1); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
                         Insert(0,  0, +1, 2); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1,  0,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(1, -1, -1, 0); // AB, ADI, ADJ, AB2
                         Insert(1, -1,  0, 2); // AB, ADI, ADJ, AB2
                         Insert(2,  0,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(2,  0,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(2,  0, -1, 0); // AB, ADI, ADJ, AB2
//Insert(2, +1,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         }

                         {} {=== NEW AND OPTIMIZED (newDI = -oldDI; newDJ = -oldDJ) ===}
                         // --- NN (7) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//Insert(0,  0,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(0, -1, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 2); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1,  0,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//Insert(1, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
//Insert(1, +1,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(2,  0,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(2,  0,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(2,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(2, -1,  0, 1); // AB, ADI, ADJ, AB2
                         {}
                       end;
                   end;
    laKagome2:     begin
                     Result := TLattice.Create(3, 7, 8); // ExpectedNumB, ExpectedNumNPerB, ATotalNeighZ
                     with Result do
                       begin
                         // --- NN (7) ---
//Insert(0,  0,  0, 1); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//Insert(0,  0,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
                         Insert(0, -1, -1, 1); // AB, ADI, ADJ, AB2
                         Insert(0,  0, -1, 2); // AB, ADI, ADJ, AB2
                         Insert(1,  0,  0, 0); // AB, ADI, ADJ, AB2
//Insert(1,  0,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (AB2>AB at ADI=ADJ=0)
//Insert(1, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
//Insert(1, +1,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(2,  0,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(2,  0,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(2,  0, +1, 0); // AB, ADI, ADJ, AB2
                         Insert(2, -1,  0, 1); // AB, ADI, ADJ, AB2
                         // --- NNN (8) ---
                         Insert(0,  0, -1, 1); // AB, ADI, ADJ, AB2
//Insert(0, +1,  0, 2); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(0, -1,  0, 1); // AB, ADI, ADJ, AB2
                         Insert(0, -1, -1, 2); // AB, ADI, ADJ, AB2
                         Insert(1,  0, +1, 0); // AB, ADI, ADJ, AB2
//Insert(1, +1, +1, 2); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
//Insert(1, +1,  0, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (i): ADI=+1 and ADJ=0)
                         Insert(1,  0, -1, 2); // AB, ADI, ADJ, AB2
//Insert(2, +1, +1, 0); // AB, ADI, ADJ, AB2 // Redundant (R2>R or L2>L at R2=R, scenario (ii): ADI=ADJ=+1)
                         Insert(2,  0, +1, 1); // AB, ADI, ADJ, AB2
                         Insert(2, -1,  0, 0); // AB, ADI, ADJ, AB2
                         Insert(2, -1, -1, 1); // AB, ADI, ADJ, AB2
                       end;
                   end;
    else
      Assert(False, 'Internal error (requested lattice not implemented)!');
  end;

  Result.FinishInserting; // Clean up. This also minimizes memory usage.
end;

procedure CreateAndOverrideLattice(ALattEnum: TLatticeEnum; var ALattice : TLattice);
begin
  if Assigned(ALattice) // Clear
    then
      FreeAndNil(ALattice);

  ALattice := CreateLattice(ALattEnum); // Create
end;

//
// --- Histogram ---
//

procedure EmptyHistogramStats;
begin
  HistMaxS := -1;
  HistTotalClusters := 0;
  HistTotalOccSites := 0;
  HistTotalAllSites := (MaxR+1)*(MaxR+1)*Latt.NumB;
  ProbWO := 0;
  ProbWE := 0;
  ProbSO := 0;
  ProbWSum := 0;
  ProbMeanClusterSize := 0;
end;

procedure InitHistogram;
var
  I: LongInt;
begin
  SetLength(Hist, MaxR); // Just an arbitrary inital size, e.g. MaxR
  for I:=Low(Hist) to High(Hist) do
    Hist[I] := 0; // Fill with zeros

  EmptyHistogramStats;
end;

procedure AddToHistogram(S: LongInt);
var
  OldLen, I: LongInt;
begin
  if S<1
    then
      Exit;

  if S>HistMaxS
    then
      HistMaxS := S;

  if S>Length(Hist)
    then
      begin
        OldLen := Length(Hist);
        SetLength(Hist, Max(2*OldLen+1, S)); // Either double the size (+1) or, if that's not enough, make it of size S
        for I:=OldLen to High(Hist) do
          Hist[I] := 0;
      end;

  Hist[S-1] := Hist[S-1] + 1;

  HistTotalClusters := HistTotalClusters + 1;
  HistTotalOccSites := HistTotalOccSites + S;
end;

procedure PostprocessHistogram;
var
  S: LongInt;
  Tmp: Extended;
  IsEven: Boolean;
begin
  ProbWO := 0;
  ProbWE := 0;
  ProbSO := 0;
  ProbMeanClusterSize := 0;

  if HistMaxS<1 // Did we find no clusters at all?
    then
      begin
        // As far as we know all occupied sites are probably isolated
        ProbWO := 1;
        ProbSO := 1;
        ProbW1 := 1;
        ProbWSum := 1;
        ProbMeanClusterSize := 1;
        Exit;
      end;

  SetLength(Prob, HistMaxS);

  IsEven := True;
  for S:=1 to HistMaxS do
    begin
      IsEven := not IsEven;

      Tmp := (Hist[S-1]/HistTotalAllSites)/P;
      Prob[S-1] := S*Tmp;

      if IsEven
        then
          ProbWE := ProbWE + Prob[S-1]
        else
          begin
            ProbWO := ProbWO + Prob[S-1];
            ProbSO := ProbSO + Tmp;
          end;

      ProbMeanClusterSize := ProbMeanClusterSize + S*Prob[S-1]; // From "_Lecture - Christensen_K_2002 - Percolation Theory - Percolation.pdf"
    end;

  ProbWSum := ProbWO + ProbWE;
  ProbW1 := Prob[0];
end;

procedure DeleteHistogram;
begin
  SetLength(Hist, 0);
  SetLength(Prob, 0);

  EmptyHistogramStats;
end;

//
// --- Read input ---
//

function ReadInput: Boolean;
var
  S: AnsiString;
  function SkipSpaces(Str: AnsiString): AnsiString;
//  function ExtractWord(Str: AnsiString): AnsiString;
  var
    I, J: LongInt;
  begin
    I := 1;
    while (I<=Length(Str)) and (Str[I]=' ') do
      I := I+1;

    J := Length(Str);
    while (J>I) and (Str[J]=' ') do
      J := J-1;
    {
    J := I+1;
    while (J<=Length(Str)) and (Str[J]<>' ') do
      J := J+1;
    if J>Length(Str)
      then
        J := Length(Str);
    }

    if J<I
      then
        Result := ''
      else
        Result := Copy(Str, I, J-I+1);
  end;
  function InterpretStr: Boolean;
  var
    I: LongInt;
  begin
    S := UpCase(SkipSpaces(S));

    PrintSparse := (Length(S)>=2) and (Copy(S, 1, 2)='-S'); // Sparse histogram output?
    if PrintSparse
      then
        S := SkipSpaces(Copy(S, 3, Length(S)-2));

    Result := False;
    for I:=Low(LatticeNames) to High(LatticeNames) do
      if S=UpCase(SkipSpaces(LatticeNames[I].Str))
        then
          begin
            LattEnum := LatticeNames[I].LattEnum;
            Result := True;
            Break;
          end;
  end;
begin
  ReadLn(MaxR, P, S);

  Result := False; // For sanity
  if not InterpretStr
    then
      WriteLn('ERROR: Unknown lattice!')
    else
      if MaxR<1
        then
          WriteLn('ERROR: Maximum radius must be at least 1!')
        else
          if (P<0) or (P>1)
            then
              WriteLn('ERROR: Site probability must be between 0 and 1!')
            else
              Result := True;
end;

//
// --- Find solution ---
//

procedure FindSolution;
var
  Layers: array of PArray2PCluster;
  R, L, I, J, B, N, DR, TmpL, TmpB: LongInt;
  Neigh: PCluster;

  function NewCluster: PCluster;
  begin
    New(Result);
    Result^.Size := 1;
    Result^.Complete := False; // By default assume it will not be complete //
  end;
  procedure InitLayers;
  var
    DR, L, B: LongInt;
  begin
    SetLength(Layers, Latt.MaxDR+1);
    for DR:=0 to Latt.MaxDR do
      begin
        New(Layers[DR]);
        SetLength(Layers[DR]^, 2*MaxR+1, Latt.NumB);
        for L:=Low(Layers[DR]^) to High(Layers[DR]^) do
          for B:=Low(Layers[DR]^[L]) to Low(Layers[DR]^[L]) do //for B:=0 to Latt.NumB-1 do
            Layers[DR]^[L, B] := nil;
      end;
  end;
  procedure DeleteLayers;
  var
    DR: LongInt;
  begin
    // WARNING: If the cluster pointers were not properly disposed of before this procedure, then they will NOT BE DISPOSED OF AT ALL, causing a memory leak!!!
    for DR:=Low(Layers) to High(Layers) do
      begin
        SetLength(Layers[DR]^, 0); //
        Dispose(Layers[DR]);
        Layers[DR] := nil; //
      end;
    SetLength(Layers, 0);
  end;
  procedure GetIJ; inline;
  begin
    I := Min(L, R);
    J := I + (R - L);
  end;
  procedure GetNeigh; inline;
  var
    I2, J2, B2, R2, L2: LongInt;
  begin
    Neigh := nil; // Default, if not in an accessible part of the lattice

    I2 := I + Latt.DI[B, N];
    if I2<0
      then
        Exit;

    J2 := J + Latt.DJ[B, N];
    if J2<0
      then
        Exit;

    R2 := Max(I2, J2);
    if R2>R //if R<R2
      then
        Exit;
    Assert((R-R2)<=Latt.MaxDR, 'Internal error (MaxDR exceeded)!)'); // Sanity

    L2 := R2 + (I2 - J2);
    if (R2=R) and (L2>L) //if (R2=R) and (L2>=L)
      then
        Exit;

    B2 := Latt.B2[B, N];
    if (R2=R) and (L2=L) and (B2>=B)
      then
        Exit;

    Neigh := Layers[R-R2]^[L2, B2]; // If we got to this line, this is safe to do and refers to an already-filled part of the lattice!
  end;
  function AdvanceLayer: Boolean; // As a side-effect it increases R! Returns "False" when the new R>MaxR+MaxDR[Latt] i.e. when we are more than MaxDR away from any actually calculated clusters.
  var
    MaxDRLayer: PArray2PCluster;
    DR, L, B, TmpL, TmpB: LongInt;
  begin
    // Cycle pointers so that layers retire by one in terms of DR
    MaxDRLayer := Layers[Latt.MaxDR]; // Remember the pointer to the oldest layer
    for DR:=Latt.MaxDR downto 1 do    // Retire each layer by one in terms of DR
      Layers[DR] := Layers[DR-1];

    // Write the completed clusters from the oldest layer to the histogram and dispose of them, clearing the oldest layer in the process
    for L:=0 to 2*(R-Latt.MaxDR) do // NB: this won't do anything if R<Latt.MaxDR (which is correct). It also won't overflow for: MaxR < R <= (MaxR+MaxDR).
      for B:=0 to Latt.NumB-1 do
        if (MaxDRLayer^[L, B]<>nil) and MaxDRLayer^[L, B]^.Complete //and (MaxDRLayer^[L, B]^.Size>0)
          then
            begin
              // Add the cluster to histogram
              AddToHistogram(MaxDRLayer^[L, B]^.Size);

              // Remove all other references to the cluster (they will only be found in the oldest layer at TmpL>L or at TmpL=L and TmpB>B)
              for TmpB:=B+1 to Latt.NumB-1 do
                if MaxDRLayer^[L, TmpB]=MaxDRLayer^[L, B]
                  then
                    MaxDRLayer^[L, TmpB] := nil;
              for TmpL:=L+1 to 2*(R-Latt.MaxDR) do
                for TmpB:=0 to Latt.NumB-1 do
                  if MaxDRLayer^[TmpL, TmpB]=MaxDRLayer^[L, B]
                    then
                      MaxDRLayer^[TmpL, TmpB] := nil;

              // Dispose of the cluster
              Dispose(MaxDRLayer^[L, B]);
              MaxDRLayer^[L, B] := nil;
            end;

    // Prepare a fresh newest (current) layer
    Layers[0] := MaxDRLayer;        // Recycle the pointer from the olders layer to use for a newest (current) layer
    for L:=0 to 2*(R-Latt.MaxDR) do // Clear what we had there before (if anything) and start afresh
      for B:=0 to Latt.NumB-1 do
        Layers[0]^[L, B] := nil;

    // Increase R
    R := R + 1;

    // Are we with the new R still in the neighbour range (MaxDR) of any actually calculated clusters (at R from 0 to MaxR)?
    Result := (R<=(MaxR+Latt.MaxDR));
    // If not, then exit, as we don't have to prepare the new oldest layer for anything (as it is already empty) and we are done.
    if not Result
      then
        Exit;

    // Find out if clusters in the new oldest layer (at DR=MaxDR) are potentially complete, by seeing if they are not present in any newer layer (DR<MaxDR).
    // The newest (current) layer is not considered yet, as it will be filled out (and any clusters marked non-complete) in the main R loop.
    if R>=Latt.MaxDR // Do we actually have an oldest layer that is MaxDR away from the new current layer (at the new R)
      then
        begin
          // Assume that there will be no neighbours in newer layers
          for L:=0 to 2*(R-Latt.MaxDR) do // NB: this won't overflow even for: MaxR < R <= (MaxR+MaxDR).
            for B:=0 to Latt.NumB-1 do
              if (Layers[Latt.MaxDR]^[L, B]<>nil)
                then
                  Layers[Latt.MaxDR]^[L, B]^.Complete := True; // Start by assuming the cluster will be complete

          // See if there actually are any neighbours in newer layers, and if there are, mark the clusters as (potentially) incomplete
          //for DR:=MaxDR[Latt]-1 downto 1 do // No need to check DR=0, as that is the new current layer and is empty
          for DR:=Latt.MaxDR-1 downto Max(1, R-MaxR) do // No need to check DR=0, as that is the new current layer and is empty. Also, no need to look at DR<(R-MaxR) [this would actually cause an overflow] when: MaxR < R <= (MaxR+MaxDR).
            for L:=0 to 2*(R-DR) do
              for B:=0 to Latt.NumB-1 do
                if (Layers[DR]^[L, B]<>nil)
                  then
                    Layers[DR]^[L, B]^.Complete := False; // Mark the cluster as incomplete
        end;
  end;
  function IsOccupied: Boolean; inline;
  begin
    Result := (Random<P);
  end;
begin
  // Init
  InitLayers;
  InitHistogram;

  R := 0;
  while R<=MaxR do
    begin
      // Implicitly starts with the assumption that clusters in the oldest layer will be finished and won't continue to this layer.
      // This assumption will be invalidated if we find any neighbours in this new current layer that will be built in this loop iteration.

      // Build up the new current layer
      for L:=0 to 2*R do
        for B:=0 to Latt.NumB-1 do
          if IsOccupied
            then
              begin
                GetIJ;

              // Scan neighboring cells for existing clusters
                for N:=0 to Latt.NumN[B]-1 do
                  begin
                    GetNeigh;

                    if Neigh<>nil // Did we find a neighbouring cluster?
                      then
                        begin
                          Neigh^.Complete := False; // The neighbouring cluster is then definitely not complete

                          // If this is a bridge between two different clusters merge them. First replace all the current cluster pointers with the neighbouring cluster pointers, and then transfer all of the current cluster's size to the neighbouring one's
                          if (Layers[0]^[L, B]<>nil) and (Layers[0]^[L, B]<>Neigh)
                            then
                              begin
                                for TmpB:=0 to B-1 do // DR=0 and TmpL=L
                                  if Layers[0]^[L, TmpB]=Layers[0]^[L, B]
                                    then
                                      Layers[0]^[L, TmpB] := Neigh;
                                for TmpL:=0 to L-1 do // DR=0
                                  for TmpB:=0 to Latt.NumB-1 do
                                    if Layers[0]^[TmpL, TmpB]=Layers[0]^[L, B]
                                      then
                                        Layers[0]^[TmpL, TmpB] := Neigh;
                                for DR:=1 to Latt.MaxDR do
                                  for TmpL:=0 to 2*(R-DR) do
                                    for TmpB:=0 to Latt.NumB-1 do
                                      if Layers[DR]^[TmpL, TmpB]=Layers[0]^[L, B]
                                        then
                                          Layers[DR]^[TmpL, TmpB] := Neigh;

                                Neigh^.Size := Neigh^.Size + Layers[0]^[L, B]^.Size; // Merge the cluster sizes

                                Dispose(Layers[0]^[L, B]); // No need for a separate pointer to the merged cluster
                              end;

                          Layers[0]^[L, B] := Neigh; // We are in the neighbouring cluster
                        end;
                  end;

                if Layers[0]^[L, B]<>nil // Are we then a part of any previous cluster?
                  then
                    Layers[0]^[L, B]^.Size := Layers[0]^[L, B]^.Size + 1 // If so, add our size (=1) to the cluster that we're in
                  else
                    Layers[0]^[L, B] := NewCluster; // Otherwise create a new cluster
              end;

      // Advance to next layer and increase R. The current layer becomes the last layer. Also, writes the completed clusters to the histogram and disposes of them.
      AdvanceLayer;
    end;

  // Finish the pattern by assuming there is nothing beyond MaxR+MaxDR. This also, implicitly, writes any clusters that were not completed before to the histogram.
  while AdvanceLayer do
    begin
      // Nothing else we have to do.
    end;

  // Done with layers
  DeleteLayers;

  // Postprocess the histogram
  PostprocessHistogram;
end;

//
// --- Write output ---
//

procedure WriteInput;
begin
  WriteLn('--- Input parameters ---');
  WriteLn;
  WriteLn('             Max R = ', MaxR);
  WriteLn('                 P = ', P:0:10);
  WriteLn('           Lattice = ', LatticeNiceName[LattEnum]);
  WriteLn('        Basis size = ', Latt.NumB);
  WriteLn('Number of neigh. Z = ', Latt.TotalNeighZ, ' (had to consider ', Latt.TotalNumN, ')');
  WriteLn('     Max neigh. DR = ', Latt.MaxDR);
  Write(  '     Sparse output = ');
  if PrintSparse
    then
      WriteLn('TRUE')
    else
      WriteLn('FALSE');
end;

procedure WriteOutput;
const
//  Delim = ' ';
  Delim = #9;
var
  S: LongInt;
begin
  WriteLn('--- Cluster histogram ---');
  WriteLn;
  WriteLn('Size', Delim, 'Number', Delim, 'SiteProbability');
  for S:=1 to HistMaxS do
    if (not PrintSparse) or (Hist[S-1]>0)
      then
        WriteLn(S, Delim, Hist[S-1], Delim, Prob[S-1]:0:10);

  WriteLn;
  WriteLn('--- Statistics ---');
  WriteLn;
  WriteLn('       Total occupied sites: ', HistTotalOccSites, ' of ', HistTotalAllSites, ' (fraction ', (HistTotalOccSites/HistTotalAllSites):0:10, ')');
  WriteLn('   Total number of clusters: ', HistTotalClusters);
  WriteLn('          Mean cluster size: ', ProbMeanClusterSize:0:10, ' (', (ProbMeanClusterSize/ProbWSum):0:10, ')');
  WriteLn('    Odd cluster probability:        W_O = ',                       ProbWO:0:10, ' (',              (ProbWO/ProbWSum):0:10, ')');
  WriteLn('   Even cluster probability:        W_E = ',                       ProbWE:0:10, ' (',              (ProbWE/ProbWSum):0:10, ')');
  WriteLn('  Total cluster probability:  W_O + W_E = ',                     ProbWSum:0:10, ' (',                          (1.0):0:10, ')');
  WriteLn(' Minimal residual half-spin:        S_O = ',                       ProbSO:0:10, ' (',              (ProbSO/ProbWSum):0:10, ')');
  WriteLn('      Lone spin probability: S_O >= W_1 = ',                       ProbW1:0:10, ' (',              (ProbW1/ProbWSum):0:10, ')');
  WriteLn('Exact lone spin probability:              ', Power(1-P, Latt.TotalNeighZ):0:10);
end;

procedure WriteTimeElapsed;
begin
  WriteLn('--- Elapsed time ---');
  WriteLn;
  WriteLn('Total elapsed time is ', ((GetTickCount64-StartAll)/1000):0:3, ' seconds (', ((StopFindSolution-StartFindSolution)/1000):0:3, ' seconds calculating percolation).');
end;

//
// --- Main program ---
//

begin
  StartAll := GetTickCount64;

  Randomize;

  if ReadInput
    then
      begin
        Latt := nil; // For safety
        CreateAndOverrideLattice(LattEnum, Latt); // Create the list of neighbours for this lattice

        WriteInput;

        StartFindSolution := GetTickCount64;
        FindSolution;
        StopFindSolution := GetTickCount64;

        WriteLn;
        WriteOutput;

        DeleteHistogram;  // Delete the histogram
        FreeAndNil(Latt); // Delete the lattice

        WriteLn;
        WriteTimeElapsed;
      end;
end.
