From 6ed3c183896d42e1ed953b2e98d79a6c5eb0cf78 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Tue, 20 Sep 2022 20:24:13 -0400 Subject: [PATCH 01/11] Implement SIMD wrapping plus-scan for Replicate --- src/builtins/slash.c | 20 ++++++++--- src/singeli/src/scan.singeli | 69 +++++++++++++++++++++++++----------- 2 files changed, 64 insertions(+), 25 deletions(-) diff --git a/src/builtins/slash.c b/src/builtins/slash.c index 44a5cd8b..dccc80c6 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -104,6 +104,16 @@ #endif #endif +#if SINGELI +extern void (*const avx2_scan_pluswrap_u8)(uint8_t* v0,uint8_t* v1,uint64_t v2,uint8_t v3); +extern void (*const avx2_scan_pluswrap_u16)(uint16_t* v0,uint16_t* v1,uint64_t v2,uint16_t v3); +extern void (*const avx2_scan_pluswrap_u32)(uint32_t* v0,uint32_t* v1,uint64_t v2,uint32_t v3); +#define avx2_scan_pluswrap_u64(V0,V1,V2,V3) for (usz i=k; i Date: Tue, 20 Sep 2022 20:33:48 -0400 Subject: [PATCH 02/11] Constant replicate using plus-scan --- src/builtins/slash.c | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/builtins/slash.c b/src/builtins/slash.c index dccc80c6..86ea9ea5 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -654,18 +654,31 @@ B slash_c2(B t, B w, B x) { if (isArr(x) && RNK(x)==1 && q_i32(w)) { usz xia = IA(x); i32 wv = o2i(w); - if (wv<=0) { + if (wv<=1) { if (wv<0) thrM("/: š•Ø cannot be negative"); - return taga(arr_shVec(TI(x,slice)(x, 0, 0))); + return wv ? x : taga(arr_shVec(TI(x,slice)(x, 0, 0))); } - if (TI(x,elType)==el_i32) { - i32* xp = i32any_ptr(x); - i32* rp; r = m_i32arrv(&rp, xia*wv); - for (usz i = 0; i < xia; i++) { - for (i64 j = 0; j < wv; j++) *rp++ = xp[i]; - } + u8 xe = TI(x,elType); + #define CONST_REP(T) \ + usz s = xia*wv; \ + T* xp = tyany_ptr(x); \ + T* rp = m_tyarrv(&r, elWidth(xe), s, el2t(xe)); \ + usz b = 1<<10; \ + T js=xp[0], px=js; \ + for (usz k=0, j=0, ij=wv; ; ) { \ + usz e = b Date: Wed, 21 Sep 2022 17:29:10 -0400 Subject: [PATCH 03/11] Use Singeli max/min scan for Indices and Sort --- src/builtins/grade.h | 58 +++++++++++++++++++++++++++++++++----------- src/builtins/slash.c | 21 +++++++++++++--- 2 files changed, 62 insertions(+), 17 deletions(-) diff --git a/src/builtins/grade.h b/src/builtins/grade.h index 5d7d634c..ee2669e4 100644 --- a/src/builtins/grade.h +++ b/src/builtins/grade.h @@ -30,21 +30,47 @@ rp[j] = xi; \ } +#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); +#define COUNT_THRESHOLD 32 +#define WRITE_SPARSE_i8 \ + for (usz i=0; i Date: Thu, 22 Sep 2022 15:40:01 -0400 Subject: [PATCH 04/11] Switch to forward inclusive sums for all radix sorting, and use Singeli Removes SWAR for 8-bit counts, since reverse sorting uses negatives --- src/builtins/grade.h | 91 +++++++++++++++++++++++++++----------------- 1 file changed, 56 insertions(+), 35 deletions(-) diff --git a/src/builtins/grade.h b/src/builtins/grade.h index ee2669e4..fd2ebeed 100644 --- a/src/builtins/grade.h +++ b/src/builtins/grade.h @@ -74,33 +74,28 @@ extern void (*const avx2_scan_min16)(int16_t* v0,int16_t* v1,uint64_t v2); TFREE(c0) // Radix sorting -#define PRE(T,K) usz p##K=s##K; s##K+=c##K[j]; c##K[j]=p##K -// 8-bit prefix sum by SWAR -#define PRE_UD(K,SL,SRE) \ - u64 p##K=s##K; s##K+=((u64*)c##K)[j]; \ - s##K+=s##K SL 8; s##K+=s##K SL 16; s##K+=s##K SL 32; \ - ((u64*)c##K)[j] = p##K|(s##K SL 8); s##K SRE 56 -#define PRE64(K) GRADE_UD(PRE_UD(K,<<,>>=), PRE_UD(K,>>,<<=)) +#define INC(P,I) GRADE_UD((P+1)[I]++,P[I]--) +#define ROFF GRADE_UD(1,0) // Radix offset #define CHOOSE_SG_SORT(S,G) S #define CHOOSE_SG_GRADE(S,G) G #define RADIX_SORT_i8(T, TYP) \ - TALLOC(T, c0, 256); T *c0o=c0+128; \ + TALLOC(T, c0, 256+ROFF); T* c0o=c0+128; \ for (usz j=0; j<256; j++) c0[j]=0; \ - for (usz i=0; i>8)]++; } \ + c1[0]=GRADE_UD(-n,c0[0]=n); \ + for (usz i=0; i>8)); } \ RADIX_SUM_2_##T; \ i16 *r0 = (i16*)(c0+2*256); \ CHOOSE_SG_##TYP( \ @@ -112,18 +107,15 @@ extern void (*const avx2_scan_min16)(int16_t* v0,int16_t* v1,uint64_t v2); for (usz i=0; i>8)]++]=g0[i]; } \ ) \ TFREE(alloc) -#define RADIX_SUM_2_u8 u64 s0=0, s1=0; FOR(j,256/8) { PRE64(0); PRE64(1); } -#define RADIX_SUM_2(T) T s0=0, s1=0; FOR(j,256) { PRE(T,0); PRE(T,1); } -#define RADIX_SUM_2_usz RADIX_SUM_2(usz) -#define RADIX_SUM_2_u32 RADIX_SUM_2(u32) #define RADIX_SORT_i32(T, TYP, I) \ - TALLOC(u8, alloc, 4*256*sizeof(T) + n*(4 + CHOOSE_SG_##TYP(0,4+sizeof(I)))); \ + TALLOC(u8, alloc, (4*256+ROFF)*sizeof(T) + n*(4 + CHOOSE_SG_##TYP(0,4+sizeof(I)))); \ T *c0=(T*)alloc, *c1=c0+256, *c2=c1+256, *c3=c2+256, *c3o=c3+128; \ for (usz j=0; j<4*256; j++) c0[j]=0; \ + c1[0]=c2[0]=c3[0]=GRADE_UD(-n,c0[0]=n); \ for (usz i=0; i> 8)]++; \ - c2 [(u8)(v>>16)]++; c3o[(i8)(v>>24)]++; } \ + INC(c0 ,(u8)v ); INC(c1 ,(u8)(v>> 8)); \ + INC(c2 ,(u8)(v>>16)); INC(c3o,(i8)(v>>24)); } \ RADIX_SUM_4_##T; \ i32 *r0 = (i32*)(c0+4*256); \ CHOOSE_SG_##TYP( \ @@ -139,10 +131,37 @@ extern void (*const avx2_scan_min16)(int16_t* v0,int16_t* v1,uint64_t v2); for (usz i=0; i>24)]++; rp[c]=g0[i]; } \ ) \ TFREE(alloc) -#define RADIX_SUM_4_u8 u64 s0=0, s1=0, s2=0, s3=0; FOR(j, 256/8) { PRE64(0); PRE64(1); PRE64(2); PRE64(3); } -#define RADIX_SUM_4(T) T s0=0, s1=0, s2=0, s3=0; FOR(j, 256) { PRE(u32,0); PRE(u32,1); PRE(u32,2); PRE(u32,3); } -#define RADIX_SUM_4_usz RADIX_SUM_4(usz) -#define RADIX_SUM_4_u32 RADIX_SUM_4(u32) + +#define PRE(K) s##K=c##K[j]+=s##K +#define RADIX_SUM_1(T) T s0=0; for(usz j=0;j<256;j++) { PRE(0); } +#define RADIX_SUM_2(T) GRADE_UD(c1[0]=0;,) T s0=0, s1=0; for(usz j=0;j<256;j++) { PRE(0); PRE(1); } +#define RADIX_SUM_4(T) GRADE_UD(c1[0]=c2[0]=c3[0]=0;,) T s0=0, s1=0, s2=0, s3=0; for(usz j=0;j<256;j++) { PRE(0); PRE(1); PRE(2); PRE(3); } + +#if SINGELI +extern void (*const avx2_scan_pluswrap_u8)(uint8_t* v0,uint8_t* v1,uint64_t v2,uint8_t v3); +extern void (*const avx2_scan_pluswrap_u32)(uint32_t* v0,uint32_t* v1,uint64_t v2,uint32_t v3); +#define RADIX_SUM_1_u8 avx2_scan_pluswrap_u8 (c0,c0, 256,0); +#define RADIX_SUM_2_u8 avx2_scan_pluswrap_u8 (c0,c0,2*256,0); +#define RADIX_SUM_2_u32 avx2_scan_pluswrap_u32(c0,c0,2*256,0); +#define RADIX_SUM_4_u8 avx2_scan_pluswrap_u8 (c0,c0,4*256,0); +#define RADIX_SUM_4_u32 avx2_scan_pluswrap_u32(c0,c0,4*256,0); +#else +#define RADIX_SUM_1_u8 RADIX_SUM_1(u8) +#define RADIX_SUM_2_u8 RADIX_SUM_2(u8) +#define RADIX_SUM_2_u32 RADIX_SUM_2(u32) +#define RADIX_SUM_4_u8 RADIX_SUM_4(u8) +#define RADIX_SUM_4_u32 RADIX_SUM_4(u32) +#endif + +#if SINGELI && !USZ_64 +#define RADIX_SUM_1_usz avx2_scan_pluswrap_u32(c0,c0, 256,0); +#define RADIX_SUM_2_usz avx2_scan_pluswrap_u32(c0,c0,2*256,0); +#define RADIX_SUM_4_usz avx2_scan_pluswrap_u32(c0,c0,4*256,0); +#else +#define RADIX_SUM_1_usz RADIX_SUM_1(usz) +#define RADIX_SUM_2_usz RADIX_SUM_2(usz) +#define RADIX_SUM_4_usz RADIX_SUM_4(usz) +#endif #define SORT_C1 CAT(GRADE_UD(and,or),c1) B SORT_C1(B t, B x) { @@ -166,9 +185,9 @@ B SORT_C1(B t, B x) { } else if (xe==el_i8) { i8* xp = i8any_ptr(x); i8* rp; r = m_i8arrv(&rp, n); - if (n<16) { + if (n < 16) { INSERTION_SORT(i8); - } else if (n<256) { + } else if (n < 256) { RADIX_SORT_i8(u8, SORT); } else { COUNTING_SORT(i8); @@ -176,7 +195,7 @@ B SORT_C1(B t, B x) { } else if (xe==el_i16) { i16* xp = i16any_ptr(x); i16* rp; r = m_i16arrv(&rp, n); - if (n < 24) { + if (n < 20) { INSERTION_SORT(i16); } else if (n < 256) { RADIX_SORT_i16(u8, SORT,); @@ -188,7 +207,7 @@ B SORT_C1(B t, B x) { } else if (xe==el_i32) { i32* xp = i32any_ptr(x); i32* rp; r = m_i32arrv(&rp, n); - if (n < 40) { + if (n < 32) { INSERTION_SORT(i32); } else if (n < 256) { RADIX_SORT_i32(u8, SORT,); @@ -386,21 +405,23 @@ done: #undef LT #undef FOR #undef PRE -#undef PRE_UD +#undef INC +#undef ROFF #undef PRE64 #undef CHOOSE_SG_SORT #undef CHOOSE_SG_GRADE #undef RADIX_SORT_i8 +#undef RADIX_SORT_i16 +#undef RADIX_SORT_i32 +#undef RADIX_SUM_1 +#undef RADIX_SUM_2 +#undef RADIX_SUM_4 #undef RADIX_SUM_1_u8 #undef RADIX_SUM_1_usz -#undef RADIX_SORT_i16 #undef RADIX_SUM_2_u8 -#undef RADIX_SUM_2 #undef RADIX_SUM_2_usz #undef RADIX_SUM_2_u32 -#undef RADIX_SORT_i32 #undef RADIX_SUM_4_u8 -#undef RADIX_SUM_4 #undef RADIX_SUM_4_usz #undef RADIX_SUM_4_u32 #undef GRADE_CAT From d647978c8f55b5b7ffad061bd348e4f68a246a72 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 22 Sep 2022 21:42:18 -0400 Subject: [PATCH 05/11] Type-independent Replicate by constant --- src/builtins/slash.c | 101 +++++++++++++++++++++++++------------------ 1 file changed, 60 insertions(+), 41 deletions(-) diff --git a/src/builtins/slash.c b/src/builtins/slash.c index 584dcd6e..d0b72109 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -527,18 +527,30 @@ B slash_c1(B t, B x) { } B slash_c2(B t, B w, B x) { - B r; - if (isArr(w) && RNK(w)==1 && depth(w)==1) { - usz wia = IA(w); + i32 wv = -1; + usz wia; + if (isArr(w)) { + if (depth(w)>1) goto base; + ur wr = RNK(w); + if (wr>1) thrF("/: Simple š•Ø must have rank 0 or 1 (%i≔=š•Ø)", wr); + if (wr<1) { B v=IGet(w, 0); decG(w); w=v; goto atom; } + wia = IA(w); if (wia==0) { decG(w); return isArr(x)? x : m_atomUnit(x); } - if (isAtm(x) || RNK(x)==0) thrM("/: š•© must have rank at least 1 for simple š•Ø"); - ur xr = RNK(x); - usz xlen = *SH(x); + } else { + atom: + if (!q_i32(w)) goto base; + wv = o2i(w); + } + if (isAtm(x) || RNK(x)==0) thrM("/: š•© must have rank at least 1 for simple š•Ø"); + ur xr = RNK(x); + usz xlen = *SH(x); + u8 xl = cellWidthLog(x); + u8 xt = arrNewType(TY(x)); + + B r; + if (wv < 0) { // Array w if (RARE(wia!=xlen)) thrF("/: Lengths of components of š•Ø must match š•© (%s ≠ %s)", wia, xlen); - u8 xl = cellWidthLog(x); - u8 xt = arrNewType(TY(x)); - u8 we = TI(w,elType); if (!elInt(we)) { w=any_squeeze(w); we=TI(w,elType); @@ -580,9 +592,9 @@ B slash_c2(B t, B w, B x) { // Make shape if needed; all cases below use it usz* rsh = NULL; if (xr > 1) { - usz* sh = rsh = m_shArr(xr)->a; - sh[0] = s; - shcpy(sh+1, SH(x)+1, xr-1); + rsh = m_shArr(xr)->a; + rsh[0] = s; + shcpy(rsh+1, SH(x)+1, xr-1); } if (xl == 0) { @@ -665,19 +677,37 @@ B slash_c2(B t, B w, B x) { } } goto decWX_ret; - } - if (isArr(x) && RNK(x)==1 && q_i32(w)) { - usz xia = IA(x); - i32 wv = o2i(w); - if (wv<=1) { - if (wv<0) thrM("/: š•Ø cannot be negative"); + } else { + if (wv <= 1) { + if (wv < 0) thrM("/: š•Ø cannot be negative"); return wv ? x : taga(arr_shVec(TI(x,slice)(x, 0, 0))); } - u8 xe = TI(x,elType); - #define CONST_REP(T) \ - usz s = xia*wv; \ - T* xp = tyany_ptr(x); \ - T* rp = m_tyarrv(&r, elWidth(xe), s, el2t(xe)); \ + if (xlen == 0) return x; + usz s = xlen * wv; + if (xl>6 || xl<3 || TI(x,elType)==el_B) { + if (xr != 1) goto base; + SLOW2("š•Ø/š•©", w, x); + B xf = getFillQ(x); + HArr_p r0 = m_harrUv(s); + SGetU(x) + for (usz i = 0; i < xlen; i++) { + B cx = incBy(GetU(x, i), wv); + for (i64 j = 0; j < wv; j++) *r0.a++ = cx; + } + r = withFill(r0.b, xf); + goto decX_ret; + } + u8 xk = xl-3; + void* rv = m_tyarrv(&r, 1< 1) { + usz* rsh = m_shArr(xr)->a; + rsh[0] = s; + shcpy(rsh+1, SH(x)+1, xr-1); + Arr* ra=a(r); SPRNK(ra,xr); PSH(ra)=rsh; PIA(ra)=s*arr_csz(x); + } + void* xv = tyany_ptr(x); + #define CONST_REP(T) { \ + T* xp = xv; T* rp = rv; \ usz b = 1<<10; \ T js=xp[0], px=js; \ for (usz k=0, j=0, ij=wv; ; ) { \ @@ -686,25 +716,14 @@ B slash_c2(B t, B w, B x) { while (ij Date: Thu, 22 Sep 2022 22:09:56 -0400 Subject: [PATCH 06/11] Merge scan-based sparse Indices/Replicate code with macros --- src/builtins/slash.c | 82 ++++++++++++++++---------------------------- 1 file changed, 29 insertions(+), 53 deletions(-) diff --git a/src/builtins/slash.c b/src/builtins/slash.c index d0b72109..cb431632 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -448,6 +448,30 @@ static B compress(B w, B x, usz wia, u8 xl, u8 xt) { return r; } +#define SCAN_CORE(WV, UPD, SET, SCAN) \ + usz b = 1<<10; \ + for (usz k=0, j=0, ij=WV; ; ) { \ + usz e = b Date: Thu, 22 Sep 2022 22:32:52 -0400 Subject: [PATCH 07/11] Constant boolean replicate --- src/builtins/slash.c | 85 ++++++++++++++++++++++++++------------------ 1 file changed, 50 insertions(+), 35 deletions(-) diff --git a/src/builtins/slash.c b/src/builtins/slash.c index cb431632..bdb2d91e 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -448,6 +448,7 @@ static B compress(B w, B x, usz wia, u8 xl, u8 xt) { return r; } +// Replicate using plus/max/xor-scan #define SCAN_CORE(WV, UPD, SET, SCAN) \ usz b = 1<<10; \ for (usz k=0, j=0, ij=WV; ; ) { \ @@ -472,6 +473,37 @@ static B compress(B w, B x, usz wia, u8 xl, u8 xt) { T js=xp[0], px=js; \ SUM_CORE(T, WV, T sx=px, (px=xp[j])-sx) +#define BOOL_REP_XOR_SCAN(WV) \ + usz b = 1<<12; \ + u64 xx=xp[0], xs=xx>>63, js=-(xx&1); xx^=xx<<1; \ + for (usz k=0, j=0, ij=WV; ; ) { \ + usz e = b>=1; j++; if (j%64==0) { u64 v=xp[j/64]; xx=v^(v<<1)^xs; xs=v>>63; } \ + rp[ij/64]^=(-(xx&1))<<(ij%64); ij+=WV; \ + } \ + for (usz i=k/64; i>63); \ + if (e==s) {break;} k=e; \ + } + +// Basic boolean loop with overwriting +#define BOOL_REP_OVER(WV, LEN) \ + u64 ri=0, rc=0, xc=0; usz j=0; \ + for (usz i = 0; i < LEN; i++) { \ + u64 v = -(u64)bitp_get(xp,i); \ + rc ^= (v^xc) << (ri%64); \ + xc = v; \ + ri += WV; usz e = ri/64; \ + if (j < e) { \ + rp[j++] = rc; \ + while (j < e) rp[j++] = v; \ + rc = v; \ + } \ + } \ + if (ri%64) rp[j] = rc; + extern B rt_slash; B slash_c1(B t, B x) { if (RARE(isAtm(x)) || RARE(RNK(x)!=1)) thrF("/: Argument must have rank 1 (%H ≔ ā‰¢š•©)", x); @@ -602,21 +634,7 @@ B slash_c2(B t, B w, B x) { u64* xp = bitarr_ptr(x); u64* rp; r = m_bitarrv(&rp, s); if (rsh) { SPRNK(a(r),xr); SH(r) = rsh; } if (s/256 <= wia) { - #define SPARSE_REP(T) \ - T* wp = T##any_ptr(w); \ - usz b = 1<<12; \ - u64 xx=xp[0], xs=xx>>63, js=-(xx&1); xx^=xx<<1; \ - for (usz k=0, j=0, ij=wp[0]; ; ) { \ - usz e = b>=1; j++; if (j%64==0) { u64 v=xp[j/64]; xx=v^(v<<1)^xs; xs=v>>63; } \ - rp[ij/64]^=(-(xx&1))<<(ij%64); ij+=wp[j]; \ - } \ - for (usz i=k/64; i>63); \ - if (e==s) {break;} k=e; \ - } + #define SPARSE_REP(T) T* wp=T##any_ptr(w); BOOL_REP_XOR_SCAN(wp[j]) if (we==el_i8 ) { SPARSE_REP(i8 ); } else if (we==el_i16) { SPARSE_REP(i16); } else { SPARSE_REP(i32); } @@ -624,19 +642,7 @@ B slash_c2(B t, B w, B x) { } else { if (we < el_i32) w = taga(cpyI32Arr(w)); i32* wp = i32any_ptr(w); - u64 ri=0, rc=0, xc=0; usz j=0; - for (usz i = 0; i < wia; i++) { - u64 v = -(u64)bitp_get(xp,i); - rc ^= (v^xc) << (ri%64); - xc = v; - ri += wp[i]; usz e = ri/64; - if (j < e) { - rp[j++] = rc; - while (j < e) rp[j++] = v; - rc = v; - } - } - if (ri%64) rp[j] = rc; + BOOL_REP_OVER(wp[i], wia) } } else { u8 xk = xl-3; @@ -675,7 +681,7 @@ B slash_c2(B t, B w, B x) { } if (xlen == 0) return x; usz s = xlen * wv; - if (xl>6 || xl<3 || TI(x,elType)==el_B) { + if (xl>6 || (xl<3 && xl!=0) || TI(x,elType)==el_B) { if (xr != 1) goto base; SLOW2("š•Ø/š•©", w, x); B xf = getFillQ(x); @@ -688,18 +694,27 @@ B slash_c2(B t, B w, B x) { r = withFill(r0.b, xf); goto decX_ret; } - u8 xk = xl-3; - void* rv = m_tyarrv(&r, 1< 1) { usz* rsh = m_shArr(xr)->a; rsh[0] = s; shcpy(rsh+1, SH(x)+1, xr-1); Arr* ra=a(r); SPRNK(ra,xr); PSH(ra)=rsh; PIA(ra)=s*arr_csz(x); } - void* xv = tyany_ptr(x); - #define CASE(L,T) case L: { REP_BY_SCAN(T, wv) goto decX_ret; } - switch (xk) { default: UD; CASE(0,u8) CASE(1,u16) CASE(2,u32) CASE(3,u64) } - #undef CASE + goto decX_ret; } base: return c2(rt_slash, w, x); From 8e24da8e8e320c2afa058a650ad01ec9889d2122 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Fri, 23 Sep 2022 08:22:41 -0400 Subject: [PATCH 08/11] Adjust Replicate thresholds --- src/builtins/slash.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/builtins/slash.c b/src/builtins/slash.c index bdb2d91e..1f1c2d83 100644 --- a/src/builtins/slash.c +++ b/src/builtins/slash.c @@ -633,7 +633,7 @@ B slash_c2(B t, B w, B x) { if (xl == 0) { u64* xp = bitarr_ptr(x); u64* rp; r = m_bitarrv(&rp, s); if (rsh) { SPRNK(a(r),xr); SH(r) = rsh; } - if (s/256 <= wia) { + if (s/1024 <= wia) { #define SPARSE_REP(T) T* wp=T##any_ptr(w); BOOL_REP_XOR_SCAN(wp[j]) if (we==el_i8 ) { SPARSE_REP(i8 ); } else if (we==el_i16) { SPARSE_REP(i16); } @@ -649,7 +649,7 @@ B slash_c2(B t, B w, B x) { void* rv = m_tyarrv(&r, 1< Date: Fri, 23 Sep 2022 15:32:27 -0400 Subject: [PATCH 09/11] Rework reshape of singleton --- src/builtins/sfns.c | 87 +++++++++++++++++++++++++-------------------- 1 file changed, 48 insertions(+), 39 deletions(-) diff --git a/src/builtins/sfns.c b/src/builtins/sfns.c index 7d7ab988..b61476ba 100644 --- a/src/builtins/sfns.c +++ b/src/builtins/sfns.c @@ -104,6 +104,13 @@ static B truncReshape(B x, usz xia, usz nia, ur nr, ShArr* sh) { // consumes all arr_shSetU(ra, nr, sh); return r; } +static void fill_words(void* rp, u64 v, u64 bytes) { + usz wds = bytes/8; + usz ext = bytes%8; + u64* p = rp; + for (usz i=0; i Date: Fri, 23 Sep 2022 16:57:41 -0400 Subject: [PATCH 10/11] Fast array reshape --- src/builtins/sfns.c | 83 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 70 insertions(+), 13 deletions(-) diff --git a/src/builtins/sfns.c b/src/builtins/sfns.c index b61476ba..3b0c9ea9 100644 --- a/src/builtins/sfns.c +++ b/src/builtins/sfns.c @@ -188,6 +188,7 @@ B shape_c2(B t, B w, B x) { decG(w); } + Arr* r; if (isArr(x)) { if (nia <= xia) { return truncReshape(x, xia, nia, nr, sh); @@ -201,21 +202,77 @@ B shape_c2(B t, B w, B x) { } if (xia <= nia/2) x = any_squeeze(x); - B xf = getFillQ(x); - MAKE_MUT(m, nia); mut_init(m, TI(x,elType)); - MUTG_INIT(m); - i64 div = nia/xia; - i64 mod = nia%xia; - for (i64 i = 0; i < div; i++) mut_copyG(m, i*xia, x, 0, xia); - mut_copyG(m, div*xia, x, 0, mod); + u8 xl = arrTypeBitsLog(TY(x)); + u8 xt = arrNewType(TY(x)); + u8* rp; + u64 bi, bf; // Bytes present, bytes wanted + if (xl == 0) { // Bits + u64* rq; r = m_bitarrp(&rq, nia); + rp = (u8*)rq; + usz nw = BIT_N(nia); + u64* xp = bitarr_ptr(x); + u64 b = xia; + if (b % 8) { + if (b < 64) { + // Need to avoid calling bit_cpy with arguments <64 bits apart + u64 v = xp[0] & (~(u64)0 >> (64-b)); + do { v |= v<64 && nia>64) rq[1] = v>>(64-b/2); + } else { + memcpy(rq, xp, (b+7)/8); + } + for (; b%8; b*=2) { + if (b>nw*32) { + if (b=bf since bf is rounded up + break; + } + bit_cpy(rq, b, rq, 0, b); + } + } else { + memcpy(rp, xp, b/8); + } + bi = b/8; + bf = 8*nw; + if (bi == 1) { memset(rp, rp[0], bf); bi=bf; } + } else { + if (TI(x,elType) == el_B) { + B xf = getFillQ(x); + MAKE_MUT(m, nia); mut_init(m, el_B); + MUTG_INIT(m); + i64 div = nia/xia; + i64 mod = nia%xia; + for (i64 i = 0; i < div; i++) mut_copyG(m, i*xia, x, 0, xia); + mut_copyG(m, div*xia, x, 0, mod); + decG(x); + Arr* ra = mut_fp(m); + arr_shSetU(ra, nr, sh); + return withFill(taga(ra), xf); + } + u8 xk = xl - 3; + rp = m_tyarrp(&r, 1<> (64-b)); + while (b<64) { v |= v<bf) l=bf; + for (; bi<=l/2; bi+=bi) memcpy(rp+bi, rp, bi); + u64 e=bi; for (; e+bi<=bf; e+=bi) memcpy(rp+e, rp, bi); + if (e Date: Sat, 24 Sep 2022 19:20:32 +0300 Subject: [PATCH 11/11] clean up generic unit reshape case --- src/builtins/sfns.c | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/builtins/sfns.c b/src/builtins/sfns.c index 3b0c9ea9..d9640efc 100644 --- a/src/builtins/sfns.c +++ b/src/builtins/sfns.c @@ -293,15 +293,14 @@ B shape_c2(B t, B w, B x) { if (c==(u8 )c) { u8* rp; r = m_c8arrp(&rp,nia); memset(rp, c, nia); } else if (c==(u16)c) { FILL(c16,u16,c*0x0001000100010001) } else { FILL(c32,u32,c*0x0000000100000001) } - #undef FILL } else { - B xf = asFill(inc(x)); + incBy(x, nia); // in addition with the existing reference, this covers the filled amount & asFill r = m_fillarrp(nia); - if (nia) incBy(x, nia-1); - else dec(x); - fill_words(fillarr_ptr(r), x.u, (u64)nia*8); - fillarr_setFill(r, xf); + if (sizeof(B)==8) fill_words(fillarr_ptr(r), x.u, (u64)nia*8); + else for (usz i = 0; i < nia; i++) fillarr_ptr(r)[i] = x; + fillarr_setFill(r, asFill(x)); } + #undef FILL } arr_shSetU(r,nr,sh); return taga(r);