代码

#define BLOCK_SIZE 16
#include <immintrin.h>

typedef unsigned int u32;
typedef unsigned long u64;
typedef unsigned long U;

#ifndef BLOCK_SIZE
#define BLOCK_SIZE 64
#endif

enum {
    CAPACITY = 100005,
    BLOCK = BLOCK_SIZE,
    WORDS = (CAPACITY + 63) / 64,
    BLOCKS = (CAPACITY + BLOCK - 1) / BLOCK + 1
};

static int sorted_id[CAPACITY], sorted_dim[3][CAPACITY];
static int count_[CAPACITY], start_[CAPACITY];
static int rank_position[3][CAPACITY];
static int lower_rank[3][CAPACITY];
static int seen[CAPACITY];
static u64 active[WORDS], current[WORDS];
static u64 prefix_bits[3][BLOCKS][WORDS];

__attribute__((target("popcnt"), always_inline))
static inline u32 intersection_count(int words, const u64 *a, const u64 *b,
                                     const u64 *c, const u64 *d) {
    u64 sum = 0;
    int i = 0;
    for (; i + 4 <= words; i += 4) {
        sum += (u64)__builtin_popcountll(a[i] & b[i] & c[i] & d[i]);
        sum += (u64)__builtin_popcountll(a[i + 1] & b[i + 1] &
                                         c[i + 1] & d[i + 1]);
        sum += (u64)__builtin_popcountll(a[i + 2] & b[i + 2] &
                                         c[i + 2] & d[i + 2]);
        sum += (u64)__builtin_popcountll(a[i + 3] & b[i + 3] &
                                         c[i + 3] & d[i + 3]);
    }
    for (; i < words; ++i)
        sum += (u64)__builtin_popcountll(a[i] & b[i] & c[i] & d[i]);
    return (u32)sum;
}

static void build_dimension(int d, int n, int words, const u32 *coordinate) {
    for (int i = 0; i < n; ++i) count_[i] = 0;
    for (int i = 0; i < n; ++i) ++count_[coordinate[i]];
    int sum = 0;
    for (int value = 0; value < n; ++value) {
        int c = count_[value];
        start_[value] = sum;
        count_[value] = sum;
        sum += c;
    }
    for (int i = 0; i < n; ++i) {
        int pos = count_[coordinate[i]]++;
        sorted_dim[d][pos] = i;
        rank_position[d][i] = pos;
        lower_rank[d][i] = start_[coordinate[i]];
    }

    for (int w = 0; w < words; ++w) current[w] = 0;
    int block_count = (n + BLOCK - 1) / BLOCK;
    for (int b = 1; b <= block_count; ++b) {
        int begin = (b - 1) * BLOCK;
        int end = b * BLOCK;
        if (end > n) end = n;
        for (int p = begin; p < end; ++p) {
            int id = sorted_dim[d][p];
            current[id >> 6] |= 1ul << (id & 63);
        }
        u64 *dst = prefix_bits[d][b];
        for (int w = 0; w < words; ++w) dst[w] = current[w];
    }
}

__attribute__((target("popcnt")))
void count_4d(int n, const u32 *x[4], u32 *out) {
    int words = (n + 63) >> 6;
    for (int i = 0; i < n; ++i) {
        out[i] = 0;
        seen[i] = 0;
    }
    for (int w = 0; w < words; ++w) active[w] = 0;

    for (int d = 0; d < 3; ++d)
        build_dimension(d, n, words, x[d + 1]);

    /* Stable counting sort by x[0]. */
    for (int i = 0; i < n; ++i) count_[i] = 0;
    for (int i = 0; i < n; ++i) ++count_[x[0][i]];
    int sum = 0;
    for (int value = 0; value < n; ++value) {
        int c = count_[value];
        count_[value] = sum;
        sum += c;
    }
    for (int i = 0; i < n; ++i) sorted_id[count_[x[0][i]]++] = i;

#ifdef PREPROCESS_ONLY
    return;
#endif

    int stamp = 0;
    for (int first = 0; first < n;) {
        int last = first + 1;
        u32 value = x[0][sorted_id[first]];
        while (last < n && x[0][sorted_id[last]] == value) ++last;

        for (int qi = first; qi < last; ++qi) {
            int id = sorted_id[qi];
            int lower0 = lower_rank[0][id];
            int lower1 = lower_rank[1][id];
            int lower2 = lower_rank[2][id];
            int b0 = (lower0 + BLOCK - 1) / BLOCK;
            int b1 = (lower1 + BLOCK - 1) / BLOCK;
            int b2 = (lower2 + BLOCK - 1) / BLOCK;
            int end0 = b0 * BLOCK; if (end0 > n) end0 = n;
            int end1 = b1 * BLOCK; if (end1 > n) end1 = n;
            int end2 = b2 * BLOCK; if (end2 > n) end2 = n;

            u32 result = intersection_count(words, active,
                                              prefix_bits[0][b0],
                                              prefix_bits[1][b1],
                                              prefix_bits[2][b2]);
#ifndef NO_TAIL
            ++stamp;
            for (int d = 0; d < 3; ++d) {
                int lower = d == 0 ? lower0 : d == 1 ? lower1 : lower2;
                int end = d == 0 ? end0 : d == 1 ? end1 : end2;
                for (int p = lower; p < end; ++p) {
                    int candidate;
                    candidate = sorted_dim[d][p];
                    if (seen[candidate] == stamp) continue;
                    seen[candidate] = stamp;
                    if (((active[candidate >> 6] >> (candidate & 63)) & 1u) &&
                        rank_position[0][candidate] < end0 &&
                        rank_position[1][candidate] < end1 &&
                        rank_position[2][candidate] < end2)
                        --result;
                }
            }
#endif
            out[id] = result;
        }

        for (int p = first; p < last; ++p) {
            int id = sorted_id[p];
            active[id >> 6] |= 1ul << (id & 63);
        }
        first = last;
    }
}

#ifndef LOCAL
static char **initial_argv;
static U initial_argc;
U getauxval(U key) {
    char **p = initial_argv + initial_argc + 1;
    while (*p) ++p;
    U *aux = (U *)(p + 1);
    while (aux[0]) {
        if (aux[0] == key) return aux[1];
        aux += 2;
    }
    return 0;
}
__attribute__((noreturn))
void __libc_start_main(int (*entry)(int, char **, char **), int argc,
                       char **argv) {
    initial_argc = (U)argc;
    initial_argv = argv;
    entry(argc, argv, (char **)0);
    __asm__ volatile("mov $60,%%eax;xor %%edi,%%edi;syscall"
                     ::: "rax", "rdi", "rcx", "r11", "memory");
    __builtin_unreachable();
}
#endif

评测结果