diff --git a/build/src/build.bqn b/build/src/build.bqn index 7ad34a16..6d915d0d 100755 --- a/build/src/build.bqn +++ b/build/src/build.bqn @@ -696,7 +696,7 @@ cachedBin‿linkerCache ← { "xa."‿"src/builtins/squeeze.c"‿"squeeze", "xa."‿"src/utils/mut.c"‿"copy", "xa."‿"src/utils/bits.c"‿"bits", "xag"‿"src/builtins/transpose.c"‿"transpose", "xag"‿"src/builtins/search.c"‿"search", "xag"‿"src/builtins/selfsearch.c"‿"selfsearch" - "xag"‿"src/builtins/scan.c"‿"scan", "xa."‿"src/builtins/fold.c"‿"fold", + "xag"‿"src/builtins/scan.c"‿"scan", "xag"‿"src/builtins/fold.c"‿"fold", "xag"‿"src/builtins/slash.c"‿"slash", "xag"‿"src/builtins/slash.c"‿"replicate", "xag"‿"src/builtins/sort.c"‿"bins", "xa."‿"src/builtins/slash.c"‿"count" diff --git a/src/builtins/cells.c b/src/builtins/cells.c index c27bc9b3..68d36f30 100644 --- a/src/builtins/cells.c +++ b/src/builtins/cells.c @@ -7,8 +7,10 @@ B fne_c1(B, B); B shape_c2(B, B, B); B transp_c2(B, B, B); -B fold_rows(Md1D* d, B x); // from fold.c -B takedrop_highrank(bool take, B w, B x); // from sfns.c +B fold_rows(Md1D* d, B x); // from fold.c +B fold_rows_bit(Md1D* d, B x, usz n, usz m); // from fold.c +B scan_rows_bit(u8, B x, usz m); // from scan.c +B takedrop_highrank(bool take, B w, B x); // from sfns.c B try_interleave_cells(B w, B x, ur xr, ur xk, usz* xsh); // from transpose.c // X - variable name; XSH - its shape; K - number of leading axes that get iterated over; SLN - number of slices that will be made; DX - additional refcount count to add to x @@ -175,6 +177,7 @@ NOINLINE B leading_axis_arith(FC2 fc2, B w, B x, usz* wsh, usz* xsh, ur mr) { // // fast special-case implementations +extern void (*const si_select_cells_bit_lt64)(uint64_t*,uint64_t*,uint32_t,uint32_t,uint32_t); // from fold.c (fold.singeli) static NOINLINE B select_cells(usz n, B x, usz cam, usz k, bool leaf) { // n {leaf? <∘⊑; ⊏}⎉¯k x; TODO probably can share some parts with takedrop_highrank and/or call ⊏? ur xr = RNK(x); assert(xr>1 && km1->flags-1; if (rtid==n_const) { f=fd->f; goto const_f; } - if ((rtid==n_fold || rtid==n_insert) && TI(x,elType)!=el_B && k==1 && xr==2 && isFun(fd->f) && isPervasiveDyExt(fd->f)) { // TODO extend to any rank x with cr==1 - usz *sh = SH(x); usz m = sh[1]; - if (m == 1) return select_cells(0, x, cam, k, false); - if (m <= 64 && m < sh[0]) return fold_rows(fd, x); + usz *sh = SH(x); + if (((rtid==n_fold && cr==1) || rtid==n_insert) && TI(x,elType)!=el_B + && isFun(fd->f) && 1==shProd(sh, k+1, xr) && sh[k] > 0) { + usz m = sh[k]; + u8 frtid = v(fd->f)->flags-1; + if (m==1 || frtid==n_ltack) return select_cells(0 , x, cam, k, false); + if ( frtid==n_rtack) return select_cells(m-1, x, cam, k, false); + if (isPervasiveDyExt(fd->f)) { + if (TI(x,elType)==el_bit) { + incG(x); // keep shape alive + B r = fold_rows_bit(fd, x, shProd(sh, 0, k), m); + if (!q_N(r)) { + if (xr > 2) { + usz* rsh = arr_shAlloc(a(r), xr-1); + shcpy(rsh, sh, k); + shcpy(rsh+k, sh+k+1, xr-1-k); + } + decG(x); return r; + } + decG(x); + } + // TODO extend to any rank + if (xr==2 && k==1 && m<=64 && mf)) goto base; + u8 frtid = v(fd->f)->flags-1; + if (frtid==n_rtack) return x; + if (TI(x,elType)==el_bit && (isPervasiveDyExt(fd->f)||frtid==n_ltack) + && 1==shProd(sh, k+1, xr)) { + B r = scan_rows_bit(frtid, x, m); if (!q_N(r)) return r; + } } } else if (TY(f) == t_md2D) { Md2D* fd = c(Md2D,f); diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 1a94a6ae..ccf88581 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -11,18 +11,52 @@ // •math.Sum: +´ with faster and more precise SIMD code for i32, f64 +// Insert with rank (˝˘ or ˝⎉k), or fold on flat array +// SHOULD optimize dyadic insert with rank +// Length 1, ⊣⊢: implemented as ⊏˘ +// SHOULD reshape identity for fast ˝˘ on empty rows +// Boolean operand, cell size 1: +// +: popcount +// Rows length ≤64: extract rows, popcount each +// TRIED scan-based version, faster for width 3 only +// COULD have bit-twiddling versions of +˝˘ on very short rows +// >64: popcount, boundary corrections (clang auto-vectorizes) +// ∨∧≠= and synonyms, rows <64: extract from scan or windowed op +// Dedicated auto-vectorizing code for sizes 2, 4 +// Extract is generic with multiplies or BMI2 (SHOULD add SIMD) +// COULD use CLMUL for faster windowed ≠ +// ∨∧ and synonyms: +// Multiples of 8 <256: comparison function, recurse if needed +// Rows 6424: read word, mask, last bit of popcount +// Rows l<64: word-based scan, correct with SWAR +// ≥64: xor words (auto-vectorizes), correct boundary, popcount +// COULD implement boolean -˝˘ with xor, +˝˘, offset +// Arithmetic on rank 2, short rows: transpose then insert, blocked +// SHOULD extend transpose-insert code to any frame and cell rank +// SHOULD have dedicated +⌊⌈ insert with rank + #include "../core.h" #include "../builtins.h" #include "../utils/mut.h" +#include "../utils/calls.h" -#if SINGELI_SIMD +#if SINGELI + extern uint64_t* const si_spaced_masks; + #define get_spaced_mask(i) si_spaced_masks[i-1] #define SINGELI_FILE fold #include "../utils/includeSingeli.h" #endif -static bool fold_ne(u64* x, u64 am) { +static u64 xor_words(u64* x, u64 l) { u64 r = 0; - for (u64 i = 0; i < (am>>6); i++) r^= x[i]; + for (u64 i = 0; i < l; i++) r^= x[i]; + return r; +} +static bool fold_ne(u64* x, u64 am) { + u64 r = xor_words(x, am>>6); if (am&63) r^= x[am>>6]<<(64-am & 63); return POPC(r) & 1; } @@ -92,8 +126,8 @@ B sum_c1(B t, B x) { } r += s; } else { - #if SINGELI_SIMD - r = simd_sum_f64(xv, ia); + #if SINGELI + r = si_sum_f64(xv, ia); #else r=0; for (usz i=0; i) } DEF_MIN_MAX(i8) DEF_MIN_MAX(i16) DEF_MIN_MAX(i32) -#if SINGELI_SIMD - static f64 min_f64(void* xv, usz ia) { return simd_fold_min_f64(xv,ia); } - static f64 max_f64(void* xv, usz ia) { return simd_fold_max_f64(xv,ia); } +#if SINGELI + static f64 min_f64(void* xv, usz ia) { return si_fold_min_f64(xv,ia); } + static f64 max_f64(void* xv, usz ia) { return si_fold_max_f64(xv,ia); } #else DEF_MIN_MAX(f64) #endif @@ -432,3 +466,102 @@ B fold_rows(Md1D* fd, B x) { return mut_fv(r); } } + +B sum_rows_bit(B x, usz n, usz m) { + u64* xp = bitarr_ptr(x); + if (m < 128) { + if (m == 2) return bi_N; // Transpose is faster + i8* rp; B r = m_i8arrv(&rp, n); + if (m <= 64) { + if (m%8 == 0) { + usz k = m/8; u64 b = (m==64? 0 : 1ull<> (j%8); + rp[i] = POPC(b & xw); + } + } else { + // Row may not fit in an aligned word + // Read a word containing the last bit, combine with saved bits + u64 b = ~(~(u64)0 >> m); + u64 prev = 0; + for (usz i=0, j=m; i> (m-sh); + } + } + } + } else { // 64> (in%64)); + rp[i] = s - o; + j = jn+1; + } + } + decG(x); return r; + } else if (m < 1<<15) { + i16* rp; B r = m_i16arrv(&rp, n); + usz l = m/64; + if (m%64==0) { + for (usz i=0; i= j+l)); + o = POPC(e >> (in%64)); + rp[i] = s - o; + j = jn+1; + } + } + decG(x); return r; + } else { + return bi_N; + } +} + +// Fold n cells of size m, stride 1 +// Return a vector regardless of argument shape, or bi_N if not handled +B fold_rows_bit(Md1D* fd, B x, usz n, usz m) { + assert(isArr(x) && TI(x,elType)==el_bit && IA(x)==n*m); + if (!v(fd->f)->flags) return bi_N; + u8 rtid = v(fd->f)->flags-1; + if (rtid==n_add) return sum_rows_bit(x, n, m); + #if SINGELI + bool is_or = rtid==n_or |rtid==n_ceil; + bool andor = rtid==n_and|rtid==n_floor|rtid==n_mul|is_or; + if (rtid==n_ne|rtid==n_eq|andor) { + if (n==0) { decG(x); return taga(arr_shVec(allZeroes(0))); } + if (andor && m < 256) while (m%8 == 0) { + usz f = CTZ(m|32); + m >>= f; usz c = m*n; + u64* yp; B y = m_bitarrv(&yp, c); + u8 e = el_i8 + f-3; + CmpASFn cmp = is_or ? CMP_AS_FN(ne, e) : CMP_AS_FN(eq, e); + CMP_AS_CALL(cmp, yp, bitarr_ptr(x), m_f64(is_or-1), c); + decG(x); if (m==1) return y; + x = y; + } + u64* xp = bitarr_ptr(x); + u64* rp; B r = m_bitarrv(&rp, n); + if (andor) si_or_rows_bit(xp, rp, n, m, !is_or); + else si_xor_rows_bit(xp, rp, n, m, rtid==n_eq); + decG(x); return r; + } + #endif + return bi_N; +} diff --git a/src/builtins/scan.c b/src/builtins/scan.c index f798881d..5e582904 100644 --- a/src/builtins/scan.c +++ b/src/builtins/scan.c @@ -1,3 +1,39 @@ +// Scan (`) +// Empty 𝕩, and length 1 if no 𝕨: return 𝕩 +// Generic operand: +// Constant: copy +// ⊢ identity, ⊣ reshape 𝕨 or first cell +// Boolean operand, rank 1: +// + AVX2 expansion (SHOULD have better generic, add SSE, NEON) +// ∨⌈ ∧×⌊ search+copy, then memset (COULD vectorize search) +// ≠ SWAR shifts, CLMUL, VPCLMUL (SHOULD add SSE, NEON) +// < SWAR +// =≤≥>- in terms of ≠<∨∧+ with adjustments +// Arithmetic operand, rank 1: +// ⌈⌊ Scalar, SSE, AVX in log(vector width) steps +// + Overflow-checked scalar or AVX2 +// Ad-hoc boolean-valued handling for ≠∨ +// SHOULD extend rank 1 special cases to cell bound 1 +// Higher-rank arithmetic, non-tiny cells: apply operand cell-wise +// SHOULD have dedicated high-rank scan optimizations + +// Scan with rank (`˘ or `⎉k) +// SHOULD optimize dyadic scan with rank +// Empty 𝕩, length 1, ⊢: return 𝕩 +// Boolean operand, cell size 1: +// ≠∨∧⊣ and synonyms, rows <64: SWAR, AVX2 (SHOULD add SSE, NEON) +// Power of two row size: autovectorized +// COULD have dedicated SIMD for CPU widths, little improvement +// ⊣ SWAR for <64, select for ≥ +// ∨⌈ ∧×⌊ SWAR with addition for small rows, search for large +// Rows 64≤l<160: SWAR specialized for ≤1 boundary +// Large rows: word-at-a-time search +// ≠ power-of-two shifts for <64, rank-1 scans and boundary corrections if ≥ +// SHOULD have a better intermediate-size (< ~256) SIMD method +// + scan in blocks, correct with mask, ⌊`, subtract +// = as ≠`⌾¬, - as (2×⊣`)-+` +// SHOULD optimize non-boolean scan with rank + #include "../core.h" #include "../utils/mut.h" #include "../builtins.h" @@ -5,11 +41,23 @@ #define F64_MIN -INFINITY #define F64_MAX INFINITY +#define CASE_N_AND case n_and: case n_mul: case n_floor +#define CASE_N_OR case n_or: case n_ceil + #if !USE_VALGRIND static u64 vg_rand(u64 x) { return x; } #endif +B slash_c1(B, B); +B shape_c2(B, B, B); +B fne_c1(B, B); +B add_c2(B, B, B); +B sub_c2(B, B, B); +B mul_c2(B, B, B); + #if SINGELI + extern uint64_t* const si_spaced_masks; + #define get_spaced_mask(i) si_spaced_masks[i-1] #define SINGELI_FILE scan #include "../utils/includeSingeli.h" #endif @@ -35,6 +83,12 @@ B scan_ne(B x, u64 p, u64 ia) { // consumes x #endif decG(x); return r; } +B scan_eq(B x, u64 ia) { // consumes x + B r = scan_ne(x, 0, ia); + u64* rp = bitarr_ptr(r); + for (usz i = 0; i < BIT_N(ia); i++) rp[i] ^= 0xAAAAAAAAAAAAAAAA; + return r; +} static B scan_or(B x, u64 ia) { // consumes x u64* xp = bitarr_ptr(x); @@ -51,7 +105,6 @@ static B scan_and(B x, u64 ia) { // consumes x decG(x); return FL_SET(r, fl_dsc|fl_squoze); } -B slash_c1(B t, B x); B scan_add_bool(B x, u64 ia) { // consumes x u64* xp = bitarr_ptr(x); u64 xs = bit_sum(xp, ia); @@ -113,7 +166,6 @@ B scan_max_num(B x, u8 xe, u64 ia) { MINMAX(max,>,MIN,or ,asc) } #undef MINMAX // Initialized: try to convert 𝕨 to type of 𝕩 // (could do better for out-of-range floats) -B shape_c2(B, B, B); #define MM2_ICASE(T,N,C,I) \ case el_##T : { \ if (wv!=(T)wv) { if (wv C 0) { r=C2(shape,m_f64(ia),w); break; } else wv=I; } \ @@ -174,9 +226,6 @@ static B scan_plus(f64 r0, B x, u8 xe, usz ia) { #endif } -B fne_c1(B, B); -B shape_c2(B, B, B); - extern B scan_arith(B f, B w, B x, usz* xsh); // from cells.c B scan_c1(Md1D* d, B x) { B f = d->f; if (isAtm(x) || RNK(x)==0) thrM("`: Argument cannot have rank 0"); @@ -204,13 +253,17 @@ B scan_c1(Md1D* d, B x) { B f = d->f; } if (!(xr==1 && xe<=el_f64)) goto base; - if (xe==el_bit) { - if (rtid==n_add ) return scan_add_bool(x, ia); // + - if (rtid==n_or | rtid==n_ceil ) return scan_or(x, ia); // ∨⌈ - if (rtid==n_and | rtid==n_mul | rtid==n_floor) return scan_and(x, ia); // ∧×⌊ - if (rtid==n_ne ) return scan_ne(x, 0, ia); // ≠ - if (rtid==n_lt) return scan_lt(x, 0, ia); // < - goto base; + if (xe==el_bit) switch (rtid) { default: goto base; + case n_add: return scan_add_bool(x, ia); // + + CASE_N_OR: return scan_or(x, ia); // ∨⌈ + CASE_N_AND: return scan_and(x, ia); // ∧×⌊ + case n_ne: return scan_ne(x, 0, ia); // ≠ + case n_eq: return scan_eq(x, ia); // = + case n_lt: return scan_lt(x, 0, ia); // < + case n_le: return bit_negate(scan_lt(bit_negate(x), 0, ia)); // ≤ + case n_gt: x=bit_negate(x); *bitarr_ptr(x)^= 1; return scan_and(x, ia); // > + case n_ge: x=bit_negate(x); *bitarr_ptr(x)^= 1; return scan_or (x, ia); // ≥ + case n_sub: return C2(sub, m_f64(2 * (*bitarr_ptr(x) & 1)), scan_add_bool(x, ia)); // - } if (rtid==n_add) return scan_plus(0, x, xe, ia); // + if (rtid==n_floor) return scan_min_num(x, xe, ia); // ⌊ @@ -251,7 +304,6 @@ B scan_c1(Md1D* d, B x) { B f = d->f; return withFill(r.b, xf); } -B add_c2(B, B, B); B scan_c2(Md1D* d, B w, B x) { B f = d->f; if (isAtm(x) || RNK(x)==0) thrM("`: 𝕩 cannot have rank 0"); ur xr = RNK(x); usz* xsh = SH(x); usz ia = IA(x); @@ -329,3 +381,52 @@ B scan_c2(Md1D* d, B w, B x) { B f = d->f; decG(x); return withFill(r.b, wf); } + +// scan cells of size m, stride 1 +B scan_rows_bit(u8 rtid, B x, usz m) { + assert(isArr(x) && TI(x,elType)==el_bit); + #if SINGELI + switch (rtid) { default: return bi_N; + case n_eq: return bit_negate(scan_rows_bit(n_ne, bit_negate(x), m)); + CASE_N_AND: CASE_N_OR: case n_ne: case n_ltack: { + usz ia = IA(x); + u64* xp = bitarr_ptr(x); + u64* rp; B r = m_bitarrc(&rp, x); + switch (rtid) { default:UD; + CASE_N_AND: si_scan_rows_and (xp, rp, ia, m); break; + CASE_N_OR: si_scan_rows_or (xp, rp, ia, m); break; + case n_ne: si_scan_rows_ne (xp, rp, ia, m); break; + case n_ltack: si_scan_rows_ltack(xp, rp, ia, m); break; + } + decG(x); return r; + } + case n_add: case n_sub: { + usz ia = IA(x); + if (m >= 128) return bi_N; + usz bl = 128; // block size + i8 buf[bl]; i8 c = 0; + u64* xp = bitarr_ptr(x); + i8* rp; B r = m_i8arrc(&rp, x); + static const u64 ms[7] = { 0x00ff00ff00ff00ff, 0x00ff0000ff0000ff, 0x000000ff000000ff, 0x0000ff00000000ff, 0x00ff0000000000ff, 0xff000000000000ff, 0 }; + u64 mm = ms[m-2>6? 6 : m-2]; usz mk = m*(POPC(mm)/8); + for (usz i = 0, j = m; i < ia; i += bl) { + usz len = ia - i; if (len > bl) len = bl; + usz e = i + len; + si_bcs8(xp + i/64, buf, len); + memset(rp+i, -c, len); + i8* bi = buf; bi-=i; // yeah this makes the pointer go out of bounds, but whatever + assert(j > i); + if (mk) while (j+mk <= e) { storeu_u64(rp+j, loadu_u64(bi+j-1) & mm); j+=mk; } + for (; j < e; j += m) rp[j] = bi[j-1]; + si_scan_max_init_i8(rp+i, rp+i, len, I8_MIN); + for (usz k = i; k < e; k++) rp[k] = bi[k] - rp[k]; + if (j == e) { j += m; c = 0; } else c = rp[e-1]; + } + if (rtid!=n_sub) { decG(x); return r; } + return C2(sub, C2(mul, m_f64(2), scan_rows_bit(n_ltack, x, m)), r); + } + } + #else + return bi_N; + #endif +} diff --git a/src/builtins/slash.c b/src/builtins/slash.c index ce4a7deb..27357e79 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -777,7 +777,11 @@ B slash_c2(B t, B w, B x) { u64* rp; r = m_bitarrv(&rp, s); #if FAST_PDEP if (wv <= 52) { - u64 m = (u64)-1 / (((u64)1< v & -promote{T,m}) infix left 35 + def load {p:*[_]_, n } = assert{0} @@ -156,7 +158,7 @@ def lvec = match { {[n]T, n, (width{T})} => 1; {T, n, w} => 0 } # base cases def { - absu,andAllZero,andnz,b_getBatch,clmul,cvt,extract,fold_addw,half, + absu,andAllZero,andnz,b_getBatch,blend,clmul,cvt,extract,fold_addw,half, homAll,homAny,homBlend,homMask,homMaskStore,homMaskStoreF,loadBatchBit, loadLow,make,maskStore,maskToHom,mulw,mulh,narrow,narrowPair,packHi,packLo,packQ,pair,pdep, pext,popcRand,sel,shl,shr,shuf,shuf16Hi,shuf16Lo,shufHalves,storeLow, diff --git a/src/singeli/src/fold.singeli b/src/singeli/src/fold.singeli index 1f262be5..f6019c71 100644 --- a/src/singeli/src/fold.singeli +++ b/src/singeli/src/fold.singeli @@ -1,5 +1,9 @@ include './base' include './mask' +if_inline (hasarch{'BMI2'}) include './bmi2' +include './spaced' +include 'util/tup' +include './scan_common' def opsh64{op}{v:([4]f64), perm} = op{v, shuf{[4]u64, v, perm}} def opsh32{op}{v:([2]f64), perm} = op{v, shuf{[4]u32, v, perm}} @@ -20,7 +24,13 @@ def reduce_pairwise{op, plog, x:*T, len, init:T} = { r } -fn fold_idem{T==f64, op}(x:*T, len:u64) : T = { +fn fold_idem{T, op}(x:*T, len:u64) : T = { + assert{len > 0} + a := load{x, 0} + @for (x over _ from 1 to len) a = op{a, x} + a +} +fn fold_idem{T==f64, op if has_simd}(x:*T, len:u64) : T = { def bulk = arch_defvw/width{T} def V = [bulk]T xv:= *V ~~ x @@ -45,10 +55,15 @@ fn fold_idem{T==f64, op}(x:*T, len:u64) : T = { extract{r, 0} } -export{'simd_fold_min_f64', fold_idem{f64,min}} -export{'simd_fold_max_f64', fold_idem{f64,max}} +export{'si_fold_min_f64', fold_idem{f64,min}} +export{'si_fold_max_f64', fold_idem{f64,max}} -fn fold_assoc_0{T==f64, op}(x:*T, len:u64) : T = { +fn fold_assoc_0{T, op}(x:*T, len:u64) : T = { + a:T = 0 + @for (x over len) a = op{a, x} + a +} +fn fold_assoc_0{T==f64, op if has_simd}(x:*T, len:u64) : T = { def bulk = arch_defvw/width{T} def V = [bulk]T xv:= *V ~~ x @@ -58,4 +73,262 @@ fn fold_assoc_0{T==f64, op}(x:*T, len:u64) : T = { if (hasarch{'AARCH64'}) vfold{op, r} else extract{mix{op, r}, 0} } -export{'simd_sum_f64', fold_assoc_0{f64,+}} +export{'si_sum_f64', fold_assoc_0{f64,+}} + + +# Short-row boolean folds: main challenge is bit packing +def fold_rows_bit_lt64{ + op, run_loop2, run_loop4, pext_res, mult_in, + off, # mask offset for generic methods + xx, rx, # input and output xor for cases not specialized to individual functions + rxs, # output xor only, where and/or are specialized + xp, rp, n, l +} = { + nw := cdiv{n*l,64} + def loop_T{T, proc_word} = { + @for (x in xp, r in *T~~rp over nw) r = cast_i{T, rxs{proc_word{x}}} + } + if (l == 2) { + def extract = if (fast_BMI2{}) pext{., u64~~64w2b01} else { + # repeated shift-or-mask + def step{x, {m, sh}} = { def a=x&m; a|a>>sh } + def ss = 1< loop_T{u32, {x} => extract{op{x, x>>1}}}} + } else if (l == 4) { + m:u64 = 64w2b0001; t := m<<3 + def extract = if (fast_BMI2{}) pext{., t} else { + {x} => { + a:= x & t + a = (a * 4r2b001) & 64w0xf000 + a = (a * 4r0x001) >> (64-16) + } + } + run_loop4{m, t, {get} => loop_T{u16, {x} => extract{get{x}}}} + } else { # generic width<64 + {e0, d} := unaligned_spaced_mask_mod{l} + el:= e0 << (l-1) # ending bit of each row + e := if (same{off,-1}) el else e0<= 32) { + store{rh, 0, cast_i{u32,rx{r}}}; ++rh; + r >>= 32; ri -= 32 + } + } + def flush{} = if (ri > 0) store{rh, 0, cast_i{u32,rx{r}}} + tup{write, flush} + } + def finish_step{...a} = { + write_bits{...a} + e = e>>d | e<<(l-d) # update end mask for next iteration + } + if (fast_BMI2{}) { + @for (x in xp over nw) finish_step{pext_res{xx{x}, e, c}, popc{e}} + } else { + # Emulate pext with 1, 2, or 3 multiply/mask steps. + # To move size-a groups spaced at distance b together, + # the multiplier has up to b/a bits spaced by b-a. + dm:= cast_i{usz, popc{el}} # minimum output bits per word + dm-= promote{usz, l&(l-1) == 0} # for divisors of 64, e0 effectively overflows; subtract 1 to correct + def loop{...par} = { + et:= e0 << clz{e0} + @for (xo in xp over nw) { + x:= xx{xo} + def {s, mod_rb} = if (length{select{par,0}} == 3) { + # fast path for l==3 + s:= op{op{x, x<<1}, op{x<<2, c}} + c = op{x>>63, x>>62} + tup{s, {rb}=>rb} + } else { + mult_in{x, e, c} + } + nb:= dm + promote{usz, e>=e0} # = popc{e} + # the multiply-mask sequence (last "mask" is the shift by (64 - nb)) + def extract = match { + { {q}, {}} => (et &~ (et - (s & e))) * q + {{...qs, q}, {...bs, b}} => (extract{qs, bs} & b) * q + } + rb:= extract{...par} >> (64 - nb) + finish_step{mod_rb{rb}, nb} + } + } + if (l == 3) { + mult0:u64 = base{1<< 2, 3**1}; top3:u64 = base{1<<9, 8**(1<<3-1)}>>2 + mult1:u64 = base{1<< 6, 3**1}; top9:u64 = base{1<<27, 3**(1<<9-1)}<<1 + mult2:u64 = base{1<<18, 3**1} + loop{tup{mult0,mult1,mult2}, tup{top3, top9}} + } else if (l < 8) { + assert{l > 4} + ld:= l-1; lld:= l*ld + {mult0, _} := unaligned_spaced_mask_mod{ld} + mult0 &= u64~~1<>l; topk|= topk<>ll + loop{tup{mult0,mult1}, tup{topk}} + } else { + {mult, _} := unaligned_spaced_mask_mod{l-1} + if (l==8) mult &= 1<<(7*8) - 1 + loop{tup{mult}, tup{}} + } + } + flush_bits{} + } +} + +fn select_rows_bit_lt64(xp:*u64, rp:*u64, n:usz, l:usz, o:usz) : void = { + assert{l < 64}; assert{o < l} # Row length, and offset within row + def run_loop2{loop} = loop{{a,b} => a>>o} + def run_loop4{m, t, loop} = loop{{x} => x<<(l-1-o)} + def pext_res{x, e, c} = pext{x, e} + def mult_in{x, e, c} = tup{x, {r}=>r} + def id{x} = x + fold_rows_bit_lt64{ + {a,b}=>a, run_loop2, run_loop4, pext_res, mult_in, o, id, id, id, + xp, rp, n, l + } +} + +fn xor_words(init:u64, x:*u64, l:usz):u64 = { + @for (x over l) init ^= x + init +} +def bit_output{rp:*T} = { + buf:u64 = 0 # Buffer for result bits + def output{i, bit, mod} = { + buf = buf>>1 | promote{u64, bit}<<63 + if ((i+1)%64==0) { store{rp, 0, mod{buf}}; ++rp } + } + def fixbuf{mod} = { buf = mod{buf} } + def flush_bits{n, mod} = { + q:=(-n)%64; if (q!=0) store{rp, 0, mod{buf} >> q} + } + def flush_bits{n} = flush_bits{n, {b}=>b} + tup{output, fixbuf, flush_bits} +} +# word and alignment of start of next row +def next_start{i, m} = { + bn := promote{u64, i+1} * promote{u64, m} + tup{bn/64, bn%64} +} +fn xor_rows_bit(xp:*u64, rp:*u64, n:usz, l:usz, eq:u1) : void = { + def p64 = promote{u64, .} + rx:= -(p64{eq} &~ p64{l}) # ne to eq conversion + if (l <= 64 and not (l%8==0 and l>16)) { + def run_loop2{loop} = loop{^} + def run_loop4{m, t, loop} = loop{{x} => { x2:= x^(x<<1); x2^(x2<<2) }} + def xor_word = prefix_byshift{^, <<} + def pext_in{xo, e, c} = { + x := xor_word{xo} + rb:= x ^ (x<>63))>>(64-l) + rb + } + def pext_res{x, e, c} = pext{pext_in{x, e, c}, e} + def mult_in{x, e, c} = tup{pext_in{x, e, c}, {r}=>r} + fold_rows_bit_lt64{ + ^, run_loop2, run_loop4, pext_res, mult_in, -1, {x}=>x, ^{rx,.}, ^{rx,.}, + xp, rp, n, l + } + } else { + def fixout = ^{rx, .} + def {add_bit, _, flush_bits} = bit_output{rp} + def add_bit{i, bit} = add_bit{i, bit, fixout} + def xor_loop{len} = { + o:u64 = 0 # Carry + j:u64 = 0; @for (i to n) { + def {jn, sh} = next_start{i, l} + s := xor_words(o, xp + j, len) + e := load{xp, jn} + s ^= e & (if (not same{len,1}) -p64{jn >= j + p64{len}} else j - jn) + o = e >> sh + add_bit{i, popc{s ^ o}} + j = jn+1 + } + } + ll := l/64 + if (l <= 64) { + bm:= u64~~2<<(l-1) - 1 + k:= l/8 + @for (i to n) add_bit{i, popc{bm & loadu{*u64~~(*u8~~xp + k*i)}}} + } else if (l < 128) xor_loop{1} + else if (l%64==0) { + @for (i to n) add_bit{i, popc{xor_words(0, xp+ll*i, ll)}} + } + else xor_loop{ll} + flush_bits{n, fixout} + } +} + +fn or_rows_bit(xp:*u64, rp:*u64, n:usz, l:usz, op_and:u1) : void = { + def {add_bit, set_out, flush_bits} = bit_output{rp} + if (l < 64) { + def run_loop2{loop} = if (op_and) loop{&} else loop{|} + def run_loop4{m, t, loop} = { + if (op_and) loop{{x} => x & ((x&~t) + m)} + else loop{{x} => x | ((x| t) - m)} + } + def xor_word = prefix_byshift{^, <<} + def pext_in{x, e, c} = { + m := e<<1 | 1 + t := e | 1<<63 + tup{x | ((x|t) - m), t} + } + def pext_res{x, e, c} = pext{...pext_in{x, e, c}} + def mult_in{x, e, c} = { + {s, t} := pext_in{x, e, c} + cs:= c; c = (s&~e) >> 63 + tup{s, {rb}=>rb|cs} + } + def xx = ^{-promote{u64, op_and}, .} + fold_rows_bit_lt64{ + |, run_loop2, run_loop4, pext_res, mult_in, -1, xx, xx, {r}=>r, + xp, rp, n, l + } + return{} + } else if (l < 128) { + c:u64 = (promote{u64, l}-1) &- op_and # a row gives 1 if its sum is >c + o:u64 = 0 + j:u64 = 0; @for (i to n) { + def {jn, sh} = next_start{i, l} + s := o + popc{load{xp,j}} + e := load{xp,jn} + s += popc{e & (j - jn)} # mask is 0 if j==jn, or -1 + o = popc{e >> sh} + add_bit{i, s > c+o, {rw}=>rw} + j = jn+1 + } + } else { + rx := -promote{u64, op_and} + def fixout = ^{rx, .} + o:u64 = 0 # Saved bits + j:u64 = 0; @for (i to n) { + def {jn, sh} = next_start{i, l} + e := load{xp,jn} ^ rx + l := ~(u64~~0) << sh + rb:u64 = 1 + if ((o | (e &~ l)) == 0) { # Search for shortcut + @for (i from j to jn) if (load{xp,i} != rx) goto{'found'} + rb = 0; setlabel{'found'} + } + o = e & l + add_bit{i, rb, fixout} + j = jn+1 + } + set_out{fixout} + } + flush_bits{n} +} +export{'si_xor_rows_bit', xor_rows_bit} +export{'si_or_rows_bit', or_rows_bit} +export{'si_select_cells_bit_lt64', select_rows_bit_lt64} diff --git a/src/singeli/src/scan.singeli b/src/singeli/src/scan.singeli index aec9976a..3d15a81b 100644 --- a/src/singeli/src/scan.singeli +++ b/src/singeli/src/scan.singeli @@ -5,6 +5,8 @@ if_inline (hasarch{'X86_64'}) { } include './mask' include './f64' +include './spaced' +include 'util/tup' include './scan_common' # Initialized scan, generic implementation @@ -295,3 +297,205 @@ export{'si_scan_plus_i32_i32', plus_scanC{i32, i32}} export{'si_scan_plus_i16_f64', plus_scanG{i16, f64}} export{'si_scan_plus_i32_f64', plus_scanG{i32, f64}} + + + +# Row-wise boolean scan +def loop_with_unaligned_mask{x, r, nw, l, step} = { + {m, d} := unaligned_spaced_mask_mod{l} + c:u64 = 0 # carry (initial value never matters) + @for (x, r over nw) { + match (step{x, c, m}) { + {{nr,nc}} => { r=nr; c=nc } + {nr} => { r=nr; c = -(r>>63) } + } + m = m>>d | m<<(l-d) + } +} +def avx2_loop_with_unaligned_mask{xp, rp, nw, l, scan_words, apply_carry} = { + {ms, d} := unaligned_spaced_mask_mod{l} + def V = [4]u64 + d4:usz = width{V} % l + m:= make{V, scan{{a,_} => a>>d | a<<(l-d), tup{ms, ...iota{3}}}} + c:= V**0 + @maskedLoop{4} (x in tup{V, xp}, + r in tup{V, rp} over promote{u64,nw}) { + s := scan_words{x, m} + pc:= c; c = shuf{V, -(s>>63), 4b2103} + r = apply_carry{s, blend{V, c, pc, 2b0001}, (m-V**1)&~m} + m = m>>d4 | m<<(l-d4) + } +} + +fn scan_rows_andor{id}(src:*u64, dst:*u64, nl:usz, l:usz) : void = { + def qand = not id + assert{l > 0} + nw := cdiv{nl, 64} + def res_m1{x,c,m} = { # result word with carry c, popc{m}<=1 + if (qand) x &~ ((x+c) & (x+m)) + else x | ((-x-c) &~ (x-m)) + } + if (l < 64) { + if ((l & (l-1)) == 0) { + if (l == 2) { + @for (r in dst, x in src over nw) { + r = (if (qand) x & (x<<1 | 64w2b01) else x | (x<<1 & 64w2b10)) + } + } else { + m:u64 = aligned_spaced_mask{l} + t := m << (l-1) + @for (r in dst, x in src over nw) { + r = (if (qand) x &~ ((t&x) ^ ((x&~t) + m)) + else x | ~((t&~x) ^ ((x|t) - m))) + } + } + # could use for l>=8; not much faster and takes up space + # def rowwise{T} = @for (r in *T~~dst, x in *T~~src over (64/width{T})*nw) r = x &~ (x+1) + } else if (hasarch{'AVX2'}) { + def scan_words{x, m:V} = { + mb:= m | V**1 + p:= if (qand) (x &~ m) >> 1 else ~(x | m) >> 1 + a:= if (qand) p + (mb & x) else p + (mb &~ x) + if (qand) p ^ a else ~(p ^ a) + } + def apply_carry{s, c, f} = { + if (qand) s & (~f | c) else s | (f & c) + } + avx2_loop_with_unaligned_mask{src, dst, nw, l, scan_words, apply_carry} + } else { + loop_with_unaligned_mask{src, dst, nw, l, {x, c, m} => { + s:= (if (qand) (x &~ m) >> 1 else ~(x | m) >> 1 ) + a:= (if (qand) s + (x&(m|c)) else s + ((m|c) &~ x)) + tup{(if (qand) s ^ a else ~(s ^ a) ), + a >> 63} # new c + }} + } + } else if (l < 160) { + q:usz = 0 # distance to next row boundary + c:u64 = id # carry + @for (r in dst, x in src over nw) { + b:= q<64 # whether there's a boundary + p:= q%64 # its position + q-= 64 - (l &- b) + r = res_m1{x, c, promote{u64, b} << p} + c = r >> 63 + } + } else { + i :usz = 0 # row bit index + wn:usz = 0 # starting word of next row + c:u64 = id # carry + def word{bit} = bit * ((1<<64) - 1) + we:= nl/64; while (wn < we) { + iw:= wn + r := res_m1{load{src, iw}, c, u64~~1 << (i%64)} + store{dst, iw, r}; ++iw + i+= l; wn = i/64 + c = r>>63 + if (c != id) while (iw < wn) { + x:= load{src, iw} + if (x != word{not id}) { + c = id + store{dst, iw, if (qand) x &~ (x+1) else x | -x}; ++iw + goto{'shortcut'} + } + store{dst, iw, x}; ++iw + } + setlabel{'shortcut'} + @for (r in dst over _ from iw to wn) r = word{id} + } + if (i%64 != 0) { + x:= load{src, wn} + store{dst, wn, if (qand) x &~ (x+c) else x | (-x-c)} + } + } +} + +fn scan_rows_neq(x:*u64, r:*u64, nl:usz, l:usz) : void = { + def scan_word = prefix_byshift{^, <<} + assert{l > 0} + nw := cdiv{nl, 64} + if (l < 64) { + if ((l & (l-1)) == 0) { + m:u64 = aligned_spaced_mask{l} + @for (r, x over nw) { + s:= scan_word{x} + b:= s<<1 & m # last bit of previous row + r = s ^ (b< pre{v, 1} + } + s:= vec_prefix_byshift{^, <<}{x} + b:= s<<1 & m # last bit of previous row + s ^ (b< { + s:= scan_word{x} + f:= (m-1)&~m # bits before first full row + b:= s<<1 & m # last bit of previous row + s ^ ((c & f) | (b<0)) >> 63) + } + } +} + +fn scan_rows_left(x:*u64, r:*u64, nl:usz, l:usz) : void = { + def scan_word = prefix_byshift{^, <<} + assert{l > 0} + nw := cdiv{nl, 64} + if (l < 64) { + if ((l & (l-1)) == 0) { + m:u64 = aligned_spaced_mask{l} + @for (r, x over nw) { b:= x & m; r = b< { + f:= (m-1)&~m # bits before first full row + b:= x & m + (c & f) | (b<>63) + i+= l; wn = i/64 + @for (r in r over _ from iw+1 to wn) r = c + } + if (i%64 != 0) store{r, wn, c} + } +} + +export{'si_scan_rows_and', scan_rows_andor{0}} +export{'si_scan_rows_or', scan_rows_andor{1}} +export{'si_scan_rows_ne', scan_rows_neq} +export{'si_scan_rows_ltack', scan_rows_left} diff --git a/src/singeli/src/search.singeli b/src/singeli/src/search.singeli index 8bc5c67c..2716a63b 100644 --- a/src/singeli/src/search.singeli +++ b/src/singeli/src/search.singeli @@ -335,8 +335,6 @@ exportT{'simd_getRangeRaw', each{getRange, tup{i8,i16,i32,f64}}} # Hash tables -oper &- ({v:T,m} => v & -promote{T,m}) infix left 35 - def rty{name} = if (to_prim{name}=='∊') i8 else i32 def ity{name} = (to_prim{name}=='⊒')**(*u32) fn hashtab{T, name}(rpi:*rty{name}, iv:*void, mi:usz, fv:*void, ni:usz, links:ity{name}) = { diff --git a/src/singeli/src/slash.singeli b/src/singeli/src/slash.singeli index 5a64a1ae..a9d0523b 100644 --- a/src/singeli/src/slash.singeli +++ b/src/singeli/src/slash.singeli @@ -400,3 +400,11 @@ fn compress_bool(w:*u64, x:*u64, r:*u64, n:u64) : void = { export{'si_compress_bool', compress_bool} export{'si_thresh_compress_bool', u64~~thresh_bool{}} + +# Spaced boolean masks, such as 2b000001000001...0000010000 for 6 +# Mask i is the smallest possible mask containing 1 every i bits: +# there would be a 1 just past the top bit. +# The number of trailing zeros is 64%i , and the popcount is 64/i . +def get_spaced_masks{i} = (1<<64 - 1<<(64%i)) / (1<>d | m<<(l-d), d} +} diff --git a/src/utils/calls.h b/src/utils/calls.h index 137b8810..a133241e 100644 --- a/src/utils/calls.h +++ b/src/utils/calls.h @@ -16,7 +16,7 @@ CMP_DEF(ge, AS); CMP_DEF(ge, AA); CMP_DEF(lt, AS); CMP_DEF(le, AS); -// will write up to 8×⌈len÷64 bytes to WHERE, i.e. whole u64-s +// will write up to 8×⌈len÷64 bytes to WHERE, i.e. whole u64-s; LEN must not be 0 #define CMP_AA_FN(FN, ELT) cmp_fns_##FN##AA[ELT] #define CMP_AS_FN(FN, ELT) cmp_fns_##FN##AS[ELT] diff --git a/test/cases/cells.bqn b/test/cases/cells.bqn index 9ca95b06..260e9c42 100644 --- a/test/cases/cells.bqn +++ b/test/cases/cells.bqn @@ -27,3 +27,31 @@ %USE tcc ⋄ ⟨⊏⟩ {𝕎 _tcc 𝕩}⌜ ⟨4‿5⥊↕20, 4‿5‿6⥊↕120, 2‿3‿4‿5⥊↕120⟩ %USE tcc ⋄ ⊏_tcc 4‿2⥊↕8 ⋄ 1⊏_tcc 4‿2⥊↕8 ⋄ ⊏_tcc ↕4‿2 ⋄ 1⊏_tcc ↕4‿2 %USE tcc ⋄ ⊑_tcc 4‿2⥊↕8 ⋄ 1⊑_tcc 4‿2⥊↕8 ⋄ ⊑_tcc ↕4‿2 ⋄ 1⊑_tcc ↕4‿2 + +( + %USE IS_HEAPVERIFY + # big ˝˘ & `˘ tester + _basicArgs ← { Test _𝕣 ranks: + {𝕊f: {{f Test 0‿1⥊˜ 𝕩⊏0‿2‿4‿8‿32}¨ ↕𝕩⥊3 ⋄ 1}¨ ranks}¨ +‿∧‿∨‿=‿≠‿⊣‿⊢ + } + + # TODO replace ○⊢ with some pureness-preserving •internal.Keep + # w F _thing x tests F _thing⎉w x + _testFoldCells ← {∨´𝕗=⊢‿⊣? 0=×´(-𝕨)↑≢𝕩? 1; heapverify? 0=≠⥊𝕩? 1 ⊣ 𝔽˝⎉𝕨 𝕩; ! (𝔽˝⎉𝕨 ≡ 𝔽˝○⊢⎉𝕨) 𝕩} + _testScanCells ← { heapverify? 0=≠⥊𝕩? 1 ⊣ 𝔽`⎉𝕨 𝕩; ! (𝔽`⎉𝕨 ≡ 𝔽`○⊢⎉𝕨) 𝕩} + + {1 𝕨 _testFoldCells 𝕩} _basicArgs 2‿3‿4 + {1 𝕨 _testScanCells 𝕩} _basicArgs 2‿3‿4 + {2 𝕨 _testScanCells 𝕩} _basicArgs 3‿4 + {2 𝕨 _testScanCells 𝕩} _basicArgs 3‿4 + + ⟨0, 0‿1, 1‿0, 2‿2, 1, 10, 10‿1, 1‿10, 30⟩ { + 𝕨‿𝕩 ⥊¨↩ + sh ← 𝕨∾𝕩 + cr ← ≠𝕩 + +‿∧‿∨‿=‿≠‿⊣‿⊢ { + cr 𝕎 _testFoldCells 𝕩 + cr 𝕎 _testScanCells 𝕩 + }⌜ {0=×´sh? 𝕩; ∾⟨𝕩, ¬⌾⊑¨ 𝕩, ¬⌾(¯1⊑⥊)¨ 𝕩⟩} sh⊸⥊¨ 0‿1 + }⌜ ⟨0, 1, 2, 8, 32, 8‿1, 4‿8, 4‿2, 59, 60, 63, 80, 81, 200, 640, 641⟩ +) diff --git a/test/run.bqn b/test/run.bqn index 783854a7..dc40b1fb 100644 --- a/test/run.bqn +++ b/test/run.bqn @@ -63,7 +63,7 @@ searchRefine ← { exitCode ← 0 Run ← { 𝕊 testname: lns ← •FLines fullfile ← "cases"•file.At filename ← testname∾".bqn" - defk‿defv ← ⟨⟩‿⟨⟩ + defk‿defv ← ⟨"IS_HEAPVERIFY"⟩‿⟨"heapverify←"∾•Repr o.heapverify⟩ i ← 0 Pos ← {∾⟨filename, ":", •Repr 1+𝕩⟩} Next ← {𝕊: (i+↩1) ⊢ i⊑lns}