代码

// Original implementation for these contribution packages, C++11 or later.
// Exact strict dominance: rank is a position; lower is the first position
// of the coordinate's equal-value group. No input array is modified.
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <vector>
#include <immintrin.h>
#pragma GCC target("avx2,popcnt")

namespace dominance {
using U = unsigned;
using W = unsigned long long;

template<int K> struct AndWords {
    static inline __m256i get(const W *const *p, U w) {
        return _mm256_and_si256(AndWords<K-1>::get(p,w),
            _mm256_loadu_si256((const __m256i*)(p[K-1]+w)));
    }
    static inline W scalar(const W *const *p, U w) {
        return AndWords<K-1>::scalar(p,w) & p[K-1][w];
    }
};
template<> struct AndWords<1> {
    static inline __m256i get(const W *const *p,U w) {
        return _mm256_loadu_si256((const __m256i*)(p[0]+w));
    }
    static inline W scalar(const W *const *p,U w) { return p[0][w]; }
};
template<int K> static U intersection(const W *const *p,U bits) {
    const __m256i lut=_mm256_setr_epi8(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);
    const __m256i mask=_mm256_set1_epi8(15);
    __m256i acc=_mm256_setzero_si256();
    U w=0,full=bits/64,ans=0;
    for(;w+4<=full;w+=4) {
        __m256i v=AndWords<K>::get(p,w);
        __m256i lo=_mm256_shuffle_epi8(lut,_mm256_and_si256(v,mask));
        __m256i hi=_mm256_shuffle_epi8(lut,_mm256_and_si256(_mm256_srli_epi16(v,4),mask));
        acc=_mm256_add_epi64(acc,_mm256_sad_epu8(_mm256_add_epi8(lo,hi),_mm256_setzero_si256()));
    }
    W sums[4];_mm256_storeu_si256((__m256i*)sums,acc);
    ans=U(sums[0]+sums[1]+sums[2]+sums[3]);
    for(;w<full;++w) ans+=__builtin_popcountll(AndWords<K>::scalar(p,w));
    if(bits%64) ans+=__builtin_popcountll(AndWords<K>::scalar(p,full)&((1ULL<<(bits%64))-1));
    return ans;
}

template<int D> void count(int n,const U *const *x,U *out) {
    if(n<=0) return;
    const U nn=U(n), words=(nn+63)/64;
    const U block=nn<=100000?64:nn<=300000?128:D<=5?256:512;
    std::vector<U> ord[D],rank[D],lower[D],hist(size_t(nn)+1),owner(nn);
    for(int d=0;d<D;++d) {
        ord[d].resize(nn);rank[d].resize(nn);lower[d].resize(nn);
        std::fill(hist.begin(),hist.end(),0);
        for(U i=0;i<nn;++i) ++hist[x[d][i]+1];
        for(U v=1;v<=nn;++v) hist[v]+=hist[v-1];
        for(U i=0;i<nn;++i) lower[d][i]=hist[x[d][i]];
        for(U i=0;i<nn;++i) ord[d][hist[x[d][i]]++]=i;
        for(U p=0;p<nn;++p) rank[d][ord[d][p]]=p;
    }
    for(U i=0;i<nn;++i) {
        U a=0;for(int d=1;d<D;++d) if(lower[d][i]<lower[a][i]) a=d;
        owner[i]=a;out[i]=0;
    }
    // Row e stores the first ceil(min(n,(e+1)*block)/64) bits.
    // Each query owned by a has lower[a] <= lower[d], so this is enough.
    const U blocks=(nn+block-1)/block;
    std::vector<size_t> off(blocks+2,0);
    for(U e=0;e<=blocks;++e) {
        size_t len=(std::min(nn,(e+1)*block)+63)/64;
        off[e+1]=off[e]+((len+7)&~size_t(7));
    }
    for(int a=0;a<D;++a) {
        bool used=false;
        for(U i=0;i<nn;++i) if(owner[i]==U(a) && lower[a][i]) {used=true;break;}
        if(!used) continue;
        const int sweep=(a+1)%D;
        int dims[D-2],nd=0;
        for(int d=0;d<D;++d) if(d!=a && d!=sweep) dims[nd++]=d;
        const int PAD=(D+3)/4*4;
        std::vector<U> records[D-2];
        std::vector<W> table[D-2],active(words,0),running(words,0);
        for(int k=0;k<D-2;++k) {
            int d=dims[k];table[k].resize(off.back());
            records[k].resize(size_t(nn)*PAD,0);
            for(U z=0;z<nn;++z) {
                U j=ord[d][z];
                U *r=records[k].data()+size_t(z)*PAD;
                r[0]=rank[a][j];r[1]=rank[sweep][j];
                for(int h=0;h<D-2;++h) r[h+2]=rank[dims[h]][j];
            }
            std::fill(running.begin(),running.end(),0);
            for(U p=0;p<=nn;++p) {
                if(p%block==0) {
                    U e=p/block;
                    size_t len=(std::min(nn,(e+1)*block)+63)/64;
                    std::memcpy(table[k].data()+off[e],running.data(),len*sizeof(W));
                }
                if(p==nn) break;
                U r=rank[a][ord[d][p]];running[r/64]|=1ULL<<(r%64);
            }
            if(nn%block) std::memcpy(table[k].data()+off[blocks],running.data(),words*sizeof(W));
        }
        U inserted=0;
        for(U pos=0;pos<nn;++pos) {
            U q=ord[sweep][pos];
            if(owner[q]!=U(a) || lower[a][q]==0) continue;
            U ta=lower[a][q],ts=lower[sweep][q],t[D-2],boundary[D-2];
            while(inserted<ts) {
                U r=rank[a][ord[sweep][inserted++]];active[r/64]|=1ULL<<(r%64);
            }
            const W *p[D-1];p[0]=active.data();
            for(int k=0;k<D-2;++k) {
                t[k]=lower[dims[k]][q];
                U e=(t[k]+block/2)/block;
                boundary[k]=std::min(nn,e*block);
                p[k+1]=table[k].data()+off[e];
            }
            long long ans=intersection<D-1>(p,ta);
            // Replace one rounded boundary at a time. Earlier dimensions
            // use exact bounds; later dimensions still use rounded bounds.
            for(int k=0;k<D-2;++k) {
                U limits[PAD];std::fill(limits,limits+PAD,nn);
                limits[0]=ta;limits[1]=ts;
                for(int h=0;h<D-2;++h) if(h!=k) limits[h+2]=h<k?t[h]:boundary[h];
                U lo=std::min(t[k],boundary[k]),hi=std::max(t[k],boundary[k]),count=0;
                if(PAD==4) {
                    __m128i v=_mm_loadu_si128((const __m128i*)limits);
                    for(U z=lo;z<hi;++z) {
                        __m128i r=_mm_loadu_si128((const __m128i*)(records[k].data()+size_t(z)*PAD));
                        count+=_mm_movemask_ps(_mm_castsi128_ps(_mm_cmpgt_epi32(v,r)))==15;
                    }
                } else {
                    __m256i v=_mm256_loadu_si256((const __m256i*)limits);
                    for(U z=lo;z<hi;++z) {
                        const U *r=records[k].data()+size_t(z)*PAD;
                        __m256i m=_mm256_cmpgt_epi32(v,_mm256_loadu_si256((const __m256i*)r));
                        if(_mm256_movemask_ps(_mm256_castsi256_ps(m))!=255) continue;
                        if(PAD>8) {
                            __m128i v2=_mm_loadu_si128((const __m128i*)(limits+8));
                            __m128i r2=_mm_loadu_si128((const __m128i*)(r+8));
                            if(_mm_movemask_ps(_mm_castsi128_ps(_mm_cmpgt_epi32(v2,r2)))!=15) continue;
                        }
                        ++count;
                    }
                }
                ans+=t[k]>=boundary[k]?static_cast<long long>(count):-static_cast<long long>(count);
            }
            out[q]=U(ans);
        }
    }
}
} // 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=5,N=1000000;
    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::count<5>(N,x,out.data());
    for(int q=0;q<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: 32 sampled answers checked");
}

评测结果