diff --git a/makefile b/makefile index cfe75915..6c5a5aba 100644 --- a/makefile +++ b/makefile @@ -247,7 +247,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 search.o selfsearch.o scan.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 fold.o scan.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 $< @@ -271,7 +271,7 @@ endif @${MAKE} i_singeli=0 singeli=0 force_build_dir=build/obj/presingeli REPLXX=0 f= lf= postmsg="singeli sources:" i_t=presingeli i_f='-O1 -DPRE_SINGELI' FFI=0 OUTPUT=build/obj/presingeli/BQN c -build_singeli: ${addprefix src/singeli/gen/, cmp.c dyarith.c copy.c equal.c squeeze.c scan.c neq.c slash.c constrep.c bits.c} +build_singeli: ${addprefix src/singeli/gen/, cmp.c dyarith.c copy.c equal.c squeeze.c fold.c scan.c neq.c slash.c constrep.c bits.c} @echo $(postmsg) src/singeli/gen/%.c: src/singeli/src/%.singeli preSingeliBin @echo $< | cut -c 17- | sed 's/^/ /' diff --git a/src/builtins/fold.c b/src/builtins/fold.c new file mode 100644 index 00000000..8937050a --- /dev/null +++ b/src/builtins/fold.c @@ -0,0 +1,244 @@ +// Fold (´) + +// Optimized operands: +// ⊣⊢ on all types +// +-∧∨=≠ and synonyms on booleans +// ≤<>≥ on booleans, monadic only, with a search +// +⌈⌊× on numbers +// Integer +: sum blocks associatively as long as sum can't exceed +-2⋆53 +// COULD implement fast numeric -´ +// ∨ on boolean-valued integers, stopping at 1 + +#include "../core.h" +#include "../builtins.h" + +#if SINGELI + #pragma GCC diagnostic push + #pragma GCC diagnostic ignored "-Wunused-variable" + #include "../singeli/gen/fold.c" + #pragma GCC diagnostic pop +#endif + +static bool fold_ne(u64* x, u64 am) { + u64 r = 0; + for (u64 i = 0; i < (am>>6); i++) r^= x[i]; + if (am&63) r^= x[am>>6]<<(64-am & 63); + return POPC(r) & 1; +} + +static i64 bit_diff(u64* x, u64 am) { + i64 r = 0; + u64 a = 0xAAAAAAAAAAAAAAAA; + for (u64 i = 0; i < (am>>6); i++) r+= POPC(x[i]^a); + if (am&63) r+= POPC((x[am>>6]^a)<<(64-am & 63)); + return r - (i64)(am/2); +} + +// It's safe to sum a block of integers as long as the current total +// is far enough from +-1ull<<53 (and integer, in dyadic fold). +static const usz sum_small_max = 1<<16; +#define DEF_INT_SUM(T,W,M,A) \ + static i64 sum_small_##T(void* xv, usz ia) { \ + i##A s=0; for (usz i=0; i0 && -lim<=c && c<=lim) { \ + usz e=i0; i0=(i0-1)&~(b-1); \ + c+=sum_small_##T(xp+i0, e-i0); \ + } \ + r = c; \ + } \ + while (i0--) r+=xp[i0]; \ + return r; \ + } +DEF_INT_SUM(i8 ,8 ,32,32) +DEF_INT_SUM(i16,16,32,32) +DEF_INT_SUM(i32,32,52,64) +#undef DEF_SUM +static f64 sum_f64(void* xv, usz i, f64 r) { + while (i--) r += ((f64*)xv)[i]; + return r; +} +static i64 (*const sum_small_fns[])(void*, usz) = { sum_small_i8, sum_small_i16, sum_small_i32 }; +static f64 (*const sum_fns[])(void*, usz, f64) = { sum_i8, sum_i16, sum_i32, sum_f64 }; + +// Try to keep to i32 product, go to f64 on overflow or non-i32 initial +#define DEF_INT_PROD(T) \ + static NOINLINE f64 prod_##T(void* xv, usz i, f64 init) { \ + while (i--) init*=((T*)xv)[i]; return init; \ + } \ + static f64 prod_int_##T(void* xv, usz ia, i32 init) { \ + T* xp = xv; \ + while (ia--) { \ + i32 i0=init; \ + if (mulOn(init,xp[ia])) return prod_##T(xv, ia+1, i0); \ + } \ + return init; \ + } +DEF_INT_PROD(i8) +DEF_INT_PROD(i16) +DEF_INT_PROD(i32) +#undef DEF_PROD +static f64 prod_f64(void* xv, usz i, f64 r) { + while (i--) r *= ((f64*)xv)[i]; + return r; +} +static f64 (*const prod_int_fns[])(void*, usz, i32) = { prod_int_i8, prod_int_i16, prod_int_i32 }; +static f64 (*const prod_fns[])(void*, usz, f64) = { prod_i8, prod_i16, prod_i32, prod_f64 }; + +#define MIN_MAX(T,C) \ + T* xp = xv; T r = xp[0]; \ + for (usz i=1; i) } +DEF_MIN_MAX(i8) DEF_MIN_MAX(i16) DEF_MIN_MAX(i32) +#if SINGELI + static f64 min_f64(void* xv, usz ia) { return avx2_fold_min_f64(xv,ia); } + static f64 max_f64(void* xv, usz ia) { return avx2_fold_max_f64(xv,ia); } +#else + DEF_MIN_MAX(f64) +#endif +#undef DEF_MIN_MAX +#undef MIN_MAX +static f64 (*const min_fns[])(void*, usz) = { min_i8, min_i16, min_i32, min_f64 }; +static f64 (*const max_fns[])(void*, usz) = { max_i8, max_i16, max_i32, max_f64 }; + +B fold_c1(Md1D* d, B x) { B f = d->f; + if (isAtm(x) || RNK(x)!=1) thrF("´: Argument must be a list (%H ≡ ≢𝕩)", x); + usz ia = IA(x); + if (ia==0) { + decG(x); + if (isFun(f)) { + B r = TI(f,identity)(f); + if (!q_N(r)) return inc(r); + } + thrM("´: No identity found"); + } + u8 xe = TI(x,elType); + if (isFun(f) && v(f)->flags) { + u8 rtid = v(f)->flags-1; + if (rtid==n_ltack) { B r = IGet(x, 0 ); decG(x); return r; } + if (rtid==n_rtack) { B r = IGet(x, ia-1); decG(x); return r; } + if (xe>el_f64) goto base; + if (xe==el_bit) { + u64* xp = bitarr_ptr(x); + f64 r; + switch (rtid) { default: goto base; + case n_add: r = bit_sum (xp, ia); break; + case n_sub: r = bit_diff(xp, ia); break; + case n_and: case n_mul: case n_floor: r = bit_has (xp, ia, 0) ^ 1; break; + case n_or: case n_ceil: r = bit_has (xp, ia, 1) ; break; + case n_ne: r = fold_ne (xp, ia) ; break; + case n_eq: r = fold_ne (xp, ia) ^ (1&~ia); break; + case n_lt: r = bit_find(xp, ia, 1) == ia-1; break; + case n_le: r = bit_find(xp, ia, 0) != ia-1; break; + case n_gt: r = bit_find(xp, ia, 0) & 1; break; + case n_ge: r =~bit_find(xp, ia, 1) & 1; break; + } + decG(x); return m_f64(r); + } + if (rtid==n_add) { // + + void *xv = tyany_ptr(x); + bool small = xe<=el_i32 & ia<=sum_small_max; + u8 sel = xe - el_i8; + f64 r = small ? sum_small_fns[sel](xv, ia) + : sum_fns[sel](xv, ia, 0); + decG(x); return m_f64(r); + } + if (rtid==n_floor) { f64 r=min_fns[xe-el_i8](tyany_ptr(x), ia); decG(x); return m_f64(r); } // ⌊ + if (rtid==n_ceil ) { f64 r=max_fns[xe-el_i8](tyany_ptr(x), ia); decG(x); return m_f64(r); } // ⌈ + if (rtid==n_mul | rtid==n_and) { // ×/∧ + void *xv = tyany_ptr(x); + u8 sel = xe - el_i8; + f64 r = xe<=el_i32 ? prod_int_fns[sel](xv, ia, 1) + : prod_fns[sel](xv, ia, 1); + decG(x); return m_f64(r); + } + if (rtid==n_or) { // ∨ + if (xe==el_i8 ) { i8* xp = i8any_ptr (x); usz i=ia; while (i--) { i8 c=xp[i]; if (c==1) break; if (c!=0) goto base; } decG(x); return m_i32(i+1 > 0); } + if (xe==el_i16) { i16* xp = i16any_ptr(x); usz i=ia; while (i--) { i16 c=xp[i]; if (c==1) break; if (c!=0) goto base; } decG(x); return m_i32(i+1 > 0); } + if (xe==el_i32) { i32* xp = i32any_ptr(x); usz i=ia; while (i--) { i32 c=xp[i]; if (c==1) break; if (c!=0) goto base; } decG(x); return m_i32(i+1 > 0); } + } + } + base:; + SLOW2("𝕎´ 𝕩", f, x); + + SGet(x) + BBB2B fc2 = c2fn(f); + B c; + if (TI(x,elType)==el_i32) { + i32* xp = i32any_ptr(x); + c = m_i32(xp[ia-1]); + for (usz i = ia-1; i>0; i--) c = fc2(f, m_i32(xp[i-1]), c); + } else { + c = Get(x, ia-1); + for (usz i = ia-1; i>0; i--) c = fc2(f, Get(x, i-1), c); + } + decG(x); + return c; +} + +B fold_c2(Md1D* d, B w, B x) { B f = d->f; + if (isAtm(x) || RNK(x)!=1) thrF("´: 𝕩 must be a list (%H ≡ ≢𝕩)", x); + usz ia = IA(x); + u8 xe = TI(x,elType); + if (isFun(f) && v(f)->flags) { + u8 rtid = v(f)->flags-1; + if (rtid==n_ltack) { + B r = w; + if (ia) { dec(w); r=IGet(x, 0); } + decG(x); return r; + } + if (rtid==n_rtack) { decG(x); return w; } + if (!isF64(w) || xe>el_f64) goto base; + f64 wf = o2fG(w); + if (xe==el_bit) { + i32 wi = wf; if (wi!=wf) goto base; + u64* xp = bitarr_ptr(x); + 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) { 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; + } + if (rtid==n_add) { // + + u8 sel = xe - el_i8; + f64 r = sum_fns[sel](tyany_ptr(x), ia, wf); + decG(x); return m_f64(r); + } + if (rtid==n_floor) { f64 r=wf; if (ia>0) { f64 m=min_fns[xe-el_i8](tyany_ptr(x), ia); if (m0) { f64 m=max_fns[xe-el_i8](tyany_ptr(x), ia); if (m>r) r=m; } decG(x); return m_f64(r); } // ⌈ + i32 wi = wf; + if (rtid==n_mul | rtid==n_and) { // ×/∧ + void *xv = tyany_ptr(x); + bool isint = xe<=el_i32 && wi==wf; + u8 sel = xe - el_i8; + f64 r = isint ? prod_int_fns[sel](xv, ia, wi) + : prod_fns[sel](xv, ia, wf); + decG(x); return m_f64(r); + } + if (rtid==n_or && (wi&1)==wf) { // ∨ + if (xe==el_i8 ) { i8* xp = i8any_ptr (x); usz i=ia; if (!wi) while (i--) { i8 c=xp[i]; if (c==1) break; if (c!=0) goto base; } decG(x); return m_i32(i+1 > 0); } + if (xe==el_i16) { i16* xp = i16any_ptr(x); usz i=ia; if (!wi) while (i--) { i16 c=xp[i]; if (c==1) break; if (c!=0) goto base; } decG(x); return m_i32(i+1 > 0); } + if (xe==el_i32) { i32* xp = i32any_ptr(x); usz i=ia; if (!wi) while (i--) { i32 c=xp[i]; if (c==1) break; if (c!=0) goto base; } decG(x); return m_i32(i+1 > 0); } + } + } + base:; + SLOW3("𝕨 F´ 𝕩", w, x, f); + + B c = w; + SGet(x) + BBB2B fc2 = c2fn(f); + for (usz i = ia; i>0; i--) c = fc2(f, Get(x, i-1), c); + decG(x); + return c; +} diff --git a/src/builtins/grade.h b/src/builtins/grade.h index 66d9e66f..48da6f19 100644 --- a/src/builtins/grade.h +++ b/src/builtins/grade.h @@ -1,3 +1,34 @@ +// Ordering functions: Sort (∧∨), Grade (⍋⍒), Bins (⍋⍒) + +// Sort and Grade +// Small ≠𝕩 sort (not grade): insertion sort +// SHOULD implement insertion grade +// Small range sort (not grade): counting sort +// Count 1s for booleans +// Large range/length ("sparse") uses scan: plus, or Singeli max/min +// COULD range check 4-byte 𝕩 to try counting sort +// 4-byte grade checks range and sum +// If consistent with a permutation, grade as one and verify after +// If small-range, grade with bucket sort +// Other 1-, 2-, 4-byte cases: radix sort +// Singeli +` used if available, speeding up shorter cases +// SHOULD skip steps where all bytes are equal +// General case: Timsort +// SHOULD check for sorted flag and scan for sortedness in all cases +// SHOULD use an adaptive quicksort for 4- and 8-byte arguments +// SHOULD widen odd cell sizes under 8 bytes in sort and grade + +// Bins +// Mixed integer and character arguments gives all 0 or ≠𝕨 +// Integers and characters: 4-byte branchless binary search +// General case: branching binary search +// SHOULD implement f64 branchless binary search +// SHOULD interleave multiple branchless binary searches +// SHOULD specialize bins on equal types at least +// SHOULD implement table-based ⍋⍒ for small-range 𝕨 +// SHOULD special-case short 𝕨 +// SHOULD partition 𝕩 when 𝕨 is large + #define GRADE_CAT(N) CAT(GRADE_UD(gradeUp,gradeDown),N) #define GRADE_NEG GRADE_UD(,-) diff --git a/src/builtins/md1.c b/src/builtins/md1.c index 01cb30de..e5f9e5ed 100644 --- a/src/builtins/md1.c +++ b/src/builtins/md1.c @@ -107,143 +107,6 @@ B each_c2(Md1D* d, B w, B x) { B f = d->f; return homFil2(f, eachd(f, w, x), wf, xf); } -static bool fold_ne(u64* x, u64 am) { - u64 r = 0; - for (u64 i = 0; i < (am>>6); i++) r^= x[i]; - if (am&63) r^= x[am>>6]<<(64-am & 63); - return POPC(r) & 1; -} -static i64 bit_diff(u64* x, u64 am) { - i64 r = 0; - u64 a = 0xAAAAAAAAAAAAAAAA; - for (u64 i = 0; i < (am>>6); i++) r+= POPC(x[i]^a); - if (am&63) r+= POPC((x[am>>6]^a)<<(64-am & 63)); - return r - (i64)(am/2); -} -B fold_c1(Md1D* d, B x) { B f = d->f; - if (isAtm(x) || RNK(x)!=1) thrF("´: Argument must be a list (%H ≡ ≢𝕩)", x); - usz ia = IA(x); - if (ia==0) { - decG(x); - if (isFun(f)) { - B r = TI(f,identity)(f); - if (!q_N(r)) return inc(r); - } - thrM("´: No identity found"); - } - u8 xe = TI(x,elType); - if (isFun(f) && v(f)->flags && xe<=el_f64) { - u8 rtid = v(f)->flags-1; - if (xe==el_bit) { - 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) { 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; - } - if (rtid==n_add) { // + - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i64 c=0; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - if (xe==el_i16) { i16* xp = i16any_ptr(x); i16 c=I16_MIN; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=I32_MIN; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - } - if (rtid==n_or) { // ∨ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool r=0; for (usz i=0; i0; i--) c = fc2(f, m_i32(xp[i-1]), c); - } else { - c = Get(x, ia-1); - for (usz i = ia-1; i>0; i--) c = fc2(f, Get(x, i-1), c); - } - decG(x); - return c; -} -B fold_c2(Md1D* d, B w, B x) { B f = d->f; - if (isAtm(x) || RNK(x)!=1) thrF("´: 𝕩 must be a list (%H ≡ ≢𝕩)", x); - usz ia = IA(x); - u8 xe = TI(x,elType); - if (q_i32(w) && isFun(f) && v(f)->flags && elInt(xe)) { - i32 wi = o2iG(w); - u8 rtid = v(f)->flags-1; - if (xe==el_bit) { - u64* xp = bitarr_ptr(x); - 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) { 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; - } - if (rtid==n_add) { // + - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i64 c=wi; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=wi; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=wi; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - } - if (rtid==n_or && (wi&1)==wi) { // ∨ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool q=wi; for (usz i=0; i0; i--) c = fc2(f, Get(x, i-1), c); - decG(x); - return c; -} - B const_c1(Md1D* d, B x) { dec(x); return inc(d->f); } B const_c2(Md1D* d, B w, B x) { dec(w); dec(x); return inc(d->f); } @@ -392,6 +255,7 @@ B cell_c2(Md1D* d, B w, B x) { B f = d->f; return bqn_merge(r); } +extern B fold_c1(Md1D* d, B x); extern B rt_insert; B insert_c1(Md1D* d, B x) { B f = d->f; if (isAtm(x) || RNK(x)==0) thrM("˝: 𝕩 must have rank at least 1"); diff --git a/src/opt/single.c b/src/opt/single.c index 91f44ee4..7eb3dd86 100644 --- a/src/opt/single.c +++ b/src/opt/single.c @@ -24,6 +24,7 @@ #include "../builtins/sort.c" #include "../builtins/search.c" #include "../builtins/selfsearch.c" +#include "../builtins/fold.c" #include "../builtins/scan.c" #include "../builtins/arithm.c" #include "../builtins/arithd.c" diff --git a/src/singeli/src/fold.singeli b/src/singeli/src/fold.singeli new file mode 100644 index 00000000..7986a361 --- /dev/null +++ b/src/singeli/src/fold.singeli @@ -0,0 +1,41 @@ +include './base' +include './sse3' +include './avx' +include './avx2' +include './mask' + +fold_idem{T==f64, op}(x:*T, len:u64) : T = { + def step = 256/width{T} + def V = [step]T + xv:= *V ~~ x + r:V = broadcast{[4]f64, 0} + if (len < step) { + # Can't overlap like the long case + assert{len > 0} + r = load{xv} + if (len > 1) { + if (len > 2) r = op{r, shuf{[4]u64, r, 4b2222}} + r = op{r, shuf{[4]u64, r, 4b1111}} + } + } else { + # Pairwise combination to shorten dependency chains + def pairwise{p, i, k} = { + def l = k-1 + op{pairwise{p, i , l}, + pairwise{p, i+(1<> pk + @for (i to f) r = op{r, pairwise{xv+(i<