Language: C++

#ifndef __INTMOD_H__
#define __INTMOD_H__

/* updated: 2019-10-03 */

#include <iostream>
#include <cassert>
#include <vector>
#include <cstdint>
#include <functional>
#include <unordered_map>

template <uint32_t Mod>
class IntMod {
	static_assert(Mod != 0, "Mod must not be 0");
	static_assert(Mod < 0x80000000, "Mod must be less than 0x80000000");

private:
	uint32_t value_m;

public:
	IntMod();
	IntMod(int32_t value);
	IntMod(uint32_t value);
	IntMod(int64_t value);
	IntMod(uint64_t value);

	const IntMod<Mod>& operator+() const;
	IntMod<Mod> operator-() const;
	IntMod<Mod>& operator++();
	IntMod<Mod>& operator--();
	IntMod<Mod> operator++(int dummy);
	IntMod<Mod> operator--(int dummy);
	IntMod<Mod>& operator+=(const IntMod<Mod>& right);
	IntMod<Mod>& operator-=(const IntMod<Mod>& right);
	IntMod<Mod>& operator*=(const IntMod<Mod>& right);
	IntMod<Mod>& operator/=(const IntMod<Mod>& right);
	IntMod<Mod> operator[](uint64_t exp) const;

	static IntMod<Mod> add(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static IntMod<Mod> subtract(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static IntMod<Mod> multiply(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static IntMod<Mod> divide(const IntMod<Mod>& a, const IntMod<Mod>& b);
	static int compare(const IntMod<Mod>& a, const IntMod<Mod>& b);

	IntMod<Mod> inverse() const;
	IntMod<Mod> power(uint64_t exp) const;
	uint32_t get() const;

	static IntMod<Mod> fact(uint32_t num);
	static IntMod<Mod> combi(uint32_t n, uint32_t r);
	static IntMod<Mod> homcombi(uint32_t n, uint32_t r);
	static IntMod<Mod> perm(uint32_t n, uint32_t r);
	static std::vector<IntMod<Mod>> calcInverseList(int max);
};

template <uint32_t Mod>
IntMod<Mod> operator+(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
IntMod<Mod> operator-(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
IntMod<Mod> operator*(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
IntMod<Mod> operator/(const IntMod<Mod>& a, const IntMod<Mod>& b);

template <uint32_t Mod, class Int>
IntMod<Mod> operator+(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator-(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator*(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator/(const IntMod<Mod>& a, const Int& b);

template <uint32_t Mod, class Int>
IntMod<Mod> operator+(const Int& a, const IntMod<Mod>& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator-(const Int& a, const IntMod<Mod>& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator*(const Int& a, const IntMod<Mod>& b);
template <uint32_t Mod, class Int>
IntMod<Mod> operator/(const Int& a, const IntMod<Mod>& b);

template <uint32_t Mod>
bool operator==(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator!=(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator<(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator<=(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator>(const IntMod<Mod>& a, const IntMod<Mod>& b);
template <uint32_t Mod>
bool operator>=(const IntMod<Mod>& a, const IntMod<Mod>& b);

template <uint32_t Mod, class Int>
bool operator==(const IntMod<Mod>& a, const Int& b);
template <uint32_t Mod, class Int>
bool operator!=(const IntMod<Mod>& a, const Int& b);

template <uint32_t Mod>
std::istream& operator>>(std::istream& ist, IntMod<Mod>& object);
template <uint32_t Mod>
std::ostream& operator<<(std::ostream& ost, const IntMod<Mod>& object);

using MInt = IntMod<1000000007>;

// for unordered_set/map
namespace std {
	template <uint32_t Mod>
	struct hash<IntMod<Mod>> {
		size_t operator() (const IntMod<Mod>& object) const {
			return std::hash<uint32_t>()(object.get());
		}
	};
}

template <uint32_t Mod>
uint32_t discreteLog(const IntMod<Mod>& a, const IntMod<Mod>& b);

/*------------------------------------*/
/*-------   Implementations   --------*/
/*------------------------------------*/

template <uint32_t Mod>
inline IntMod<Mod>::IntMod() 
	: value_m(0) {
}

template <uint32_t Mod>
inline IntMod<Mod>::IntMod(int32_t value) {
	int32_t tmp = value % int32_t(Mod);
	value_m = tmp >= 0 ? tmp : Mod - uint32_t(-tmp);
}

template <uint32_t Mod>
inline IntMod<Mod>::IntMod(uint32_t value)
	: value_m(value % Mod) {
}

template <uint32_t Mod>
inline IntMod<Mod>::IntMod(int64_t value) {
	int32_t tmp = value % int64_t(Mod);
	value_m = tmp >= 0 ? tmp : Mod - uint32_t(-tmp);
}

template <uint32_t Mod>
inline IntMod<Mod>::IntMod(uint64_t value)
	: value_m(value % Mod) {
}

template <uint32_t Mod>
inline const IntMod<Mod>& IntMod<Mod>::operator+() const {
	return *this;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::operator-() const {
	IntMod<Mod> ret;
	if (value_m == 0) {
		ret.value_m = 0;
	} else {
		ret.value_m = Mod - value_m;
	}
	return ret;
}

template <uint32_t Mod>
inline IntMod<Mod>& IntMod<Mod>::operator++() {
	++value_m;
	if (value_m == Mod) {
		value_m = 0;
	}
	return *this;
}

template <uint32_t Mod>
inline IntMod<Mod>& IntMod<Mod>::operator--() {
	if (value_m == 0) {
		value_m = Mod;
	}
	--value_m;
	return *this;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::operator++(int dummy) {
	IntMod ret(*this);
	++*this;
	return ret;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::operator--(int dummy) {
	IntMod ret(*this);
	--*this;
	return ret;
}

template <uint32_t Mod>
inline IntMod<Mod>& IntMod<Mod>::operator+=(const IntMod<Mod>& right) {
	value_m += right.value_m;		// note: value_m < 0x80000000
	if (value_m >= Mod) {
		value_m -= Mod;
	}
	return *this;
}

template <uint32_t Mod>
inline IntMod<Mod>& IntMod<Mod>::operator-=(const IntMod<Mod>& right) {
	if (value_m < right.value_m) {
		value_m += Mod;				// note: value_m < 0x80000000
	}
	value_m -= right.value_m;
	return *this;
}

template <uint32_t Mod>
inline IntMod<Mod>& IntMod<Mod>::operator*=(const IntMod<Mod>& right) {
	value_m = (uint64_t(value_m) * right.value_m) % Mod;
	return *this;
}

template <uint32_t Mod>
inline IntMod<Mod>& IntMod<Mod>::operator/=(const IntMod<Mod>& right) {
	*this *= right.inverse();
	return *this;
}

// for power
template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::operator[](uint64_t exp) const {
	return power(exp);
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::add(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret += b;
	return ret;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::subtract(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret -= b;
	return ret;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::multiply(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret *= b;
	return ret;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::divide(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	IntMod<Mod> ret(a);
	ret /= b;
	return ret;
}

template <uint32_t Mod>
inline int IntMod<Mod>::compare(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	if (a.value_m < b.value_m) return -1;
	if (a.value_m > b.value_m) return 1;
	return 0;
}

/* Mod must be a prime. */
template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::inverse() const {
	assert(value_m != 0);
	return power(Mod - 2);
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::power(uint64_t exp) const {
	IntMod<Mod> product(1);
	IntMod<Mod> factor(*this);
	while (exp > 0) {
		if (exp & 1) {
			product *= factor;
		}
		factor *= factor;
		exp >>= 1;
	}
	return product;
}

template <uint32_t Mod>
inline uint32_t IntMod<Mod>::get() const {
	return value_m;
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::fact(uint32_t num) {
	assert(num <= 0x80000000);

	static std::vector<IntMod<Mod>> table(1, 1);
	if (table.size() > num) return table[num];

	int old_size = table.size();
	table.resize(num + 1);
	for (int i = old_size; i <= num; ++i) {
		table[i] = table[i - 1] * i;
	}
	return table[num];
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::combi(uint32_t n, uint32_t r) {
	assert(n >= r);
	return fact(n) / (fact(n - r) * fact(r));
}

template<uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::homcombi(uint32_t n, uint32_t r) {
	assert(n != 0 || r == 0);

	if (n == 0 && r == 0) return IntMod<Mod>(1);
	return combi(n + r - 1, r);
}

template <uint32_t Mod>
inline IntMod<Mod> IntMod<Mod>::perm(uint32_t n, uint32_t r) {
	assert(n >= r);
	return fact(n) / fact(n - r);
}

template <uint32_t Mod>
inline std::vector<IntMod<Mod>> IntMod<Mod>::calcInverseList(int max) {
	assert(max < Mod);

	std::vector<IntMod<Mod>> table(max + 1);
	table[1] = 1;
	for (int i = 2; i <= max; ++i) {
		table[i] = (uint64_t(Mod - Mod / i) * table[Mod % i].Get()) % Mod;
	}
	return table;
}

template <uint32_t Mod>
inline IntMod<Mod> operator+(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::add(a, b);
}

template <uint32_t Mod>
inline IntMod<Mod> operator-(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::subtract(a, b);
}

template <uint32_t Mod>
inline IntMod<Mod> operator*(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::multiply(a, b);
}

template <uint32_t Mod>
inline IntMod<Mod> operator/(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::divide(a, b);
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator+(const IntMod<Mod>& a, const Int& b) {
	return a + IntMod<Mod>(b);
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator-(const IntMod<Mod>& a, const Int& b) {
	return a - IntMod<Mod>(b);
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator*(const IntMod<Mod>& a, const Int& b) {
	return a * IntMod<Mod>(b);
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator/(const IntMod<Mod>& a, const Int& b) {
	return a / IntMod<Mod>(b);
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator+(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) + b;
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator-(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) - b;
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator*(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) * b;
}

template<uint32_t Mod, class Int>
inline IntMod<Mod> operator/(const Int& a, const IntMod<Mod>& b) {
	return IntMod<Mod>(a) / b;
}

template <uint32_t Mod>
inline bool operator==(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) == 0;
}

template <uint32_t Mod>
inline bool operator!=(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) != 0;
}

template <uint32_t Mod>
inline bool operator<(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) < 0;
}

template <uint32_t Mod>
inline bool operator<=(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) <= 0;
}

template <uint32_t Mod>
inline bool operator>(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) > 0;
}

template <uint32_t Mod>
inline bool operator>=(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	return IntMod<Mod>::compare(a, b) >= 0;
}

template <uint32_t Mod, class Int>
bool operator==(const IntMod<Mod>& a, const Int& b) {
	return a == IntMod<Mod>(b);
}

template <uint32_t Mod, class Int>
bool operator!=(const IntMod<Mod>& a, const Int& b) {
	return a != IntMod<Mod>(b);
}

template <uint32_t Mod>
inline std::istream& operator>>(std::istream& ist, IntMod<Mod>& object) {
	uint64_t value;
	ist >> value;
	object = IntMod<Mod>(value);
	return ist;
}

template <uint32_t Mod>
inline std::ostream& operator<<(std::ostream& ost, const IntMod<Mod>& object) {
	ost << object.get();
	return ost;
}

/* 離散対数を求める O(sqrt(Mod)) */
/* 存在すれば a^x == b となる x を、存在しなければ -1 を返す */
template <uint32_t Mod>
uint32_t discreteLog(const IntMod<Mod>& a, const IntMod<Mod>& b) {
	auto sqrtInt = [](uint32_t num) {
		uint32_t lo = 0;
		uint32_t hi = num;
		while (hi - lo > 1) {
			uint64_t mid = (lo + hi) / 2;
			if (mid * mid >= num) {
				hi = mid;
			} else {
				lo = mid;
			}
		}
		return hi;
	};

	using T = IntMod<Mod>;
	const uint32_t middle = sqrtInt(Mod);
	std::unordered_map<T, uint32_t> dict;
	dict.reserve(middle);

	const T baby_base = a;
	T baby_temp = 1;
	for (int i = 0; i < middle; ++i) {
		if (dict.count(baby_temp)) continue;
		dict[baby_temp] = i;
		baby_temp *= baby_base;
	}

	const T giant_base = a.power(middle).inverse();
	T giant_temp = b;
	for (int k = 0; k < middle; ++k) {
		auto it = dict.find(giant_temp);
		if (it != dict.end()) {
			return k * middle + it->second;
		}
		giant_temp *= giant_base;
	}

	return 0xFFFFFFFF;
}

#endif

#ifndef __NUMBERS_H__
#define __NUMBERS_H__

/* updated: 2019-10-10 */

#include <vector>
#include <cstdint>
#include <map>

uint64_t gcd(uint64_t a, uint64_t b);
uint64_t lcm(uint64_t a, uint64_t b);

uint64_t gcd(const std::vector<uint64_t>& values);
uint64_t lcm(const std::vector<uint64_t>& values);

int64_t extendedGcd(int64_t a, int64_t b, int64_t& x, int64_t& y);
std::pair<uint64_t, uint64_t> extendedGcdUunsigned(uint64_t a, uint64_t b);

/* ax - by = c を解く関数を別に作ってもいいかも */

bool isPrime(uint64_t n);
std::vector<bool> primeTable(uint64_t max);
std::vector<uint64_t> primeList(uint64_t max);

std::map<uint64_t, int> factorize(uint64_t n);
std::vector<uint64_t> divisors(uint64_t n);

uint64_t totient(uint64_t n);
std::vector<uint64_t> totientTable(uint64_t max, const std::vector<uint64_t>& primes_ref = {});

std::map<uint64_t, int> moebius(uint64_t n);
std::vector<int> moebiusTable(uint64_t max, const std::vector<uint64_t>& primes_ref = {});


template <class T>
T power(T base, uint64_t exp);

template <class T>
void transformDivisorZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref = {});
template <class T>
void transformInverseDivisorZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref = {});
template <class T>
void transformMultiplierZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref = {});
template <class T>
void transformInverseMultiplierZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref = {});
template <class T>
std::vector<T> convolveOnLcm(const std::vector<T>& a, const std::vector<T>& b);
template <class T>
std::vector<T> convolveOnGcd(const std::vector<T>& a, const std::vector<T>& b);



/* 累乗; 置く場所がないのでここに置いておきます O(logN) */
template <class T>
T power(T base, uint64_t exp) {
	T product = T(1);
	while (exp > 0) {
		if (exp & 1) {
			product *= base;
		}
		base *= base;
		exp >>= 1;
	}
	return product;
}

/* 約数ゼータ変換(破壊的) O(NloglogN); 第0要素は使用しない */
/* g[n] = sum [k divides n] f[k] */
template <class T>
void transformDivisorZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref) {
	const uint64_t n = arr.size() - 1;
	const std::vector<uint64_t>& primes = primes_ref.empty() ? primeList(n) : primes_ref;

	for (uint64_t p : primes) {
		for (uint64_t i = 1; i <= n / p; ++i) {
			arr[i * p] += arr[i];
		}
	}
}

/* 約数ゼータ変換の逆変換(破壊的) O(NloglogN) */
template <class T>
void transformInverseDivisorZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref) {
	const uint64_t n = arr.size() - 1;
	const std::vector<uint64_t>& primes = primes_ref.empty() ? primeList(n) : primes_ref;

	for (uint64_t p : primes) {
		for (uint64_t i = n / p; i >= 1; --i) {
			arr[i * p] -= arr[i];
		}
	}
}

/* 倍数ゼータ変換(破壊的) O(NloglogN) */
/* g[n] = sum [k is divided by n] f[k] */
template<class T>
void transformMultiplierZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref) {
	const uint64_t n = arr.size() - 1;
	const std::vector<uint64_t>& primes = primes_ref.empty() ? primeList(n) : primes_ref;

	for (uint64_t p : primes) {
		for (uint64_t i = n / p; i >= 1; --i) {
			arr[i] += arr[i * p];
		}
	}
}


/* 倍数ゼータ変換の逆変換(破壊的) O(NloglogN) */
template<class T>
void transformInverseMultiplierZeta(std::vector<T>& arr, const std::vector<uint64_t>& primes_ref) {
	const uint64_t n = arr.size() - 1;
	const std::vector<uint64_t>& primes = primes_ref.empty() ? primeList(n) : primes_ref;

	for (uint64_t p : primes) {
		for (uint64_t i = 1; i <= n / p; ++i) {
			arr[i] -= arr[i * p];
		}
	}
}

/* LCM畳み込み O(NloglogN); 第0要素は使用しない */
/* ret[n] = sum [lcm(i, j) == n] (a[i] * b[j]) */
template<class T>
std::vector<T> convolveOnLcm(const std::vector<T>& a, const std::vector<T>& b) {
	assert(a.size() == b.size());
	const int size = a.size();

	std::vector<uint64_t> primes = primeList(size - 1);
	std::vector<T> a_copied = a;
	std::vector<T> b_copied = b;
	transformDivisorZeta(a_copied, primes);
	transformDivisorZeta(b_copied, primes);

	std::vector<T> ret(size);
	for (int i = 1; i < size; ++i) {
		ret[i] = a_copied[i] * b_copied[i];
	}

	transformInverseDivisorZeta(ret, primes);
	return ret;
}

/* GCD畳み込み O(NloglogN); 第0要素は使用しない */
/* ret[n] = sum [gcd(i, j) == n] (a[i] * b[j]) */
template<class T>
inline std::vector<T> convolveOnGcd(const std::vector<T>& a, const std::vector<T>& b) {
	assert(a.size() == b.size());
	const int size = a.size();

	std::vector<uint64_t> primes = primeList(size - 1);
	std::vector<T> a_copied = a;
	std::vector<T> b_copied = b;
	transformMultiplierZeta(a_copied, primes);
	transformMultiplierZeta(b_copied, primes);

	std::vector<T> ret(size);
	for (int i = 1; i < size; ++i) {
		ret[i] = a_copied[i] * b_copied[i];
	}

	transformInverseMultiplierZeta(ret, primes);
	return ret;
}

#endif

#ifdef ONLY_MY_ENVIR
#include "Numbers.h"
#endif

#include <numeric>
#include <cassert>
#include <algorithm>

uint64_t gcd(uint64_t a, uint64_t b) {
	if (b == 0) return a;
	return gcd(b, a % b);
}

uint64_t lcm(uint64_t a, uint64_t b) {
	return a / gcd(a, b) * b;
}

uint64_t gcd(const std::vector<uint64_t>& values) {
	using FuncType = uint64_t(*)(uint64_t, uint64_t);
	return std::accumulate(values.begin(), values.end(), 0ULL, (FuncType)gcd);
}

uint64_t lcm(const std::vector<uint64_t>& values) {
	using FuncType = uint64_t(*)(uint64_t, uint64_t);
	return std::accumulate(values.begin(), values.end(), 1ULL, (FuncType)lcm);
}

/* ax + by = gcd(a, b) の整数解(絶対値最小?)を求める
 * a > 1 || b > 1 のとき abs(x) <= b / (2g) && abs(y) <= a / (2g) を満たす (g = gcd(a, b))
 * 返り値は gcd(a, b) */
int64_t extendedGcd(int64_t a, int64_t b, int64_t& x, int64_t& y) {
	if (b == 0) {
		x = 1;
		y = 0;
		return a;
	}
	uint64_t ret = extendedGcd(b, a % b, y, x);
	y -= (a / b) * x;
	return ret;
}

/* ax - by = gcd(a, b) の非負整数解を求める (x < b / g && y < a / g を満たす, g = gcd(a, b)) */
std::pair<uint64_t, uint64_t> extendedGcdUunsigned(uint64_t a, uint64_t b) {
	assert(a != 0 && b != 0);

	/* もともとの方程式の解 (x, y) が、
	 * 現在の (a, b) についての ap - bq = 1 の解 (p, q) の一次式で
	 * (x, y) = mat * (p, q) と書ける(ように mat を定める) */
	uint64_t mat[2][2] = {
		{ 1, 0 },
		{ 0, 1 }
	};

	while (true) {
		if (a > b) {
			if (a % b == 0) return std::make_pair(
				mat[0][0] + (a / b - 1) * mat[0][1],
				mat[1][0] + (a / b - 1) * mat[1][1]
			);
			mat[0][0] += (a / b) * mat[0][1];
			mat[1][0] += (a / b) * mat[1][1];
			a %= b;
		} else {
			if (b % a == 0) return std::make_pair(mat[0][0], mat[1][0]);
			mat[0][1] += (b / a) * mat[0][0];
			mat[1][1] += (b / a) * mat[1][0];
			b %= a;
		}
	}
}

/* 素数判定 O(sqrt(N)) */
bool isPrime(uint64_t n) {
	if (n <= 1) return false;
	for (uint64_t i = 2; i * i <= n; ++i) {
		if (n % i == 0) return false;
	}
	return true;
}

/* エラトステネスの篩 O(NloglogN) */
std::vector<bool> primeTable(uint64_t max) {
	std::vector<bool> ret(max + 1, true);
	ret[0] = false;
	ret[1] = false;

	auto nextPrime = [max, &ret](uint64_t p) -> uint64_t {
		do {
			++p;
		} while (p <= max && !ret[p]);
		return p;
	};

	for (uint64_t p = 2; p * p <= max; p = nextPrime(p)) {
		for (uint64_t i = p * p; i <= max; i += p) {
			ret[i] = false;
		}
	}

	return ret;
}

/* 素数の昇順リストを生成 O(N * loglogN) */
std::vector<uint64_t> primeList(uint64_t max) {
	std::vector<bool> is_prime = primeTable(max);

	std::vector<uint64_t> list;
	for (uint64_t i = 2; i <= max; ++i) {
		if (is_prime[i]) {
			list.push_back(i);
		}
	}

	return list;
}

/* 素因数分解 O(sqrt(N)) */
std::map<uint64_t, int> factorize(uint64_t n) {
	std::map<uint64_t, int> ret;
	
	for (uint64_t i = 2; i * i <= n; ++i) {
		while (n % i == 0) {
			++ret[i];
			n /= i;
		}
	}

	if (n > 1) {
		++ret[n];
	}

	return ret;
}

/* 約数列挙 O(sqrt(N)) */
std::vector<uint64_t> divisors(uint64_t n) {
	std::map<uint64_t, int> factors = factorize(n);
	std::vector<uint64_t> ret;
	ret.push_back(1);

	for (const auto& pair : factors) {
		uint64_t prime = pair.first;
		int exp = pair.second;

		const uint64_t old_size = ret.size();
		for (uint64_t i = 0; i < old_size; ++i) {
			uint64_t divisor = ret[i];
			for (int k = 1; k <= exp; ++k) {
				divisor *= prime;
				ret.push_back(divisor);
			}
		}
	}
	std::sort(ret.begin(), ret.end());
	
	return ret;
}

/* オイラーのφ関数 O(sqrt(N)) */
uint64_t totient(uint64_t n) {
	assert(n > 0);

	std::map<uint64_t, int> factors = factorize(n);

	uint64_t ret = n;
	for (const auto& pair : factors) {
		uint64_t prime = pair.first;
		ret = ret / prime * (prime - 1);
	}
	return ret;
}

/* [1, n] に関してオイラーのφ関数の表を生成 O(NloglogN) */
std::vector<uint64_t> totientTable(uint64_t max, const std::vector<uint64_t>& primes_ref) {
	const std::vector<uint64_t>& primes = primes_ref.empty() ? primeList(max) : primes_ref;
	
	std::vector<uint64_t> ret(max + 1);
	std::iota(ret.begin(), ret.end(), uint64_t());
	for (uint64_t p : primes) {
		for (uint64_t i = p; i <= max; i += p) {
			ret[i] = ret[i] / p * (p - 1);
		}
	}

	return ret;
}

/* n の約数それぞれに関してメビウス関数の値を求める O(sqrt(N)) */
std::map<uint64_t, int> moebius(uint64_t n) {
	assert(n > 0);

	std::vector<uint64_t> primes;
	for (const auto& pair : factorize(n)) {
		primes.push_back(pair.first);
	}

	int prime_count = primes.size();
	std::map<uint64_t, int> ret;
	for (uint64_t s = 0; s < (uint64_t(1) << prime_count); ++s) {
		uint64_t product = 1;
		int value = 1;
		for (int i = 0; i < prime_count; ++i) {
			if (s & (uint64_t(1) << i)) {
				product *= primes[i];
				value = -value;
			}
		}
		ret[product] = value;
	}

	return ret;
}

/* [1, n] に関してメビウス関数の表を生成 O(NloglogN) */
std::vector<int> moebiusTable(uint64_t max, const std::vector<uint64_t>& primes_ref) {
	const std::vector<uint64_t>& primes = primes_ref.empty() ? primeList(max) : primes_ref;

	std::vector<int> ret(max + 1, 1);
	ret[0] = 0;

	/* p を約数に持つなら -1 倍 */
	for (uint64_t p : primes) {
		for (uint64_t i = p; i <= max; i += p) {
			ret[i] *= -1;
		}
	}

	/* p^2 を約数に持つなら 0 */
	for (uint64_t p : primes) {
		uint64_t p_sqr = p * p;		// オーバーフロー対策は必要ないと思うのでしない
		if (p_sqr > max) break;
		for (uint64_t i = p_sqr; i <= max; i += p_sqr) {
			ret[i] = 0;
		}
	}
	return ret;
}



Compiler output:
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/14/../../../x86_64-linux-gnu/Scrt1.o: in function `_start':
(.text+0x17): undefined reference to `main'
collect2: error: ld returned 1 exit status