代码

// Exact 2D strict-dominance counting for 0 <= coordinates < n <= 100000000.
// Original implementation: 9-bit in-place radix sorting, then an AVX2 16-ary
// prefix tree. All equal-x queries precede updates. Digit-zero prefixes are
// always zero, so their unused slots hold overflow-record indices. Other slots
// store low 16 bits; sparse records store the high 16 bits without approximation.
// At a fixed level, an overflowed node contains at least 65536 inserted points;
// nodes are disjoint, so levels*(n/65536+1)+2 records suffice. Metadata indices
// remain below 65536 for the supported bounds. The packed sort allocation is
// reused for index order plus the tree, never modifying either input array.
#include <immintrin.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC push_options
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
namespace duck_count2_impl {
// Non-stable radix partition by x. Eight independent cycles hide cache-miss latency.
static void radix(uint64_t*a,unsigned n,int shift){
 if(n<=32){for(unsigned i=1;i<n;++i){uint64_t v=a[i];unsigned j=i;while(j&&a[j-1]>v){a[j]=a[j-1];--j;}a[j]=v;}return;}
 constexpr unsigned B=512,L=8;
 unsigned counts[B]={},start[B+1],next[B];
 for(unsigned i=0;i<n;++i)++counts[(a[i]>>(32+shift))&(B-1)];
 unsigned total=0;for(unsigned d=0;d<B;++d){start[d]=next[d]=total;total+=counts[d];}start[B]=total;
 for(unsigned d=0;d<B;++d){
  uint64_t values[L];unsigned holes[L],dest[L];bool active[L];
  unsigned live=0;
  for(unsigned lane=0;lane<L;++lane){active[lane]=next[d]<start[d+1];if(active[lane]){holes[lane]=next[d]++;values[lane]=a[holes[lane]];++live;}}
  while(live){
#pragma GCC unroll 8
   for(unsigned lane=0;lane<L;++lane){
    dest[lane]=~0u;if(!active[lane])continue;
    unsigned target=(values[lane]>>(32+shift))&(B-1);
    if(target==d){a[holes[lane]]=values[lane];if(next[d]<start[d+1]){holes[lane]=next[d]++;values[lane]=a[holes[lane]];}else{active[lane]=false;--live;}}
    else dest[lane]=next[target]++;
   }
#pragma GCC unroll 8
   for(unsigned lane=0;lane<L;++lane)if(dest[lane]!=~0u){uint64_t v=a[dest[lane]];a[dest[lane]]=values[lane];values[lane]=v;}
  }
 }
 if(shift)for(unsigned d=0;d<B;++d)if(counts[d])radix(a+start[d],counts[d],shift-9);
}
void solve(int n,const unsigned*x,const unsigned*y,unsigned*out){
 if(n<=0)return;
 size_t lengths[10],offsets[10],all=0;unsigned levels=0,v=n;
 do{lengths[levels]=(v+15u)&~15u;offsets[levels]=all;all+=lengths[levels++];v=(v+15)/16;}while(v>1);
 size_t padded=(size_t(n)+7)&~size_t(7);
 void *storage=nullptr;if(posix_memalign(&storage,4096,padded*8))abort();
 uint64_t*packed=(uint64_t*)storage;
 for(int i=0;i<n;++i)packed[i]=((uint64_t)x[i]<<32)|unsigned(i);
 int shift=0;while((uint64_t(1)<<(shift+9))<unsigned(n))shift+=9;
 radix(packed,unsigned(n),shift);
 unsigned*order=(unsigned*)packed,*tree=order+padded;
 unsigned previous=~0u;
 for(int i=0;i<n;++i){uint64_t v;memcpy(&v,packed+i,sizeof(v));unsigned xv=unsigned(v>>32),id=unsigned(v)|(xv!=previous?0x80000000u:0);memcpy(order+i,&id,sizeof(id));previous=xv;}
 unsigned short *compact=(unsigned short*)tree;
 memset(compact,0,all*2);
 unsigned cap=levels*(unsigned(n)/65536+1)+2,used=0;
 if(cap>=65536)abort();
 void *high_storage=nullptr;if(posix_memalign(&high_storage,64,(size_t)cap*32))abort();
 unsigned short *highs=(unsigned short*)high_storage;memset(highs,0,(size_t)cap*32);
 alignas(32) __m256i masks[16];
 for(unsigned d=0;d<16;++d){unsigned short temp[16];for(unsigned j=0;j<16;++j)temp[j]=j>d;masks[d]=_mm256_loadu_si256((const __m256i*)temp);}
 const __m256i ones=_mm256_set1_epi16(-1);
 for(unsigned left=0;left<(unsigned)n;){
  unsigned right=left+1;while(right<(unsigned)n&&!(order[right]&0x80000000u))++right;

  for(unsigned q=left;q<right;++q){
   if(q+128<(unsigned)n){unsigned id=order[q+128]&0x7fffffffu;__builtin_prefetch(y+id);__builtin_prefetch(out+id,1);}
   if(q+32<(unsigned)n){unsigned yy=y[order[q+32]&0x7fffffffu];for(unsigned l=0;l<levels;++l){__builtin_prefetch(compact+offsets[l]+yy,1);yy>>=4;}}
   unsigned id=order[q]&0x7fffffffu,yy=y[id],ans=0;
   for(unsigned l=0;l<levels;++l){unsigned d=yy&15;if(d){const unsigned short*p=compact+offsets[l]+(yy&~15u);ans+=p[d];if(p[0])ans+=(unsigned)highs[(size_t)p[0]*16+d]<<16;}yy>>=4;}
   out[id]=ans;
  }
  for(unsigned q=left;q<right;++q){
   unsigned id=order[q]&0x7fffffffu,yy=y[id];
   for(unsigned l=0;l<levels;++l){
    unsigned d=yy&15;unsigned short*p=compact+offsets[l]+(yy&~15u);
    __m256i a=_mm256_load_si256((const __m256i*)p);
    __m256i carry=_mm256_and_si256(_mm256_cmpeq_epi16(a,ones),masks[d]);
    _mm256_store_si256((__m256i*)p,_mm256_add_epi16(a,masks[d]));
    if(!_mm256_testz_si256(carry,carry)){
     if(!p[0]){if(++used>=cap)abort();p[0]=(unsigned short)used;}
     unsigned short*h=highs+(size_t)p[0]*16;
     _mm256_store_si256((__m256i*)h,_mm256_add_epi16(_mm256_load_si256((const __m256i*)h),carry));
    }
    yy>>=4;
   }
  }
  left=right;
 }
 free(highs);free(packed);
}
}
namespace dominance {void count2(int n,const unsigned*x,const unsigned*y,unsigned*out){duck_count2_impl::solve(n,x,y,out);}}
#pragma GCC pop_options

void count_2d(int n, const unsigned *x, const unsigned *y, unsigned *out) {
    dominance::count2(n,x,y,out);
}

评测结果