unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics,unit1, Controls, Forms, Dialogs;

Function GASDEV:real;
procedure KSTWO(DATA1:array of real; N1:integer;
                    DATA2:array of real; N2:integer;var D, PROB:real);
Function ERF(X:real):real;

implementation
Function GAMMLN(xx:real):real;
const
  STP=2.50662827465; HALF=0.5; ONE=1.0; FPF=5.5;
var
  x,tmp,ser:double;
  j:integer;
  cof:array[1..6] of double;
begin
  COF[1]:=76.18009173;      COF[2]:=-86.50532033;
  COF[3]:=24.01409822;      COF[4]:=-1.231739516;
  COF[5]:= 0.120858003e-2;  COF[6]:=-0.536382e-5;
  X:=XX-ONE;
  TMP:=X+FPF;
  TMP:=(X+HALF)*Ln(TMP)-TMP;
  SER:=ONE;
  For J:=1 To 6 do
  begin
      X:=X+ONE;
      SER:=SER+COF[J]/X
  end;
  GAMMLN:=TMP+Ln(STP*SER);
end;

procedure GCF(var GAMMCF:real; A:real; X:real; var GLN:real);
label 1;
const
  ITMAX=100;    EPS=0.0000003;
var
  N:integer; GOLD,G,FAC,B1,B0,ANF,ANA,AN,A1,A0:real;
begin
  GLN:=GAMMLN(A);
  GOLD:=0;
  A0:=1;
  A1:=X;
  B0:=0;
  B1:=1;
  FAC:=1;
  For N:=1 To ITMAX do
  begin
    AN:=N;
    ANA:=AN - A;
    A0:=(A1 + A0 * ANA)* FAC;
    B0:=(B1 + B0 * ANA)* FAC;
    ANF:=AN * FAC;
    A1:=X * A0 + ANF * A1;
    B1:=X * B0 + ANF * B1;
    If A1 <> 0 Then
    begin
      FAC:=1 / A1;
      G:=B1 * FAC;
      If Abs((G - GOLD)/ G)< EPS Then GoTo 1;
      GOLD:=G;
    end
  end;
  ShowMessage('A too large, ITMAX too small');
1: GAMMCF:=Exp(-X + A * Ln(X)- GLN)* G;
end;

procedure GSER(var GAMSER:REAL;A:real;X:REAL;var GLN:REAL);
label 1;
const
  ITMAX=100;    EPS=0.3E-6;
var
  N:Integer;
  SUM,DEL,AP:Real;
begin
  GLN:=GAMMLN(A);
  If X <= 0 Then
    If X < 0 Then
      begin
        ShowMessage('警告：X<0，退出程序。');
        goto 1;
      end
  else
    begin
      GAMSER:=0;
      goto 1;
    end;
  AP:=A;
  Sum:=1 / A;
  DEL:=Sum;
  For N:=1 To ITMAX do
    begin
      AP:=AP + 1;
      DEL:=DEL * X / AP;
      Sum:=Sum + DEL;
      If Abs(DEL) < Abs(Sum) * EPS Then GoTo 1;
    end;
  ShowMessage('A too large, ITMAX too small');
1: GAMSER:=Sum * Exp(-X + A * Ln(X) - GLN);
end;

Function GAMMP(A,X:real):real;
var
  GAMMCF,GLN:REAL;
begin
    If (X < 0) Or (A <= 0) Then  ShowMessage('PAUSE');
    If X < A + 1 Then
    begin
        GSER(GAMMCF, A, X, GLN);
        GAMMP:=GAMMCF;
    end
    Else
    begin
        GCF(GAMMCF, A, X, GLN);
        GAMMP:= 1 - GAMMCF;
    end;
end;

Function ERF(X:real):real;
begin
    If X < 0 Then
        ERF:=-GAMMP(0.5, sqr(X))
    Else
        ERF:=GAMMP(0.5, sqr(X));
end;

Function GASDEV:real;
var
   V1,V2,FAC,R:real;
begin
    If ISET^= 0 Then
    begin
      repeat
        V1:=2 * Random - 1;
        V2:=2 * Random - 1;
        R:=Sqr(V1) + Sqr(V2);
      until (R < 1);
      FAC:=Sqrt(-2 * Ln(R) / R);
      GSET^:=V1 * FAC;
      GASDEV:=V2 * FAC;
      ISET^:=1;
    end
    Else
    begin
      GASDEV:=GSET^;
      ISET^:=0;
    end;
end;

Procedure SORT(N:integer; var RA:array of real);
Label 99;
var
    I,J,L,IR:integer;   RRA:real;
begin
    L:= N div 2 + 1;
    IR:=N;
    While true do
    begin
        If L > 1 Then
        begin
            L:=L - 1;
            RRA:=RA[L];
        end
        Else
        begin
            RRA:=RA[IR];
            RA[IR]:=RA[1];
            IR:=IR - 1;
            If IR = 1 Then
            begin
              RA[1]:=RRA;
              goto 99;
            end;
        end;
        I:=L;
        J:=L + L;
        While J <= IR do
        begin
            If J < IR Then
                If RA[J] < RA[J + 1] Then J:=J + 1;
            If RRA < RA[J] Then
            begin
                RA[I]:=RA[J];
                I:=J;
                J:=J + J;
            end
            Else
                J:=IR + 1;
        end;
        RA[I]:=RRA;
    end;
99: end;

Function PROBKS(ALAM:real):real;
label 99;
const
    EPS1 = 0.001;    EPS2 = 0.00000001;
var 
    A2,FAC,TERMBF,TERM,AAA,PRO:real;  J:integer;
begin
    A2:=-2 * Sqr(ALAM);
    FAC:=2; 
    PRO:=0; 
    TERMBF:=0; 
    For J:=1 To 100 do
    begin
        TERM:=FAC * Exp(A2 * J*J);
        PRO:=PRO + TERM;
        AAA:=Abs(TERM);
        If (AAA < EPS1 * TERMBF) Or (AAA < EPS2 * PRO) Then
            goto 99;
        FAC:=-FAC;
        TERMBF:=Abs(TERM);
    end;
    PRO:=1;
99: PROBKS:=PRO;
end;

procedure KSTWO(DATA1:array of real; N1:integer;
                    DATA2:array of real; N2:integer;var D, PROB:real);
label 1;
var
    J1,J2:integer;  FO1,FO2,FN1,FN2,EN1,EN2,AM,DT:real;
begin
    SORT(N1, DATA1);
    SORT(N2, DATA2);
    EN1:=N1;
    EN2:=N2;
    J1:=1;
    J2:=1;
    FO1:=0;
    FO2:=0;
    D:=0; 
1:  If (J1 <= N1) And (J2 <= N2) Then
    begin
        If DATA1[J1] < DATA2[J2] Then
        begin
            FN1:=J1 / EN1;
            If Abs(FN1 - FO2) > Abs(FO1 - FO2) Then
                AM:=Abs(FN1 - FO2)
            Else
                AM:=Abs(FO1 - FO2);
            DT:=AM;
            If DT > D Then D:=DT;
            FO1:=FN1;
            J1:=J1 + 1;
        end
        Else
        begin
            FN2:=J2 / EN2;
            If Abs(FN2 - FO1) > Abs(FO2 - FO1) Then
                AM:=Abs(FN2 - FO1)
            Else
                AM:=Abs(FO2 - FO1);
            DT:=AM;
            If DT > D Then D:=DT;
            FO2:=FN2;
            J2:=J2 + 1;
        end;
        GoTo 1;
    end;
    PROB:=PROBKS(Sqrt(EN1 * EN2 / (EN1 + EN2)) * D);
end;
end.