mdian2.txt

《Delphi常用数值算法集》的配书源码

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;