代码

/*
References:
- Duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/104819: copied public grouped AVX2 CSA bitset implementation and its inline memory-operand XOR; no separate license notice appears in that source.
- Our account saffah_codex_6s_agg2, https://duck.ac/submission/104874: copied our accepted high-plane 16-bit accumulation variant and adapted its result storage.
- Wojciech Mula, Nathan Kurz and Daniel Lemire, "Faster Population Counts Using AVX2 Instructions", https://arxiv.org/abs/1611.07612: carry-save population-count idea, attributed under the paper's CC BY 4.0 terms; no text or software code copied.
Approach:
- Store answers by contiguous grouped query position during bitset processing, then scatter them into original order after all groups. This keeps random result writes away from the hot bitset read phase.
Purpose:
- Experimental official speed and correctness measurement of deferred result scattering on wc2017b2.
*/
// wc2017b2: 2-plane (4 bitset array) encoding + query grouping by relative bit offset.
//
// win(a,b) (a=s1[i], b=s2[i+delta]) holds iff a == ((b-1) mod 3).
// Encode a char v by an injective 2-bit code:
//   s1 side:  e0 = (a>=1), e1 = (a>=2)          (00->0, 10->1, 11->2)
//   s2 side:  f0 = (g!=0), f1 = (g>=2) with g=(b-1)mod3   (00->0, 10->1, 11->2)
// then win  <=>  (e0==f0) && (e1==f1)  <=>  ~((e0^f0)|(e1^f1)).
// Count over a window = l - popcount(diff), 4 arrays instead of 6.
// delta = y-x with |delta| < n; rho = |delta| & 63 selects a pre-rotated copy, so the
// two streams are bit-aligned. Queries are grouped by rho so that only one rotated copy
// is hot at a time (the working set then fits in L2 instead of streaming from L3).
#include <cstring>
#include <immintrin.h>

#define MAXN 300005
#define MAXW (((MAXN + 63) / 64 + 12 + 3) & ~3)   /* multiple of 4 words = 32B rows */

static unsigned long long E0[MAXW + 2] __attribute__((aligned(64))), E1[MAXW + 2] __attribute__((aligned(64))), F0[MAXW + 2] __attribute__((aligned(64))), F1[MAXW + 2] __attribute__((aligned(64)));
// The 256 pre-rotated copies (4 x 64 x 4694 x 8 B = 9388 KB) are replaced by ONE two-array
// scratch: the query loop is group-outer and uses only one rotation per group, so the rotation
// is generated for that group immediately before its inner loop.  9388 KB -> 75 KB.
static unsigned long long R0[MAXW + 2] __attribute__((aligned(64))), R1[MAXW + 2] __attribute__((aligned(64)));
static int ORD[MAXN];
static __m256i FMASK[256] __attribute__((aligned(64)));
static __m256i LMASK[256] __attribute__((aligned(64)));
static unsigned long long RQ[MAXN];
static unsigned OUT[MAXN];
static int GHEAD[520], GCUR[520];
static int g_nw;
static unsigned long long WMASK[65];

#define SOLVE_ARGS int n, int q, char *s1, char *s2, int *q_x, int *q_y, int *q_len, unsigned *ans

#pragma GCC push_options
#pragma GCC target("avx2,popcnt")

// judge's g++-9 with no -march lowers _mm256_loadu_* into vmovdqu+{vinserti128}/{vextracti128}
// (the AVX-256 split-load tax).  This one asm instruction emits the real 32-byte load.
__attribute__((target("avx2"))) static inline __m256i ldu256(const void *p) {
  __m256i r; __asm__("vmovdqu %1, %0" : "=x"(r) : "m"(*(const __m256i *)p)); return r; }

static const signed char LUTA[32] __attribute__((aligned(32))) =
  {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4, 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
static const signed char LOW4A[32] __attribute__((aligned(32))) =
  {15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15};
#define VLUT (*((const __m256i *)LUTA))
#define VLOW4 (*((const __m256i *)LOW4A))

static inline unsigned pc_avx2_diff(const unsigned long long *A0, const unsigned long long *A1,
                                    const unsigned long long *B0, const unsigned long long *B1,
                                    int from, int to) {
  unsigned acc = 0;
  __m256i accv = _mm256_setzero_si256();
  int w = from;
  for (; w + 4 <= to; w += 4) {
    __m256i z = _mm256_or_si256(
        _mm256_xor_si256(ldu256(A0 + w),
                         ldu256(B0 + w)),
        _mm256_xor_si256(ldu256(A1 + w),
                         ldu256(B1 + w)));
    __m256i c8 = _mm256_add_epi8(
        _mm256_shuffle_epi8(VLUT, _mm256_and_si256(z, VLOW4)),
        _mm256_shuffle_epi8(VLUT, _mm256_and_si256(_mm256_srli_epi16(z, 4), VLOW4)));
    accv = _mm256_add_epi64(accv, _mm256_sad_epu8(c8, _mm256_setzero_si256()));
  }
  unsigned long long t[4];
  _mm256_storeu_si256((__m256i *)t, accv);
  acc = (unsigned)(t[0] + t[1] + t[2] + t[3]);
  for (; w < to; w++)
    acc += (unsigned)__builtin_popcountll((A0[w] ^ B0[w]) | (A1[w] ^ B1[w]));
  return acc;
}
#pragma GCC pop_options

#pragma GCC push_options
#pragma GCC target("avx2,popcnt")

__attribute__((target("avx2"),always_inline)) static inline __m256i csa_diff(
    const unsigned long long *a0,const unsigned long long *a1,
    const unsigned long long *b0,const unsigned long long *b1,int k) {
  __m256i r,t;
  __asm__("vmovdqu %2, %0\n\t"
          "vpxor   %3, %0, %0\n\t"
          "vmovdqu %4, %1\n\t"
          "vpxor   %5, %1, %1\n\t"
          "vpor    %1, %0, %0"
          : "=&x"(r),"=&x"(t)
          : "m"(*(const __m256i *)(a0+k)),"m"(*(const __m256i *)(b0+k)),
            "m"(*(const __m256i *)(a1+k)),"m"(*(const __m256i *)(b1+k)));
  return r;
}
__attribute__((target("avx2"),always_inline)) static inline void csa_add3(
    __m256i &hi,__m256i &lo,__m256i a,__m256i b,__m256i c) {
  __m256i u=_mm256_xor_si256(a,b);
  hi=_mm256_or_si256(_mm256_and_si256(a,b),_mm256_and_si256(u,c));
  lo=_mm256_xor_si256(u,c);
}
__attribute__((target("avx2"),always_inline)) static inline __m256i csa_count_bytes(__m256i x) {
  return _mm256_add_epi8(
    _mm256_shuffle_epi8(VLUT,_mm256_and_si256(x,VLOW4)),
    _mm256_shuffle_epi8(VLUT,_mm256_and_si256(_mm256_srli_epi16(x,4),VLOW4)));
}
__attribute__((target("avx2"),always_inline)) static inline __m256i csa_count(__m256i x) {
  return _mm256_sad_epu8(csa_count_bytes(x),_mm256_setzero_si256());
}

static void run_avx2(int n, int q, int *q_x, int *q_y, int *q_len, unsigned *ans) {
  (void)n;
  for (int gi = 0; gi < 512; gi++) {
    int g = (gi >> 2) & 63;
    int ph = gi & 3;                                  /* this sub-pass's phase class */
    if (GHEAD[gi] == GHEAD[gi + 1]) continue;
    const unsigned long long *S0, *S1;
    if ((gi >> 2) < 64) { S0 = F0; S1 = F1; } else { S0 = E0; S1 = E1; }
    int corr = 0;                                     /* B pointer correction (words) */
    if (g) {              /* rotation, built AT THIS PHASE: dst[m] = src_rot[m + ph] */
      unsigned long long sh = 64 - g;
      corr = ph;
      {                                                /* vectorised rotation: 4 words per step */
        __m128i vg = _mm_cvtsi32_si128((long long)g), vsh = _mm_cvtsi32_si128((long long)sh);
        int i = ph, lim = g_nw + 2;
        for (; i + 4 <= lim; i += 4) {
          __m256i a0 = ldu256(S0 + i);
          __m256i b0 = ldu256(S0 + i + 1);
          __m256i a1 = ldu256(S1 + i);
          __m256i b1 = ldu256(S1 + i + 1);
          _mm256_store_si256((__m256i *)(R0 + i - ph), _mm256_or_si256(_mm256_srl_epi64(a0, vg), _mm256_sll_epi64(b0, vsh)));
          _mm256_store_si256((__m256i *)(R1 + i - ph), _mm256_or_si256(_mm256_srl_epi64(a1, vg), _mm256_sll_epi64(b1, vsh)));
        }
        for (; i <= g_nw + 1; i++) {
          R0[i - ph] = (S0[i] >> g) | (S0[i + 1] << sh);
          R1[i - ph] = (S1[i] >> g) | (S1[i + 1] << sh);
        }
      }
      S0 = R0; S1 = R1;
    }
    for (int t = GHEAD[gi]; t < GHEAD[gi + 1]; t++) {
      unsigned long long rq = RQ[t];
      int k = ORD[t];
      int x = (int)(rq & 0x1FFFFFu), y = (int)((rq >> 21) & 0x1FFFFFu), l = (int)(rq >> 42);
      int delta = y - x;
      int dpos = delta >= 0;
      int dw = dpos ? (delta >> 6) : ((-delta) >> 6);
      /* branchless: the query start is min(x,y) and the A side is (E0,E1) when
         y>=x else (F0,F1).  The 50/50 branch here mispredicts on half the queries. */
      int s = dpos ? x : y;
      const unsigned long long *AB = dpos ? E0 : F0, *AC = dpos ? E1 : F1;
      int w0, wl, d, lo, hi;
      const unsigned long long *A0, *A1, *B0, *B1;
      w0 = s >> 6; wl = (s + l - 1) >> 6; d = wl - w0; lo = s & 63; hi = ((s + l - 1) & 63) + 1;
      A0 = AB + w0; A1 = AC + w0; B0 = S0 + w0 + dw - corr; B1 = S1 + w0 + dw - corr;
      int wA = w0 & ~3, wZ = wl & ~3;
      int o = w0 - wA, oL = wl - wZ;
      const unsigned long long *pA0 = A0 - o, *pA1 = A1 - o, *pB0 = B0 - o, *pB1 = B1 - o;
      int bend = wZ - wA;
      __m256i accv = _mm256_setzero_si256();
      if (bend == 0) {
        __m256i z = _mm256_or_si256(
            _mm256_xor_si256(_mm256_load_si256((const __m256i *)(const void *)pA0), ldu256(pB0)),
            _mm256_xor_si256(_mm256_load_si256((const __m256i *)(const void *)pA1), ldu256(pB1)));
        z = _mm256_and_si256(z, _mm256_and_si256(FMASK[o * 64 + lo], LMASK[oL * 64 + hi - 1]));
        __m256i c8 = _mm256_add_epi8(
            _mm256_shuffle_epi8(VLUT, _mm256_and_si256(z, VLOW4)),
            _mm256_shuffle_epi8(VLUT, _mm256_and_si256(_mm256_srli_epi16(z, 4), VLOW4)));
        accv = _mm256_add_epi64(accv, _mm256_sad_epu8(c8, _mm256_setzero_si256()));
      } else {
        {
          __m256i z = _mm256_or_si256(
              _mm256_xor_si256(_mm256_load_si256((const __m256i *)(const void *)pA0), ldu256(pB0)),
              _mm256_xor_si256(_mm256_load_si256((const __m256i *)(const void *)pA1), ldu256(pB1)));
          z = _mm256_and_si256(z, FMASK[o * 64 + lo]);
          __m256i c8 = _mm256_add_epi8(
              _mm256_shuffle_epi8(VLUT, _mm256_and_si256(z, VLOW4)),
              _mm256_shuffle_epi8(VLUT, _mm256_and_si256(_mm256_srli_epi16(z, 4), VLOW4)));
          accv = _mm256_add_epi64(accv, _mm256_sad_epu8(c8, _mm256_setzero_si256()));
        }
        int b = 4;
        __m256i z=_mm256_setzero_si256(), a1=z, a2=z, a4=z, a8=z, t0=z, t1=z, t2=z, t3=z, t4=z, t5=z, t6=z, t7=z, t8=z, t9=z, t10=z, t11=z, t12=z, t13=z, fl=z;
        __m256i high16=z;
        for (; b + 64 <= bend; b += 64) {
#define CD(K) csa_diff(pA0,pA1,pB0,pB1,b+(K))
          csa_add3(t0,a1,a1,CD(0),CD(4));
          csa_add3(t1,a1,a1,CD(8),CD(12));
          csa_add3(t2,a2,a2,t0,t1);
          csa_add3(t3,a1,a1,CD(16),CD(20));
          csa_add3(t4,a1,a1,CD(24),CD(28));
          csa_add3(t5,a2,a2,t3,t4);
          csa_add3(t6,a4,a4,t2,t5);
          csa_add3(t7,a1,a1,CD(32),CD(36));
          csa_add3(t8,a1,a1,CD(40),CD(44));
          csa_add3(t9,a2,a2,t7,t8);
          csa_add3(t10,a1,a1,CD(48),CD(52));
          csa_add3(t11,a1,a1,CD(56),CD(60));
          csa_add3(t12,a2,a2,t10,t11);
          csa_add3(t13,a4,a4,t9,t12);
          csa_add3(fl,a8,a8,t6,t13);
          high16=_mm256_add_epi16(high16,
              _mm256_maddubs_epi16(csa_count_bytes(fl),_mm256_set1_epi8(1)));
#undef CD
        }
        if (b > 4) {
          __m256i hs=_mm256_madd_epi16(high16,_mm256_set1_epi16(1));
          __m256i h64=_mm256_add_epi64(_mm256_and_si256(hs,_mm256_set1_epi64x(0xffffffffULL)),
                                      _mm256_srli_epi64(hs,32));
          accv=_mm256_add_epi64(accv,_mm256_slli_epi64(h64,4));
          accv=_mm256_add_epi64(accv,csa_count(a1));
          accv=_mm256_add_epi64(accv,_mm256_slli_epi64(csa_count(a2),1));
          accv=_mm256_add_epi64(accv,_mm256_slli_epi64(csa_count(a4),2));
          accv=_mm256_add_epi64(accv,_mm256_slli_epi64(csa_count(a8),3));
        }
        for (; b < bend; b += 4) {
          __m256i z = csa_diff(pA0,pA1,pB0,pB1,b);
          __m256i c8 = _mm256_add_epi8(
              _mm256_shuffle_epi8(VLUT, _mm256_and_si256(z, VLOW4)),
              _mm256_shuffle_epi8(VLUT, _mm256_and_si256(_mm256_srli_epi16(z, 4), VLOW4)));
          accv = _mm256_add_epi64(accv, _mm256_sad_epu8(c8, _mm256_setzero_si256()));
        }
        {
          __m256i z = _mm256_or_si256(
              _mm256_xor_si256(_mm256_load_si256((const __m256i *)(const void *)(pA0 + bend)), ldu256(pB0 + bend)),
              _mm256_xor_si256(_mm256_load_si256((const __m256i *)(const void *)(pA1 + bend)), ldu256(pB1 + bend)));
          z = _mm256_and_si256(z, LMASK[oL * 64 + hi - 1]);
          __m256i c8 = _mm256_add_epi8(
              _mm256_shuffle_epi8(VLUT, _mm256_and_si256(z, VLOW4)),
              _mm256_shuffle_epi8(VLUT, _mm256_and_si256(_mm256_srli_epi16(z, 4), VLOW4)));
          accv = _mm256_add_epi64(accv, _mm256_sad_epu8(c8, _mm256_setzero_si256()));
        }
      }
      { unsigned long long sumv[4]; _mm256_storeu_si256((__m256i *)sumv, accv);
        OUT[t] = (unsigned)l - (unsigned)(sumv[0] + sumv[1] + sumv[2] + sumv[3]); }
    }
  }
}

#pragma GCC pop_options

static void run_scalar(int n, int q, int *q_x, int *q_y, int *q_len, unsigned *ans) {
  (void)n;
  for (int gi = 0; gi < 512; gi++) {
    int g = (gi >> 2) & 63;
    int ph = gi & 3;
    if (GHEAD[gi] == GHEAD[gi + 1]) continue;
    const unsigned long long *S0, *S1;
    if ((gi >> 2) < 64) { S0 = F0; S1 = F1; } else { S0 = E0; S1 = E1; }
    int corr = 0;
    if (g) {                                          /* build THIS sub-pass's rotation in place */
      unsigned long long sh = 64 - g;
      corr = ph;
      for (int i = ph; i <= g_nw + 1; i++) {
        R0[i - ph] = (S0[i] >> g) | (S0[i + 1] << sh);
        R1[i - ph] = (S1[i] >> g) | (S1[i + 1] << sh);
      }
      S0 = R0; S1 = R1;
    }
    for (int t = GHEAD[gi]; t < GHEAD[gi + 1]; t++) {
      unsigned long long rq = RQ[t];
      int k = ORD[t];
      int x = (int)(rq & 0x1FFFFFu), y = (int)((rq >> 21) & 0x1FFFFFu), l = (int)(rq >> 42);
      int delta = y - x;
      int dpos = delta >= 0;
      int dw = dpos ? (delta >> 6) : ((-delta) >> 6);
      /* branchless: the query start is min(x,y) and the A side is (E0,E1) when
         y>=x else (F0,F1).  The 50/50 branch here mispredicts on half the queries. */
      int s = dpos ? x : y;
      const unsigned long long *AB = dpos ? E0 : F0, *AC = dpos ? E1 : F1;
      int w0, wl, d, lo, hi;
      const unsigned long long *A0, *A1, *B0, *B1;
      w0 = s >> 6; wl = (s + l - 1) >> 6; d = wl - w0; lo = s & 63; hi = ((s + l - 1) & 63) + 1;
      A0 = AB + w0; A1 = AC + w0; B0 = S0 + w0 + dw - corr; B1 = S1 + w0 + dw - corr;
      unsigned dd = 0;
      for (int w = 0; w <= d; w++) {
        unsigned long long m = (A0[w] ^ B0[w]) | (A1[w] ^ B1[w]);
        if (w == 0) m &= ~WMASK[lo];
        if (w == d) m &= WMASK[hi];
        dd += (unsigned)__builtin_popcountll(m);
      }
      ans[k] = (unsigned)l - dd;
    }
  }
}

__attribute__((target("avx2,popcnt"))) static void build_all(int n, const char *s1, const char *s2) {
  int nw = (n + 63) >> 6;
  for (int i = 0; i <= nw + 1; i++) { E0[i] = E1[i] = F0[i] = F1[i] = 0; }
  {
    /* Vectorised plane build.  E0 = (a != '0'), E1 = (a == '2'),
       F0 = (c != '1'), F1 = (c == '0')   -- algebraically identical to the
       branchy version, but with no data-dependent branches at all
       (the scalar form mispredicts ~4x per 3 characters on random input). */
    const __m256i Z0 = _mm256_set1_epi8('0'), Z1 = _mm256_set1_epi8('1'), Z2 = _mm256_set1_epi8('2');
    int i = 0;
    for (; i + 32 <= n; i += 32) {
      int w = i >> 6, b = i & 63;
      __m256i v1 = _mm256_loadu_si256((const __m256i *)(s1 + i));
      __m256i v2 = _mm256_loadu_si256((const __m256i *)(s2 + i));
      unsigned e0 = ~(unsigned)_mm256_movemask_epi8(_mm256_cmpeq_epi8(v1, Z0));
      unsigned e1 = (unsigned)_mm256_movemask_epi8(_mm256_cmpeq_epi8(v1, Z2));
      unsigned f0 = ~(unsigned)_mm256_movemask_epi8(_mm256_cmpeq_epi8(v2, Z1));
      unsigned f1 = (unsigned)_mm256_movemask_epi8(_mm256_cmpeq_epi8(v2, Z0));
      if (b == 0) { E0[w] = e0; E1[w] = e1; F0[w] = f0; F1[w] = f1; }
      else { E0[w] |= (unsigned long long)e0 << 32; E1[w] |= (unsigned long long)e1 << 32;
             F0[w] |= (unsigned long long)f0 << 32; F1[w] |= (unsigned long long)f1 << 32; }
    }
    for (; i < n; i++) {
      int w = i >> 6, b = i & 63;
      int a = s1[i] - '0', c = s2[i] - '0';
      if (a >= 1) E0[w] |= 1ull << b;
      if (a >= 2) E1[w] |= 1ull << b;
      int g = (c + 2) % 3;
      if (g != 0) F0[w] |= 1ull << b;
      if (g >= 2) F1[w] |= 1ull << b;
    }
  }
  g_nw = nw;
  for (int i = 0; i < 64; i++) WMASK[i] = (1ull << i) - 1ull;
  WMASK[64] = ~0ull;
  for (int o = 0; o < 4; o++)
    for (int lo = 0; lo < 64; lo++) {
      unsigned long long *p = (unsigned long long *)&FMASK[o * 64 + lo];
      for (int j = 0; j < 4; j++) p[j] = (j < o) ? 0ull : (j == o ? ~WMASK[lo] : ~0ull);
    }
  for (int o = 0; o < 4; o++)
    for (int hh = 1; hh <= 64; hh++) {
      unsigned long long *p = (unsigned long long *)&LMASK[o * 64 + hh - 1];
      for (int j = 0; j < 4; j++) p[j] = (j < o) ? ~0ull : (j == o ? WMASK[hh] : 0ull);
    }
}

static inline int have_avx2(void) {
  unsigned a, b, c, d;
  __asm__ volatile("cpuid" : "=a"(a), "=b"(b), "=c"(c), "=d"(d) : "a"(1), "c"(0));
  int osxsave = (int)((c >> 27) & 1u), avx = (int)((c >> 28) & 1u);
  if (!(osxsave && avx)) return 0;
  unsigned lo, hi;
  __asm__ volatile("xgetbv" : "=a"(lo), "=d"(hi) : "c"(0));
  if ((lo & 6u) != 6u) return 0;
  __asm__ volatile("cpuid" : "=a"(a), "=b"(b), "=c"(c), "=d"(d) : "a"(7), "c"(0));
  return (int)((b >> 5) & 1u);
}

void solve(SOLVE_ARGS) {
  build_all(n, s1, s2);
  // group queries by the relative bit offset (sign aware), counting sort
  for (int i = 0; i <= 512; i++) GHEAD[i] = 0;
  for (int k = 0; k < q; k++) {
    int d = q_y[k] - q_x[k];
    int ad = (d >= 0 ? d : -d);
    int g = ((d >= 0) ? 0 : 64) + (ad & 63);
    GHEAD[g * 4 + ((ad >> 6) & 3) + 1]++;
  }
  for (int i = 0; i < 512; i++) GHEAD[i + 1] += GHEAD[i];
  for (int i = 0; i < 512; i++) GCUR[i] = GHEAD[i];
  for (int k = 0; k < q; k++) {
    int x = q_x[k], y = q_y[k], l = q_len[k];
    int d = y - x;
    int ad = (d >= 0 ? d : -d);
    int g = ((d >= 0) ? 0 : 64) + (ad & 63);
    int p = GCUR[g * 4 + ((ad >> 6) & 3)]++;
    ORD[p] = k;
    RQ[p] = (unsigned long long)(unsigned)x | ((unsigned long long)(unsigned)y << 21)
          | ((unsigned long long)(unsigned)l << 42);
  }
#ifdef TESTSCALAR
  run_scalar(n, q, q_x, q_y, q_len, ans);
#else
  if (have_avx2()) {
    run_avx2(n, q, q_x, q_y, q_len, ans);
    for (int t=0; t<q; ++t) ans[ORD[t]]=OUT[t];
  } else run_scalar(n, q, q_x, q_y, q_len, ans);
#endif
}

评测结果