1 #include "data.table.h"
2
3 // from good ol' Numerical Recipes in C
4
iswap(int * a,int * b)5 static inline void iswap(int *a, int *b) {int tmp=*a; *a=*b; *b=tmp;}
dswap(double * a,double * b)6 static inline void dswap(double *a, double *b) {double tmp=*a; *a=*b; *b=tmp;}
i64swap(int64_t * a,int64_t * b)7 static inline void i64swap(int64_t *a, int64_t *b) {int64_t tmp=*a; *a=*b; *b=tmp;}
8
9 #undef BODY
10 #define BODY(SWAP) \
11 if (n==0) return NA_REAL; \
12 int med = n/2 - (n%2==0); \
13 unsigned long ir=n-1, l=0; \
14 for(;;) { \
15 if (ir <= l+1) { \
16 if (ir == l+1 && x[ir] < x[l]) { \
17 SWAP(x+l, x+ir); \
18 } \
19 break; \
20 } else { \
21 unsigned long mid=(l+ir) >> 1; \
22 SWAP(x+mid, x+l+1); \
23 if (x[l] > x[ir]) { \
24 SWAP(x+l, x+ir); \
25 } \
26 if (x[l+1] > x[ir]) { \
27 SWAP(x+l+1, x+ir); \
28 } \
29 if (x[l] > x[l+1]) { \
30 SWAP(x+l, x+l+1); \
31 } \
32 unsigned long i=l+1, j=ir; \
33 a=x[l+1]; \
34 for (;;) { \
35 do i++; while (x[i] < a); \
36 do j--; while (x[j] > a); \
37 if (j < i) break; \
38 SWAP(x+i, x+j); \
39 } \
40 x[l+1]=x[j]; \
41 x[j]=a; \
42 if (j >= med) ir=j-1; \
43 if (j <= med) l=i; \
44 } \
45 } \
46 a = x[med]; \
47 if (n%2 == 1) { \
48 return (double)a; \
49 } else { \
50 b = x[med+1]; \
51 for (int i=med+2; i<n; ++i) { \
52 if (x[i]<b) b=x[i]; \
53 } \
54 return ((double)a+(double)b)/2.0; \
55 }
56
dquickselect(double * x,int n)57 double dquickselect(double *x, int n) {
58 double a, b;
59 BODY(dswap);
60 }
61
iquickselect(int * x,int n)62 double iquickselect(int *x, int n) {
63 int a, b;
64 BODY(iswap);
65 }
66
i64quickselect(int64_t * x,int n)67 double i64quickselect(int64_t *x, int n) {
68 int64_t a, b;
69 BODY(i64swap);
70 }
71
72