unit Unit2;

interface
uses
  unit1;
procedure MDIAN1(var X:array of real; N:integer;var XMED:real);
Function GASDEV:real;
procedure MDIAN2(X:array of real; N:integer;var XMED:real);

implementation
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;

procedure MDIAN1(var X:array of real; N:integer;var XMED:real);
var
    N2:integer;
begin
    SORT(N, X);
    N2:=N div 2;
    If 2 * N2 = N Then
        XMED:=0.5 * (X[N2] + X[N2 + 1])
    Else
        XMED:=X[N2 + 1];
end;

procedure MDIAN2(X:array of real; N:integer;var XMED:real);
label 1;
var
    BIG,AMP,AFAC,A,AP,AM,SUM,SUMX,EPS,XP,XM,XX,DUM,BBB,AA:real;
    J,NP,NM:integer;
begin
    BIG:=1E+30;
    AFAC:=1.5;
    AMP:=1.5;
    A:=0.5 * (X[1] + X[N]);
    EPS:=Abs(X[N] - X[1]);
    AP:=BIG;
    AM:=-BIG;
1:  Sum:=0; 
    SUMX:=0; 
    NP:=0;
    NM:=0;
    XP:=BIG;
    XM:=-BIG;
    For J:=1 To N do
    begin
        XX:=X[J];
        If XX <> A Then
        begin
            If XX > A Then
            begin
                NP:=NP + 1;
                If XX < XP Then XP:=XX;
            end
            Else If XX < A Then
            begin
                NM:=NM + 1;
                If XX > XM Then XM:=XX;
            end;
            DUM:=1  / (EPS + Abs(XX - A));
            Sum:=Sum + DUM;
            SUMX:=SUMX + XX * DUM;
        end;
    end; 
    If NP - NM >= 2 Then
    begin
        AM:=A;
        If SUMX / Sum - A < 0 Then
            BBB:=0
        Else
            BBB:=SUMX / Sum - A;
        AA:=XP + BBB * AMP;
        If AA > AP Then AA:=0.5 * (A + AP);
        EPS:=AFAC * Abs(AA - A);
        A:=AA;
        GoTo 1;
    end
    Else If NM - NP >= 2 Then
    begin
        AP:=A;
        If SUMX / Sum - A < 0 Then
            BBB:=SUMX / Sum - A
        Else
            BBB:=0;
        AA:=XM + BBB * AMP;
        If AA < AM Then AA:=0.5 * (A + AM);
        EPS:=AFAC * Abs(AA - A);
        A:=AA;
        GoTo 1;
    end
    Else
    begin
        If (N Mod 2) = 0 Then
        begin
            If NP = NM Then
                XMED:=0.5 * (XP + XM)
            Else If NP > NM Then
                XMED:=0.5 * (A + XP)
            Else
                XMED:=0.5 * (XM + A);
        end
        Else
        begin
            If NP = NM Then
                XMED:=A
            Else If NP > NM Then
                XMED:=XP
            Else
                XMED:=XM;
        end;
    end;
end;

end.