diff --git a/src/builtins/selfsearch.c b/src/builtins/selfsearch.c index 2a6fade8..bc5aa5ef 100644 --- a/src/builtins/selfsearch.c +++ b/src/builtins/selfsearch.c @@ -4,6 +4,40 @@ B not_c1(B t, B x); B shape_c1(B t, B x); +B slash_c2(B t, B w, B x); + +#if defined(__SSE4_2__) +static inline u32 hash32(u32 x) { return _mm_crc32_u32(0x973afb51, x); } +#else +// Murmur3 +static inline u32 hash32(u32 x) { + x ^= x >> 16; x *= 0x85ebca6b; + x ^= x >> 13; x *= 0xc2b2ae35; + x ^= x >> 16; + return x; +} +#endif +static inline u64 hash64(u64 x) { + x ^= x >> 33; x *= 0xff51afd7ed558ccd; + x ^= x >> 33; x *= 0xc4ceb9fe1a85ec53; + x ^= x >> 33; + return x; +} + +static bool canCompare64_norm(B x, usz n) { + u8 e = TI(x,elType); + if (e == el_B) return 0; + if (e == el_f64) { + f64* pf = f64any_ptr(x); + u64* pu = (u64*)pf; + for (usz i = 0; i < n; i++) { + f64 f = pf[i]; + if (f!=f) return 0; + if (pu[i] == m_f64(-0.0).u) return 0; + } + } + return 1; +} #define GRADE_UD(U,D) U #include "radix.h" @@ -54,9 +88,76 @@ u8 radix_offsets_2_u32(usz* c0, u32* v0, usz n) { } \ decG(x); TFREE(alloc); +static NOINLINE void memset32(u32* p, u32 v, usz l) { for (usz i=0; i20) msl=20; \ + usz sh = W - (msl<14? msl : 12+(msl&1)); /* Shift to fit to table */ \ + usz sz = 1 << (W - sh); /* Initial size */ \ + usz msz = 1ull << msl; /* Max sz */ \ + usz b = 64; /* Block size */ \ + /* Resize or abort if more than 1/2^thresh collisions/element */ \ + usz thresh = THRESH; \ + /* Filling e slots past the end requires e*(e+1)/2 collisions, so */ \ + /* n entries with <2 each can fill b? i+b : n; \ + for (; i < e; i++) { \ + T h = hash##W(xp[i]); T j0 = h>>sh, j = j0; T k; \ + while (k=hash[j], k!=h & k!=x0) j++; \ + cc += j-j0; \ + RESWRITE \ + } \ + if (i == n) break; \ + i64 dc = (i64)cc - ((THRESHMUL*i)>>thresh); \ + if (dc >= 0) { \ + if (sz == msz || (RAD && i < n/2 && sz >= 1<<18)) break; /*Abort*/ \ + /* Avoid resizing if close to the end */ \ + if (cc>(5+log-(W-sh))) continue;\ + /* Resize hash, factor of 4 */ \ + usz m = 2; \ + usz dif = sz*((1<>sh, k = k0; while (hash[k]!=x0) k++; \ + cc += k-k0; \ + hash[k] = h; AUXMOVE \ + } \ + if (cc >= STOP) break; \ + thresh = THRESH; \ + } \ + } \ + TFREE(halloc); \ + if (i==n) { decG(x); return RESULT; } +#define SELFHASHTAB_VAL(T, W, RAD, STOP, RES0, RESULT, RESWRITE, THRESHMUL, THRESH, INIT) \ + SELFHASHTAB(T, W, RAD, STOP, RES0, RESULT, RESWRITE, THRESHMUL, THRESH, \ + /*AUXSIZE*/sizeof(u32), \ + /* AUXINIT */ \ + u32* val = (u32*)(hash+sz+ext) + msz-sz; \ + memset32(val, 0, sz+ext); \ + INIT , \ + /*AUXEXTEND*/val -= dif; memset32(val, 0, dif); , \ + /*AUXMOVE*/u32 v = val[j]; val[j] = 0; val[k] = v;) + B memberOf_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("∊: Argument cannot have rank 0"); - usz n = *SH(x); + u64 n = *SH(x); if (n<=1) { decG(x); return n ? taga(arr_shVec(allOnes(1))) : emptyIVec(); } u8 lw = cellWidthLog(x); @@ -69,14 +170,15 @@ B memberOf_c1(B t, B x) { return r; } #define BRUTE(T) \ - i##T* xp = xv; \ - u64* rp; B r = m_bitarrv(&rp, n); bitp_set(rp, 0, 1); \ - for (usz i=1; i=(1<<15)? 3 : 5) + + // Radix-assisted lookup when hash table gives up usz rx = 256, tn = 1<<16; // Radix; table length u32* v0 = (u32*)xv; - i8* r0; B r = m_i8arrv(&r0, n); + i8* r0 = rp; TALLOC(u8, alloc, 6*n+(4+(tn>3*n?tn:3*n)+(2*rx+1)*sizeof(usz))); // timeline @@ -111,8 +220,15 @@ B memberOf_c1(B t, B x) { u8 *tab= (u8 *)(r1); // tn [+] tab tn ##### RADIX_LOOKUP_32(1, =0) - return num_squeeze(r); + return taga(cpyBitArr(r)); } + if (lw == 6 && canCompare64_norm(x, n)) { + if (n<20) { BRUTE(64); } + i8* rp; B r = m_i8arrv(&rp, n); + HASHTAB(u64, 64, 0, n, sz==msz? 0 : sz>=(1<<18)? 0 : sz>=(1<<14)? 3 : 5) + decG(r); // Fall through + } + #undef HASHTAB #undef BRUTE if (RNK(x)>1) x = toCells(x); @@ -126,7 +242,7 @@ B memberOf_c1(B t, B x) { B count_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⊒: Argument cannot have rank 0"); - usz n = *SH(x); + u64 n = *SH(x); if (n<=1) { decG(x); return n ? taga(arr_shVec(allZeroes(1))) : emptyIVec(); } if (n>(usz)I32_MAX+1) thrM("⊒: Argument length >2⋆31 not supported"); @@ -134,14 +250,14 @@ B count_c1(B t, B x) { if (lw==0) { x = toI8Any(x); lw = cellWidthLog(x); } void* xv = tyany_ptr(x); #define BRUTE(T) \ - i##T* xp = xv; \ - i8* rp; B r = m_i8arrv(&rp, n); rp[0]=0; \ - for (usz i=1; i=(1<<14)? 3 : 5) // Radix-assisted lookup usz rx = 256, tn = 1<<16; // Radix; table length u32* v0 = (u32*)xv; - i32* r0; B r = m_i32arrv(&r0, n); + i32* r0 = rp; TALLOC(u8, alloc, 6*n+(4+4*(tn>n?tn:n)+(2*rx+1)*sizeof(usz))); // timeline @@ -178,6 +304,13 @@ B count_c1(B t, B x) { RADIX_LOOKUP_32(0, ++) return num_squeeze(r); } + if (lw == 6 && canCompare64_norm(x, n)) { + if (n<20) { BRUTE(64); } + i32* rp; B r = m_i32arrv(&rp, n); + HASHTAB(u64, 64, 0, n, sz==msz? 0 : sz>=(1<<18)? 0 : sz>=(1<<14)? 3 : 5) + decG(r); // Fall through + } + #undef HASHTAB #undef BRUTE if (RNK(x)>1) x = toCells(x); @@ -192,10 +325,9 @@ B count_c1(B t, B x) { return r; } -extern B rt_indexOf; B indexOf_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⊐: 𝕩 cannot have rank 0"); - usz n = *SH(x); + u64 n = *SH(x); if (n<=1) { decG(x); return n ? taga(arr_shVec(allZeroes(1))) : emptyIVec(); } if (n>(usz)I32_MAX+1) thrM("⊐: Argument length >2⋆31 not supported"); @@ -206,15 +338,22 @@ B indexOf_c1(B t, B x) { return shape_c1(m_f64(0), r); } #define BRUTE(T) \ - i##T* xp = xv; \ - i8* rp; B r = m_i8arrv(&rp, n); rp[0]=0; \ - TALLOC(i##T, uniq, n); uniq[0]=xp[0]; \ - for (usz i=1, u=1; i8 && nmax) max = c; } i64 dst = 1 + (max-(i64)min); - if ((dst=(1<<18)? 1 : sz>=(1<<14)? 4 : 6) + decG(r); // Fall through } + if (lw==6 && canCompare64_norm(x, n)) { + if (n<16) { BRUTE(64); } + i32* rp; B r = m_i32arrv(&rp, n); + u64* xp = tyany_ptr(x); + HASHTAB(u64, 64, sz==msz? 0 : sz>=(1<<17)? 1 : sz>=(1<<13)? 4 : 6) + decG(r); // Fall through + } + #undef HASHTAB #undef BRUTE + #undef DOTAB if (RNK(x)>1) x = toCells(x); i32* rp; B r = m_i32arrv(&rp, n); @@ -274,8 +422,6 @@ B indexOf_c1(B t, B x) { return r; } -B slash_c2(B t, B w, B x); -extern B rt_find; B find_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⍷: Argument cannot have rank 0"); usz n = *SH(x);