/* p番目の要素を選択する平均時間O(n)の確率アルゴリズム、
LAZYSELECTのプログラム例 */

#include <stdio.h>
#include <stdlib.h> 
#include <math.h> 
#define N 10000       /* 配列Aの最大サイズ */
#define SEED 13       /* 乱数の種 */ 
/* 関数の宣言 */
int lazyselect(int i, int j, int p, int *A); 
int rand_from(int min, int max); 
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要素の選択のテストプログラム -- LAZYSELECT -- */
{
 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=lazyselect(0, n-1, p, A);          /* LAZYSELECTの実行 */
 printf("%d-th smallest element = %d\n", p, ap); /* 結果の出力 */
}

int lazyselect(int i, int j, int p, int *A) 
/* 配列A[i],A[i+1},...,A[j]内のp番目の大きさの要素を選択する
確率アルゴリズムLAZYSELECT */
{
 int k, nB, rdm, L, H, kL, kH, l, h, n;
 double q;

 n=j-i+1;                        /* nは配列のサイズ */
 if(n<50)                        /* 配列サイズの小さい場合 */
   {
    quicksort(i, j, A);
    return(A[i+p-1]);
   }
 else                            /* 配列サイズの大きい場合 */
   {
    nB=(int)(pow((double)n, 0.75)+0.5);   
                                 /* サイズnBの配列Bの生成 */
    for(k=0; k<nB; k++)          /* Aの前半にBを移動 */
      {
       rdm=rand_from(i+k, j);    /* i+kとjの間の乱整数を生成 */
       swap(i+k, rdm, A);
      }
    q=p*pow((double)n, -0.25);   /* p番目のBにおける推定位置 */
    l=(int)(q-sqrt((double)n));      /* 幅の左端 */
    if(l<1) l=1;
    h=(int) ceil(q+sqrt((double)n)); /* 幅の右端 */
    if(h>nB) h=nB;
    L=lazyselect(i, i+nB-1, l, A);   /* Bの大きさl番目の要素 */
    H=lazyselect(i, i+nB-1, h, A);   /* Bの大きさh番目の要素 */
    kH=partition(i, j, H, A);    /* A1とA2の分割 */
    if(kH<i+p)                   /* p番目はA2にある */
      return(lazyselect(kH, j, p-kH+i, A));
    else                         /* p番目はA1にある */
      {
       kL=partition(i, kH-1, L, A);  /* A11とA12の分割 */
       if(kL>i+p-1)              /* p番目はA11にある */
         return(lazyselect(i, kL-1, p, A));
       else                      /* p番目はA12にある */
         return(lazyselect(kL, kH-1, p-kL+i, A));
      }
   }
}

int rand_from(int min, int max)
/* 区間[min, max]からランダムに整数を選ぶ */
/* randは0からRAND_MAXまでの乱整数を生成する標準関数 */
{
 return((int)(((double)rand()/((unsigned)RAND_MAX+1))
   *(max-min+1))+min);
}

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を満たす */
/* 上記のkを出力 */
{
 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;
}