#include "k.h"

double median(double *x,int *n)
{
    /* System generated locals */
    int i_1;

    /* Local variables */
    static double hold,xmed,*y;
    static int nmid;
    static int i, iflag;
    static int nmidp1, iupper;

    /* Parameter adjustments */
    --x;

    /* Function Body */

/*     PURPOSE--THIS SUBROUTINE COMPUTES THE */
/*              SAMPLE MEDIAN */
/*              OF THE DATA IN THE INPUT VECTOR X. */
/*              THE SAMPLE MEDIAN = THAT VALUE SUCH THAT HALF THE */
/*              DATA SET IS BELOW IT AND HALF ABOVE IT. */
/*     INPUT  ARGUMENTS--X      = THE VECTOR OF */
/*                                (UNSORTED OR SORTED) OBSERVATIONS. */
/*                     --N      = THE INTEGER NUMBER OF OBSERVATIONS */
/*                                IN THE VECTOR X. */
/*     OUTPUT--THE COMPUTED VALUE OF THE */
/*             SAMPLE MEDIAN. */
/*     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N */
/*                   FOR THIS SUBROUTINE IS 15000. */
/*     OTHER DATAPAC   SUBROUTINES NEEDED--SORT. */
/*     REFERENCES--KENDALL AND STUART, THE ADVANCED THEORY OF */
/*                 STATISTICS, VOLUME 1, EDITION 2, 1963, PAGE 326. */
/*               --KENDALL AND STUART, THE ADVANCED THEORY OF */
/*                 STATISTICS, VOLUME 2, EDITION 1, 1961, PAGE 49. */
/*               --DAVID, ORDER STATISTICS, 1970, PAGE 139. */
/*               --SNEDECOR AND COCHRAN, STATISTICAL METHODS, */
/*                 EDITION 6, 1967, PAGE 123. */
/*               --DIXON AND MASSEY, INTRODUCTION TO STATISTICAL */
/*                 ANALYSIS, EDITION 2, 1957, PAGE 70. */
/*     WRITTEN BY--JAMES J. FILLIBEN */
/*                 STATISTICAL ENGINEERING LABORATORY (205.03) */
/*                 NATIONAL BUREAU OF STANDARDS */
/*                 WASHINGTON, D. C. 20234 */
/*                 PHONE:  301-921-2315 */
/*     ORIGINAL VERSION--JUNE      1972. */
/*     UPDATED         --SEPTEMBER 1975. */
/*     UPDATED         --NOVEMBER  1975. */
/*     UPDATED         --FEBRUARY  1976. */

/* --------------------------------------------------------------------- 
*/
/*     CHECK THE INPUT ARGUMENTS FOR ERRORS */
    if (*n < 1) {
	goto L50;
    }
    if (*n == 1) {
	goto L55;
    }
    hold = x[1];
    i_1 = *n;
    for (i = 2; i <= i_1; ++i) {
	if (x[i] != hold) {
	    goto L90;
	}
/* L60: */
    }
    xmed = x[1];
    goto L101;
L50:
	printf("Invalid vector length in median routine");
	exit(1);
L55:
    xmed = x[1];
    goto L101;
L90:

/* -----START POINT----------------------------------------------------- 
*/
	if ((y=(double *)malloc((*n)*sizeof(double)))==NULL) {
  		printf("I can't allocate memory: median routine");
 		exit(1);
	}
    sort(&x[1], n, y);
    iflag = *n - (*n / 2 << 1);
    nmid = *n / 2;
    nmidp1 = nmid + 1;
    if (iflag == 0) {
	xmed = (y[nmid - 1] + y[nmidp1 - 1]) / 2.;
    }
    if (iflag == 1) {
	xmed = y[nmidp1 - 1];
    }

L101:
	return xmed;
}

void sort(double *x,int *n,double *y)
{
    /* Local variables */
    static double amed, hold, med;
    static int i, j, k, l, m, il[36], iu[36],i_1;
    static double tt;
    static int ip1, nm1, mid, jmi, lmi, jmk, ipr;

    /* Parameter adjustments */
    --x;
    --y;

    /* Function Body */

/*     PURPOSE--THIS SUBROUTINE SORTS (IN ASCENDING ORDER) */
/*              THE N ELEMENTS OF THE VECTOR X */
/*              AND PUTS THE RESULTING N SORTED VALUES INTO THE */
/*              VECTOR Y. */
/*     INPUT  ARGUMENTS--X      = THE VECTOR OF */
/*                                OBSERVATIONS TO BE SORTED. */
/*                     --N      = THE INTEGER NUMBER OF OBSERVATIONS */
/*                                IN THE VECTOR X. */
/*     OUTPUT ARGUMENTS--Y      = THE VECTOR */
/*                                INTO WHICH THE SORTED DATA VALUES */
/*                                FROM X WILL BE PLACED. */
/*     OUTPUT--THE VECTOR Y */
/*             CONTAINING THE SORTED */
/*             (IN ASCENDING ORDER) VALUES */
/*             OF THE VECTOR X. */
/*     RESTRICTIONS--THE DIMENSIONS OF THE VECTORS IL AND IU */
/*                   (DEFINED AND USED INTERNALLY WITHIN */
/*                   THIS SUBROUTINE) DICTATE THE MAXIMUM */
/*                   ALLOWABLE VALUE OF N FOR THIS SUBROUTINE. */
/*                   IF IL AND IU EACH HAVE DIMENSION K, */
/*                   THEN N MAY NOT EXCEED 2**(K+1) - 1. */
/*                   FOR THIS SUBROUTINE AS WRITTEN, THE DIMENSIONS */
/*                   OF IL AND IU HAVE BEEN SET TO 36, */
/*                   THUS THE MAXIMUM ALLOWABLE VALUE OF N IS */
/*                   APPROXIMATELY 137 BILLION. */
/*                   SINCE THIS EXCEEDS THE MAXIMUM ALLOWABLE */
/*                   VALUE FOR AN INTEGER VARIABLE IN MANY COMPUTERS, */
/*                   AND SINCE A SORT OF 137 BILLION ELEMENTS */
/*                   IS PRESENTLY IMPRACTICAL AND UNLIKELY, */
/*                   THEN THERE IS NO PRACTICAL RESTRICTION */
/*                   ON THE MAXIMUM VALUE OF N FOR THIS SUBROUTINE. */
/*                   (IN LIGHT OF THE ABOVE, NO CHECK OF THE */
/*                   UPPER LIMIT OF N HAS BEEN INCORPORATED */
/*                   INTO THIS SUBROUTINE.) */
/*     COMMENT--IF THE ANALYST DESIRES A SORT 'IN PLACE', */
/*              THIS IS DONE BY HAVING THE SAME */
/*              OUTPUT VECTOR AS INPUT VECTOR IN THE CALLING SEQUENCE. */
/*              THUS, FOR EXAMPLE, THE CALLING SEQUENCE */
/*              CALL SORT(X,N,X) */
/*              IS ALLOWABLE AND WILL RESULT IN */
/*              THE DESIRED 'IN-PLACE' SORT. */
/*     COMMENT--THE SORTING ALGORTHM USED HEREIN */
/*              IS THE BINARY SORT. */
/*              THIS ALGORTHIM IS EXTREMELY FAST AS THE */
/*              FOLLOWING TIME TRIALS INDICATE. */
/*              THESE TIME TRIALS WERE CARRIED OUT ON THE */
/*              UNIVAC 1108 EXEC 8 SYSTEM AT NBS */
/*              IN AUGUST OF 1974. */
/*              BY WAY OF COMPARISON, THE TIME TRIAL VALUES */
/*              FOR THE EASY-TO-PROGRAM BUT EXTREMELY */
/*              INEFFICIENT BUBBLE SORT ALGORITHM HAVE */
/*              ALSO BEEN INCLUDED-- */
/*              NUMBER OF RANDOM        BINARY SORT       BUBBLE SORT */
/*               NUMBERS SORTED */
/*                N = 10                 .002 SEC          .002 SEC */
/*                N = 100                .011 SEC          .045 SEC */
/*                N = 1000               .141 SEC         4.332 SEC */
/*                N = 3000               .476 SEC        37.683 SEC */
/*                N = 10000             1.887 SEC      NOT COMPUTED */
/*     REFERENCES--CACM MARCH 1969, PAGE 186 (BINARY SORT ALGORITHM */
/*                 BY RICHARD C. SINGLETON). */
/*               --CACM JANUARY 1970, PAGE 54. */
/*               --CACM OCTOBER 1970, PAGE 624. */
/*               --JACM JANUARY 1961, PAGE 41. */
/*     WRITTEN BY--JAMES J. FILLIBEN */
/*                 STATISTICAL ENGINEERING LABORATORY (205.03) */
/*                 NATIONAL BUREAU OF STANDARDS */
/*                 WASHINGTON, D. C. 20234 */
/*                 PHONE--301-921-2315 */
/*     ORIGINAL VERSION--JUNE      1972. */
/*     UPDATED         --NOVEMBER  1975. */

/* --------------------------------------------------------------------- 
*/
/*     COPY THE VECTOR X INTO THE VECTOR Y */
    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {
	y[i] = x[i];
    }

/*     CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED */

    nm1 = *n - 1;
    i_1 = nm1;
    for (i = 1; i <= i_1; ++i) {
	ip1 = i + 1;
	if (y[i] <= y[ip1]) {
	    goto L200;
	}
	goto L250;
L200:
    ;}
    return;
L250:
    m = 1;
    i = 1;
    j = *n;
L305:
    if (i >= j) {
	goto L370;
    }
L310:
    k = i;
    mid = (i + j) / 2;
    amed = y[mid];
    if (y[i] <= amed) {
	goto L320;
    }
    y[mid] = y[i];
    y[i] = amed;
    amed = y[mid];
L320:
    l = j;
    if (y[j] >= amed) {
	goto L340;
    }
    y[mid] = y[j];
    y[j] = amed;
    amed = y[mid];
    if (y[i] <= amed) {
	goto L340;
    }
    y[mid] = y[i];
    y[i] = amed;
    amed = y[mid];
    goto L340;
L330:
    y[l] = y[k];
    y[k] = tt;
L340:
    --l;
    if (y[l] > amed) {
	goto L340;
    }
    tt = y[l];
L350:
    ++k;
    if (y[k] < amed) {
	goto L350;
    }
    if (k <= l) {
	goto L330;
    }
    lmi = l - i;
    jmk = j - k;
    if (lmi <= jmk) {
	goto L360;
    }
    il[m - 1] = i;
    iu[m - 1] = l;
    i = k;
    ++m;
    goto L380;
L360:
    il[m - 1] = k;
    iu[m - 1] = j;
    j = l;
    ++m;
    goto L380;
L370:
    --m;
    if (m == 0) {
	return;
    }
    i = il[m - 1];
    j = iu[m - 1];
L380:
    jmi = j - i;
    if (jmi >= 11) {
	goto L310;
    }
    if (i == 1) {
	goto L305;
    }
    --i;
L390:
    ++i;
    if (i == j) {
	goto L370;
    }
    amed = y[i + 1];
    if (y[i] <= amed) {
	goto L390;
    }
    k = i;
L395:
    y[k + 1] = y[k];
    --k;
    if (amed < y[k]) {
	goto L395;
    }
    y[k + 1] = amed;
    goto L390;
}