home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / MDIAN2.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  1.6 KB  |  63 lines
  1. PROCEDURE mdian2(x: narray; n: integer; VAR xmed: real);
  2. (* Programs using routine MDIAN2 must define the type
  3. TYPE
  4.    narray = ARRAY [1..n] OF real;
  5. in the calling routine *)
  6. LABEL 1;
  7. CONST
  8.    big=1.0e30;
  9.    afac=1.5;
  10.    amp=1.5;
  11. VAR
  12.    np,nm,j: integer;
  13.    xx,xp,xm,sumx,sum,eps: real;
  14.    stemp,dum,ap,am,aa,a: real;
  15. BEGIN
  16.    a := 0.5*(x[1]+x[n]);
  17.    eps := abs(x[n]-x[1]);
  18.    ap := big;
  19.    am := -big;
  20. 1:   sum := 0.0; sumx := 0.0; np := 0; nm := 0; xp := big; xm := -big;
  21.    FOR j := 1 TO n DO BEGIN
  22.       xx := x[j];
  23.       IF (xx <> a) THEN BEGIN
  24.          IF (xx > a) THEN BEGIN
  25.             np := np+1;
  26.             IF (xx < xp) THEN xp := xx END
  27.          ELSE IF  (xx < a) THEN BEGIN
  28.             nm := nm+1;
  29.             IF (xx > xm) THEN xm := xx
  30.          END;
  31.          dum := 1.0/(eps+abs(xx-a));
  32.          sum := sum+dum;
  33.          sumx := sumx+xx*dum
  34.       END
  35.    END;
  36.    stemp := (sumx/sum)-a;
  37.    IF ((np-nm) >= 2) THEN BEGIN
  38.       am := a;
  39.       IF (stemp < 0.0) THEN aa := xp
  40.       ELSE aa := xp+stemp*amp;
  41.       IF (aa > ap) THEN aa := 0.5*(a+ap);
  42.       eps := afac*abs(aa-a);
  43.       a := aa;
  44.       GOTO 1 END
  45.    ELSE IF ((nm-np) >= 2) THEN BEGIN
  46.       ap := a;
  47.       IF (stemp > 0.0) THEN aa := xm
  48.       ELSE aa := xm+stemp*amp;
  49.       IF (aa < am) THEN aa := 0.5*(a+am);
  50.       eps := afac*abs(aa-a);
  51.       a := aa;
  52.       GOTO 1 END
  53.    ELSE IF (n MOD 2) = 0 THEN BEGIN
  54.       IF (np = nm) THEN xmed := 0.5*(xp+xm)
  55.       ELSE IF (np > nm) THEN xmed := 0.5*(a+xp)
  56.       ELSE xmed := 0.5*(xm+a) END
  57.    ELSE BEGIN
  58.       IF (np = nm) THEN xmed := a
  59.       ELSE IF (np > nm) THEN xmed := xp
  60.       ELSE xmed := xm
  61.    END
  62. END;
  63.