#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_AVX2 #define SINGELI_FILE scan #include "../utils/includeSingeli.h" #if __PCLMUL__ #define SINGELI_FILE neq #include "../utils/includeSingeli.h" #endif #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_AVX2 && __PCLMUL__ clmul_scan_ne(p, xp, rp, BIT_N(ia)); #if USE_VALGRIND if (ia&63) rp[ia>>6] = vg_def_u64(rp[ia>>6]); #endif #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); SHOULD_INLINE B scan2_min_num(B w, B x, u8 xe, usz ia) { MINMAX2(min,<,MAX,and,1,dsc) } SHOULD_INLINE 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 static B scan_lt(B x, u64 p, usz ia) { u64* xp = bitarr_ptr(x); u64* rp; B r=m_bitarrv(&rp,ia); usz n=BIT_N(ia); u64 m10 = 0x5555555555555555; for (usz i=0; i>63); rp[i] = p = x & (m10 ^ (x + c)); } decG(x); return r; } static B scan_plus(f64 r0, B x, u8 xe, usz ia) { assert(xe!=el_bit && elNum(xe)); B r; void* rp = m_tyarrv(&r, xe==el_f64? sizeof(f64) : sizeof(i32), ia, xe==el_f64? t_f64arr : t_i32arr); #if SINGELI_AVX2 switch(xe) { default:UD; case el_i8: { if (!q_fi32(r0) || simd_scan_plus_i8_i32 (i8any_ptr(x), r0, rp, ia)!=ia) goto cs_i8_f64; decG(x); return r; } case el_i16: { if (!q_fi32(r0) || simd_scan_plus_i16_i32(i16any_ptr(x), r0, rp, ia)!=ia) goto cs_i16_f64; decG(x); return r; } case el_i32: { if (!q_fi32(r0) || simd_scan_plus_i32_i32(i32any_ptr(x), r0, rp, ia)!=ia) goto cs_i32_f64; decG(x); return r; } case el_f64: { f64* xp=f64any_ptr(x); f64 c=r0; for (usz i=0; if; if (isAtm(x) || RNK(x)==0) thrM("`: Argument cannot have rank 0"); ur xr = RNK(x); usz ia = IA(x); if (*SH(x)<=1 || ia==0) return x; if (RARE(!isFun(f))) { if (isMd(f)) thrM("Calling a modifier"); B xf = getFillR(x); MAKE_MUT(rm, ia); usz csz = arr_csz(x); mut_copy(rm, 0, x, 0, csz); mut_fill(rm, csz, f, ia-csz); return withFill(mut_fcd(rm, x), xf); } u8 xe = TI(x,elType); if (v(f)->flags) { u8 rtid = v(f)->flags-1; if (rtid==n_rtack) return x; if (rtid==n_ltack) { usz csz = arr_csz(x); B s = C1(fne, incG(x)); Arr* r = TI(x,slice)(x, 0, csz); return C2(shape, s, taga(r)); } if (!(xr==1 && xe<=el_f64)) goto base; if (xe==el_bit) { 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) return scan_lt(x, 0, ia); // < goto base; } if (rtid==n_add) return scan_plus(0, x, xe, ia); // + if (rtid==n_floor) return scan_min_num(x, xe, ia); // โŒŠ if (rtid==n_ceil ) return scan_max_num(x, xe, ia); // โŒˆ if (rtid==n_ne) { // โ‰  if (!elInt(xe)) goto base; f64 x0 = IGetU(x,0).f; if (!q_fbit(x0)) goto base; u64* rp; B r = m_bitarrv(&rp, ia); bool c = x0; rp[0] = c; if (xe==el_i8 ) { i8* xp=i8any_ptr (x); 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); if (isArr(w)? !ptr_eqShape(SH(w), RNK(w), xsh+1, xr-1) : xr!=1) thrF("`: Shape of ๐•จ must match the cell of ๐•ฉ (%H โ‰ก โ‰ข๐•จ, %H โ‰ก โ‰ข๐•ฉ)", w, x); if (ia==0) { dec(w); return x; } if (RARE(!isFun(f))) { if (isMd(f)) thrM("Calling a modifier"); B xf = getFillR(x); MAKE_MUT(rm, ia); mut_fill(rm, 0, f, ia); return withFill(mut_fcd(rm, x), xf); } u8 xe = TI(x,elType); if (v(f)->flags) { u8 rtid = v(f)->flags-1; if (rtid==n_rtack) { dec(w); return x; } if (rtid==n_ltack) return C2(shape, C1(fne, x), w); if (!(xr==1 && elNum(xe) && xe<=el_f64)) goto base; if (rtid==n_floor && isF64(w)) return scan2_min_num(w, x, xe, ia); // โŒŠ if (rtid==n_ceil && isF64(w)) return scan2_max_num(w, x, xe, ia); // โŒˆ if (rtid==n_add) { // + if (xe==el_bit) { if (!q_i64(w)) goto base; i64 wv = o2i64G(w); if (wv<=(-(1LL<<53)) || wv>=(1LL<<53) || wv+(i64)ia >= (1LL<<53)) goto base; B t = scan_add_bool(x, ia); return wv==0? t : C2(add, w, t); } if (isF64(w) && elInt(xe)) return scan_plus(o2fG(w), x, xe, ia); } if (rtid==n_ne) { // โ‰  bool wBit = q_bit(w); if (xe==el_bit) return scan_ne(x, -(u64)(wBit? o2bG(w) : 1&~*bitarr_ptr(x)), ia); if (!wBit || !elInt(xe)) goto base; bool c = o2bG(w); u64* rp; B r = m_bitarrv(&rp, ia); if (xe==el_i8 ) { i8* xp=i8any_ptr (x); for (usz i=0; i