From b9a3a7886dddb0442a22ce6a9dd8645d9d3f88b2 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sun, 20 Nov 2022 19:57:00 -0500 Subject: [PATCH 01/13] Ordering function implementation comments --- src/builtins/grade.h | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) 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(,-) From 5bfdc26d00c894e347fdc3e4202e8a1df65405e3 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sun, 20 Nov 2022 20:06:35 -0500 Subject: [PATCH 02/13] Move fold to its own file --- makefile | 2 +- src/builtins/fold.c | 142 ++++++++++++++++++++++++++++++++++++++++++++ src/builtins/md1.c | 138 +----------------------------------------- src/opt/single.c | 1 + 4 files changed, 145 insertions(+), 138 deletions(-) create mode 100644 src/builtins/fold.c diff --git a/makefile b/makefile index cfe75915..b8649cf9 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 $< diff --git a/src/builtins/fold.c b/src/builtins/fold.c new file mode 100644 index 00000000..25a203f9 --- /dev/null +++ b/src/builtins/fold.c @@ -0,0 +1,142 @@ +#include "../core.h" +#include "../builtins.h" + +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; +} 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" From ebfd002793cf2006ddd93c51aa65e23805070961 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sun, 20 Nov 2022 21:06:05 -0500 Subject: [PATCH 03/13] =?UTF-8?q?Blocked=20+=C2=B4=20and=20f64=20=E2=8C=88?= =?UTF-8?q?=C2=B4=20=E2=8C=8A=C2=B4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/fold.c | 62 ++++++++++++++++++++++++++++++--------------- 1 file changed, 41 insertions(+), 21 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 25a203f9..9a5a3aef 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -32,19 +32,37 @@ B fold_c1(Md1D* d, B x) { B f = d->f; 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; + 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; + } + decG(x); return m_f64(r); } if (rtid==n_add) { // + - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i64 c=0; for (usz i=0; i0 && -lim<=c && c<=lim) { \ + usz e=i0; i0=(i0-1)&~(b-1); \ + i##A s=0; for (usz i=i0; if; if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } if (xe==el_f64) { f64* xp = f64any_ptr(x); f64 c=1; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } } - if (rtid==n_floor) { // ⌊ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i8 c=I8_MAX ; 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); } - } + #define CASE(T,C) case el_##T: { \ + T* xp = xv; T c = xp[0]; \ + for (usz i=0; i) } // ⌈ + #undef FC + #undef CASE if (rtid==n_or) { // ∨ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool r=0; for (usz i=0; i Date: Sun, 20 Nov 2022 21:51:30 -0500 Subject: [PATCH 04/13] =?UTF-8?q?Move=20+=E2=8C=88=E2=8C=8A=20folds=20to?= =?UTF-8?q?=20function=20tables,=20unifying=20monadic=20and=20dyadic=20cas?= =?UTF-8?q?es?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/fold.c | 98 ++++++++++++++++++++++++++------------------- 1 file changed, 56 insertions(+), 42 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 9a5a3aef..1adb97f4 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -16,6 +16,51 @@ static i64 bit_diff(u64* x, u64 am) { 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 }; + +#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) DEF_MIN_MAX(f64) +#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); @@ -45,23 +90,10 @@ B fold_c1(Md1D* d, B x) { B f = d->f; } if (rtid==n_add) { // + void *xv = tyany_ptr(x); - f64 r; - #define CASE_INT(T,M,A) case el_i##T: { \ - usz b=1<<(M-T); i64 lim = (1ull<<53) - (1ull<0 && -lim<=c && c<=lim) { \ - usz e=i0; i0=(i0-1)&~(b-1); \ - i##A s=0; for (usz i=i0; if; if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } if (xe==el_f64) { f64* xp = f64any_ptr(x); f64 c=1; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } } - #define CASE(T,C) case el_##T: { \ - T* xp = xv; T c = xp[0]; \ - for (usz i=0; i) } // ⌈ - #undef FC - #undef CASE + 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_or) { // ∨ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool r=0; for (usz i=0; if; 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_floor) { f64 r=wi; 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); } // ⌈ if (rtid==n_or && (wi&1)==wi) { // ∨ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool q=wi; for (usz i=0; i Date: Sun, 20 Nov 2022 22:00:40 -0500 Subject: [PATCH 05/13] Handle some f64 cases in dyadic fold --- src/builtins/fold.c | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 1adb97f4..d3e23275 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -96,14 +96,14 @@ B fold_c1(Md1D* d, B x) { B f = d->f; : 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) { // ×/∧ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_f64(c); } if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } if (xe==el_f64) { f64* xp = f64any_ptr(x); f64 c=1; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } } - 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_or) { // ∨ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool r=0; for (usz i=0; if; 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); + if (isF64(w) && isFun(f) && v(f)->flags && xe<=el_f64) { + f64 wf = w.f; u8 rtid = v(f)->flags-1; 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; } @@ -148,17 +149,19 @@ B fold_c2(Md1D* d, B w, B x) { B f = d->f; } if (rtid==n_add) { // + u8 sel = xe - el_i8; - f64 r = sum_fns[sel](tyany_ptr(x), ia, wi); + 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) { // ×/∧ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } + if (xe==el_i8 && wi==wf) { i8* xp = i8any_ptr (x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } + if (xe==el_i16 && wi==wf) { i16* xp = i16any_ptr(x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } + if (xe==el_i32 && wi==wf) { i32* xp = i32any_ptr(x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } + if (xe==el_f64 ) { f64* xp = f64any_ptr(x); f64 c=wf; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } } - if (rtid==n_floor) { f64 r=wi; 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); } // ⌈ - if (rtid==n_or && (wi&1)==wi) { // ∨ + if (rtid==n_or && (wi&1)==wf) { // ∨ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool q=wi; for (usz i=0; i Date: Mon, 21 Nov 2022 08:07:38 -0500 Subject: [PATCH 06/13] Do tack-fold with selection --- src/builtins/fold.c | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index d3e23275..4c55020a 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -73,8 +73,11 @@ B fold_c1(Md1D* d, B x) { B f = d->f; thrM("´: No identity found"); } u8 xe = TI(x,elType); - if (isFun(f) && v(f)->flags && xe<=el_f64) { + 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; @@ -132,9 +135,16 @@ 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 (isF64(w) && isFun(f) && v(f)->flags && xe<=el_f64) { - f64 wf = w.f; + 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 = w.f; if (xe==el_bit) { i32 wi = wf; if (wi!=wf) goto base; u64* xp = bitarr_ptr(x); From 52608ba6ce1a2346c64a9b7e4a19c324aa4694e3 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Mon, 21 Nov 2022 11:47:50 -0500 Subject: [PATCH 07/13] Singeli min/max-fold --- makefile | 2 +- src/builtins/fold.c | 15 ++++++++++++- src/singeli/src/fold.singeli | 41 ++++++++++++++++++++++++++++++++++++ 3 files changed, 56 insertions(+), 2 deletions(-) create mode 100644 src/singeli/src/fold.singeli diff --git a/makefile b/makefile index b8649cf9..6c5a5aba 100644 --- a/makefile +++ b/makefile @@ -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 index 4c55020a..fcced21e 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -1,6 +1,13 @@ #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]; @@ -55,7 +62,13 @@ static f64 (*const sum_fns[])(void*, usz, f64) = { sum_i8, sum_i16, sum_i32, sum #define DEF_MIN_MAX(T) \ static f64 min_##T(void* xv, usz ia) { MIN_MAX(T,<) } \ static f64 max_##T(void* xv, usz ia) { MIN_MAX(T,>) } -DEF_MIN_MAX(i8) DEF_MIN_MAX(i16) DEF_MIN_MAX(i32) DEF_MIN_MAX(f64) +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 }; 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< Date: Mon, 21 Nov 2022 15:19:27 -0500 Subject: [PATCH 08/13] =?UTF-8?q?Implement=20=C3=97=C2=B4=20for=20all=20nu?= =?UTF-8?q?meric=20cases,=20using=20f64=20product=20when=20i32=20fails?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/fold.c | 43 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index fcced21e..86c83ca3 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -55,6 +55,30 @@ static f64 sum_f64(void* xv, usz i, f64 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 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; if; 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) { // ×/∧ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_f64(c); } - if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_f64) { f64* xp = f64any_ptr(x); f64 c=1; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } + 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); bool r=0; for (usz i=0; if; if (rtid==n_ceil ) { f64 r=wf; if (ia>0) { 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) { // ×/∧ - if (xe==el_i8 && wi==wf) { i8* xp = i8any_ptr (x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_i16 && wi==wf) { i16* xp = i16any_ptr(x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_i32 && wi==wf) { i32* xp = i32any_ptr(x); i32 c=wi; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } - if (xe==el_f64 ) { f64* xp = f64any_ptr(x); f64 c=wf; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } + 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); bool q=wi; for (usz i=0; i Date: Mon, 21 Nov 2022 15:36:06 -0500 Subject: [PATCH 09/13] Boolean comparison folds, monadic only --- src/builtins/fold.c | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 86c83ca3..e687498b 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -119,12 +119,16 @@ B fold_c1(Md1D* d, B x) { B f = d->f; 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_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); } From bda653e86542bf45b65cbfd3901a5cf23d505c27 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Mon, 21 Nov 2022 15:51:03 -0500 Subject: [PATCH 10/13] =?UTF-8?q?Shortcutting=20for=20boolean-valued=20int?= =?UTF-8?q?eger=20=E2=88=A8=C2=B4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/fold.c | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index e687498b..553ea173 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -150,9 +150,9 @@ B fold_c1(Md1D* d, B x) { B f = d->f; decG(x); return m_f64(r); } if (rtid==n_or) { // ∨ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool r=0; for (usz i=0; i 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:; @@ -216,9 +216,9 @@ B fold_c2(Md1D* d, B w, B x) { B f = d->f; decG(x); return m_f64(r); } if (rtid==n_or && (wi&1)==wf) { // ∨ - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool q=wi; for (usz i=0; i 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:; From 31106c730c62c6fa37f9547e9e48f1622324ba3b Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Mon, 21 Nov 2022 15:56:15 -0500 Subject: [PATCH 11/13] Fold implementation comments --- src/builtins/fold.c | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 553ea173..ac56d18a 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -1,3 +1,14 @@ +// 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" From e7cc8ec444e900740531fad92beff438ae338e5d Mon Sep 17 00:00:00 2001 From: dzaima Date: Mon, 21 Nov 2022 23:26:25 +0200 Subject: [PATCH 12/13] use o2fG instead of union field --- src/builtins/fold.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index ac56d18a..66e7aaca 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -197,7 +197,7 @@ B fold_c2(Md1D* d, B w, B x) { B f = d->f; } if (rtid==n_rtack) { decG(x); return w; } if (!isF64(w) || xe>el_f64) goto base; - f64 wf = w.f; + f64 wf = o2fG(w); if (xe==el_bit) { i32 wi = wf; if (wi!=wf) goto base; u64* xp = bitarr_ptr(x); From 6f5c329e22f73419556170614e218826108eaeb6 Mon Sep 17 00:00:00 2001 From: dzaima Date: Tue, 22 Nov 2022 14:24:43 +0200 Subject: [PATCH 13/13] noinline float product --- src/builtins/fold.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 66e7aaca..8937050a 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -68,7 +68,7 @@ static f64 (*const sum_fns[])(void*, usz, f64) = { sum_i8, sum_i16, sum_i32, sum // Try to keep to i32 product, go to f64 on overflow or non-i32 initial #define DEF_INT_PROD(T) \ - static f64 prod_##T(void* xv, usz i, f64 init) { \ + 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) { \