/* p番目の要素を選択する最悪時間O(n)のアルゴリズム、
SELECTのプログラム例 */

#include <stdio.h>
#define N 500                          /* 配列Aの最大サイズ */
/* 関数の宣言 */
int select(int i, int j, int p, int *A);
void quicksort(int i, int j, int *A);
int partition(int i, int j, int a, int *A);
int pivot(int i, int j, int *A);
void swap(int i, int j, int *A);

main()
/* p番目の要素の選択のテストプログラム -- SELECT -- */
{
 int A[N], n, i, p, ap;
 FILE *file;

 file=fopen("selectdata", "r");        /* データファイルの読込 */
 fscanf(file, "%d", &n);
 fscanf(file, "%d", &p); 
 if(n>N)                               /* 配列サイズのチェック */
   {
    printf("Illegal array size n = %d for N = %d\n", n, N); 
    exit(1);
   }
 for(i=0; i<n; i++) fscanf(file, "%d", &A[i]); 
 printf("n = %d\nA = ", n);            /* Aの出力 */
 for(i=0; i<n; i++) printf("%d ", A[i]);
 printf("\n"); 
 ap=select(0, n-1, p, A);              /* SELECTの実行 */
 printf("%d-th smallest element = %d\n", p, ap); /* 結果の出力 */
}

int select(int i, int j, int p, int *A)
/* 配列A[i],...,A[j]からp番目の大きさの要素を出力 */
{
 int m, a, h, k;

 if(j-i+1 < 50)      /* Aのサイズ小の場合 */
   {
    quicksort(i, j, A);
    return(A[i+p-1]);
   }
 else                /* Aのサイズ大の場合 */
   {
    m=(j-i+1)/5;
    for(h=0; h<m; h++)                  /* 5要素ずつmグループに分割 */
      {
       quicksort(i+5*h, i+5*h+4, A);    /* 各グループにクイックソート */
       swap(i+h, i+5*h+2, A);           /* 各グループの中央値を先頭へ */
      }
    a=select(i, i+m-1, (j-i+6)/10, A);  /* a: m個の中央値の中央値 */
    k=partition(i, j, a, A);            /* aによって分割 */
    if(p<=k-i) return(select(i, k-1, p, A)); /* 正しい部分へ再帰的適用 */
    else return(select(k, j, p-k+i, A));
   }
}

void quicksort(int i, int j, int *A)
/* 配列A[i],...,A[j]をクイックソートにより整列 */
{
 int a, pv, k;
 pv=pivot(i, j, A);
 if(pv!=-1)                    /* 軸要素が見つかった場合 */
   {
    a=A[pv];                   /* 軸要素 */
    k=partition(i, j, a, A);   /* 軸要素A[pv]により分割 */
    quicksort(i, k-1, A);      /* 前半の整列 */
    quicksort(k, j, A);        /* 後半の整列 */
   }
 return;
}

int partition(int i, int j, int a, int *A)
/* A[i],...,A[j]を軸要素aによって分割 */
/* 前半はA[i],...,A[k-1]<a, 後半はA[k],...,A[j]>=aを満たす */
{
 int l, r, k;
 l=i; r=j;                     /* lは左から、rは右から */
 while(A[l]<a) l=l+1; while(A[r]>=a) r=r-1;
 while(l<=r)                   /* A[l]>=aとA[r]<aが成立 */
   {
    swap(l, r, A); l=l+1; r=r-1;                /* A[l]とA[r]の交換 */
    while(A[l]<a) l=l+1; while(A[r]>=a) r=r-1;  /* つぎのlとrを計算 */
   }
 k=l;
 return(k);
}

int pivot(int i, int j, int *A)
/* A[i],...,A[j]から軸要素A[pv]を選びpvを出力 */
/* A[pv]はA[i]と最初に異なるA[k]の大きい方；すべて同じならpv=-1 */
{
 int pv, k, h, s0, s1, s2;
 s0=j-i; s1=(s0+1)/2; 
 if(s0%2==0) s2=s1; else s2=s1-1; 
 h=1;
 k=i+s1;                      /* s1とs2の交互のステップ幅でkを増加 */
 while(A[i]==A[k] && h<s0)
   {
    if(h%2==1) k=k+s2; else k=k+s1;  /* 次のk */
    if(k>j) k=k-j+i-1;               /* kがjを越える場合の修正 */
    h=h+1;
   }
 if(A[i]>A[k]) pv=i; else if(A[i]<A[k]) pv=k; else pv=-1;
 return(pv);
}

void swap(int i, int j, int *A)
/* A[i]とA[j]の位置の交換 */
{
 int temp;
 temp=A[i]; A[i]=A[j]; A[j]=temp;
 return;
}