4-byte sorting: insertion and radix

This commit is contained in:
Marshall Lochbaum 2022-08-26 20:54:12 -04:00
parent f946a3de41
commit 9d59dff522

View File

@ -2,24 +2,18 @@
#define LT SORT_UD(<,>)
#define TIM_B SORT_UD(bA,bD)
#define TIM_I SORT_UD(iA,iD)
#define SORT_CMP(W, X) SORT_UD(compare(W, X), compare(X, W))
#define SORT_NAME TIM_B
#define SORT_TYPE B
#include "sortTemplate.h"
#define SORT_CMP(W, X) (SORT_UD((W) - (i64)(X), (X) - (i64)(W)))
#define SORT_NAME TIM_I
#define SORT_TYPE i32
#include "sortTemplate.h"
#define FOR(I,MAX) SORT_UD(for (usz I=0; I<MAX; I++), for (usz I=MAX; I--; ))
#define PRE_UD(K,SL,SRE) \
u64 p##K=s##K; s##K+=cw##K[j]; \
u64 p##K=s##K; s##K+=((u64*)c##K)[j]; \
s##K+=s##K SL 8; s##K+=s##K SL 16; s##K+=s##K SL 32; \
cw##K[j] = p##K|(s##K SL 8); s##K SRE 56
((u64*)c##K)[j] = p##K|(s##K SL 8); s##K SRE 56
#define PRE64(K) SORT_UD(PRE_UD(K,<<,>>=), PRE_UD(K,>>,<<=))
#define INSERTION_SORT(T) \
@ -53,12 +47,12 @@
B SORT_C1(B t, B x) {
if (isAtm(x) || RNK(x)==0) thrM(SORT_UD("","")": Argument cannot have rank 0");
if (RNK(x)!=1) return bqn_merge(SORT_C1(t, toCells(x)));
usz xia = IA(x);
if (xia <= 1) return x;
usz n = IA(x);
if (n <= 1) return x;
u8 xe = TI(x,elType);
B r;
if (xe==el_i8) {
i8* xp = i8any_ptr(x); usz n=xia;
i8* xp = i8any_ptr(x);
i8* rp; r = m_i8arrv(&rp, n);
if (n<16) {
INSERTION_SORT(i8);
@ -66,15 +60,14 @@ B SORT_C1(B t, B x) {
TALLOC(u8, c0, 256); u8 *c0o=c0+128; // Offset for signedness
for (usz j=0; j<256; j++) c0[j]=0;
for (usz i=0; i<n; i++) c0o[xp[i]]++;
u64 s=0; u64 *cw=(u64*)c0;
FOR(j, 256/8) { PRE64(); } // Prefix sum
u64 s0=0; FOR(j, 256/8) { PRE64(0); } // Prefix sum
for (usz i=0; i<n; i++) { i8 xi=xp[i]; u8 c=c0o[xi]++; rp[c]=xi; }
TFREE(c0);
} else {
COUNTING_SORT(i8);
}
} else if (xe==el_i16) {
i16* xp = i16any_ptr(x); usz n=xia;
i16* xp = i16any_ptr(x);
i16* rp; r = m_i16arrv(&rp, n);
if (n < 24) {
INSERTION_SORT(i16);
@ -90,7 +83,7 @@ B SORT_C1(B t, B x) {
for (usz i=0; i<n; i++) { i16 v=r0[i]; T c=c1o[(i8)(v>>8)]++; rp[c]=v; } \
TFREE(alloc)
RADIX2(u8,
u64 s0=0; u64 s1=0; u64 *cw0=(u64*)c0; u64 *cw1=(u64*)c1;
u64 s0=0; u64 s1=0;
FOR(j, 256/8) { PRE64(0); PRE64(1); }
);
} else if (n < 1<<15) { // Radix sort
@ -107,15 +100,46 @@ B SORT_C1(B t, B x) {
}
} else if (xe==el_i32) {
i32* xp = i32any_ptr(x);
i32* rp; r = m_i32arrv(&rp, xia);
memcpy(rp, xp, xia*4);
CAT(TIM_I,tim_sort)(rp, xia);
i32* rp; r = m_i32arrv(&rp, n);
if (n < 40) {
INSERTION_SORT(i32);
} else if (n <= 256) {
#define RADIX4(T, PRE_SUM) \
TALLOC(u8, alloc, 4*256*sizeof(T) + n*4); \
T *c0=(T*)alloc, *c1=c0+256, *c2=c1+256, *c3=c2+256, *c3o=c3+128; \
for (usz j=0; j<4*256; j++) c0[j]=0; \
for (usz i=0; i<n; i++) { i32 v=xp[i]; \
c0 [(u8)v ]++; c1 [(u8)(v>> 8)]++; \
c2 [(u8)(v>>16)]++; c3o[(i8)(v>>24)]++; } \
PRE_SUM; \
i32 *r0 = (i32*)(c0+4*256); \
for (usz i=0; i<n; i++) { i32 v=xp[i]; T c=c0 [(u8)v ]++; r0[c]=v; } \
for (usz i=0; i<n; i++) { i32 v=r0[i]; T c=c1 [(u8)(v>> 8)]++; rp[c]=v; } \
for (usz i=0; i<n; i++) { i32 v=rp[i]; T c=c2 [(u8)(v>>16)]++; r0[c]=v; } \
for (usz i=0; i<n; i++) { i32 v=r0[i]; T c=c3o[(i8)(v>>24)]++; rp[c]=v; } \
TFREE(alloc)
RADIX4(u8,
u64 s0=0; u64 s1=0; u64 s2=0; u64 s3=0;
FOR(j, 256/8) { PRE64(0); PRE64(1); PRE64(2); PRE64(3); }
);
} else {
RADIX4(usz,
usz s0=0; usz s1=0; usz s2=0; usz s3=0;
FOR(j, 256) {
u32 p0=s0; s0+=c0[j]; c0[j]=p0;
u32 p1=s1; s1+=c1[j]; c1[j]=p1;
u32 p2=s2; s2+=c2[j]; c2[j]=p2;
u32 p3=s3; s3+=c3[j]; c3[j]=p3;
}
);
#undef RADIX4
}
} else {
B xf = getFillQ(x);
HArr_p ra = m_harrUv(xia);
HArr_p ra = m_harrUv(n);
SGet(x)
for (usz i = 0; i < xia; i++) ra.a[i] = Get(x,i);
CAT(TIM_B,tim_sort)(ra.a, xia);
for (usz i = 0; i < n; i++) ra.a[i] = Get(x,i);
CAT(TIM_B,tim_sort)(ra.a, n);
r = withFill(ra.b,xf);
}
decG(x);
@ -124,7 +148,6 @@ B SORT_C1(B t, B x) {
#undef SORT_C1
#undef LT
#undef TIM_B
#undef TIM_I
#undef FOR
#undef PRE_UD
#undef PRE64