代码

/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/125840: retained accepted ten-million-point engine and SIMD leaf accumulator.
 * saffah_codex_6a_agg3, https://duck.ac/submission/121175: cumulative-counter concept from the million-point solver.
 * saffah_codex_maker, https://duck.ac/submission/125716: inherited CDQ, 24-bit records and exact ties; all source credits retained.
 * Idea: Store inclusive prefixes inside each16-child counter group. Insertions
 * and rollbacks add/subtract a SIMD suffix mask; each query reads one prefix
 * instead of horizontally reducing16 counters at every level. The final byte/
 * word entry can wrap exactly256/65536 but is never queried (partial counts
 * only read positions0..14); every queried prefix fits its unsigned storage.
 * Purpose: Official experiment trading wider updates for far cheaper hierarchy queries.
 */
/* References:
 * saffah_codex_maker, https://duck.ac/submission/125716: retained the accepted
 * ten-million-point 24-bit record CDQ implementation and its tie handling.
 * saffah_cc_v41_260924, https://duck.ac/submission/86372: inherited original
 * million-point CDQ design as credited by the public reference.
 * Idea: Accumulate comparison masks in four 64-bit SIMD counters inside the
 * small-node prefix scan, reducing to scalar only once after the scan. This
 * removes a movemask, popcount and scalar dependence on every four candidates.
 * Purpose: Official experiment improving the repeated small-node SIMD scan.
 */
// Adapted from public Judge Duck submission 86372 by saffah_cc_v41_260924:
// https://duck.ac/submission/86372 (original 1,000,000-point solver).
// Extension to n <= 10,000,000: 24-bit ranks/positions, split answer fields,
// one extra z-counter level, and direct x-order record construction.
// Only the input coordinates are read; no generator, seed, or fixed data.
// Coordinates may repeat. Every dimension is compared strictly (<).
#include <string.h>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2,popcnt")
#pragma GCC optimize("O3")
namespace duck3d_fast_public {
typedef unsigned u32;
typedef unsigned long long u64;
typedef unsigned char u8;
typedef unsigned short u16;

static const int MAXN = 10000001;

// W0: y[24] | x-order position[24] | answer low[16].
// W1: equal-z group start T[24] | distinct z rank R[24] | answer high[16].
// Splitting the answer retains a 16-byte record while supporting 24-bit ranks.
#define YP(v) ((u32)((v) >> 40))
#define PP(v) ((u32)(((v) >> 16) & 0xFFFFFFu))
#define RW(v) ((u32)(((v) >> 16) & 0xFFFFFFu))
#define TW(v) ((u32)((v) >> 40))
static inline void add_answer(u64 &w0, u64 &w1, u32 c) {
  const u32 low = u32(w0 & 0xFFFFu) + c;
  w0 = (w0 & ~u64(0xFFFFu)) | (low & 0xFFFFu);
  w1 += u64(low >> 16);
}
static inline u32 answer(u64 w0, u64 w1) {
  return u32(w0 & 0xFFFFu) | (u32(w1 & 0xFFFFu) << 16);
}

static u64 W0[2][MAXN];
static u64 W1[2][MAXN];
static u32 XS[MAXN];                 // x value at x-order position
static u32 XPM[MAXN];                // orig index at x-order position
static u32 HIST[MAXN + 1], HIST2[MAXN + 1], CNT[MAXN + 1];

// ---- z counter structure: rank domain [0, n) ----
static u64 ZB[(MAXN + 63) / 64 + 16];          // bit per rank
static u8  ZW1[(MAXN + 15) / 16 + 32];         // per 16 ranks
static u16 ZW2[(MAXN + 255) / 256 + 32];       // per 256 ranks
static u16 ZW3[(MAXN + 4095) / 4096 + 32];     // per 4096 ranks
static u32 ZW4[(MAXN + 65535) / 65536 + 32];   // per 65536 ranks
static u32 ZW5[(MAXN + 1048575) / 1048576 + 32]; // per 1048576 ranks

static inline u32 sum_u8(const u8 *p, u32 k) {  // k <= 16
  __m128i v = _mm_loadu_si128((const __m128i *)p);
  __m128i m = _mm_cmpgt_epi8(_mm_set1_epi8((char)k),
                             _mm_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15));
  __m128i s = _mm_sad_epu8(_mm_and_si128(v, m), _mm_setzero_si128());
  return (u32)_mm_cvtsi128_si32(s) + (u32)(_mm_extract_epi16(s, 4) & 0xFFFF);
}
static inline u32 sum_u16(const u16 *p, u32 k) {  // k <= 16
  __m256i v = _mm256_loadu_si256((const __m256i *)p);
  __m256i m = _mm256_cmpgt_epi16(_mm256_set1_epi16((short)k),
                                 _mm256_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15));
  __m256i x = _mm256_and_si256(v, m);
  __m256i s = _mm256_madd_epi16(x, _mm256_set1_epi16(1));
  __m128i lo = _mm256_castsi256_si128(s), hi = _mm256_extracti128_si256(s, 1);
  __m128i t = _mm_add_epi32(lo, hi);
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
  return (u32)_mm_cvtsi128_si32(t);
}
static inline u32 sum_u32(const u32 *p, u32 k) {  // k <= 16
  __m256i v0 = _mm256_loadu_si256((const __m256i *)p);
  __m256i v1 = _mm256_loadu_si256((const __m256i *)(p + 8));
  __m256i m0 = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)k), _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
  __m256i m1 = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)k), _mm256_setr_epi32(8, 9, 10, 11, 12, 13, 14, 15));
  __m256i x = _mm256_add_epi32(_mm256_and_si256(v0, m0), _mm256_and_si256(v1, m1));
  __m128i lo = _mm256_castsi256_si128(x), hi = _mm256_extracti128_si256(x, 1);
  __m128i t = _mm_add_epi32(lo, hi);
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
  return (u32)_mm_cvtsi128_si32(t);
}

alignas(32) static u8 P8[16][16];
alignas(32) static u16 P16[16][16];
alignas(32) static u32 P32[16][16];
static bool prefix_ready=false;
static void prefix_init(){
 if(prefix_ready)return;
 for(int i=0;i<16;i++)for(int j=0;j<16;j++)P8[i][j]=P16[i][j]=P32[i][j]=(j>=i);
 prefix_ready=true;
}
template<bool SUB> static inline void zchange(u32 r){
 if(SUB)ZB[r>>6]&=~(1ull<<(r&63));else ZB[r>>6]|=1ull<<(r&63);
 u32 a=r>>4,b=r>>8,c=r>>12,d=r>>16,e=r>>20;
 __m128i q8=_mm_loadu_si128((const __m128i*)(ZW1+(a&~15u)));
 __m128i m8=_mm_load_si128((const __m128i*)P8[a&15]);
 q8=SUB?_mm_sub_epi8(q8,m8):_mm_add_epi8(q8,m8);
 _mm_storeu_si128((__m128i*)(ZW1+(a&~15u)),q8);
 __m256i q2=_mm256_loadu_si256((const __m256i*)(ZW2+(b&~15u))),m2=_mm256_load_si256((const __m256i*)P16[b&15]);
 __m256i q3=_mm256_loadu_si256((const __m256i*)(ZW3+(c&~15u))),m3=_mm256_load_si256((const __m256i*)P16[c&15]);
 _mm256_storeu_si256((__m256i*)(ZW2+(b&~15u)),SUB?_mm256_sub_epi16(q2,m2):_mm256_add_epi16(q2,m2));
 _mm256_storeu_si256((__m256i*)(ZW3+(c&~15u)),SUB?_mm256_sub_epi16(q3,m3):_mm256_add_epi16(q3,m3));
 for(int h=0;h<16;h+=8){
  __m256i q4=_mm256_loadu_si256((const __m256i*)(ZW4+(d&~15u)+h)),m4=_mm256_load_si256((const __m256i*)(P32[d&15]+h));
  __m256i q5=_mm256_loadu_si256((const __m256i*)(ZW5+h)),m5=_mm256_load_si256((const __m256i*)(P32[e]+h));
  _mm256_storeu_si256((__m256i*)(ZW4+(d&~15u)+h),SUB?_mm256_sub_epi32(q4,m4):_mm256_add_epi32(q4,m4));
  _mm256_storeu_si256((__m256i*)(ZW5+h),SUB?_mm256_sub_epi32(q5,m5):_mm256_add_epi32(q5,m5));
 }
}
static inline void zadd(u32 r){zchange<false>(r);}
static inline void zsub(u32 r){zchange<true>(r);}
static inline u32 zquery(u32 T){
 u32 s=0,k5=T>>20,k4=(T>>16)&15,k3=(T>>12)&15,k2=(T>>8)&15,k1=(T>>4)&15;
 if(k5)s+=ZW5[k5-1];
 if(k4)s+=ZW4[((T>>20)<<4)+k4-1];
 if(k3)s+=ZW3[((T>>16)<<4)+k3-1];
 if(k2)s+=ZW2[((T>>12)<<4)+k2-1];
 if(k1)s+=ZW1[((T>>8)<<4)+k1-1];
 u32 rem=T&63;
 if(rem){u32 lo=rem&~15u;u64 mask=((1ull<<rem)-1)^((1ull<<lo)-1);s+=(u32)__builtin_popcountll(ZB[T>>6]&mask);}
 return s;
}


#ifndef SIMD_LIM
#define SIMD_LIM 1024
#endif

// count k in [l,i) with (w1[k] & mask) < T   (SIMD, 4 records per iteration)
static inline u32 count_lt(const u64 *w1, int l, int i, u32 T) {
  u32 c = 0;
  __m256i acc=_mm256_setzero_si256();
  int k = l;
  const __m256i m24 = _mm256_set1_epi64x(0xFFFFFF);
  const __m256i tv = _mm256_set1_epi64x((long long)T);
  for (; k + 4 <= i; k += 4) {
    __m256i v = _mm256_loadu_si256((const __m256i *)(w1 + k));
    __m256i r = _mm256_and_si256(_mm256_srli_epi64(v, 16), m24);
    __m256i cm = _mm256_cmpgt_epi64(tv, r);
    acc = _mm256_sub_epi64(acc,cm);
  }
  __m128i sums=_mm_add_epi64(_mm256_castsi256_si128(acc),_mm256_extracti128_si256(acc,1));
  c=(u32)((unsigned long long)_mm_cvtsi128_si64(sums)+(unsigned long long)_mm_extract_epi64(sums,1));
  for (; k < i; k++) if (RW(w1[k]) < T) c++;
  return c;
}

// common case: no x-ties span the midpoint
static void cross_clean(int l, int m, int r, int src) {
  u64 *w0 = W0[src];
  u64 *w1 = W1[src];
  int i = l;
  if (r - l < SIMD_LIM) {
    for (int j = m; j < r; j++) {
      u32 yj = YP(w0[j]);
      while (i < m && YP(w0[i]) < yj) i++;
      u32 c = count_lt(w1, l, i, TW(w1[j]));
      if (c) add_answer(w0[j], w1[j], c);
    }
  } else {
    for (int j = m; j < r; j++) {
      u32 yj = YP(w0[j]);
      while (i < m && YP(w0[i]) < yj) { zadd(RW(w1[i])); i++; }
      u32 c = zquery(TW(w1[j]));
      if (c) add_answer(w0[j], w1[j], c);
    }
    for (int k = l; k < i; k++) zsub(RW(w1[k]));
  }
}

// x-tie spans the midpoint: restrict queries/inserts by x-order position
static void cross_filter(int l, int m, int r, int src, u32 plo, u32 phi, u32 qlo, u32 qhi) {
  u64 *w0 = W0[src];
  u64 *w1 = W1[src];
  int i = l;
  if (r - l < SIMD_LIM) {
    for (int j = m; j < r; j++) {
      u32 pj = PP(w0[j]);
      if (pj < plo || pj >= phi) continue;
      u32 yj = YP(w0[j]);
      while (i < m && YP(w0[i]) < yj) i++;
      u32 T = TW(w1[j]), c = 0;
      for (int k = l; k < i; k++) {
        u32 pk = PP(w0[k]);
        if (pk >= qlo && pk < qhi && RW(w1[k]) < T) c++;
      }
      if (c) add_answer(w0[j], w1[j], c);
    }
  } else {
    for (int j = m; j < r; j++) {
      u32 pj = PP(w0[j]);
      if (pj < plo || pj >= phi) continue;
      u32 yj = YP(w0[j]);
      while (i < m && YP(w0[i]) < yj) {
        u32 pi = PP(w0[i]);
        if (pi >= qlo && pi < qhi) zadd(RW(w1[i]));
        i++;
      }
      u32 c = zquery(TW(w1[j]));
      if (c) add_answer(w0[j], w1[j], c);
    }
    for (int k = l; k < i; k++) {
      u32 pk = PP(w0[k]);
      if (pk >= qlo && pk < qhi) zsub(RW(w1[k]));
    }
  }
}

static void merge_run(int l, int m, int r, int src, int dst) {
  const u64 *s0 = W0[src], *s1 = W1[src];
  u64 *d0 = W0[dst], *d1 = W1[dst];
  int a = l, b = m, t = l;
  while (a < m && b < r) {
    int c = (YP(s0[a]) <= YP(s0[b]));
    int s = c ? a : b;
    d0[t] = s0[s]; d1[t] = s1[s];
    a += c; b += 1 - c; t++;
  }
  while (a < m) { d0[t] = s0[a]; d1[t] = s1[a]; a++; t++; }
  while (b < r) { d0[t] = s0[b]; d1[t] = s1[b]; b++; t++; }
}

// returns the buffer holding the y-sorted range [l,r)
static int cdq(int l, int r, int in) {
  if (r - l <= 1) return in;
  int m = (l + r) >> 1;
  int b1 = cdq(l, m, in);
  int b2 = cdq(m, r, in);
  if (b1 != b2) {
    memcpy(W0[b1] + m, W0[b2] + m, (size_t)(r - m) * 8);
    memcpy(W1[b1] + m, W1[b2] + m, (size_t)(r - m) * 8);
    b2 = b1;
  }
  if (XS[m - 1] < XS[m]) {
    cross_clean(l, m, r, b1);
  } else {
    u32 g = XS[m];
    int gs = m - 1; while (gs > l && XS[gs - 1] == g) gs--;
    int ge = m;     while (ge < r && XS[ge] == g) ge++;
    if (ge < r) cross_filter(l, m, r, b1, (u32)ge, 0xFFFFFFFFu, 0, 0xFFFFFFFFu);
    if (gs > l) cross_filter(l, m, r, b1, (u32)m, (u32)ge, 0, (u32)gs);
  }
  int out = 1 - b1;
  merge_run(l, m, r, b1, out);
  return out;
}

static void solve(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
  if (n <= 1) { if (n == 1) out[0] = 0; return; }

  // ---------- z compression: R (unique rank), T (group start rank) ----------
  memset(HIST, 0, (size_t)(n + 1) * 4);
  for (int i = 0; i < n; i++) HIST[z[i]]++;
  {
    u32 acc = 0;
    for (int v = 0; v < n; v++) { u32 c = HIST[v]; HIST[v] = acc; acc += c; }
    HIST[n] = acc;
  }
  memcpy(HIST2, HIST, (size_t)n * 4);

  // ---------- x-order scatter, no packed pre-record scratch ----------
  memset(CNT, 0, (size_t)n * sizeof(u32));
  for (int i = 0; i < n; ++i) ++CNT[x[i]];
  {
    u32 acc = 0;
    for (int v = 0; v < n; ++v) { const u32 c = CNT[v]; CNT[v] = acc; acc += c; }
  }
  for (int i = 0; i < n; ++i) {
    const u32 xv = x[i], zv = z[i], p = CNT[xv]++;
    const u32 T = HIST[zv], R = HIST2[zv]++;
    W0[0][p] = (u64(y[i]) << 40) | (u64(p) << 16);
    W1[0][p] = (u64(T) << 40) | (u64(R) << 16);
    XS[p] = xv;
    XPM[p] = u32(i);
  }
  // ---------- CDQ ----------
  int b = cdq(0, n, 0);
  const u64 *s0 = W0[b];
  for (int p = 0; p < n; p++) out[XPM[PP(s0[p])]] = answer(s0[p], W1[b][p]);
}

} // namespace duck3d_fast_public
void count_3d(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
  // Statically provisioned for the requested 10,000,000-point upper bound.
  if (n > 10000000) __builtin_trap();
  duck3d_fast_public::prefix_init();
  duck3d_fast_public::solve(n,x,y,z,out);
}
#pragma GCC pop_options

评测结果