{$MODE OBJFPC} { -*- delphi -*- }
{$INCLUDE settings.inc}
program test_arrayutils_shuffle;
uses
   arrayutils, hashtable, hashfunctions, genericutils;

const
   ListLength = 5;
   ShuffleAttempts = 10000;

type
   TPermutationName = String[ListLength];
   TPermutationNameUtils = specialize DefaultUtils <TPermutationName>;
   TPermutationHashTable = specialize THashTable <TPermutationName, Cardinal, TPermutationNameUtils>;

function Hash(const Key: TPermutationName): Cardinal;
begin
   Result := AnsiStringHash32(Key);
end;

function Factorial(N: Cardinal): Cardinal;
begin
   Result := N;
   while (N > 2) do
   begin
      Dec(N);
      Result := Result * N;
   end;
end;

function CalculatePValue(ChiSquare: Double; DegreesOfFreedom: Cardinal): Double;
var
   Q, P, T: Double;
   K, A: Cardinal;
begin
   // Having failed to find any articles online explaining how to calculate the p-value that
   // were not aimed at people with PhDs, I found this page:
   //    http://www.easycalculation.com/statistics/chi-squared-p-value.php
   // ...which implemented something in JavaScript which I have ported here.
   if ((DegreesOfFreedom = 1) and (ChiSquare > 1000.0)) then
   begin
      Result := 0;
      Exit;
   end;
   if ((ChiSquare > 1000) or (DegreesOfFreedom > 1000)) then
   begin
      Q := CalculatePValue((ChiSquare - DegreesOfFreedom) * (ChiSquare - DegreesOfFreedom) / (2.0 * ChiSquare), 1) / 2.0;
      if (ChiSquare > DegreesOfFreedom) then
         Result := Q
      else
         Result := 1.0 - Q;
      Exit;
   end;
   P := Exp(-0.5 * ChiSquare);
   if (DegreesOfFreedom mod 2 = 1) then
      P := P * Sqrt(2 * ChiSquare / Pi);
   K := DegreesOfFreedom;
   while (K > 2) do
   begin
      P := P * ChiSquare / K;
      Dec(K, 2);
   end;
   T := P;
   A := DegreesOfFreedom;
   while (T > 0.0000000001*P) do
   begin
      Inc(A, 2);
      T := T * ChiSquare / A;
      P := P + T;
   end;
   Result := 1.0 - P;
end;

var
   Base, Test: TPermutationName;
   Permutations: TPermutationHashTable;
   PermutationCount, Index: Cardinal;
   ChiSquare, Ideal, Delta, P: Double;
begin
   //Randomize(); // leave it at default seed so we can reproduce results
   SetLength(Base, ListLength);
   for Index := 1 to Length(Base) do
      Base[Index] := Chr(Ord('A')+Index-1);
   PermutationCount := Factorial(ListLength);
   Ideal := ShuffleAttempts / PermutationCount;
   Assert(Ideal >= 5);
   Permutations := TPermutationHashTable.Create(@Hash, PermutationCount);
   for Index := 1 to ShuffleAttempts do
   begin
      Test := Base;
      FisherYatesShuffle(Test[1], Length(Test), SizeOf(Test[1]));
      Permutations[Test] := Permutations[Test]+1;
   end;
   ChiSquare := 0;
   for Test in Permutations do
   begin
      Delta := (Permutations[Test] - Ideal);
      ChiSquare := ChiSquare + Delta*Delta/Ideal;
   end;
   P := CalculatePValue(ChiSquare, PermutationCount-1);
   Permutations.Free();
   if (P < 0.05) then
   begin
      Writeln('FisherYatesShuffle does not seem random. ChiSquared test of ', ShuffleAttempts, ' shuffles of ', ListLength, '-item arrays gave probability of this extreme a distribution happening if the shuffle truly is random at p=', P:0:5);
      Halt(1);
   end;
end.