<|start_filename|>benchmark/linux-perf-events-wrapper.h<|end_filename|>
#pragma once
#include <memory>
#include <memory>
#include <unordered_map>
#ifdef __linux__
#include "linux-perf-events.h"
#endif

#ifdef __linux__
typedef LinuxEvents<PERF_TYPE_HARDWARE> EventClass;
#endif

class LinuxEventsWrapper {
  public:
    LinuxEventsWrapper(const std::vector<int> event_codes) {
#ifdef __linux__
      for(int ecode: event_codes) {
        event_obj.emplace(ecode, std::shared_ptr<EventClass>(new EventClass(ecode)));
        event_res.emplace(ecode, 0);
      }
#endif
    }
    void start() {
#ifdef __linux__
      for (const auto& [ecode, ptr]: event_obj) {
        ptr->start();  
      }
#endif
    }
    void end() {
#ifdef __linux__
      for (const auto& [ecode, ptr]: event_obj) {
        event_res[ecode] = ptr->end();  
      }
#endif
    }
    // Throws an exception if the code is not present
    unsigned long get_result(int ecode) {
#ifdef __linux__
      return event_res.at(ecode);
#else
      return 0;
#endif
    }
  private:
#ifdef __linux__
    std::unordered_map<int, std::shared_ptr<EventClass>> event_obj;
    std::unordered_map<int, unsigned long> event_res;
#endif
};


<|start_filename|>benchmark/ztimer.h<|end_filename|>
/**
 * This code is released under the
 * Apache License Version 2.0 http://www.apache.org/licenses/.
 *
 * (c) <NAME>, http://lemire.me/en/
 */

#ifndef ZTIMER_2014
#define ZTIMER_2014

#include <sys/time.h>
#include <sys/resource.h>


/**
 *  author: <NAME>
 */
class WallClockTimer {
public:
    struct timeval t1, t2;
    WallClockTimer() :
        t1(), t2() {
        gettimeofday(&t1, 0);
        t2 = t1;
    }
    void reset() {
        gettimeofday(&t1, 0);
        t2 = t1;
    }
    uint64_t elapsed() {
        return ((t2.tv_sec - t1.tv_sec) * 1000ULL * 1000ULL) + ((t2.tv_usec - t1. tv_usec));
    }
    uint64_t split() {
        gettimeofday(&t2, 0);
        return elapsed();
    }
};

/**
 *  author: <NAME>
 */
class CPUTimer {
public:
    //clock_t t1, t2;
    struct rusage t1,t2;

    CPUTimer() :
        t1(), t2() {
        getrusage(RUSAGE_SELF, &t1);
        //t1 = clock();
        t2 = t1;
    }
    void reset() {
        getrusage(RUSAGE_SELF, &t1);
        t2 = t1;
    }
    // proxy for userelapsed
    uint64_t elapsed() {
        return totalelapsed();
    }

    uint64_t totalelapsed() {
        return userelapsed() + systemelapsed();
    }
    // returns the *user* CPU time in micro seconds (mu s)
    uint64_t userelapsed() {
        return ((t2.ru_utime.tv_sec - t1.ru_utime.tv_sec) * 1000ULL * 1000ULL) + ((t2.ru_utime.tv_usec - t1.ru_utime.tv_usec)
                );
    }

    // returns the *system* CPU time in micro seconds (mu s)
    uint64_t systemelapsed() {
        return ((t2.ru_stime.tv_sec - t1.ru_stime.tv_sec) * 1000ULL * 1000ULL) + ((t2.ru_stime.tv_usec - t1.ru_stime.tv_usec)
                );
    }

    uint64_t split() {
        getrusage(RUSAGE_SELF, &t2);
        return elapsed();
    }
};

#endif



<|start_filename|>benchmark/linux-perf-events.h<|end_filename|>
#pragma once

#include <unistd.h>             // for syscall
#include <sys/ioctl.h>          // for ioctl
#include <asm/unistd.h>         // for __NR_perf_event_open
#include <linux/perf_event.h>   // for perf event constants

#include <cerrno>               // for errno
#include <cstring>              // for memset
#include <stdexcept>


template <int TYPE = PERF_TYPE_HARDWARE>
class LinuxEvents {

    int fd;
    perf_event_attr attribs;

public:
    LinuxEvents(int config) : fd(0) {
        memset(&attribs, 0, sizeof(attribs));
        attribs.type        = TYPE;
        attribs.size        = sizeof(attribs);
        attribs.config      = config;
        attribs.disabled        = 1;
        attribs.exclude_kernel  = 1;
        attribs.exclude_hv      = 1;

        const int pid = 0;    // the current process
        const int cpu = -1;   // all CPUs
        const int group = -1; // no group
        const unsigned long flags = 0;
        fd = syscall(__NR_perf_event_open, &attribs, pid, cpu, group, flags);
        if (fd == -1) {
            report_error("perf_event_open");
        }
    }

    ~LinuxEvents() {
        close(fd);
    }

    void start() {
        if (ioctl(fd, PERF_EVENT_IOC_RESET, 0) == -1) {
            report_error("ioctl(PERF_EVENT_IOC_RESET)");
        }

        if (ioctl(fd, PERF_EVENT_IOC_ENABLE, 0) == -1) {
            report_error("ioctl(PERF_EVENT_IOC_ENABLE)");
        }
    }

    unsigned long end() {
        if (ioctl(fd, PERF_EVENT_IOC_DISABLE, 0) == -1) {
            report_error("ioctl(PERF_EVENT_IOC_DISABLE)");
        }

        unsigned long result;
        if (read(fd, &result, sizeof(result)) == -1) {
            report_error("read");
        }

        return result;
    }

private:
    void report_error(const std::string& context) {
        throw std::runtime_error(context + ": " + std::string(strerror(errno)));
    }

};



<|start_filename|>include/fastscancount_avx2.h<|end_filename|>
#ifndef FASTSCANCOUNT_AVX2_H
#define FASTSCANCOUNT_AVX2_H

// this code expects an x64 processor with AVX2

#ifdef _MSC_VER
#include <intrin.h>
#else
#include <x86intrin.h>
#endif

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <vector>

namespace fastscancount {
namespace {
// credit: implementation and design by <NAME>
static inline size_t find_next_gt(uint8_t *array, const size_t size,
                                  const uint8_t threshold) {
  size_t vsize = size / 32;
  __m256i *varray = (__m256i *)array;
  const __m256i comprand = _mm256_set1_epi8(threshold);
  int bits = 0;

  for (size_t i = 0; i < vsize; i++) {
    __m256i v = _mm256_loadu_si256(varray + i);
    __m256i cmp = _mm256_cmpgt_epi8(v, comprand);
    if ((bits = _mm256_movemask_epi8(cmp))) {
      return i * 32 + __builtin_ctz(bits);
    }
  }

  // tail handling
  for (size_t i = vsize * 32; i < size; i++) {
    auto v = array[i];
    if (v > threshold)
      return i;
  }

  return SIZE_MAX;
}

void populate_hits_avx(std::vector<uint8_t> &counters, size_t range,
                       size_t threshold, size_t start,
                       std::vector<uint32_t> &out) {
  uint8_t *array = counters.data();

  size_t ro = range;
  while (true) {
    size_t next = find_next_gt(array, range, (uint8_t)threshold);
    if (next == SIZE_MAX)
      break;
    out.push_back(start + next);
    range -= (next + 1);
    array += (next + 1);
    start += (next + 1);
  }
}

void update_counters(const uint32_t *&it_, uint8_t *counters,
                     uint32_t range_end) {
  const uint32_t *it = it_;
  for (uint32_t e; (e = *it) < range_end; ++it) {
    counters[e]++;
  }
  it_ = it;
}

void update_counters_final(const uint32_t *&it_, const uint32_t *end,
                           uint8_t *counters) {
  uint64_t e;
  const uint32_t *it = it_;
  for (; it != end; it++) {
    counters[*it]++;
  }
  it_ = end;
}
} // namespace

void fastscancount_avx2(const std::vector<const std::vector<uint32_t>*> &data,
                        std::vector<uint32_t> &out, uint8_t threshold) {
  const size_t cache_size = 40000;
  std::vector<uint8_t> counters(cache_size);
  out.clear();
  const size_t dsize = data.size();

  struct data_info {
    const uint32_t *cur; // current pointer into data
    const uint32_t *end; // pointer to end
    uint32_t last;       // value of last element
    data_info(const uint32_t *cur, const uint32_t *end, uint32_t last)
        : cur{cur}, end{end}, last{last} {}
  };

  std::vector<data_info> iter_data;
  iter_data.reserve(dsize);
  for (auto &d : data) {
    iter_data.emplace_back(d->data(), d->data() + d->size(), d->back());
  }

  uint32_t largest = 0;
  for (size_t c = 0; c < data.size(); c++) {
    if (largest < (*data[c])[data[c]->size() - 1])
      largest = (*data[c])[data[c]->size() - 1];
  }
  auto cdata = counters.data();
  for (uint32_t start = 0; start < largest; start += cache_size) {
    memset(cdata, 0, cache_size * sizeof(counters[0]));
    for (auto &id : iter_data) {
      // determine if the loop will end because we get to the end of
      // data, or because we get to the end of the range
      if (__builtin_expect(id.last >= start + cache_size, 1)) {
        // the iteration is guaranteed to end because an element becomes >=
        // range_end, so we don't need to check for end of data
        update_counters(id.cur, cdata - start, start + cache_size);
      } else {
        // the iteration is guaranteed to end because we get to the end of the
        // data
        update_counters_final(id.cur, id.end, cdata - start);
      }
    }

    populate_hits_avx(counters, cache_size, threshold, start, out);
  }
}

} // namespace fastscancount
#endif


<|start_filename|>benchmark/counters.cpp<|end_filename|>
// Fine-grained statistics is available only on Linux
#include "fastscancount.h"
#include "ztimer.h"
#ifdef __AVX2__
#include "fastscancount_avx2.h"
#endif
#ifdef __AVX512F__
#include "fastscancount_avx512.h"
#endif
#include "linux-perf-events-wrapper.h"
#include "maropuparser.h"
#include <algorithm>
#include <cstdint>
#include <cstdio>
#include <immintrin.h>
#include <iostream>
#include <vector>
#include <stdexcept>

#define REPEATS 10
#define RUNNINGTESTS

void scancount(const std::vector<const std::vector<uint32_t>*> &data,
               std::vector<uint32_t> &out, size_t threshold) {
  uint64_t largest = 0;
  for(auto z : data) {
    const std::vector<uint32_t> & v = *z;
    if(v[v.size() - 1] > largest) largest = v[v.size() - 1];
  }
  std::vector<uint8_t> counters(largest+1);
  out.clear();
  for (size_t c = 0; c < data.size(); c++) {
    const std::vector<uint32_t> &v = *data[c];
    for (size_t i = 0; i < v.size(); i++) {
      counters[v[i]]++;
    }
  }
  for (uint32_t i = 0; i < counters.size(); i++) {
    if (counters[i] > threshold)
      out.push_back(i);
  }
}

void calc_boundaries(uint32_t largest, uint32_t range_size, 
                    const std::vector<uint32_t>& data, 
                    std::vector<uint32_t>& range_ends) {
  if (!range_size) {
    throw std::runtime_error("range_size must be > 0");
  }
  uint32_t end = 0;
  range_ends.clear();
  
  for (uint32_t start = 0; start <= largest; start += range_size) {
    uint32_t curr_max = std::min(largest, start + range_size - 1);
    while (end < data.size() && data[end] <= curr_max) {
      end++;
    }
    range_ends.push_back(end);
  }
} 

const uint32_t range_size_avx512 = 40000;

void calc_alldata_boundaries(const std::vector<std::vector<uint32_t>>& data,
                             std::vector<std::vector<uint32_t>>& range_ends,
                             size_t range_size) {
  uint32_t largest = 0;
  range_ends.clear();
  range_ends.resize(data.size());
  for(const auto& v : data) {
    if (!v.empty() && v[v.size() - 1] > largest) largest = v[v.size() - 1];
  }
  for (unsigned i = 0; i < data.size(); ++i) {
    calc_boundaries(largest, range_size, data[i], range_ends[i]); 
  }
}

template <typename F>
void test(F f, const std::vector<const std::vector<uint32_t>*>& data_ptrs,
          std::vector<uint32_t>& answer, unsigned threshold, const std::string &name) {
  scancount(data_ptrs, answer, threshold);
  size_t s1 = answer.size();
  auto a1 (answer);
  std::sort(a1.begin(), a1.end());
  answer.clear();
  f();
  size_t s2 = answer.size();
  auto a2 (answer);
  std::sort(a2.begin(), a2.end());
  if (a1 != a2) {
    std::cout << "s1: " << s1 << " s2: " << s2 << std::endl;
    for(size_t j = 0; j < std::min(s1, s2); j++) {
      std::cout << j << " " << a1[j] << " vs " << a2[j] ;

      if(a1[j] != a2[j]) std::cout << " oh oh ";
      std::cout << std::endl;
    }
    throw std::runtime_error("bug: " + name);
  }
}

template <typename F>
void bench(F f, const std::string &name,
           LinuxEventsWrapper &unified,
           float& elapsed,
           std::vector<uint32_t> &answer, size_t sum, size_t expected,
           bool print) {
  WallClockTimer tm;
  unified.start();
  f();
  unified.end();
  elapsed += tm.split();
  if (answer.size() != expected)
    std::cerr << "bug: expected " << expected << " but got " << answer.size()
              << "\n";
#ifdef __linux__
  if (print) {
    double cycles = unified.get_result(PERF_COUNT_HW_CPU_CYCLES);
    double instructions = unified.get_result(PERF_COUNT_HW_INSTRUCTIONS);
    double misses = unified.get_result(PERF_COUNT_HW_BRANCH_MISSES);
    std::cout << name << std::endl;
    std::cout << cycles / sum << " cycles/element " << std::endl;
    std::cout << instructions / cycles << " instructions/cycles " << std::endl;
    std::cout << misses / sum << " miss/element " << std::endl;
  }
#endif
}

void demo_data(const std::vector<std::vector<uint32_t>>& data,
              const std::vector<std::vector<uint32_t>>& queries,
              size_t threshold) {
  size_t N = 0;
  for (const auto& data_elem : data) {
    size_t sz = data_elem.size();
    if (sz) {
      N = std::max(N, (size_t)data_elem[sz-1] + 1);
    }
  }

  std::vector<uint32_t> answer;
  answer.reserve(N);

  std::vector<int> evts = {
#ifdef __linux__
                           PERF_COUNT_HW_CPU_CYCLES,
                           PERF_COUNT_HW_INSTRUCTIONS,
                           PERF_COUNT_HW_BRANCH_MISSES,
                           PERF_COUNT_HW_CACHE_REFERENCES,
                           PERF_COUNT_HW_CACHE_MISSES
#endif
                          };
  LinuxEventsWrapper unified(evts);

  std::vector<std::vector<uint32_t>> range_boundaries;
  calc_alldata_boundaries(data, range_boundaries, range_size_avx512);

  std::vector<const std::vector<uint32_t>*> data_ptrs;
  std::vector<const std::vector<uint32_t>*> range_ptrs;

  float elapsed = 0, elapsed_fast = 0, elapsed_avx = 0, elapsed_avx512 = 0;

  size_t sum_total = 0;

  for (size_t qid = 0; qid < queries.size(); ++qid) {
    const auto& query_elem = queries[qid];
    data_ptrs.clear();
    range_ptrs.clear();
    size_t sum = 0;
    for (uint32_t idx : query_elem) {
      if (idx >= data.size()) {
        std::stringstream err;
        err << "Inconsistent data, posting " << idx << 
               " is >= # of postings " << data.size() << " query id " << qid;
        throw std::runtime_error(err.str());
      }
      sum += data[idx].size();
      data_ptrs.push_back(&data[idx]);
      range_ptrs.push_back(&range_boundaries[idx]);
    }
    sum_total += sum;

    scancount(data_ptrs, answer, threshold);
    const size_t expected = answer.size();

#ifdef RUNNINGTESTS
    test(
      [&](){
        fastscancount::fastscancount(data_ptrs, answer, threshold);
      }, data_ptrs, answer, threshold, "fastscancount"
    );
#ifdef __AVX2__
    test(
      [&](){
        fastscancount::fastscancount_avx2(data_ptrs, answer, threshold);
      }, data_ptrs, answer, threshold, "fastscancount_avx2"
    );
#endif

#ifdef __AVX512F__
    test(
      [&](){
        fastscancount::fastscancount_avx512(range_size_avx512, data_ptrs, range_ptrs, answer, threshold);
      }, data_ptrs, answer, threshold, "fastscancount_avx512"
    );
#endif

#endif
    std::cout << "Qid: " << qid << " got " << expected << " hits\n";

    bool last = (qid == queries.size() - 1);

    bench(
        [&]() {
          scancount(data_ptrs, answer, threshold);
        },
        "baseline scancount", unified, elapsed, answer, sum,
        expected, last);

    bench(
        [&]() {
          fastscancount::fastscancount(data_ptrs, answer, threshold);
        },
        "optimized cache-sensitive scancount", unified, elapsed_fast, answer, sum,
        expected, last);
#ifdef __AVX2__
    bench(
        [&]() {
          fastscancount::fastscancount_avx2(data_ptrs, answer, threshold);
        },
        "AVX2-based scancount", unified, elapsed_avx, answer, sum, expected, last);
#endif
#ifdef __AVX512F__
    bench(
        [&]() {
          fastscancount::fastscancount_avx512(range_size_avx512, data_ptrs, range_ptrs, answer, threshold);
        },
        "AVX512-based scancount", unified, elapsed_avx512, answer, sum, expected, last);
#endif
  }
  std::cout << "Elems per millisecond:" << std::endl;
  std::cout << "scancount: " << (sum_total/(elapsed/1e3)) << std::endl; 
  std::cout << "fastscancount: " << (sum_total/(elapsed_fast/1e3)) << std::endl; 
#ifdef __AVX2__
  std::cout << "fastscancount_avx2: " << (sum_total/(elapsed_avx/1e3)) << std::endl; 
#endif
#ifdef __AVX512F__
  std::cout << "fastscancount_avx512: " << (sum_total/(elapsed_avx512/1e3)) << std::endl; 
#endif

}


void demo_random(size_t N, size_t length, size_t array_count, size_t threshold) {
  std::vector<std::vector<uint32_t>> data(array_count);

  std::vector<const std::vector<uint32_t>*> data_ptrs;
  std::vector<uint32_t> answer;
  answer.reserve(N);

  size_t sum = 0;
  for (size_t c = 0; c < array_count; c++) {
    std::vector<uint32_t> &v = data[c];
    for (size_t i = 0; i < length; i++) {
      v.push_back(rand() % N);
    }
    std::sort(v.begin(), v.end());
    v.resize(std::distance(v.begin(), unique(v.begin(), v.end())));
    sum += v.size();
    data_ptrs.push_back(&data[c]);
  }


  std::vector<std::vector<uint32_t>> range_boundaries;
  calc_alldata_boundaries(data, range_boundaries, range_size_avx512);
  std::vector<const std::vector<uint32_t>*> range_ptrs;
  for (size_t c = 0; c < array_count; c++) {
    range_ptrs.push_back(&range_boundaries[c]);
  }

  std::vector<int> evts = {
#ifdef __linux__
                           PERF_COUNT_HW_CPU_CYCLES,
                           PERF_COUNT_HW_INSTRUCTIONS,
                           PERF_COUNT_HW_BRANCH_MISSES,
                           PERF_COUNT_HW_CACHE_REFERENCES,
                           PERF_COUNT_HW_CACHE_MISSES
#endif
                          };
  LinuxEventsWrapper unified(evts);
  float elapsed = 0, elapsed_fast = 0, elapsed_avx = 0, elapsed_avx512 = 0;
  scancount(data_ptrs, answer, threshold);
  const size_t expected = answer.size();
  std::cout << "Got " << expected << " hits\n";
  size_t sum_total = sum * REPEATS;
  for (size_t t = 0; t < REPEATS; t++) {
    bool last = (t == REPEATS - 1);

    bench(
        [&]() {
          scancount(data_ptrs, answer, threshold);
        },
        "baseline scancount", unified, elapsed, answer, sum,
        expected, last);
  }

  for (size_t t = 0; t < REPEATS; t++) {
    bool last = (t == REPEATS - 1);

#ifdef RUNNINGTESTS
    test(
      [&](){
        fastscancount::fastscancount(data_ptrs, answer, threshold);
      }, data_ptrs, answer, threshold, "fastscancount"
    );
#endif

    bench(
        [&]() {
          fastscancount::fastscancount(data_ptrs, answer, threshold);
        },
        "optimized cache-sensitive scancount", unified, elapsed_fast, answer, sum,
        expected, last);
  }

  for (size_t t = 0; t < REPEATS; t++) {
    bool last = (t == REPEATS - 1);

#ifdef __AVX2__
#ifdef RUNNINGTESTS
    test(
      [&](){
        fastscancount::fastscancount_avx2(data_ptrs, answer, threshold);
      }, data_ptrs, answer, threshold, "fastscancount_avx2"
    );
#endif
    bench(
        [&]() {
          fastscancount::fastscancount_avx2(data_ptrs, answer, threshold);
        },
        "AVX2-based scancount", unified, elapsed_avx, answer, sum, expected, last);
#endif
  }

  for (size_t t = 0; t < REPEATS; t++) {
    bool last = (t == REPEATS - 1);
#ifdef __AVX512F__
#ifdef RUNNINGTESTS
    test(
      [&](){
        fastscancount::fastscancount_avx512(range_size_avx512, data_ptrs, range_ptrs, answer, threshold);
      }, data_ptrs, answer, threshold, "fastscancount_avx512"
    );
#endif

    bench(
        [&]() {
          fastscancount::fastscancount_avx512(range_size_avx512, data_ptrs, range_ptrs, answer, threshold);
        },
        "AVX512-based scancount", unified, elapsed_avx512, answer, sum, expected, last);
#endif
  }

  std::cout << "Elems per millisecond:" << std::endl;
  std::cout << "scancount: " << (sum_total/(elapsed/1e3)) << std::endl; 
  std::cout << "fastscancount: " << (sum_total/(elapsed_fast/1e3)) << std::endl; 
#ifdef __AVX2__
  std::cout << "fastscancount_avx2: " << (sum_total/(elapsed_avx/1e3)) << std::endl; 
#endif
#ifdef __AVX512F__
  std::cout << "fastscancount_avx512: " << (sum_total/(elapsed_avx512/1e3)) << std::endl; 
#endif
}

void usage(const std::string& err="") {
  if (!err.empty()) {
    std::cerr << err << std::endl;
  }
  std::cerr << "usage: --postings <postings file> --queries <queries file> --threshold <threshold>" << std::endl;
}

int main(int argc, char *argv[]) {
  // A very naive way to process arguments, 
  // but it's ok unless we need to extend it substantially.
  if (argc != 1) {
    if (argc != 7) {
      usage("");
      return EXIT_FAILURE; 
    }
    std::string postings_file, queries_file;
    int threshold = -1;
    for (int i = 1; i < argc; ++i) {
      if (std::string(argv[i]) == "--postings") {
        postings_file = argv[++i];
      } else if (std::string(argv[i]) == "--queries") {
        queries_file = argv[++i];
      } else if (std::string(argv[i]) == "--threshold") {
        threshold = std::atoi(argv[++i]);
      }
    }
    if (postings_file.empty() || queries_file.empty() || threshold < 0) {
      usage("Specify queries, postings, and the threshold!");
      return EXIT_FAILURE; 
    }
    std::vector<uint32_t> tmp; 
    std::vector<std::vector<uint32_t>> data;
    {
      MaropuGapReader drdr(postings_file);
      if (!drdr.open()) {
        usage("Cannot open: " + postings_file);
        return EXIT_FAILURE; 
      }
      while (drdr.loadIntegers(tmp)) {
        data.push_back(tmp);
      }
    }
    std::vector<std::vector<uint32_t>> queries;
    {
      MaropuGapReader qrdr(queries_file);
      if (!qrdr.open()) {
        usage("Cannot open: " + queries_file);
        return EXIT_FAILURE; 
      }
      while (qrdr.loadIntegers(tmp)) {
        queries.push_back(tmp);
      }
    }
              
    try { 
      demo_data(data, queries, threshold);
    } catch (const std::exception& e) {
      std::cerr << "Exception: " << e.what() << std::endl;
      return EXIT_FAILURE;
    }
  } else {
    try {
      // Previous demo with threshold 3
      //demo_random(20000000, 50000, 100, 3);
      for (unsigned k = 1; k < 10; ++k) {
        std::cout << "Demo threshold:" << k << std::endl;
        demo_random(20000000, 50000, 100, k);
        std::cout << "=======================" << std::endl;
      }
    } catch (const std::exception& e) {
      std::cerr << "Exception: " << e.what() << std::endl;
      return EXIT_FAILURE;
    }
  }
  return EXIT_SUCCESS;
}


<|start_filename|>include/fastscancount.h<|end_filename|>
#ifndef FASTSCANCOUNT_H
#define FASTSCANCOUNT_H

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <vector>

// credit: implementation and design by <NAME> and <NAME>

namespace fastscancount {

namespace {

// used by natefastscancount
uint32_t *natefastscancount_maincheck(uint8_t *counters, size_t &it,
                                      const uint32_t *d, size_t start,
                                      size_t range, uint8_t threshold,
                                      uint32_t *out) {
  range += start;
  counters -= start;
  size_t i = it;
  for (uint32_t val = d[i]; val < range; val = d[++i]) {
    uint8_t c = counters[val];
    if (c == threshold) *out++ = val;
    counters[val] = c + 1;
  }
  it = i;
  return out;
}

// used by natefastscancount
uint32_t *natefastscancount_finalcheck(uint8_t *counters, size_t &it,
                                       const uint32_t *d, size_t start,
                                       size_t itend, uint8_t threshold,
                                       uint32_t *out) {
  uint8_t *const deccounters = counters - start;
  size_t i = it;
  for (; i < itend; i++) {
    uint32_t val = d[i];
    uint8_t *location = deccounters + val;
    uint8_t c = *location;
    if (c == threshold) {
      *out++ = val;
    }
    *location = c + 1;
  }
  it = i;
  return out;
}
} // namespace

void fastscancount(const std::vector<const std::vector<uint32_t>*> &data,
                   std::vector<uint32_t> &out, uint8_t threshold) {
  size_t cache_size = 65536;
  size_t range = cache_size;
  std::vector<uint8_t> counters(cache_size);
  size_t ds = data.size();
  out.resize( 4 * range); // let us add lots of capacity
  uint32_t *output = out.data();
  uint32_t *initout = out.data();
  std::vector<size_t> iters(ds);
  size_t countsofar = 0;
  uint32_t largest = 0;
  for (size_t c = 0; c < ds; c++) {
    if (largest < (*data[c])[data[c]->size() - 1])
      largest = (*data[c])[data[c]->size() - 1];
  }
  // we are assuming that all vectors in data are non-empty
  for (size_t start = 0; start < largest; start += range) {
    // make sure that the capacity is sufficient
    countsofar = output - initout;
    if (out.size() - countsofar < range) {
      out.resize(out.size() + 4 * range);
      initout = out.data();
      output = out.data() + countsofar;
    }
    memset(counters.data(), 0, range);
    for (size_t c = 0; c < ds; c++) {
      size_t it = iters[c]; // recover where we were
      const std::vector<uint32_t> &d = *data[c];
      const size_t itend = d.size();
      if (it == itend) // check that there is data to be processed
        continue;      // exhausted
      // check if we need to be careful:
      bool near_the_end = (d[itend - 1] < start + range);
      if (near_the_end) {
        output = natefastscancount_finalcheck(counters.data(), it, d.data(),
                                              start, itend, threshold, output);
      } else {
        output = natefastscancount_maincheck(counters.data(), it, d.data(),
                                             start, range, threshold, output);
      }
      iters[c] = it; // store it for next round
    }
  }
  countsofar = output - initout;
  out.resize(countsofar);
}
} // namespace fastscancount

#endif


<|start_filename|>Makefile<|end_filename|>
OPT := -O3
# Leo really doubts -mavx2 helps anything, but one can
# disable avx512 tests by enforcing -mavx2
#CXXFLAGS := -std=c++17 $(OPT) -mavx2
CXXFLAGS := -std=c++17 $(OPT) -march=native

counter: benchmark/counters.cpp include/*.h Makefile
	$(CXX) $(CXXFLAGS) $(CXXEXTRA) -o counter benchmark/counters.cpp -Ibenchmark -Iinclude

clean:
	rm -f counter

