diff --git a/src/builtins/sysfn.c b/src/builtins/sysfn.c index b94f292d..59ba601a 100644 --- a/src/builtins/sysfn.c +++ b/src/builtins/sysfn.c @@ -260,34 +260,56 @@ B rand_deal_c2(B t, B w, B x) { if (RARE(xi<0)) thrM("(rand).Deal: 𝕩 cannot be negative"); if (RARE(wi>xi)) thrM("(rand).Deal: 𝕨 cannot exceed 𝕩"); if (wi==0) return emptyIVec(); + B r; RAND_START; - TALLOC(i32,s,xi); - for (i64 i = 0; i < xi; i++) s[i] = i; - for (i64 i = 0; i < wi; i++) { - i32 j = wy2u0k(wyrand(&seed), xi-i) + i; - i32 c = s[j]; - s[j] = s[i]; - s[i] = c; + if (wi > xi/64) { + // Dense shuffle + TALLOC(i32,s,xi); + for (i64 i = 0; i < xi; i++) s[i] = i; + for (i64 i = 0; i < wi; i++) { + i32 j = wy2u0k(wyrand(&seed), xi-i) + i; + i32 c = s[j]; + s[j] = s[i]; + s[i] = c; + } + i32* rp; r = m_i32arrv(&rp, wi); + memcpy(rp, s, wi*4); + TFREE(s); + } else { + // Hash-based shuffle + i32* rp; r = m_i32arrv(&rp, wi); + i64 sz = 1; + while (sz < wi*2) sz*= 2; + TALLOC(i32, hash, 2*sz); i32* val = hash+1; + for (u64 i = 0; i < 2*sz; i++) hash[i] = 0; + for (i64 i = 0; i < wi; i++) rp[i] = i; + u64 mask = 2*(sz-1); + for (i64 i = 0; i < wi; i++) { + u64 j = wy2u0k(wyrand(&seed), xi-i) + i; + if (jxi)) thrM("(rand).Subset: 𝕨 cannot exceed 𝕩"); if (wi==0) return emptyIVec(); + B r; + if (wi==xi) { // Only one complete subset; will hang without this + i32* rp; r = m_i32arrv(&rp, wi); + for (u64 i = 0; i < wi; i++) rp[i] = i; + return r; + } + bool invert = wi > xi/2; + i32 wn = invert ? xi-wi : wi; RAND_START; - i32* rp; B r = m_i32arrv(&rp, wi); - i64 sz = 1; - while (sz < wi*2) sz*= 2; - sz*= 2; - H_i2i* map = m_i2i(sz); - for (i64 i = 0; i < wi; i++) rp[i] = i; - for (i64 i = 0; i < wi; i++) { - i32 j = wy2u0k(wyrand(&seed), xi-i) + i; - if (j xi/8) { + // Bit set (as bytes) + TALLOC(u8, set, xi); + for (u64 i = 0; i < xi; i++) set[i] = 0; + for (i32 i = xi-wn; i < xi; i++) { + i32 j = wy2u0k(wyrand(&seed), i+1); + if (set[j]) j=i; + set[j] = 1; + } + i32* rp; r = m_i32arrv(&rp, wi); + if (!invert) { for (u64 i = 0; i < xi; i++) if ( set[i]) *rp++=i; } + else { for (u64 i = 0; i < xi; i++) if (!set[i]) *rp++=i; } + TFREE(set); + } else { + // Sorted "hash" set + u64 sh = 0; + for (u64 xt=xi/4; xt>=wn; xt>>=1) sh++; + u64 sz = ((xi-1)>>sh)+1 + wn; + TALLOC(i32, hash, sz); + for (u64 i = 0; i < sz; i++) hash[i] = xi; + for (i32 i = xi-wn; i < xi; i++) { + i32 j = wy2u0k(wyrand(&seed), i+1); + u64 p = (u64)j >> sh; + while (true) { + i32 h = hash[p]; + if (LIKELY(j>sh; continue; } + p++; + } + } + i32* rp; r = m_i32arrv(&rp, wi); + if (!invert) { + for (u64 i = 0; i < sz; i++) if (hash[i]!=xi) *rp++=hash[i]; } else { - rp[i] = swap_i2i(&map, j, rp[i], j); + i32 e = 0; + for (u64 i = 0; i < sz; i++) { + i32 f = hash[i]; + if (f!=xi) { while (e