|
The Tiny HTM Library
|
00001 00009 #include "tinyhtm/select.h" 00010 00011 #ifdef __cplusplus 00012 extern "C" { 00013 #endif 00014 00015 /* Lookup tables for the 4 and 5 element median finding algorithms. 00016 Computed with the following python 2.6+ script (with n = 4, 5): 00017 00018 @code 00019 import itertools 00020 00021 def computeLut(n): 00022 nbits = (n * (n - 1)) / 2 00023 lut = [-1] * 2**nbits 00024 array = range(n) 00025 median = array[len(array) >> 1] 00026 for p in itertools.permutations(array): 00027 res = [] 00028 for i in xrange(n - 1): 00029 for j in xrange(i + 1, n): 00030 res.append(1 if p[i] < p[j] else 0) 00031 index = 0 00032 for i in xrange(len(res)): 00033 index += res[i] << (len(res) - i - 1) 00034 lut[index] = p.index(median) 00035 return lut 00036 @endcode 00037 */ 00038 static const signed char _htm_lut4[64] = { 00039 1, 1,-1, 3, 2,-1, 2, 3,-1,-1,-1, 0,-1,-1,-1, 0, 00040 -1,-1,-1,-1, 0,-1, 0,-1,-1,-1,-1,-1,-1,-1, 3, 2, 00041 0, 0,-1,-1,-1,-1,-1,-1,-1, 3,-1, 1,-1,-1,-1,-1, 00042 2,-1,-1,-1, 1,-1,-1,-1, 2, 3,-1, 1, 1,-1, 3, 2 00043 }; 00044 00045 static const signed char _htm_lut5[1024] = { 00046 2, 2,-1, 4, 3,-1, 3, 4,-1,-1,-1, 1,-1,-1,-1, 1,-1,-1,-1,-1, 1,-1, 1,-1,-1,-1,-1,-1,-1,-1, 4, 3, 00047 1, 1,-1,-1,-1,-1,-1,-1,-1, 4,-1, 2,-1,-1,-1,-1, 3,-1,-1,-1, 2,-1,-1,-1, 3, 4,-1, 2, 2,-1, 4, 3, 00048 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1,-1,-1,-1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 3, 00049 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 2,-1,-1,-1, 3, 00050 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1,-1, 1,-1,-1,-1,-1,-1,-1,-1, 4,-1, 00051 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1, 2,-1, 4,-1, 00052 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0, 0, 00053 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0, 0, 00054 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00055 1, 1,-1,-1,-1,-1,-1,-1,-1, 4,-1,-1,-1,-1,-1,-1, 3,-1,-1,-1,-1,-1,-1,-1, 3, 4,-1,-1,-1,-1,-1,-1, 00056 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00057 -1,-1,-1,-1,-1,-1,-1,-1,-1, 0,-1, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0,-1, 0,-1,-1,-1,-1, 00058 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00059 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0,-1,-1,-1, 0,-1,-1,-1, 0,-1,-1,-1, 0,-1,-1,-1, 00060 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00061 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 4, 3,-1, 3, 4,-1, 2, 2, 00062 2, 2,-1, 4, 3,-1, 3, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00063 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00064 -1,-1,-1, 0,-1,-1,-1, 0,-1,-1,-1, 0,-1,-1,-1, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00065 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00066 -1,-1,-1,-1, 0,-1, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 0,-1, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00067 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00068 -1,-1,-1,-1,-1,-1, 4, 3,-1,-1,-1,-1,-1,-1,-1, 3,-1,-1,-1,-1,-1,-1, 4,-1,-1,-1,-1,-1,-1,-1, 1, 1, 00069 -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00070 0, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00071 0, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00072 -1, 4,-1, 2,-1,-1,-1,-1,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00073 -1, 4,-1,-1,-1,-1,-1,-1,-1, 1,-1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00074 3,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00075 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1,-1,-1,-1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 00076 3, 4,-1, 2, 2,-1, 4, 3,-1,-1,-1, 2,-1,-1,-1, 3,-1,-1,-1,-1, 2,-1, 4,-1,-1,-1,-1,-1,-1,-1, 1, 1, 00077 3, 4,-1,-1,-1,-1,-1,-1,-1, 1,-1, 1,-1,-1,-1,-1, 1,-1,-1,-1, 1,-1,-1,-1, 4, 3,-1, 3, 4,-1, 2, 2 00078 }; 00079 00080 /* Returns the index of the median of 2 doubles. 00081 */ 00082 HTM_INLINE size_t _htm_median2(const double *array) 00083 { 00084 return (array[0] < array[1]) ? 0 : 1; 00085 } 00086 00087 /* Returns the index of the median of 3 doubles. 00088 */ 00089 HTM_INLINE size_t _htm_median3(const double *array) 00090 { 00091 double x0 = array[0]; 00092 double x1 = array[1]; 00093 double x2 = array[2]; 00094 if (x0 < x1) { 00095 return x1 < x2 ? 1 : (x0 < x2 ? 2 : 0); 00096 } else { 00097 return x1 < x2 ? (x0 < x2 ? 0 : 2) : 1; 00098 } 00099 } 00100 00101 /* Returns the index of the median of 4 doubles using a branchless algorithm. 00102 */ 00103 HTM_INLINE size_t _htm_median4(const double *array) 00104 { 00105 double a = array[0]; 00106 double b = array[1]; 00107 double c = array[2]; 00108 double d = array[3]; 00109 /* avoids branches, but always requires 6 comparisons. */ 00110 int i = (((int) (a < b)) << 5) | 00111 (((int) (a < c)) << 4) | 00112 (((int) (a < d)) << 3) | 00113 (((int) (b < c)) << 2) | 00114 (((int) (b < d)) << 1) | 00115 ((int) (c < d)); 00116 return (size_t) _htm_lut4[i]; 00117 } 00118 00119 /* Returns the index of the median of 5 doubles using a branchless algorithm. 00120 */ 00121 HTM_INLINE size_t _htm_median5(const double *array) 00122 { 00123 double a = array[0]; 00124 double b = array[1]; 00125 double c = array[2]; 00126 double d = array[3]; 00127 double e = array[4]; 00128 /* avoids branches, but always performs 10 comparisons */ 00129 int i = (((int) (a < b)) << 9) | 00130 (((int) (a < c)) << 8) | 00131 (((int) (a < d)) << 7) | 00132 (((int) (a < e)) << 6) | 00133 (((int) (b < c)) << 5) | 00134 (((int) (b < d)) << 4) | 00135 (((int) (b < e)) << 3) | 00136 (((int) (c < d)) << 2) | 00137 (((int) (c < e)) << 1) | 00138 ((int) (d < e)); 00139 return (size_t) _htm_lut5[i]; 00140 } 00141 00142 /* Returns the index of the median of medians for an array. 00143 00144 The following pre-conditions are assumed: 00145 - array != 0 00146 - n > 0 00147 */ 00148 static size_t _htm_mm(double *array, size_t n) 00149 { 00150 size_t i, j, m = 0; 00151 while (1) { 00152 if (n <= 5) { 00153 switch (n) { 00154 case 1: m = 0; break; 00155 case 2: m = _htm_median2(array); break; 00156 case 3: m = _htm_median3(array); break; 00157 case 4: m = _htm_median4(array); break; 00158 case 5: m = _htm_median5(array); break; 00159 } 00160 break; 00161 } 00162 for (i = 0, j = 0; i < n - 4; i += 5, ++j) { 00163 size_t m5 = _htm_median5(array + i) + i; 00164 double tmp = array[j]; 00165 array[j] = array[m5]; 00166 array[m5] = tmp; 00167 } 00168 n = j; 00169 } 00170 return m; 00171 } 00172 00173 /* Partitions the given array around the value of the i-th element, and 00174 returns the index of the pivot value after partitioning. 00175 00176 If the array contains many values identical to the pivot value, 00177 lop-sided partitions can be generated by a naive algorithm. In the 00178 worst case (e.g. all elements are identical), an array of n elements will 00179 be partitioned into two sub-arrays of size 1 and n - 1. This leads to 00180 quadratic run-time for selection/sorting algorithms that use the 00181 partitioning primitive. 00182 00183 This implementation therefore counts values identical to the pivot 00184 while partitioning. If the partitions are overly lop-sided, a second 00185 pass over the larger partition is performed. This pass assigns pivot 00186 values to the smaller partition until the sizes are balanced or there 00187 are no more pivot values left. 00188 00189 The run-time of this function is O(n). 00190 00191 The following pre-conditions are assumed: 00192 - array != 0 00193 - n > 0 00194 - i < n 00195 */ 00196 static size_t _htm_wcpart(double *array, size_t n, size_t i) 00197 { 00198 size_t u, v, neq; 00199 double pivot = array[i]; 00200 array[i] = array[n - 1]; 00201 /* partition around pivot */ 00202 for (u = 0, v = 0, neq = 0; v < n - 1; ++v) { 00203 if (array[v] < pivot) { 00204 double tmp = array[u]; 00205 array[u] = array[v]; 00206 array[v] = tmp; 00207 ++u; 00208 } else if (array[v] == pivot) { 00209 ++neq; 00210 } 00211 } 00212 array[n - 1] = array[u]; 00213 array[u] = pivot; 00214 if (neq > 0 && u < (n >> 2)) { 00215 /* lop-sided partition - use values identical 00216 to the pivot value to increase u */ 00217 if (u + neq > (n >> 1)) { 00218 neq = (n >> 1) - u; 00219 } 00220 for (v = u + 1; neq != 0; ++v) { 00221 if (array[v] == pivot) { 00222 ++u; 00223 array[v] = array[u]; 00224 array[u] = pivot; 00225 --neq; 00226 } 00227 } 00228 } 00229 return u; 00230 } 00231 00232 /* Chooses a pivot value using the median-of-3 strategy. 00233 00234 The following pre-conditions are assumed: 00235 - array != 0 00236 - n > 0 00237 */ 00238 static size_t _htm_med3pivot(double *array, size_t n) { 00239 double a, b, c; 00240 size_t m; 00241 if (n <= 5) { 00242 switch (n) { 00243 case 1: return 0; 00244 case 2: return _htm_median2(array); 00245 case 3: return _htm_median3(array); 00246 case 4: return _htm_median4(array); 00247 case 5: return _htm_median5(array); 00248 } 00249 } 00250 m = n >> 1; 00251 a = array[0]; 00252 b = array[m]; 00253 c = array[n - 1]; 00254 if (a < b) { 00255 return b < c ? m : (a < c ? n - 1 : 0); 00256 } else { 00257 return b < c ? (a < c ? 0 : n - 1) : m; 00258 } 00259 } 00260 00261 /* Partitions the given array around the value of the i-th element, and 00262 returns the index of the pivot value after partitioning. 00263 00264 The following pre-conditions are assumed: 00265 - array != 0 00266 - n > 0 00267 - i < n 00268 */ 00269 static size_t _htm_part(double *array, size_t n, size_t i) 00270 { 00271 size_t u, v; 00272 const double pivot = array[i]; 00273 array[i] = array[n - 1]; 00274 for (u = 0, v = 0; v < n - 1; ++v) { 00275 if (array[v] < pivot) { 00276 double tmp = array[u]; 00277 array[u] = array[v]; 00278 array[v] = tmp; 00279 ++u; 00280 } 00281 } 00282 array[n - 1] = array[u]; 00283 array[u] = pivot; 00284 return u; 00285 } 00286 00287 00288 double htm_select(double *array, size_t n, size_t k) 00289 { 00290 size_t tot = 0; 00291 const size_t thresh = n*3; 00292 00293 while (1) { 00294 size_t i = _htm_med3pivot(array, n); 00295 i = _htm_part(array, n, i); 00296 if (k == i) { 00297 break; 00298 } else if (k < i) { 00299 n = i; 00300 } else { 00301 array += i + 1; 00302 n -= i + 1; 00303 k -= i + 1; 00304 } 00305 tot += n; 00306 if (tot > thresh) { 00307 return htm_selectmm(array, n, k); 00308 } 00309 } 00310 return array[k]; 00311 } 00312 00313 00314 double htm_selectmm(double *array, size_t n, size_t k) 00315 { 00316 while (1) { 00317 size_t i = _htm_mm(array, n); 00318 i = _htm_wcpart(array, n, i); 00319 if (k == i) { 00320 break; 00321 } else if (k < i) { 00322 n = i; 00323 } else { 00324 array += i + 1; 00325 n -= i + 1; 00326 k -= i + 1; 00327 } 00328 } 00329 return array[k]; 00330 } 00331 00332 00333 double htm_min(const double *array, size_t n) 00334 { 00335 double m; 00336 size_t i; 00337 00338 m = array[0]; 00339 for (i = 1; i < n; ++i) { 00340 if (array[i] < m) { 00341 m = array[i]; 00342 } 00343 } 00344 return m; 00345 } 00346 00347 00348 #ifdef __cplusplus 00349 } 00350 #endif 00351
1.7.6.1