'****************************************************************
'*    Calculate the median value of an array without sorting    *
'* ------------------------------------------------------------ *
'* REFERENCE:                                                   *
'*            "NUMERICAL RECIPES By W.H. Press, B.P. Flannery,  *
'*             S.A. Teukolsky and W.T. Vetterling, Cambridge    *
'*             University Press, 1986" [BIBLI 08].              *
'*                                                              *
'*                                Basic Version By J-P Moreau   *
'*                                     (www.jpmoreau.fr)        *
'* ------------------------------------------------------------ *
'* SAMPLE RUN:  (see input file tmdian1.dat).                   *
'*                                                              *
'* N= 80                                                        *
'*                                                              *
'* Input Table:                                                 *
'*                                                              *
'*  407.8 192.8 851.1 604.4 932.3 799.4 914.5 965.8 453.7 295.1 *
'*  154.5 977.4 410.2 916.2 934.7 504.8 823.2 225.2 456.6  49.0 *
'*  933.5 663.0 335.3 346.6 568.7 956.1 654.7 300.7 379.6 591.9 *
'*  992.9 689.6 644.7 305.4 148.2 257.2 664.6 612.1 713.0  99.7 *
'*   46.5 167.6 984.6 847.2  55.4  82.7 999.0  10.7 877.7 929.4 *
'*  398.1 972.8 874.1 755.1 472.1 122.8 671.4  35.5 128.8  76.8 *
'*  454.2 959.2 510.1 791.3 122.8 176.6 237.9 995.8 548.3 309.8 *
'*  162.6 996.5 750.0 250.6 577.7 761.1 101.9 797.1 539.0 723.5 *
'*                                                              *
'* Median Value: 558.5                                          *
'*                                                              *
'*                                                              *
'****************************************************************
DEFINT I-N

  CLS
  OPEN "c:\temp1\tmdian1.dat" FOR INPUT AS #1

  'read size of table
  INPUT #1, N
  DIM A(N)

  'read table of numbers
  FOR i = 1 TO N
    INPUT #1, A(i)
  NEXT i

  CLOSE #1

  PRINT
  PRINT " N="; N
  PRINT
  PRINT " Input Table:"
  GOSUB 2000 'call TWRIT(N,A)

  'call MDIAN1 subroutine
  GOSUB 1000 'call MDIAN1(A,N,AMED)

  PRINT
  PRINT " Median Value: "; AMED
  PRINT
  INPUT "", R$

END


'*******************************************************
'* Given an array X of N numbers, returns their median *
'* value XMED. The array X is not modified, and is     *
'* accessed sequentially in each consecutive pass.     *
'*******************************************************
1000 'Subroutine MDIAN1(A, N, AMED)
  BIG = 1.2E+16: AFAC = 1.5: AMP = 1.5
' Here, AMP is an overconvergence factor: on each iteration,
' we move the guess by this factor. AFAC is a factor used to
' optimize the size of the "smoothing constant" EPS at each
' iteration.
' A1,AA,AM,AP,DUM,EPS,SUM,SUMX,XM,XP,XX:Real
' J,NM,NP:Integer

  A1 = .5 * (A(1) + A(N))
  EPS = ABS(A(N) - A(1))
  AP = BIG
  AM = -BIG
1 SUM = 0#
  SUMX = 0#
  NP = 0
  NM = 0
  XP = BIG
  XM = -BIG
  FOR J = 1 TO N
    XX = A(J)
    IF XX <> A1 THEN
      IF XX > A1 THEN
	NP = NP + 1
	IF XX < XP THEN XP = XX
      ELSEIF XX < A1 THEN
	NM = NM + 1
	IF XX > XM THEN XM = XX
      END IF
      DUM = 1 / (EPS + ABS(XX - A1))
      SUM = SUM + DUM
      SUMX = SUMX + XX * DUM
    END IF
  NEXT J
  IF NP - NM >= 2 THEN     'guess is too low, make another pass
    AM = A1
    AA = XP + MAX(0!, SUMX / SUM - A1) * AMP
    IF AA > AP THEN AA = .5 * (A1 + AP)
    EPS = AFAC * ABS(AA - A1)
    A1 = AA
    GOTO 1
  ELSEIF NM - NP >= 2 THEN 'guess is too high
    AM = A1
    AA = XM + MIN(0!, SUMX / SUM - A1) * AMP
    IF AA < AM THEN AA = .5 * (A1 + AM)
    EPS = AFAC * ABS(AA - A1)
    A1 = AA
    GOTO 1
  ELSE                     'guess is ok
    IF (N MOD 2) = 0 THEN  'for even N median is an average
      IF NP = NM THEN
	AMED = .5 * (XP + XM)
      ELSEIF NP > NM THEN
	AMED = .5 * (A1 + XP)
      ELSE
	AMED = .5 * (XM + A1)
      END IF
    ELSE
      IF NP = NM THEN
	AMED = A1
      ELSEIF NP > NM THEN
	AMED = XP
      ELSE
	AMED = XM
      END IF
    END IF
  END IF
RETURN

'write table of size N to standard output
2000 'SUBROUTINE TWRIT(N,A)
  F$ = " ####.#"
  PRINT
  FOR i = 1 TO N
    PRINT USING F$; A(i);
    IF INT(i / 10) = i / 10 THEN PRINT
  NEXT i
RETURN

'end of file tmdian1.bas