math.hpp

#pragma once
 
#include <atomic>
#include <chrono>
#include <concepts>
#include <numeric>
 
namespace tf {

template <typename T>
requires (std::is_unsigned_v<std::decay_t<T>> && sizeof(T) == 8)
constexpr T next_pow2(T x) {
  if(x == 0) return 1;
  x--;
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  x |= x >> 32;
  x++;
  return x;
}

template <typename T>
requires (std::is_unsigned_v<std::decay_t<T>> && sizeof(T) == 4)
constexpr T next_pow2(T y) {
  if(y == 0) return 1;
  y--;
  y |= y >> 1;
  y |= y >> 2;
  y |= y >> 4;
  y |= y >> 8;
  y |= y >> 16;
  y++;
  return y;
}

template <std::integral T>
constexpr bool is_pow2(const T& x) {
  return x && (!(x&(x-1)));
}

template <size_t N>
constexpr size_t static_floor_log2() {
  return (N < 2) ? 0 : 1 + static_floor_log2<N / 2>();
  //auto log = 0;
  //while (N >>= 1) {
  //  ++log;
  //}
  //return log;
}
 

template <typename RandItr, typename C>
RandItr median_of_three(RandItr l, RandItr m, RandItr r, C cmp) {
  return cmp(*l, *m) ? (cmp(*m, *r) ? m : (cmp(*l, *r) ? r : l ))
                     : (cmp(*r, *m) ? m : (cmp(*r, *l) ? r : l ));
}

template <typename RandItr, typename C>
RandItr pseudo_median_of_nine(RandItr beg, RandItr end, C cmp) {
  size_t N = std::distance(beg, end);
  size_t offset = N >> 3;
  return median_of_three(
    median_of_three(beg, beg+offset, beg+(offset*2), cmp),
    median_of_three(beg+(offset*3), beg+(offset*4), beg+(offset*5), cmp),
    median_of_three(beg+(offset*6), beg+(offset*7), end-1, cmp),
    cmp
  );
}

template<typename Iter, typename Compare>
void sort2(Iter a, Iter b, Compare comp) {
  if (comp(*b, *a)) std::iter_swap(a, b);
}

template<typename Iter, typename Compare>
void sort3(Iter a, Iter b, Iter c, Compare comp) {
  sort2(a, b, comp);
  sort2(b, c, comp);
  sort2(a, b, comp);
}

template <std::integral T>
T unique_id() {
  static std::atomic<T> counter{0};
  return counter.fetch_add(1, std::memory_order_relaxed);
}

template <typename T>
inline void atomic_max(std::atomic<T>& v, const T& max_v) noexcept {
  T prev = v.load(std::memory_order_relaxed);
  while(prev < max_v && 
        !v.compare_exchange_weak(prev, max_v, std::memory_order_relaxed,
                                              std::memory_order_relaxed)) {
  }
}

template <typename T>
inline void atomic_min(std::atomic<T>& v, const T& min_v) noexcept {
  T prev = v.load(std::memory_order_relaxed);
  while(prev > min_v && 
        !v.compare_exchange_weak(prev, min_v, std::memory_order_relaxed,
                                              std::memory_order_relaxed)) {
  }
}

template <typename T>
inline T seed() noexcept {
  return std::chrono::system_clock::now().time_since_epoch().count();
}
 
// ------------------------------------------------------------------------------------------------
// coprime
// ------------------------------------------------------------------------------------------------

constexpr size_t coprime(size_t N) {
  if(N < 3) {
    return 1;
  }
  for (size_t x = N; --x > 0;) {
    if (std::gcd(x, N) == 1) {
      return x;
    }
  }
  return 1;
}

template <size_t N>
constexpr std::array<size_t, N> make_coprime_lut() {
  static_assert(N>0, "N must be greater than 0");
  std::array<size_t, N> coprimes{};
  for (size_t n = 0; n < N; ++n) {
    coprimes[n] = coprime(n);
  }
  return coprimes;
}
 
//template <typename T>
//constexpr T lemire_range(T x, T range) {
//  return (uint32_t)(((uint64_t)x * (uint64_t)range) >> 32);
//}
 

template <typename T>
class Xorshift {
  static_assert(std::is_unsigned<T>::value, "Xorshift requires an unsigned integral type.");
 
  public:

  Xorshift() = default;

  Xorshift(T value) : _state(value) {}

  void seed(T value) {
    _state = value;
  }
  
  T operator()() {
    if constexpr (sizeof(T) == 8) {
      // Xorshift64 constants
      _state ^= _state << 13;
      _state ^= _state >> 7;
      _state ^= _state << 17;
      return _state * 0x2545F4914F6CDD1DULL;
    }
    else if constexpr (sizeof(T) == 4) {
      // Xorshift32 constants
      _state ^= _state << 13;
      _state ^= _state >> 17;
      _state ^= _state << 5;
      return _state;
    }
    else {
      static_assert(sizeof(T) == 0, "Unsupported bit-width for Xorshift. Use uint32_t or uint64_t.");
    }
  }
 
  private:
 
  T _state;  // must be initialized with non-zero
};
 
}  // end of namespace tf -----------------------------------------------------