diff --git a/makefile b/makefile index bf26179b..34d354ab 100644 --- a/makefile +++ b/makefile @@ -212,7 +212,7 @@ ${bd}/%.o: src/jit/%.c @echo $< | cut -c 5- @$(CC_INC) $@.d -o $@ -c $< -builtins: ${addprefix ${bd}/, arithm.o arithd.o cmp.o sfns.o squeeze.o select.o slash.o group.o sort.o selfsearch.o md1.o md2.o fns.o sysfn.o internal.o inverse.o} +builtins: ${addprefix ${bd}/, arithm.o arithd.o cmp.o sfns.o squeeze.o select.o slash.o group.o sort.o search.o selfsearch.o md1.o md2.o fns.o sysfn.o internal.o inverse.o} ${bd}/%.o: src/builtins/%.c @echo $< | cut -c 5- @$(CC_INC) $@.d -o $@ -c $< diff --git a/src/builtins/fns.c b/src/builtins/fns.c index a6916e2e..ccdf0d28 100644 --- a/src/builtins/fns.c +++ b/src/builtins/fns.c @@ -198,120 +198,11 @@ B fne_c2(B t, B w, B x) { } -extern B rt_indexOf; -B indexOf_c2(B t, B w, B x) { - if (!isArr(w) || RNK(w)==0) thrM("⊐: 𝕨 must have rank at least 1"); - if (RNK(w)==1) { - if (!isArr(x) || RNK(x)==0) { - usz wia = IA(w); - B el = isArr(x)? IGetU(x,0) : x; - i32 res = wia; - if (TI(w,elType)==el_i32) { - if (q_i32(el)) { - i32* wp = i32any_ptr(w); - i32 v = o2iG(el); - for (usz i = 0; i < wia; i++) { - if (wp[i] == v) { res = i; break; } - } - } - } else { - SGetU(w) - for (usz i = 0; i < wia; i++) { - if (equal(GetU(w,i), el)) { res = i; break; } - } - } - decG(w); dec(x); - i32* rp; Arr* r = m_i32arrp(&rp, 1); - arr_shAlloc(r, 0); - rp[0] = res; - return taga(r); - } else { - usz wia = IA(w); - usz xia = IA(x); - // TODO O(wia×xia) for small wia or xia - i32* rp; B r = m_i32arrc(&rp, x); - H_b2i* map = m_b2i(64); - SGetU(x) - SGetU(w) - for (usz i = 0; i < wia; i++) { - bool had; u64 p = mk_b2i(&map, GetU(w,i), &had); - if (!had) map->a[p].val = i; - } - for (usz i = 0; i < xia; i++) rp[i] = getD_b2i(map, GetU(x,i), wia); - free_b2i(map); decG(w); decG(x); - return wia<=I8_MAX? taga(cpyI8Arr(r)) : wia<=I16_MAX? taga(cpyI16Arr(r)) : r; - } - } - return c2(rt_indexOf, w, x); -} - -B enclosed_0; -B enclosed_1; -extern B rt_memberOf; -extern B eq_c2(B,B,B); -extern B or_c2(B,B,B); -B memberOf_c2(B t, B w, B x) { - if (isAtm(x) || RNK(x)!=1) goto bad; - if (isAtm(w)) goto single; - ur wr = RNK(w); - if (wr==0) { - B w0 = IGet(w, 0); - dec(w); - w = w0; - goto single; - } - if (wr==1) goto many; - goto bad; - - bad: return c2(rt_memberOf, w, x); - - B r; - single: { - usz xia = IA(x); - SGetU(x) - for (usz i = 0; i < xia; i++) if (equal(GetU(x, i), w)) { r = inc(enclosed_1); goto dec_wx; } - r = inc(enclosed_0); - dec_wx:; dec(w); - goto dec_x; - } - - - many: { - u8 we = TI(w,elType); usz wia = IA(w); - u8 xe = TI(x,elType); usz xia = IA(x); - if (xia<=16 && wia>16 && wef; u64* xp = bitarr_ptr(x); if (rtid==n_add) { B r = m_f64(bit_sum (xp, ia)); decG(x); return r; } if (rtid==n_sub) { B r = m_f64(bit_diff(xp, ia)); decG(x); return r; } - if (rtid==n_and | rtid==n_mul | rtid==n_floor) { bool r=1; for (usz i=0; i<(ia>>6); i++) if (~xp[i]){r=0;break;} if(~bitp_l1(xp,ia))r=0; decG(x); return m_i32(r); } - if (rtid==n_or | rtid==n_ceil ) { bool r=0; for (usz i=0; i<(ia>>6); i++) if ( xp[i]){r=1;break;} if( bitp_l0(xp,ia))r=1; decG(x); return m_i32(r); } + if (rtid==n_and | rtid==n_mul | rtid==n_floor) { B r = m_i32(!bit_has(xp, ia, 0)); decG(x); return r; } + if (rtid==n_or | rtid==n_ceil ) { B r = m_i32( bit_has(xp, ia, 1)); decG(x); return r; } if (rtid==n_ne) { bool r=fold_ne(xp, ia) ; decG(x); return m_i32(r); } if (rtid==n_eq) { bool r=fold_ne(xp, ia) ^ (1&~ia); decG(x); return m_i32(r); } goto base; @@ -383,8 +383,8 @@ B fold_c2(Md1D* d, B w, B x) { B f = d->f; if (rtid==n_add) { B r = m_f64(wi + bit_sum (xp, ia)); decG(x); return r; } if (rtid==n_sub) { B r = m_f64((ia&1?-wi:wi) + bit_diff(xp, ia)); decG(x); return r; } if (wi!=(wi&1)) goto base; - if (rtid==n_and | rtid==n_mul | rtid==n_floor) { bool r=wi; if ( r) { for (usz i=0; i<(ia>>6); i++) if (~xp[i]){r=0;break;} if(~bitp_l1(xp,ia))r=0; } decG(x); return m_i32(r); } - if (rtid==n_or | rtid==n_ceil ) { bool r=wi; if (!r) { for (usz i=0; i<(ia>>6); i++) if ( xp[i]){r=1;break;} if( bitp_l0(xp,ia))r=1; } decG(x); return m_i32(r); } + if (rtid==n_and | rtid==n_mul | rtid==n_floor) { B r = m_i32(wi && !bit_has(xp, ia, 0)); decG(x); return r; } + if (rtid==n_or | rtid==n_ceil ) { B r = m_i32(wi || bit_has(xp, ia, 1)); decG(x); return r; } if (rtid==n_ne) { bool r=wi^fold_ne(xp, ia) ; decG(x); return m_i32(r); } if (rtid==n_eq) { bool r=wi^fold_ne(xp, ia) ^ (1&ia); decG(x); return m_i32(r); } goto base; diff --git a/src/builtins/radix.h b/src/builtins/radix.h new file mode 100644 index 00000000..bcd68f7f --- /dev/null +++ b/src/builtins/radix.h @@ -0,0 +1,42 @@ +#pragma once + +// Radix sorting utilities +// These are leaky macros and assume counts are c0, c1,... +// which must be adjacent in memory. Search functions will also need +// #define GRADE_UD(U,D) U +// to do the appropriate sums for an ascending radix sort. + +#define RDX_PRE(K) s##K=c##K[j]+=s##K +#define RDX_SUM_1(T) T s0=0; for(usz j=0;j<256;j++) { RDX_PRE(0); } +#define RDX_SUM_2(T) GRADE_UD(c1[0]=0;,) T s0=0, s1=0; for(usz j=0;j<256;j++) { RDX_PRE(0); RDX_PRE(1); } +#define RDX_SUM_4(T) GRADE_UD(c1[0]=c2[0]=c3[0]=0;,) T s0=0, s1=0, s2=0, s3=0; for(usz j=0;j<256;j++) { RDX_PRE(0); RDX_PRE(1); RDX_PRE(2); RDX_PRE(3); } + +#if SINGELI +extern void (*const avx2_scan_pluswrap_u8)(uint8_t* v0,uint8_t* v1,uint64_t v2,uint8_t v3); +extern void (*const avx2_scan_pluswrap_u32)(uint32_t* v0,uint32_t* v1,uint64_t v2,uint32_t v3); +#define RADIX_SUM_1_u8 avx2_scan_pluswrap_u8 (c0,c0, 256,0); +#define RADIX_SUM_1_u32 avx2_scan_pluswrap_u32(c0,c0, 256,0); +#define RADIX_SUM_2_u8 avx2_scan_pluswrap_u8 (c0,c0,2*256,0); +#define RADIX_SUM_2_u32 avx2_scan_pluswrap_u32(c0,c0,2*256,0); +#define RADIX_SUM_4_u8 avx2_scan_pluswrap_u8 (c0,c0,4*256,0); +#define RADIX_SUM_4_u32 avx2_scan_pluswrap_u32(c0,c0,4*256,0); +#else +#define RADIX_SUM_1_u8 RDX_SUM_1(u8) +#define RADIX_SUM_1_u32 RDX_SUM_1(u32) +#define RADIX_SUM_2_u8 RDX_SUM_2(u8) +#define RADIX_SUM_2_u32 RDX_SUM_2(u32) +#define RADIX_SUM_4_u8 RDX_SUM_4(u8) +#define RADIX_SUM_4_u32 RDX_SUM_4(u32) +#endif + +#if SINGELI && !USZ_64 +#define RADIX_SUM_1_usz avx2_scan_pluswrap_u32(c0,c0, 256,0); +#define RADIX_SUM_2_usz avx2_scan_pluswrap_u32(c0,c0,2*256,0); +#define RADIX_SUM_4_usz avx2_scan_pluswrap_u32(c0,c0,4*256,0); +#else +#define RADIX_SUM_1_usz RDX_SUM_1(usz) +#define RADIX_SUM_2_usz RDX_SUM_2(usz) +#define RADIX_SUM_4_usz RDX_SUM_4(usz) +#endif + +u8 radix_offsets_2_u32(usz* c0, u32* v0, usz n); // selfsearch.c diff --git a/src/builtins/search.c b/src/builtins/search.c new file mode 100644 index 00000000..bfa805f8 --- /dev/null +++ b/src/builtins/search.c @@ -0,0 +1,245 @@ +#include "../core.h" +#include "../utils/hash.h" +#include "../utils/mut.h" +#include "../utils/talloc.h" + +#define C2(F,X,W) F##_c2(m_f64(0),X,W) +#define C2i(F, W, X) C2(F, m_i32(W), X) +extern B eq_c2(B,B,B); +extern B ne_c2(B,B,B); +extern B or_c2(B,B,B); +extern B add_c2(B,B,B); +extern B sub_c2(B,B,B); +extern B mul_c2(B,B,B); + +static u64 elRange(u8 eltype) { return 1ull<<(1<ft? it : ft; /* Table allocation width */ \ + TALLOC(TY, tab0, t); TY* tab = tab0 + t/2; \ + usz m=IN##ia, n=FOR##ia; \ + void* ip = tyany_ptr(IN); \ + void* fp = tyany_ptr(FOR); \ + /* Initialize */ \ + if (FOR##e==el_i16 && n16 && we20 && we<=el_i16 && xe<=el_i16) { + B r; + TABLE(w, x, i32, wia, i) + return r; + } + i32* rp; B r = m_i32arrc(&rp, x); + H_b2i* map = m_b2i(64); + SGetU(x) + SGetU(w) + for (usz i = 0; i < wia; i++) { + bool had; u64 p = mk_b2i(&map, GetU(w,i), &had); + if (!had) map->a[p].val = i; + } + for (usz i = 0; i < xia; i++) rp[i] = getD_b2i(map, GetU(x,i), wia); + free_b2i(map); decG(w); decG(x); + return wia<=I8_MAX? taga(cpyI8Arr(r)) : wia<=I16_MAX? taga(cpyI16Arr(r)) : r; + } + } + return c2(rt_indexOf, w, x); +} + +B enclosed_0; +B enclosed_1; +extern B rt_memberOf; +B memberOf_c2(B t, B w, B x) { + if (isAtm(x) || RNK(x)!=1) goto bad; + if (isAtm(w)) goto single; + ur wr = RNK(w); + if (wr==0) { + B w0 = IGet(w, 0); + dec(w); + w = w0; + goto single; + } + if (wr==1) goto many; + goto bad; + + bad: return c2(rt_memberOf, w, x); + + B r; + single: { + usz xia = IA(x); + SGetU(x) + for (usz i = 0; i < xia; i++) if (equal(GetU(x, i), w)) { r = incG(enclosed_1); goto dec_wx; } + r = incG(enclosed_0); + dec_wx:; dec(w); + goto dec_x; + } + + + many: { + u8 we = TI(w,elType); usz wia = IA(w); + u8 xe = TI(x,elType); usz xia = IA(x); + if (xia == 0) { r=taga(arr_shVec(allZeroes(wia))); decG(w); goto dec_x; } + #define WEQ(V) C2(eq, incG(w), V) + if (xe==el_bit) { + u64* xp = bitarr_ptr(x); + u64 x0 = 1 & xp[0]; + r = WEQ(m_usz(x0)); + if (bit_has(xp, xia, !x0)) r = C2(or, r, WEQ(m_usz(!x0))); + decG(w); goto dec_x; + } + if (xia<=(xe==el_i16?8:16) && wia>16 && we20 && we<=el_i16 && xe<=el_i16) { + B r; + TABLE(x, w, i8, 0, 1) + return num_squeeze(r); + } + H_Sb* set = m_Sb(64); + SGetU(x) SGetU(w) + bool had; + for (usz i = 0; i < xia; i++) mk_Sb(&set, GetU(x,i), &had); + u64* rp; r = m_bitarrv(&rp, wia); + for (usz i = 0; i < wia; i++) bitp_set(rp, i, has_Sb(set, GetU(w,i))); + free_Sb(set); decG(w); + goto dec_x; + } + + dec_x:; + decG(x); + return r; +} + +extern B rt_count; +B count_c2(B t, B w, B x) { + if (!isArr(w) || RNK(w)==0) thrM("⊒: 𝕨 must have rank at least 1"); + if (RNK(w)!=1) return c2(rt_count, w, x); + if (!isArr(x) || IA(x)<=1) return indexOf_c2(m_f64(0), w, x); + u8 we = TI(w,elType); usz wia = IA(w); + u8 xe = TI(x,elType); usz xia = IA(x); + i32* rp; B r = m_i32arrc(&rp, x); + TALLOC(usz, wnext, wia+1); + wnext[wia] = wia; + if (we<=el_i16 && xe<=el_i16) { + if (we==el_bit) { w = toI8Any(w); we = TI(w,elType); } + if (xe==el_bit) { x = toI8Any(x); xe = TI(x,elType); } + usz it = elRange(we); // Range of writes + usz ft = elRange(xe); // Range of lookups + usz t = it>ft? it : ft; // Table allocation width + TALLOC(i32, tab0, t); i32* tab = tab0 + t/2; + usz m=wia, n=xia; + void* ip = tyany_ptr(w); + void* fp = tyany_ptr(x); + // Initialize + if (xe==el_i16 && na[p].val : wia; + map->a[p].val = i; + } + for (usz i = 0; i < xia; i++) { + bool had; u64 p = getQ_b2i(map, GetU(x,i), &had); + usz j = wia; + if (had) { j = map->a[p].val; map->a[p].val = wnext[j]; } + rp[i] = j; + } + free_b2i(map); + } + TFREE(wnext); decG(w); decG(x); + return wia<=I8_MAX? taga(cpyI8Arr(r)) : wia<=I16_MAX? taga(cpyI16Arr(r)) : r; +} + + +void search_gcFn() { + mm_visit(enclosed_0); + mm_visit(enclosed_1); +} + +void search_init() { + gc_addFn(search_gcFn); + { u64* p; Arr* a=m_bitarrp(&p, 1); arr_shAlloc(a,0); *p= 0; enclosed_0=taga(a); } + { u64* p; Arr* a=m_bitarrp(&p, 1); arr_shAlloc(a,0); *p=~0ULL; enclosed_1=taga(a); } +} diff --git a/src/builtins/selfsearch.c b/src/builtins/selfsearch.c index 37e647fd..5bf2ada2 100644 --- a/src/builtins/selfsearch.c +++ b/src/builtins/selfsearch.c @@ -3,6 +3,58 @@ #include "../utils/mut.h" #include "../utils/talloc.h" +B not_c1(B t, B x); +B shape_c1(B t, B x); + +#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 (; i>24)]++; c1[(u8)(v>>16)]++; } - // Exclusive prefix sum - usz s0=0, s1=0; - for (usz j=0; j>24); usz c=c0[k]++; v1[c]=v; } - for (usz i=0; i>16); 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]=1; } - for (; i(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; \ @@ -127,7 +164,6 @@ B count_c1(B t, B x) { i32* r0; B r = m_i32arrv(&r0, n); 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 # @@ -140,30 +176,7 @@ B count_c1(B t, B x) { u32 *r1 = (u32*)v1; // n [+..] r1 4*n ######## u32 *tab= (u32*)v1; // tn [+] tab 4*tn ########### - // Count keys - for (usz j=0; j<2*rx; j++) c0[j] = 0; - for (usz i=0; i>24)]++; c1[(u8)(v>>16)]++; } - // Exclusive prefix sum - usz s0=0, s1=0; - for (usz j=0; j>24); usz c=c0[k]++; v1[c]=v; } - for (usz i=0; i>16); 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]=0; } - for (; ia[p].val : def; } + static inline u64 N(getQ) (Map* m, KT k, bool* had) { + u64 h1 = H1(k); u64 h2 = H2(k, h1); + return N(find)(m, k, h1, h2, had); + } #endif