diff --git a/makefile b/makefile index 34d354ab..d83de3b5 100644 --- a/makefile +++ b/makefile @@ -212,7 +212,7 @@ ${bd}/%.o: src/jit/%.c @echo $< | cut -c 5- @$(CC_INC) $@.d -o $@ -c $< -builtins: ${addprefix ${bd}/, arithm.o arithd.o cmp.o sfns.o squeeze.o select.o slash.o group.o sort.o search.o selfsearch.o md1.o md2.o fns.o sysfn.o internal.o inverse.o} +builtins: ${addprefix ${bd}/, arithm.o arithd.o cmp.o sfns.o squeeze.o select.o slash.o group.o sort.o search.o selfsearch.o 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/grade.h b/src/builtins/grade.h index f59f73ec..66d9e66f 100644 --- a/src/builtins/grade.h +++ b/src/builtins/grade.h @@ -31,22 +31,22 @@ } #if SINGELI -extern void (*const avx2_scan_max8)(int8_t* v0,int8_t* v1,uint64_t v2); -extern void (*const avx2_scan_min8)(int8_t* v0,int8_t* v1,uint64_t v2); -extern void (*const avx2_scan_max16)(int16_t* v0,int16_t* v1,uint64_t v2); -extern void (*const avx2_scan_min16)(int16_t* v0,int16_t* v1,uint64_t v2); +extern void (*const avx2_scan_max_i8)(int8_t* v0,int8_t* v1,uint64_t v2); +extern void (*const avx2_scan_min_i8)(int8_t* v0,int8_t* v1,uint64_t v2); +extern void (*const avx2_scan_max_i16)(int16_t* v0,int16_t* v1,uint64_t v2); +extern void (*const avx2_scan_min_i16)(int16_t* v0,int16_t* v1,uint64_t v2); #define COUNT_THRESHOLD 32 #define WRITE_SPARSE_i8 \ for (usz i=0; ixn) w = take_c2(m_f64(0), m_f64(xn), w); - rp[1] = slash_c2(m_f64(0), inc(w), inc(x)); - rp[0] = slash_c2(m_f64(0), not_c1(m_f64(0), w), x); + if (wia>xn) w = C2(take, m_f64(xn), w); + rp[1] = C2(slash, inc(w), inc(x)); + rp[0] = C2(slash, bit_negate(w), x); return taga(r); } // Needed to make sure wia>0 for ip[wia-1] below @@ -128,9 +127,6 @@ static B group_simple(B w, B x, ur xr, usz wia, usz xn, usz* xsh, u8 we) { // Few changes in 𝕨: move in chunks if (xn>64 && notB && change<(xn*width)/32) { - #define C1(F,X ) F##_c1(m_f64(0),X ) - #define C2(F,X,W) F##_c2(m_f64(0),X,W) - u64* mp; B m = m_bitarrv(&mp, xn); u8* wp0 = tyany_ptr(w); we = TI(w,elType); @@ -139,8 +135,6 @@ static B group_simple(B w, B x, ur xr, usz wia, usz xn, usz* xsh, u8 we) { B ind = C1(slash, m); w = C2(select, inc(ind), w); - #undef C1 - #undef C2 if (TI(ind,elType)!=el_i32) ind = taga(cpyI32Arr(ind)); if (TI(w ,elType)!=el_i32) w = taga(cpyI32Arr(w )); wia = IA(ind); @@ -179,10 +173,10 @@ static B group_simple(B w, B x, ur xr, usz wia, usz xn, usz* xsh, u8 we) { // Many Β―1s: filter out, then continue if (xn>32 && neg>(bits?0:xn/4)+xn/8) { - if (wia>xn) w = take_c2(m_f64(0), m_f64(xn), w); - B m = ne_c2(m_f64(0), m_f64(-1), inc(w)); - w = slash_c2(m_f64(0), inc(m), w); - x = slash_c2(m_f64(0), m, x); xn = *SH(x); + if (wia>xn) w = C2(take, m_f64(xn), w); + B m = C2(ne, m_f64(-1), inc(w)); + w = C2(slash, inc(m), w); + x = C2(slash, m, x); xn = *SH(x); neg = 0; } if (TI(w,elType)!=el_i32) w = taga(cpyI32Arr(w)); @@ -331,8 +325,8 @@ B group_c2(B t, B w, B x) { B group_c1(B t, B x) { if (isArr(x) && RNK(x)==1 && TI(x,arrD1)) { usz ia = IA(x); - B range = ud_c1(t, m_f64(ia)); - return group_c2(m_f64(0), x, range); + B range = C1(ud, m_f64(ia)); + return C2(group, x, range); } return c1(rt_group, x); } diff --git a/src/builtins/md1.c b/src/builtins/md1.c index 4b3c4d9b..01cb30de 100644 --- a/src/builtins/md1.c +++ b/src/builtins/md1.c @@ -107,29 +107,6 @@ B each_c2(Md1D* d, B w, B x) { B f = d->f; return homFil2(f, eachd(f, w, x), wf, xf); } -#if SINGELI && __PCLMUL__ - #pragma GCC diagnostic push - #pragma GCC diagnostic ignored "-Wunused-variable" - #include "../singeli/gen/neq.c" - #pragma GCC diagnostic pop -#endif -B scan_ne(u64 p, u64* xp, u64 ia) { - u64* rp; B r=m_bitarrv(&rp,ia); -#if SINGELI && __PCLMUL__ - clmul_scan_ne(p, xp, rp, BIT_N(ia)); -#else - for (usz i = 0; i < BIT_N(ia); i++) { - u64 c = xp[i]; - u64 r = c ^ (c<<1); - r^= r<< 2; r^= r<< 4; r^= r<<8; - r^= r<<16; r^= r<<32; r^= p; - rp[i] = r; - p = -(r>>63); // repeat sign bit - } -#endif - return r; -} - static bool fold_ne(u64* x, u64 am) { u64 r = 0; for (u64 i = 0; i < (am>>6); i++) r^= x[i]; @@ -143,204 +120,6 @@ static i64 bit_diff(u64* x, u64 am) { if (am&63) r+= POPC((x[am>>6]^a)<<(64-am & 63)); return r - (i64)(am/2); } - - -#if SINGELI - #pragma GCC diagnostic push - #pragma GCC diagnostic ignored "-Wunused-variable" - #include "../singeli/gen/scan.c" - #pragma GCC diagnostic pop -#endif -B slash_c1(B f, B x); -B scan_bit_sum(B x, u64* xp, u64 ia, u64 xs) { // consumes x - u8 re = xs<=I8_MAX? el_i8 : xs<=I16_MAX? el_i16 : el_i32; - if (xs < ia/128) { - B ones = slash_c1(m_f64(0), x); - MAKE_MUT(r0, ia) mut_init(r0, re); MUTG_INIT(r0); - SGetU(ones) - usz ri = 0; - for (usz i = 0; i < xs; i++) { - usz e = o2s(GetU(ones, i)); - mut_fillG(r0, ri, m_i32(i), e-ri); - ri = e; - } - if (rif; - if (isAtm(x) || RNK(x)==0) thrM("`: Argument cannot have rank 0"); - ur xr = RNK(x); - usz ia = IA(x); - if (ia==0) return x; - B xf = getFillQ(x); - u8 xe = TI(x,elType); - if (xr==1 && xe<=el_f64 && isFun(f) && v(f)->flags) { - u8 rtid = v(f)->flags-1; - if (xe==el_bit) { - u64* xp=bitarr_ptr(x); - if (rtid==n_add) { - u64 xs = bit_sum(xp, ia); - if (xs>I32_MAX) goto base; - if (xs<=1) { if (xs==0) return x; goto bit_or; } - return FL_SET(scan_bit_sum(x, xp, ia, xs), fl_asc|fl_squoze); - } - if (rtid==n_or | rtid==n_ceil ) { bit_or:; u64* rp; B r=m_bitarrv(&rp,ia); usz n=BIT_N(ia); u64 xi; usz i=0; while(i>63); - rp[i] = p = x & (m10 ^ (x + c)); - } - decG(x); return r; - } - goto base; - } - if (rtid==n_add) { // + - if (xe==el_i8 ) { i8* xp=i8any_ptr (x); i32* rp; B r=m_i32arrv(&rp, ia); i32 c=0; for (usz i=0; ic)c=xp[i]; rp[i]=c; } decG(x); return r; } - if (xe==el_i16) { i16* xp=i16any_ptr(x); i16* rp; B r=m_i16arrv(&rp, ia); i16 c=I16_MIN; for (usz i=0; ic)c=xp[i]; rp[i]=c; } decG(x); return r; } - if (xe==el_i32) { i32* xp=i32any_ptr(x); i32* rp; B r=m_i32arrv(&rp, ia); i32 c=I32_MIN; for (usz i=0; ic)c=xp[i]; rp[i]=c; } decG(x); return r; } - #endif - } - if (rtid==n_ne) { // β‰  - f64 x0 = IGetU(x,0).f; if (x0!=0 && x0!=1) goto base; - if (xe==el_i8 ) { i8* xp=i8any_ptr (x); u64* rp; B r=m_bitarrv(&rp,ia); bool c=x0; rp[0]=c; for (usz i=1; if; - if (isAtm(x) || RNK(x)==0) thrM("`: 𝕩 cannot have rank 0"); - ur xr = RNK(x); usz* xsh = SH(x); usz ia = IA(x); - B wf = getFillQ(w); - u8 xe = TI(x,elType); - if (xr==1 && q_i32(w) && elInt(xe) && isFun(f) && v(f)->flags) { - u8 rtid = v(f)->flags-1; - i32 wv = o2iG(w); - if (xe==el_bit) { - u64* xp=bitarr_ptr(x); - if (rtid==n_add) { i32* rp; B r=m_i32arrv(&rp, ia); i64 c=wv; for (usz i=0; ic)c=xp[i]; rp[i]=c; } decG(x); return r; } - if (xe==el_i16 && wv==(i16)wv) { i16* xp=i16any_ptr(x); i16* rp; B r=m_i16arrv(&rp, ia); i16 c=wv; for (usz i=0; ic)c=xp[i]; rp[i]=c; } decG(x); return r; } - if (xe==el_i32 && wv==(i32)wv) { i32* xp=i32any_ptr(x); i32* rp; B r=m_i32arrv(&rp, ia); i32 c=wv; for (usz i=0; ic)c=xp[i]; rp[i]=c; } decG(x); return r; } - } - if (rtid==n_ne) { // β‰  - if (!q_ibit(wv)) { goto base; } bool c=wv; - if (xe==el_i8 ) { i8* xp=i8any_ptr (x); u64* rp; B r=m_bitarrv(&rp, ia); for (usz i=0; if; if (isAtm(x) || RNK(x)!=1) thrF("Β΄: Argument must be a list (%H ≑ ≒𝕩)", x); usz ia = IA(x); diff --git a/src/builtins/scan.c b/src/builtins/scan.c new file mode 100644 index 00000000..af61f243 --- /dev/null +++ b/src/builtins/scan.c @@ -0,0 +1,261 @@ +#include "../core.h" +#include "../utils/mut.h" +#include "../builtins.h" +#include +#define F64_MIN -INFINITY +#define F64_MAX INFINITY + +#if !USE_VALGRIND +static u64 vg_rand(u64 x) { return x; } +#endif + +#if SINGELI + #pragma GCC diagnostic push + #pragma GCC diagnostic ignored "-Wunused-variable" + #include "../singeli/gen/scan.c" + #pragma GCC diagnostic pop +#endif + +#if SINGELI && __PCLMUL__ + #pragma GCC diagnostic push + #pragma GCC diagnostic ignored "-Wunused-variable" + #include "../singeli/gen/neq.c" + #pragma GCC diagnostic pop +#endif +B scan_ne(B x, u64 p, u64 ia) { // consumes x + u64* xp = bitarr_ptr(x); + u64* rp; B r=m_bitarrv(&rp,ia); +#if SINGELI && __PCLMUL__ + clmul_scan_ne(p, xp, rp, BIT_N(ia)); +#else + for (usz i = 0; i < BIT_N(ia); i++) { + u64 c = xp[i]; + u64 r = c ^ (c<<1); + r^= r<< 2; r^= r<< 4; r^= r<<8; + r^= r<<16; r^= r<<32; r^= p; + rp[i] = r; + p = -(r>>63); // repeat sign bit + } +#endif + decG(x); return r; +} + +static B scan_or(B x, u64 ia) { // consumes x + u64* xp = bitarr_ptr(x); + u64* rp; B r=m_bitarrv(&rp,ia); + usz n=BIT_N(ia); u64 xi; usz i=0; + while (i,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; } \ + T* xp=T##any_ptr(x); T* rp; r=m_##T##arrv(&rp, ia); MINMAX_SCAN(T,N,C,wv); \ + break; } +#define MINMAX2(NAME,C,INIT,BIT,BI,ORD) \ + i32 wv=0; if (q_i32(w)) { wv=w.f; } else { x=taga(cpyF64Arr(x)); xe=el_f64; } \ + B r; switch (xe) { default:UD; \ + case el_bit: if (wv C BI) r=C2(shape,m_f64(ia),w); else return scan_##BIT(x, ia); break; \ + MM2_ICASE(i8 ,NAME,C,I8_##INIT ) \ + MM2_ICASE(i16,NAME,C,I16_##INIT) \ + MM_CASE(i32,NAME,C,wv) \ + MM_CASE(f64,NAME,C,w.f) \ + } \ + decG(x); return FL_SET(r, fl_##ORD); +static B scan2_min_num(B w, B x, u8 xe, usz ia) { MINMAX2(min,<,MAX,and,1,dsc) } +static B scan2_max_num(B w, B x, u8 xe, usz ia) { MINMAX2(max,>,MIN,or ,0,asc) } +#undef MINMAX2 +#undef MM_CASE +#undef MM_CASE2 +#undef MINMAX_SCAN + +B scan_c1(Md1D* d, B x) { B f = d->f; + if (isAtm(x) || RNK(x)==0) thrM("`: Argument cannot have rank 0"); + ur xr = RNK(x); + usz ia = IA(x); + if (ia==0) return x; + B xf = getFillQ(x); + u8 xe = TI(x,elType); + if (xr==1 && xe<=el_f64 && isFun(f) && v(f)->flags) { + u8 rtid = v(f)->flags-1; + if (xe==el_bit) { + u64* xp=bitarr_ptr(x); + 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) { + u64* rp; B r=m_bitarrv(&rp,ia); usz n=BIT_N(ia); + u64 m10 = 0x5555555555555555; + u64 p = 0; + for (usz i=0; i>63); + rp[i] = p = x & (m10 ^ (x + c)); + } + decG(x); return r; + } + goto base; + } + if (rtid==n_add) { // + + if (xe==el_i8 ) { i8* xp=i8any_ptr (x); i32* rp; B r=m_i32arrv(&rp, ia); i32 c=0; for (usz i=0; if; + if (isAtm(x) || RNK(x)==0) thrM("`: 𝕩 cannot have rank 0"); + ur xr = RNK(x); usz* xsh = SH(x); usz ia = IA(x); + B wf = getFillQ(w); + u8 xe = TI(x,elType); + if (xr==1 && xe<=el_f64 && isFun(f) && v(f)->flags) { + u8 rtid = v(f)->flags-1; + if (rtid==n_floor) return scan2_min_num(w, x, xe, ia); // ⌊ + if (rtid==n_ceil ) return scan2_max_num(w, x, xe, ia); // ⌈ + if (!q_i32(w)) goto base; + i32 wv = o2iG(w); + if (xe==el_bit) { + u64* xp=bitarr_ptr(x); + if (rtid==n_add) { i32* rp; B r=m_i32arrv(&rp, ia); i64 c=wv; for (usz i=0; i(usz)I32_MAX+1) thrM("βŠ’: Argument length >2⋆31 not supported"); u8 lw = cellWidthLog(x); - if (lw==0) { x = toI8Any(x); lw = cellWidthLog(x); } + if (lw==0) { + u64* xp = bitarr_ptr(x); + B r; + #define COUNT_BOOL(T) \ + T* rp; r = m_##T##arrv(&rp, n); \ + usz n1 = 0; \ + for (usz i=0; i>=1; \ + rp[i] = b? n1 : i-n1; \ + n1 += b; \ + } \ + } + if (n <= 128) { COUNT_BOOL(i8) } + else if (n <= 1<<15) { COUNT_BOOL(i16) } + else { COUNT_BOOL(i32) } + decG(x); return r; + #undef COUNT_BOOL + } + if (use_sorted(x, lw) && n>16 && (lw>4 || n<1<<16)) { // β†•βˆ˜β‰ (⊣-⌈`βˆ˜Γ—)∊ + B c = shift_ne(x, n, lw, 1); + B i = C1(ud, m_f64(n)); + B m = C2(mul, c, inc(i)); + return C2(sub, i, scan_max_num(m, TI(m,elType), n)); + } void* xv = tyany_ptr(x); #define BRUTE(T) \ i##T* xp = xv; \ @@ -280,7 +360,7 @@ B count_c1(B t, B x) { /*THRESHMUL*/1, THRESH, \ /*INIT*/u32 ctr0 = 1;) if (lw==5) { - if (n<12) { BRUTE(32); } + if (n<20) { BRUTE(32); } i32* rp; B r = m_i32arrv(&rp, n); HASHTAB(u32, 32, 1, n/2, sz==msz? 1 : sz>=(1<<14)? 3 : 5) // Radix-assisted lookup @@ -334,8 +414,11 @@ B indexOf_c1(B t, B x) { u8 lw = cellWidthLog(x); void* xv = tyany_ptr(x); if (lw == 0) { - B r = 1&*(u64*)xv ? not_c1(m_f64(0), x) : x; - return shape_c1(m_f64(0), r); + B r = 1&*(u64*)xv ? bit_negate(x) : x; + return C1(shape, r); + } + if (use_sorted(x, lw) && n>8) { + return scan_add_bool(shift_ne(x, n, lw, 0), n); } #define BRUTE(T) \ i##T* xp = xv; \ @@ -426,5 +509,13 @@ B find_c1(B t, B x) { if (isAtm(x) || RNK(x)==0) thrM("⍷: Argument cannot have rank 0"); usz n = *SH(x); if (n<=1) return x; - return slash_c2(m_f64(0), memberOf_c1(m_f64(0), inc(x)), x); + if (TI(x,elType)==el_bit && RNK(x)==1) { + u64* xp = bitarr_ptr(x); + u64 x0 = 1 & *xp; + usz i = bit_find(xp, n, !x0); decG(x); + u64* rp; B r = m_bitarrv(&rp, 1 + (ish{4b1110,sh{4b2200,a}}} +} +'avx2_scan_min_init_i8' = avx2_scan_idem{i8 , min}; 'avx2_scan_max_init_i8' = avx2_scan_idem{i8 , max} +'avx2_scan_min_init_i16' = avx2_scan_idem{i16, min}; 'avx2_scan_max_init_i16' = avx2_scan_idem{i16, max} +'avx2_scan_min_init_i32' = avx2_scan_idem{i32, min}; 'avx2_scan_max_init_i32' = avx2_scan_idem{i32, max} +'avx2_scan_min_init_f64' = avx2_scan_idem{f64, min}; 'avx2_scan_max_init_f64' = avx2_scan_idem{f64, max} + +avx2_scan_idem_id{T, op}(x:*T, r:*T, len:u64) : void = { def m = 1 << (width{T}-1) - avx2_scan_idem{T, op, (if (match{op,min}) m-1; else -m)} + def id = (if (match{op,min}) m-1; else -m) + avx2_scan_idem{T, op}(x, r, len, id) } -'avx2_scan_min8' = avx2_scan_idem{i8 , min} -'avx2_scan_max8' = avx2_scan_idem{i8 , max} -'avx2_scan_min16' = avx2_scan_idem{i16, min} -'avx2_scan_max16' = avx2_scan_idem{i16, max} -'avx2_scan_min32' = avx2_scan_idem{i32, min} -'avx2_scan_max32' = avx2_scan_idem{i32, max} +'avx2_scan_min_i8' = avx2_scan_idem_id{i8 , min}; 'avx2_scan_max_i8' = avx2_scan_idem_id{i8 , max} +'avx2_scan_min_i16' = avx2_scan_idem_id{i16, min}; 'avx2_scan_max_i16' = avx2_scan_idem_id{i16, max} +'avx2_scan_min_i32' = avx2_scan_idem_id{i32, min}; 'avx2_scan_max_i32' = avx2_scan_idem_id{i32, max} # Associative scan avx2_scan_assoc_0{T, op}(x:*T, r:*T, len:u64, init:T) : void = { diff --git a/src/utils/calls.h b/src/utils/calls.h index c9cc4d1b..b0bd741d 100644 --- a/src/utils/calls.h +++ b/src/utils/calls.h @@ -1,5 +1,12 @@ #pragma once +#define C1(F, X) F##_c1(m_f64(0), X) +#define C2(F,W,X) F##_c2(m_f64(0),W,X) +#define M1C1(M,F, X) m1c1_unsafe(M##_c1, bi_##F, X) +#define M1C2(M,F,W,X) m1c2_unsafe(M##_c2, bi_##F, W, X) +static inline B m1c1_unsafe(D1C1 m, B f, B x) { Md1D d; d.f=f; return m(&d, x); } +static inline B m1c2_unsafe(D1C2 m, B f, B w, B x) { Md1D d; d.f=f; return m(&d, w, x); } + typedef void (*M_CopyF)(void*, usz, B, usz, usz); typedef void (*M_FillF)(void*, usz, B, usz); extern M_CopyF copyFns[el_MAX]; @@ -28,4 +35,4 @@ CMP_DEF(le, AS); #define CMP_AA_IMM(FN, ELT, WHERE, WP, XP, LEN) CMP_AA_CALL(CMP_AA_FN(FN, ELT), WHERE, WP, XP, LEN) #define CMP_AS_IMM(FN, ELT, WHERE, WP, X, LEN) CMP_AS_CALL(CMP_AS_FN(FN, ELT), WHERE, WP, X, LEN) -void bit_negatePtr(u64* rp, u64* xp, usz count); // count is number of u64-s \ No newline at end of file +void bit_negatePtr(u64* rp, u64* xp, usz count); // count is number of u64-s