代码

// 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 <memory>
#include <immintrin.h>
#pragma GCC target("avx2,popcnt")
#pragma GCC optimize("O3")

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:256;
    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;
    }
    const int PAD=(D+3)/4*4;
    std::vector<U> packedRanks(size_t(nn)*PAD,0);
    for(U i=0;i<nn;++i) for(int d=0;d<D;++d) packedRanks[size_t(i)*PAD+d]=rank[d][i];
    // Only store the prefix actually needed by a query using each row.
    const U blocks=(nn+block-1)/block;
    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;
        std::vector<U> lead[D-2],needed[D-2];
        std::vector<size_t> off[D-2];
        std::unique_ptr<W[]> table[D-2];
        std::vector<W> active(words,0),running(words,0);
        for(int k=0;k<D-2;++k) {
            int d=dims[k];
            needed[k].resize(blocks+1,0);off[k].resize(blocks+2,0);
            for(U q=0;q<nn;++q) if(owner[q]==U(a)) {
                U e=(lower[d][q]+block/2)/block;
                needed[k][e]=std::max(needed[k][e],(lower[a][q]+63)/64);
            }
            for(U e=0;e<=blocks;++e) off[k][e+1]=off[k][e]+((needed[k][e]+7)&~size_t(7));
            table[k].reset(new W[off[k].back()]);
            lead[k].resize(nn);
            for(U z=0;z<nn;++z) lead[k][z]=rank[a][ord[d][z]];
            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=needed[k][e];
                    if(len) std::memcpy(table[k].get()+off[k][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 && needed[k][blocks]) std::memcpy(table[k].get()+off[k][blocks],running.data(),needed[k][blocks]*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].get()+off[k][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[a]=ta;limits[sweep]=ts;
                for(int h=0;h<D-2;++h) if(h!=k) limits[dims[h]]=h<k?t[h]:boundary[h];
                U lo=std::min(t[k],boundary[k]),hi=std::max(t[k],boundary[k]),count=0;
                __m256i va=_mm256_set1_epi32(ta);
                const U *ld=lead[k].data();
                U z=lo;
                for(;z+8<=hi;z+=8) {
                    __m256i l=_mm256_loadu_si256((const __m256i*)(ld+z));
                    unsigned mask=unsigned(_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(va,l))));
                    while(mask) {
                        U bit=__builtin_ctz(mask);mask&=mask-1;
                        const U *r=packedRanks.data()+size_t(ord[dims[k]][z+bit])*PAD;
                        if(PAD==4) {
                            __m128i v=_mm_loadu_si128((const __m128i*)limits);
                            __m128i w=_mm_loadu_si128((const __m128i*)r);
                            count+=_mm_movemask_ps(_mm_castsi128_ps(_mm_cmpgt_epi32(v,w)))==15;
                        } else {
                            __m256i v=_mm256_loadu_si256((const __m256i*)limits);
                            __m256i w=_mm256_loadu_si256((const __m256i*)r);
                            if(_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(v,w)))!=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;
                        }
                    }
                }
                for(;z<hi;++z) {
                    if(ld[z]>=ta) continue;
                    const U *r=packedRanks.data()+size_t(ord[dims[k]][z])*PAD;
                    bool ok=true;
                    for(int h=0;h<D;++h) if(r[h]>=limits[h]) {ok=false;break;}
                    count+=ok;
                }
                ans+=t[k]>=boundary[k]?static_cast<long long>(count):-static_cast<long long>(count);
            }
            out[q]=U(ans);
        }
    }
}
} // namespace dominance

void count_9d(int n, const unsigned *x[9], unsigned *out) {
    dominance::count<9>(n,x,out);
}

评测结果