hyper fibonacci sequence

#include <cassert>

#include <iostream>

#include <type_traits>

#include <vector>

​

template <auto& Mod>

struct modint {

public:

  using value_type = std::decay_t<decltype(Mod)>;

​

  modint() : val() {}

  modint(const value_type v) : val(v) { clamp_large(); }

​

  [[nodiscard]] modint operator+(const modint rhs) const {

    return modint(val + rhs.val, true);

  }

  [[nodiscard]] modint operator-(const modint rhs) const {

    return modint(val - rhs.val, true);

  }

  [[nodiscard]] modint operator*(const modint rhs) const {

    return modint(val * rhs.val);

  }

​

  modint& operator=(const modint rhs) {

    val = rhs.val;

    return *this;

  }

  modint& operator+=(const modint rhs) {

    val += rhs.val;

    clamp_small();

    return *this;

  }

  modint& operator-=(const modint rhs) {

    val -= rhs.val;

    clamp_small();

    return *this;

  }

  modint& operator*=(const modint rhs) {

    val *= rhs.val;

    clamp_large();

    return *this;

  }

​

  [[nodiscard]] modint pow(int index) const {

    modint res = 1, base = *this;

​

    while (index) {

      if (index & 1) {

        res *= base;

      }

      base *= base;

      index >>= 1;

    }

​

    return res;

  }

​

  [[nodiscard]] operator value_type() const {

    return val;

  }

​

  friend std::istream& operator>>(std::istream& is, modint& x) {

    return is >> x.val;

  }

  friend std::ostream& operator<<(std::ostream& os, const modint x) {

    return os << x.val;

  }

​

private:

  value_type val;

​

  void clamp_small() {

    if (val < 0) {

      val += Mod;

    } else if (val >= Mod) {

      val -= Mod;

    }

  }

​

  void clamp_large() {

    val %= Mod;

  }

​

  modint(const value_type v, bool) : val(v) { clamp_small(); }

};

​

int main() {

  static int M;

  using mint = modint<M>;

​

  mint a, b;

  int c, k;

  std::cin >> M >> a >> b >> c >> k;

  --k;

​

  std::vector index(M, std::vector<int>(M, -1));

​

  std::vector<mint> terms;

  terms.reserve(M * M + 1);

  terms.emplace_back(a);

  terms.emplace_back(b);

​

  int loop_starts_at, loop_length;

  int calculated = 2;

​

  while (true) {

    if (index[terms[calculated - 2]][terms[calculated - 1]] != -1) {

      loop_starts_at = index[terms[calculated - 2]][terms[calculated - 1]] - 2;

      loop_length    = calculated - loop_starts_at - 2;

      terms.pop_back();

      terms.pop_back();

      break;

    } else {

      index[terms[calculated - 2]][terms[calculated - 1]] = calculated;

      terms.emplace_back(terms[calculated - 2].pow(c) + terms[calculated - 1].pow(c));

      ++calculated;

    }

  }

​

  if (k > loop_starts_at) {

    k -= loop_starts_at;

    k %= loop_length;

    k += loop_starts_at;

  }

​

  std::cout << terms[k] << '\n';

}