pe_parallel_algo

#ifndef PE_PARALLEL_ALGO_
#define PE_PARALLEL_ALGO_

#include "pe_base"
#include "pe_int128"
#include "pe_type_traits"
#include "pe_mod"

namespace pe {
SL void ParallelExecute(const std::vector<std::function<void()>>& functions,
                        int TN = 8) {
  const int size = static_cast<int>(std::size(functions));
#if ENABLE_OPENMP
#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
#endif
  for (int i = 0; i < size; ++i) {
    functions[i]();
  }
}

template <int TN, typename T>
SL void ParallelSort(T* s, T* e) {
  static_assert(TN > 0, "TN > 0");

  const int64 n = e - s;
  const int64 each = n / TN;

  int64 pos[TN + 1];
  for (int i = 0; i < TN; ++i) pos[i] = i * each;
  pos[TN] = n;

#if ENABLE_OPENMP
#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
#endif
  for (int i = 0; i < TN; ++i) {
    std::sort(s + pos[i], s + pos[i + 1]);
  }

  for (int remain = TN; remain > 1;) {
#if ENABLE_OPENMP
#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
#endif
    for (int i = 0; i <= remain - 2; i += 2) {
      std::inplace_merge(s + pos[i], s + pos[i + 1], s + pos[i + 2]);
    }

    int top = 0;
    for (int i = 0; i <= remain; ++i) {
      if (!(i % 2 == 1 && i + 1 <= remain)) {
        pos[top++] = pos[i];
      }
    }
    remain = top - 1;
  }
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    FindFirst(T first, T last, const std::function<bool(T)>& f) {
  const T end = last + 1;
  if (first > last) return end;
  while (first <= last && !f(first)) ++first;
  return first;
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    FindFirst(T first, const std::function<bool(T)>& f) {
  while (!f(first)) ++first;
  return first;
}

template <typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    FindFirst(T first, const std::function<T(T, T)>& f) {
  for (T j = first;; j += B) {
    const T y = j + B - 1;
    T x = f(j, y);
    if (x <= y) {
      return x;
    }
  }
}

// |block  0|block  1|...|block  TN-1|block TN|block TN+1|...
// |thread 0|thread 1|...|thread TN-1|thead  0|thread   1|...
template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindFirst(T first, T last, const std::function<T(T, T)>& f) {
#if !ENABLE_OPENMP
  return f(first, last);
#else
  const T end = last + 1;
  if (first > last) return end;
  if (TN <= 1 || last - first + 1 <= B) return f(first, last);

  omp_lock_t locker;
  omp_init_lock(&locker);

  const T STEP = static_cast<T>(TN) * B;
  T best = last + 1;

#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
  for (int id = 0; id < TN; ++id) {
    for (T j = first + static_cast<T>(id) * B; j <= last; j += STEP) {
      const T y = std::min(last, j + B - 1);
      T x = f(j, y);
      if (x <= y) {
        omp_set_lock(&locker);
        if (x < best) best = x;
        omp_unset_lock(&locker);
        break;
      } else {
        int should_exit = 0;
        omp_set_lock(&locker);
        if (best < j) should_exit = 1;
        omp_unset_lock(&locker);
        if (should_exit) break;
      }
    }
  }

  omp_destroy_lock(&locker);

  return best;

#endif
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindFirst(T first, T last, const std::function<bool(T)>& f) {
  auto g = [&](T a, T b) -> T {
    while (a <= b && !f(a)) ++a;
    return a;
  };
  return ParallelFindFirst<TN, T, B>(first, last, g);
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindFirst(T first, const std::function<T(T, T)>& f) {
#if !ENABLE_OPENMP
  return FindFirst<T, B>(first, f);
#else
  if (TN <= 1) {
    return FindFirst<T, B>(first, f);
  }
  omp_lock_t locker;
  omp_init_lock(&locker);

  const T guard = first - 1;
  const T STEP = static_cast<T>(TN) * B;
  T best = guard;

#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
  for (int id = 0; id < TN; ++id) {
    for (T j = first + static_cast<T>(id) * B;; j += STEP) {
      const T y = j + B - 1;
      T x = f(j, y);
      if (x <= y) {
        omp_set_lock(&locker);
        if (best == guard || x < best) best = x;
        omp_unset_lock(&locker);
        break;
      } else {
        int should_exit = 0;
        omp_set_lock(&locker);
        if (best != guard && best < j) should_exit = 1;
        omp_unset_lock(&locker);
        if (should_exit) break;
      }
    }
  }

  omp_destroy_lock(&locker);

  return best;

#endif
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindFirst(T first, const std::function<bool(T)>& f) {
  auto g = [&](T a, T b) -> T {
    while (a <= b && !f(a)) ++a;
    return a;
  };
  return ParallelFindFirst<TN, T, B>(first, g);
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    FindLast(T first, T last, const std::function<bool(T)>& f) {
  const T end = first - 1;
  if (first > last) return end;
  while (last >= first && !f(last)) --last;
  return last;
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    FindLast(T last, const std::function<bool(T)>& f) {
  while (!f(last)) --last;
  return last;
}

template <typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    FindLast(T last, const std::function<T(T, T)>& f) {
  for (T j = last;; j -= B) {
    const T y = j - B + 1;
    T x = f(y, j);
    if (x >= y) {
      return x;
    }
  }
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindLast(T first, T last, const std::function<T(T, T)>& f) {
#if !ENABLE_OPENMP
  return f(first, last);
#else
  const T end = first - 1;
  if (first > last) return end;
  if (TN <= 1 || last - first + 1 <= B) return f(first, last);

  omp_lock_t locker;
  omp_init_lock(&locker);

  const T STEP = static_cast<T>(TN) * B;
  T best = first - 1;

#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
  for (int id = 0; id < TN; ++id) {
    for (T j = last - static_cast<T>(id) * B; j >= first; j -= STEP) {
      const T y = std::max(first, j - B + 1);
      T x = f(y, j);
      if (x >= y) {
        omp_set_lock(&locker);
        if (x > best) best = x;
        omp_unset_lock(&locker);
        break;
      } else {
        int should_exit = 0;
        omp_set_lock(&locker);
        if (best > j) should_exit = 1;
        omp_unset_lock(&locker);
        if (should_exit) break;
      }
    }
  }

  omp_destroy_lock(&locker);

  return best;

#endif
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindLast(T first, T last, const std::function<bool(T)>& f) {
  auto g = [&](T a, T b) -> T {
    while (b >= a && !f(b)) --b;
    return b;
  };
  return ParallelFindLast<TN, T, B>(first, last, g);
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindLast(T last, const std::function<T(T, T)>& f) {
#if !ENABLE_OPENMP
  return FindLast(last, f);
#else
  if (TN <= 1) {
    return FindLast(last, f);
  }
  omp_lock_t locker;
  omp_init_lock(&locker);

  const T guard = last + 1;
  const T STEP = static_cast<T>(TN) * B;
  T best = guard;

#pragma omp parallel for schedule(dynamic, 1) num_threads(TN)
  for (int id = 0; id < TN; ++id) {
    for (T j = last - static_cast<T>(id) * B;; j -= STEP) {
      const T y = j - B + 1;
      T x = f(j, y);
      if (x >= y) {
        omp_set_lock(&locker);
        if (best == guard || x > best) best = x;
        omp_unset_lock(&locker);
        break;
      } else {
        int should_exit = 0;
        omp_set_lock(&locker);
        if (best != guard && best > j) should_exit = 1;
        omp_unset_lock(&locker);
        if (should_exit) break;
      }
    }
  }

  omp_destroy_lock(&locker);

  return best;

#endif
}

template <int TN, typename T, int B = 10000>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T)
    ParallelFindLast(T last, const std::function<bool(T)>& f) {
  auto g = [&](T a, T b) -> T {
    while (b >= a && !f(b)) --b;
    return b;
  };
  return ParallelFindLast<TN, T, B>(last, g);
}

#if ENABLE_OPENMP
template <typename T>
struct PSum {
  PSum() { Reset(); }

  PSum& Reset() {
    std::fill(data, data + 128, 0);
    return *this;
  }

  PSum& operator+=(T v) {
    data[omp_get_thread_num()] += ExtractValue(v);
    return *this;
  }

  PSum& Add(T v) {
    *this += v;
    return *this;
  }

  PSum& operator++() { return Add(1); }

  PSum& operator++(int) { return Add(1); }

  PSum& operator-=(T v) {
    data[omp_get_thread_num()] -= ExtractValue(v);
    return *this;
  }

  PSum& Sub(T v) {
    *this -= v;
    return *this;
  }

  PSum& operator--() { return Sub(1); }

  PSum& operator--(int) { return Sub(1); }

  T CalSum() const {
    T r = 0;
    for (int i = 0; i < 128; ++i) {
      r += data[i];
    }
    return r;
  }

  T Cal() const { return CalSum(); }
  T value() const { return CalSum(); }

  operator T() const { return CalSum(); }
  T operator()() const { return CalSum(); }

  // Note: false sharing issue
  T data[128];
};

template <typename T>
struct PSumMod {
  PSumMod(T mod) : mod(mod) { Reset(); }

  PSumMod& Reset() {
    std::fill(data, data + 128, 0);
    return *this;
  }

  PSumMod& operator+=(T v) {
    const int id = omp_get_thread_num();
    data[id] = AddMod<T>(data[id], ExtractValue(v), mod);
    return *this;
  }

  PSumMod& Add(T v) {
    *this += v;
    return *this;
  }

  PSumMod& operator++() { return Add(1); }

  PSumMod& operator++(int) { return Add(1); }

  PSumMod& operator-=(T v) {
    const int id = omp_get_thread_num();
    data[id] = SubMod<T>(data[id], ExtractValue(v), mod);
    return *this;
  }

  PSumMod& Sub(T v) {
    *this -= v;
    return *this;
  }

  PSumMod& operator--() { return Sub(1); }

  PSumMod& operator--(int) { return Sub(1); }

  T CalSum() const {
    T r = 0;
    for (int i = 0; i < 128; ++i) {
      r = AddMod(r, data[i], mod);
    }
    return r;
  }

  T Cal() const { return CalSum(); }
  T value() const { return CalSum(); }

  operator T() const { return CalSum(); }
  T operator()() const { return CalSum(); }

  // Note: false sharing issue
  T data[128];
  const T mod;
};

template <typename T, typename C = std::less<T>>
struct PMin {
  PMin() {
    omp_init_lock(&locker);
    enable_log = true;
    invert_cmp = false;
    Reset();
  }

  ~PMin() { omp_destroy_lock(&locker); }

  PMin& Reset() {
    for (int i = 0; i < 128; ++i) {
      data[i].reset();
    }
    result.reset();
    return *this;
  }

  PMin& SetEnableLog(bool enable_log = true) {
    this->enable_log = enable_log;
    return *this;
  }

  PMin& SetInvertCmp(bool invert_cmp = true) {
    this->invert_cmp = invert_cmp;
    return *this;
  }

  PMin& CheckMin(const T& v) {
    const int id = omp_get_thread_num();
    if (!data[id].has_value() || Compare(v, data[id].value())) {
      data[id].emplace(v);
      CheckGlobalMin(data[id]);
    }
    return *this;
  }

  PMin& CheckGlobalMin(std::optional<T> o) {
    omp_set_lock(&locker);
    if (!result.has_value() || Compare(o.value(), result.value())) {
      result = o;
      if (enable_log) {
        std::cerr << "PMin: " << result.value() << std::endl;
      }
    }
    omp_unset_lock(&locker);
    return *this;
  }

  int Compare(const T& a, const T& b) {
    int t = cmp(a, b);
    return invert_cmp ? !t : t;
  }

  T value() const { return result.value(); }

  operator T() const { return result.value(); }
  T operator()() const { return result.value(); }

  // Note: false sharing issue
  std::optional<T> data[128];
  std::optional<T> result;
  omp_lock_t locker;
  bool enable_log;
  bool invert_cmp;
  C cmp;
};
#endif
}  // namespace pe
#endif