-
Notifications
You must be signed in to change notification settings - Fork 1k
Expand file tree
/
Copy pathshellsort.c
More file actions
71 lines (68 loc) · 2.51 KB
/
shellsort.c
File metadata and controls
71 lines (68 loc) · 2.51 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include "data.table.h"
/* shellsort's gap sequences
# Sedgewick 1982
sedgewick1982 = function(k) 4^k+3*2^(k-1)+1
# Sedgewick 1986
sedgewick1986 = function(k) if (k%%2) 8*2^k-6*2^((k+1)/2)+1 else 9*(2^k-2^(k/2))+1
# Tokuda 1992
tokuda = function(k, y=2.25) ceiling((y^k-1)/(y-1))
# Ciura 2001
ciura2001 = function(k) if (k==1) 1 else if (k==2) 4 else if (k==3) 10 else if (k==4) 23 else if (k==5) 57 else if (k==6) 132 else if (k==7) 301 else if (k==8) 701 else floor(2.25 * ciura(k-1))
# Lee 2021
lee = function(k, y=2.243609061420001) tokuda(k, y=y)
*/
static const int sedgewick1982[17] = {
1073790977, 268460033, 67121153, 16783361, 4197377, 1050113,
262913, 65921, 16577, 4193, 1073, 281, 77, 23, 8, 1, 0
};
#define NGAPS 16
/* shellsort
* uses sedgewick1982 gap sequence formula as it turned out to be the most efficient - tested various n, k (for rnorm only!)
* currently used in frollmedian algo="fast" when no NAs
*/
void shellsort(const double *x, int n, int *o) {
for (int i=0; i < n; i++) o[i] = i;
int gap = 0;
while (sedgewick1982[gap] > n) gap++;
for (int h=sedgewick1982[gap]; gap<NGAPS; h=sedgewick1982[++gap]) {
for (int i=h; i<n; i++) {
int io = o[i];
int j = i;
while (j >= h && (x[o[j-h]] > x[io] || (x[o[j-h]] == x[io] && /*this makes stable sort*/ o[j-h] > io))) {
o[j] = o[j-h];
j -= h;
}
o[j] = io;
}
}
}
/* shellsort with NA support
* we support NAs by replacing them by Inf, so they are at the end
* to still support Inf in the data we use isna mask to check if Inf was really Inf or it was NA
*
* not used anymore as frollmedian algo="exact" switched from shellshort to quickselect
* could be used when algo="fast" will have NA support
*
* NAs has to be decoded to Inf inside x
* isna mask must be computed
*/
// void shellsortna(double *x, int n, int *o, bool *isna) {
// for (int i=0; i<n; i++) o[i] = i;
// int gap = 0;
// while (sedgewick1982[gap] > n) gap++;
// for (int h=sedgewick1982[gap]; gap<NGAPS; h=sedgewick1982[++gap]) {
// for (int i=h; i<n; i++) {
// int io = o[i];
// int j = i;
// while (j >= h && (
// (x[o[j-h]] > x[io]) || (
// (x[o[j-h]] == x[io] && (/*both not NA*/(!isna[o[j-h]] && !isna[io]) || /*both NA*/(isna[o[j-h]] && isna[io])) && /*this makes stable sort*/ o[j-h] > io)
// || (/*NA and real R_PosInf*/x[o[j-h]] == x[io] && isna[o[j-h]] && !isna[io])
// ))) {
// o[j] = o[j-h];
// j -= h;
// }
// o[j] = io;
// }
// }
// }