Merge pull request #61 from mlochbaum/fold

Various Fold improvements
This commit is contained in:
dzaima 2022-11-22 15:15:26 +02:00 committed by GitHub
commit 68a525b8ba
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 320 additions and 139 deletions

View File

@ -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/^/ /'

244
src/builtins/fold.c Normal file
View File

@ -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; i<ia; i++) s+=((T*)xv)[i]; return s; \
} \
static f64 sum_##T(void* xv, usz ia, f64 init) { \
usz b=1<<(M-W); i64 lim = (1ull<<53) - (1ull<<M); \
T* xp = xv; \
f64 r=init; i64 c=init; usz i0=ia; \
if (c == init) { \
while (i0>0 && -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<ia; i++) if (xp[i] C r) r=xp[i]; \
return r;
#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)
#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 (m<r) r=m; } decG(x); return m_f64(r); } // ⌊
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) { // ×/∧
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;
}

View File

@ -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(,-)

View File

@ -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; i<ia; i++) c+= xp[i]; decG(x); return m_f64(c); } // won't worry about 64TB array sum float inaccuracy for now
if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=0; for (usz i=0; i<ia; i++) if (addOn(c,xp[i]))goto base; decG(x); return m_i32(c); }
if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=0; for (usz i=0; i<ia; i++) if (addOn(c,xp[i]))goto base; decG(x); return m_i32(c); }
if (xe==el_f64) { f64* xp = f64any_ptr(x); f64 c=0; for (usz i=ia; i--; ) c+= xp[i]; decG(x); return m_f64(c); }
}
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) { // ⌊
if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i8 c=I8_MAX ; for (usz i=0; i<ia; i++) if (xp[i]<c) c=xp[i]; decG(x); return m_i32(c); }
if (xe==el_i16) { i16* xp = i16any_ptr(x); i16 c=I16_MAX; for (usz i=0; i<ia; i++) if (xp[i]<c) c=xp[i]; decG(x); return m_i32(c); }
if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=I32_MAX; for (usz i=0; i<ia; i++) if (xp[i]<c) c=xp[i]; decG(x); return m_i32(c); }
}
if (rtid==n_ceil) { // ⌈
if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i8 c=I8_MIN ; for (usz i=0; i<ia; i++) if (xp[i]>c) 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; i<ia; i++) if (xp[i]>c) 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; i<ia; i++) if (xp[i]>c) 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; i<ia; i++) { i8 c=xp[i]; if (c!=0&&c!=1)goto base; r|=c; } decG(x); return m_i32(r); }
if (xe==el_i16) { i16* xp = i16any_ptr(x); bool r=0; for (usz i=0; i<ia; i++) { i16 c=xp[i]; if (c!=0&&c!=1)goto base; r|=c; } decG(x); return m_i32(r); }
if (xe==el_i32) { i32* xp = i32any_ptr(x); bool r=0; for (usz i=0; i<ia; i++) { i32 c=xp[i]; if (c!=0&&c!=1)goto base; r|=c; } decG(x); return m_i32(r); }
}
}
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 (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; i<ia; i++) c+=xp[i]; decG(x); return m_f64(c); }
if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=wi; for (usz i=0; i<ia; i++) if (addOn(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=0; i<ia; i++) if (addOn(c,xp[i]))goto base; decG(x); return m_i32(c); }
}
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 (rtid==n_floor) { // ⌊
if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i32 c=wi; for (usz i=0; i<ia; i++) if (xp[i]<c) 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; i<ia; i++) if (xp[i]<c) 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; i<ia; i++) if (xp[i]<c) c=xp[i]; decG(x); return m_i32(c); }
}
if (rtid==n_ceil) { // ⌈
if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i32 c=wi; for (usz i=0; i<ia; i++) if (xp[i]>c) 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; i<ia; i++) if (xp[i]>c) 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; i<ia; i++) if (xp[i]>c) 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; i<ia; i++) { i8 c=xp[i]; if (c!=0&&c!=1)goto base; q|=c; } decG(x); return m_i32(q); }
if (xe==el_i16) { i16* xp = i16any_ptr(x); bool q=wi; for (usz i=0; i<ia; i++) { i16 c=xp[i]; if (c!=0&&c!=1)goto base; q|=c; } decG(x); return m_i32(q); }
if (xe==el_i32) { i32* xp = i32any_ptr(x); bool q=wi; for (usz i=0; i<ia; i++) { i32 c=xp[i]; if (c!=0&&c!=1)goto base; q|=c; } decG(x); return m_i32(q); }
}
}
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;
}
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");

View File

@ -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"

View File

@ -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<<l), l}}
}
def pairwise{p, i, k==0} = load{p, i}
def pk = 2 # Combine 1<<pk values in a step
r = load{*V ~~ (x+len-step)}
e:= (len-1)/step
f:= e >> pk
@for (i to f) r = op{r, pairwise{xv+(i<<pk), 0, pk}}
@for (xv over i from f<<pk to e) r = op{r, xv}
r = op{r, shuf{[4]u64, r, 4b2301}}
r = op{r, shuf{[4]u64, r, 4b1032}}
}
extract{r, 0}
}
'avx2_fold_min_f64' = fold_idem{f64,min}
'avx2_fold_max_f64' = fold_idem{f64,max}