代码

// Original implementation: radix-sort indices by x, then a Fenwick sweep.
// Queries precede updates for each equal-x group; prefix excludes equal y.
// Working storage: 8*n+O(1) bytes beyond the caller's input/output buffers.
#include <algorithm>
#include <cstddef>
#include <vector>
#include <cstdint>
#include <cstring>
#pragma GCC target("popcnt")
namespace dominance {
template<class Emit> void count2_each(int n,const unsigned *x,const unsigned *y,Emit emit) {
    if(n<=0) return;
    std::vector<unsigned> index(n),scratch(n);
    for(int i=0;i<n;++i) index[i]=unsigned(i);
    for(unsigned shift=0;shift<32;shift+=8) {
        unsigned hist[256]={};
        for(int i=0;i<n;++i) ++hist[(x[index[i]]>>shift)&255];
        unsigned sum=0;
        for(unsigned &v:hist) {unsigned c=v;v=sum;sum+=c;}
        for(int i=0;i<n;++i) {unsigned j=index[i];scratch[hist[(x[j]>>shift)&255]++]=j;}
        index.swap(scratch);
    }
    std::fill(scratch.begin(),scratch.end(),0);
    for(int left=0;left<n;) {
        int right=left+1;
        while(right<n && x[index[left]]==x[index[right]]) ++right;
        for(int z=left;z<right;++z) {
            unsigned j=index[z],sum=0;
            for(unsigned k=y[j];k;k&=k-1) sum+=scratch[k-1];
            emit(j,sum);
        }
        for(int z=left;z<right;++z)
            for(unsigned k=y[index[z]];k<unsigned(n);k|=k+1) ++scratch[k];
        left=right;
    }
}
// Stable MSD radix counting. The two index buffers use 8*n bytes.
// Equal-x points are sorted by descending y, so none can count each other.
void count2(int n,const unsigned *x,const unsigned *y,unsigned *out) {
    if(n<=0) return;
    std::vector<unsigned> a(n),b(n,0);
    for(int i=0;i<n;++i) ++b[x[i]];
    unsigned total=0;
    for(unsigned &v:b) {unsigned c=v;v=total;total+=c;}
    for(int i=0;i<n;++i) a[b[x[i]]++]=unsigned(i);
    for(int l=0;l<n;) {
        int r=l+1;while(r<n && x[a[r]]==x[a[l]]) ++r;
        if(r>l+1) std::sort(a.begin()+l,a.begin()+r,[y](unsigned i,unsigned j){return y[i]>y[j];});
        l=r;
    }
    std::fill(out,out+n,0);
    std::vector<unsigned> groups(1,0),next;
    groups.push_back(unsigned(n));
    int shift=0;while((uint64_t(1)<<(shift+6))<unsigned(n)) shift+=6;
    for(;shift>=0;shift-=6) {
        next.clear();
        for(size_t g=0;g+1<groups.size();++g) {
            unsigned l=groups[g],r=groups[g+1],cnt[64]={},dest[64];
            if(shift) {
                for(unsigned z=l;z<r;++z) ++cnt[(y[a[z]]>>shift)&63];
                unsigned at=l;
                for(unsigned d=0;d<64;++d) {dest[d]=at;if(cnt[d])next.push_back(at);at+=cnt[d];cnt[d]=0;}
            }
            for(unsigned base=l;base<r;base+=64) {
                unsigned len=std::min(64u,r-base),dig[64],prefix[64];
                uint64_t pos[64]={},lt[64],mask=0;
                for(unsigned i=0;i<len;++i) {
                    unsigned d=(y[a[base+i]]>>shift)&63;
                    dig[i]=d;pos[d]|=uint64_t(1)<<i;
                }
                unsigned sum=0;
                for(unsigned d=0;d<64;++d) {lt[d]=mask;mask|=pos[d];prefix[d]=sum;sum+=cnt[d];}
                for(unsigned i=0;i<len;++i) {
                    unsigned d=dig[i],id=a[base+i];
                    out[id]+=prefix[d]+__builtin_popcountll(lt[d]&((uint64_t(1)<<i)-1));
                    ++cnt[d];if(shift)b[dest[d]++]=id;
                }
            }
        }
        if(shift) {next.push_back(unsigned(n));groups.swap(next);a.swap(b);}
    }
}
} // namespace dominance

#include <cassert>
#include <cstdio>
static unsigned long long state=0x7145aa819275cdefULL;
static unsigned rnd() {state+=0x9e3779b97f4a7c15ULL;unsigned long long z=state;z=(z^(z>>30))*0xbf58476d1ce4e5b9ULL;z=(z^(z>>27))*0x94d049bb133111ebULL;return unsigned(z^(z>>31));}
int main() {
    const int D=2,N=100000000;
    std::vector<unsigned> coords[D],out(N);
    const unsigned *x[D];
    for(int d=0;d<D;++d) {coords[d].resize(N);x[d]=coords[d].data();for(int i=0;i<N;++i) coords[d][i]=rnd()%N;}
    dominance::count2(N,x[0],x[1],out.data());
    for(int q=0;q<(N>10000000?2:32);++q) {unsigned i=rnd()%N,want=0;
        for(int j=0;j<N;++j) {bool ok=true;for(int d=0;d<D;++d) if(x[d][j]>=x[d][i]) {ok=false;break;} want+=ok;}
        assert(out[i]==want);
    }
    std::puts("dominance benchmark: sampled answers checked");
}

评测结果