#include "../core.h" #include "../utils/hash.h" #include "../utils/talloc.h" 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" u8 radix_offsets_2_u32(usz* c0, u32* v0, usz n) { usz rx = 256; usz* c1 = c0 + rx; // Count keys for (usz j=0; j<2*rx+1; j++) c0[j] = 0; for (usz i=0; i>16)]++; (c1+1)[(u8)(v>>24)]++; } u32 v=v0[0]; // Inclusive prefix sum; note c offsets above if ((c1+1)[(u8)(v>>24)] < n) { c1[0]-=n; RADIX_SUM_2_u32; return 2; } if ((c0+1)[(u8)(v>>16)] < n) { RADIX_SUM_1_u32; return 1; } return 0; } #undef GRADE_UD #define RADIX_LOOKUP_32(INIT, SETTAB) \ u8 bytes = radix_offsets_2_u32(c0, v0, n); \ usz tim = tn/(64/sizeof(*tab)); /* sparse table init max */ \ if (bytes==0) { \ if (n>16); usz c=c0[k]++; v1[c]=v; } \ if (bytes==1) { \ /* Table lookup, getting radix boundaries from c0 */ \ for (usz i=0; i>16)]; \ if (l-i < tim) for (usz ii=i; ii>24); usz c=c1[k]++; v2[c]=v; } \ /* Table lookup */ \ u32 tv=v2[0]>>16; v2[n]=~v2[n-1]; \ for (usz l=0, i=0; l>16; if (tv!=t0) break; tab[(u16)v]=INIT; } \ for (; 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"); u64 n = *SH(x); if (n<=1) { decG(x); return n ? taga(arr_shVec(allOnes(1))) : emptyIVec(); } u8 lw = cellWidthLog(x); void* xv = tyany_ptr(x); if (lw == 0) { usz i = bit_find(xv, n, 1 &~ *(u64*)xv); decG(x); B r = taga(arr_shVec(allZeroes(n))); u64* rp = tyany_ptr(r); rp[0]=1; if (i8 && n=(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 = rp; TALLOC(u8, alloc, 6*n+(4+(tn>3*n?tn:3*n)+(2*rx+1)*sizeof(usz))); // timeline // Allocations len count radix hash deradix bytes layout: usz *c0 = (usz*)(alloc)+1; // rx [+++................] c0 rx # usz *c1 = (usz*)(c0+rx); // rx [++................] c1 rx # u8 *k0 = (u8 *)(c1+rx); // n [+.............] k0 n ## u32 *v2 = (u32*)(k0+n); // n+1 [+.......] v2 4*n+4 ######## u8 *k1 = (u8 *)(v2+n+1); // n [+............] k1 n ## u32 *v1 = (u32*)(k1); // n [+-] v1 4*n ######## u8 *r2 = (u8 *)(v2); // n [+.....] r2 n ## u8 *r1 = (u8 *)(k1+n); // n [+..] r1 n ## u8 *tab= (u8 *)(r1); // tn [+] tab tn ##### RADIX_LOOKUP_32(1, =0) 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); u64* rp; B r = m_bitarrv(&rp, n); H_Sb* set = m_Sb(64); SGetU(x) for (usz i = 0; i < n; i++) bitp_set(rp, i, !ins_Sb(&set, GetU(x,i))); free_Sb(set); decG(x); return r; } B count_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⊒: Argument cannot have rank 0"); 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"); u8 lw = cellWidthLog(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; i8 && n=(1<<14)? 3 : 5) // Radix-assisted lookup usz rx = 256, tn = 1<<16; // Radix; table length u32* v0 = (u32*)xv; i32* r0 = rp; TALLOC(u8, alloc, 6*n+(4+4*(tn>n?tn:n)+(2*rx+1)*sizeof(usz))); // timeline // Allocations len count radix hash deradix bytes layout: usz *c0 = (usz*)(alloc)+1; // rx [+++................] c0 rx # usz *c1 = (usz*)(c0+rx); // rx [++................] c1 rx # u8 *k0 = (u8 *)(c1+rx); // n [+.............] k0 n ## u8 *k1 = (u8 *)(k0+n); // n [+............] k1 n ## u32 *v2 = (u32*)(k1+n); // n+1 [+....-] v2 4*n ######## u32 *v1 = (u32*)(v2+n+1); // n [+..] v1 4*n ######## u32 *r2 = (u32*)v2; // n [+.....] r2 4*n ######## u32 *r1 = (u32*)v1; // n [+..] r1 4*n ######## u32 *tab= (u32*)v1; // tn [+] tab 4*tn ########### 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); i32* rp; B r = m_i32arrv(&rp, n); H_b2i* map = m_b2i(64); SGetU(x) for (usz i = 0; i < n; i++) { bool had; u64 p = mk_b2i(&map, GetU(x,i), &had); rp[i] = had? ++map->a[p].val : (map->a[p].val = 0); } decG(x); free_b2i(map); return r; } B indexOf_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⊐: 𝕩 cannot have rank 0"); 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"); u8 lw = cellWidthLog(x); void* xv = tyany_ptr(x); if (lw == 0) { B r = 1&*(u64*)xv ? not_c1(m_f64(0), x) : 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); H_b2i* map = m_b2i(64); SGetU(x) i32 ctr = 0; for (usz i = 0; i < n; i++) { bool had; u64 p = mk_b2i(&map, GetU(x,i), &had); if (had) rp[i] = map->a[p].val; else rp[i] = map->a[p].val = ctr++; } free_b2i(map); decG(x); return r; } B find_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⍷: Argument cannot have rank 0"); usz n = *SH(x); if (n<=1) return x; return slash_c2(m_f64(0), memberOf_c1(m_f64(0), inc(x)), x); }