Language: C++

#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> extendedGcdUnsigned(uint64_t a, uint64_t b);
uint64_t inverse(uint64_t a, uint64_t mod);


/* 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>
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(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>
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> extendedGcdUnsigned(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;
		}
	}
}

/* mod における a の逆数(a と mod は互いに素) */
uint64_t inverse(uint64_t a, uint64_t mod) {
	assert(a != 0 && mod != 0);
	assert(gcd(a, mod) == 1);

	std::pair<uint64_t, uint64_t> res = extendedGcdUnsigned(a, mod);
	return res.first;
}

/* 素数判定 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;
}

#ifndef __DYNAMIC_MOD_H__
#define __DYNAMIC_MOD_H__

/* updated: 2020-06-26 */

// required: Number

#include <vector>
#include <cstdint>
#include <iostream>

/* 複数の mod を扱うためのクラス (64 bit 対応) */
/* 演算の結果をどう返すかなど、まだ考えることはありそう */
class DynamicMod {
	uint64_t value_m;
	uint64_t mod_m;

	static uint64_t default_mod_m;

public:
	DynamicMod();
	explicit DynamicMod(int32_t value, uint64_t mod = default_mod_m);
	explicit DynamicMod(uint32_t value, uint64_t mod = default_mod_m);
	explicit DynamicMod(int64_t value, uint64_t mod = default_mod_m);
	explicit DynamicMod(uint64_t value, uint64_t mod = default_mod_m);
	DynamicMod(const DynamicMod&) = default;
	DynamicMod& operator=(const DynamicMod&) = default;

	const DynamicMod& operator+() const;
	DynamicMod operator-() const;
	DynamicMod& operator++();
	DynamicMod& operator--();
	DynamicMod operator++(int dummy);
	DynamicMod operator--(int dummy);
	DynamicMod& operator+=(const DynamicMod& right);
	DynamicMod& operator-=(const DynamicMod& right);
	DynamicMod& operator*=(const DynamicMod& right);

	template <class Int>
	DynamicMod& operator+=(const Int& right);
	template <class Int>
	DynamicMod& operator-=(const Int& right);
	template <class Int>
	DynamicMod& operator*=(const Int& right);

	static bool hasSameMod(const DynamicMod& a, const DynamicMod& b);

	static DynamicMod add(const DynamicMod& a, const DynamicMod& b);
	static DynamicMod subtract(const DynamicMod& a, const DynamicMod& b);
	static DynamicMod multiply(const DynamicMod& a, const DynamicMod& b);
	static DynamicMod divide(const DynamicMod& a, const DynamicMod& b);
	static int compare(const DynamicMod& a, const DynamicMod& b);

	DynamicMod inverse() const;
	DynamicMod power(uint64_t exp) const;

	uint64_t getValue() const;
	uint64_t getMod() const;

	static DynamicMod conjunction(const DynamicMod& a, const DynamicMod& b);
	static DynamicMod conjunction(const std::vector<DynamicMod>& values);

	static void setDefaultMod(uint64_t mod);
};

DynamicMod operator+(const DynamicMod& left, const DynamicMod& right);
DynamicMod operator-(const DynamicMod& left, const DynamicMod& right);
DynamicMod operator*(const DynamicMod& left, const DynamicMod& right);
template <class Int>
DynamicMod operator+(const Int& left, const DynamicMod& right);
template <class Int>
DynamicMod operator-(const Int& left, const DynamicMod& right);
template <class Int>
DynamicMod operator*(const Int& left, const DynamicMod& right);
template <class Int>
DynamicMod operator+(const DynamicMod& left, const Int& right);
template <class Int>
DynamicMod operator-(const DynamicMod& left, const Int& right);
template <class Int>
DynamicMod operator*(const DynamicMod& left, const Int& right);

bool operator==(const DynamicMod& left, const DynamicMod& right);
bool operator!=(const DynamicMod& left, const DynamicMod& right);
bool operator<(const DynamicMod& left, const DynamicMod& right);
bool operator<=(const DynamicMod& left, const DynamicMod& right);
bool operator>(const DynamicMod& left, const DynamicMod& right);
bool operator>=(const DynamicMod& left, const DynamicMod& right);

template<class Int>
DynamicMod& DynamicMod::operator+=(const Int& right) {
	*this += DynamicMod(right, getMod());
	return *this;
}

template<class Int>
DynamicMod& DynamicMod::operator-=(const Int& right) {
	*this -= DynamicMod(right, getMod());
	return *this;
}

template<class Int>
DynamicMod& DynamicMod::operator*=(const Int& right) {
	*this *= DynamicMod(right, getMod());
	return *this;
}

template<class Int>
DynamicMod operator+(const Int& left, const DynamicMod& right) {
	return DynamicMod(left, right.getMod()) + right;
}

template<class Int>
DynamicMod operator-(const Int& left, const DynamicMod& right) {
	return DynamicMod(left, right.getMod()) - right;
}

template<class Int>
DynamicMod operator*(const Int& left, const DynamicMod& right) {
	return DynamicMod(left, right.getMod()) * right;
}

template<class Int>
DynamicMod operator+(const DynamicMod& left, const Int& right) {
	return left + DynamicMod(right, left.getMod());
}

template<class Int>
DynamicMod operator-(const DynamicMod& left, const Int& right) {
	return left - DynamicMod(right, left.getMod());
}

template<class Int>
DynamicMod operator*(const DynamicMod& left, const Int& right) {
	return left * DynamicMod(right, left.getMod());
}

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

#include <algorithm>
#include <map>
#include <cassert>
#include <stdexcept>
#include <limits>

uint64_t DynamicMod::default_mod_m = 1;

DynamicMod::DynamicMod()
	: value_m(0), mod_m(default_mod_m) {
}

DynamicMod::DynamicMod(int32_t value, uint64_t mod) 
	: DynamicMod(int64_t(value), mod) {
}

DynamicMod::DynamicMod(uint32_t value, uint64_t mod)
	: DynamicMod(uint64_t(value), mod) {
}

DynamicMod::DynamicMod(int64_t value, uint64_t mod) {
	assert(mod != 0);

	const uint64_t thres = uint64_t(1) << 63;
	if (mod >= thres) {
		// mod >= 2^63 のとき
		if (mod == thres && uint64_t(value) == thres) {
			// mod == 2^63, value = -2^63 のときは例外として 0
			value_m = 0;
		} else {
			// それ以外の場合は正ならそのまま、負なら mod を足す (キャスト注意, 規格にはあっているはず)
			value_m = value >= 0 ? uint64_t(value) : uint64_t(value) + mod;
		}
	} else {
		// mod を符号付きにキャストできるならキャストして剰余を計算
		int64_t tmp = value % int64_t(mod);
		value_m = tmp >= 0 ? tmp : mod - uint64_t(-tmp);
	}
	mod_m = mod;
}

DynamicMod::DynamicMod(uint64_t value, uint64_t mod)
	: value_m(value % mod), mod_m(mod) {
	assert(mod != 0);
}

const DynamicMod& DynamicMod::operator+() const {
	return *this;
}

DynamicMod DynamicMod::operator-() const {
	DynamicMod ret;
	ret.mod_m = mod_m;
	if (value_m == 0) {
		ret.value_m = 0;
	} else {
		ret.value_m = mod_m - value_m;
	}
	return ret;
}

DynamicMod& DynamicMod::operator++() {
	++value_m;
	if (value_m == mod_m) {
		value_m = 0;
	}
	return *this;
}

DynamicMod& DynamicMod::operator--() {
	if (value_m == 0) {
		value_m = mod_m;
	}
	--value_m;
	return *this;
}

DynamicMod DynamicMod::operator++(int dummy) {
	DynamicMod ret(*this);
	++*this;
	return ret;
}

DynamicMod DynamicMod::operator--(int dummy) {
	DynamicMod ret(*this);
	--*this;
	return ret;
}

DynamicMod& DynamicMod::operator+=(const DynamicMod& right) {
	if (!hasSameMod(*this, right)) throw std::runtime_error("not compatible with different mods");

	if (mod_m - value_m <= right.value_m) {
		value_m -= mod_m;
	}
	value_m += right.value_m;
	return *this;
}

DynamicMod& DynamicMod::operator-=(const DynamicMod& right) {
	if (!hasSameMod(*this, right)) throw std::runtime_error("not compatible with different mods");

	if (value_m < right.value_m) {
		value_m += mod_m;
	}
	value_m -= right.value_m;
	return *this;
}

DynamicMod& DynamicMod::operator*=(const DynamicMod& right) {
	if (!hasSameMod(*this, right)) throw std::runtime_error("not compatible with different mods");
	*this = multiply(*this, right);
	return *this;
}

bool DynamicMod::hasSameMod(const DynamicMod& a, const DynamicMod& b) {
	return a.mod_m == b.mod_m;
}

DynamicMod DynamicMod::add(const DynamicMod& a, const DynamicMod& b) {
	if (!hasSameMod(a, b)) throw std::runtime_error("not compatible with different mods");

	DynamicMod ret(a);
	ret += b;
	return ret;
}

DynamicMod DynamicMod::subtract(const DynamicMod& a, const DynamicMod& b) {
	if (!hasSameMod(a, b)) throw std::runtime_error("not compatible with different mods");

	DynamicMod ret(a);
	ret -= b;
	return ret;
}

DynamicMod DynamicMod::multiply(const DynamicMod& a, const DynamicMod& b) {
	if (!hasSameMod(a, b)) throw std::runtime_error("not compatible with different mods");

	auto multiply32 = [](uint64_t a, uint64_t b, uint64_t mod) -> uint64_t {
		return a * b % mod;
	};

	auto multiply64 = [](uint64_t a, uint64_t b, uint64_t mod) -> uint64_t {
		auto assignAdd = [mod](uint64_t& a, uint64_t b) -> void {
			if (mod - a <= b) {
				a -= mod;
			}
			a += b;
		};

		uint64_t product = 0;
		uint64_t addend = a;
		uint64_t multiplier = b;
		while (multiplier > 0) {
			if (multiplier & 1) {
				assignAdd(product, addend);
			}
			assignAdd(addend, addend);
			multiplier >>= 1;
		}
		return product;
	};

	const uint64_t threshold = uint64_t(1) << 32;
	const uint64_t mod = a.mod_m;
	DynamicMod ret;
	if (mod < threshold) {
		ret.value_m = multiply32(a.value_m, b.value_m, mod);
	} else {
		ret.value_m = multiply64(a.value_m, b.value_m, mod);
	}
	ret.mod_m = mod;
	return ret;
}

/* bx = a (mod m) を満たす x の条件を与える (引数と戻り値の mod は一致するとは限らない) */
DynamicMod DynamicMod::divide(const DynamicMod& a, const DynamicMod& b) {
	if (!hasSameMod(a, b)) throw std::runtime_error("not compatible with different mods");
	if (b.getValue() == 0) throw std::runtime_error("division by zero");

	const uint64_t va = a.value_m;
	const uint64_t vb = b.value_m;
	const uint64_t mod = a.mod_m;

	const uint64_t gcd = ::gcd(vb, mod);
	if (va % gcd != 0) throw std::runtime_error("no solution");

	auto pair = extendedGcdUnsigned(vb, mod);
	return (va / gcd) * DynamicMod(pair.first, mod / gcd);
}

int DynamicMod::compare(const DynamicMod& a, const DynamicMod& b) {
	if (a.mod_m < b.mod_m) return -1;
	if (a.mod_m > b.mod_m) return 1;
	if (a.value_m < b.value_m) return -1;
	if (a.value_m > b.value_m) return 1;
	return 0;
}

/* mod_m must be a prime. */
DynamicMod DynamicMod::inverse() const {
	assert(value_m != 0);
	return power(mod_m - 2);
}

DynamicMod DynamicMod::power(uint64_t exp) const {
	DynamicMod product(1, mod_m);
	DynamicMod factor(*this);
	while (exp > 0) {
		if (exp & 1) {
			product *= factor;
		}
		factor *= factor;
		exp >>= 1;
	}
	return product;
}

uint64_t DynamicMod::getValue() const {
	return value_m;
}

uint64_t DynamicMod::getMod() const {
	return mod_m;
}

/* 中国剰余定理により a と b の条件に共通して当てはまる条件を求める O(log(mod)) */
/* mod がオーバーフローしたとき mod = 0 にすることも考えている */
DynamicMod DynamicMod::conjunction(const DynamicMod& a, const DynamicMod& b) {
	/* 後の計算の都合上 a.value_m <= b.value_m としておく */
	if (a.value_m > b.value_m) return conjunction(b, a);

	const uint64_t ma = a.mod_m;
	const uint64_t mb = b.mod_m;
	const uint64_t va = a.value_m;
	const uint64_t vb = b.value_m;

	const uint64_t gcd = ::gcd(ma, mb);
	const uint64_t lcm = ma / gcd * mb;
	const uint64_t diff = vb - va;

	/* 例外ではなく mod_m == 0 のオブジェクトを返すなども考えている */
	if (diff % gcd != 0) throw std::runtime_error("no solution");

	/* 方程式 ma * x - mb * y = gcd の解を求める */
	auto pair = extendedGcdUnsigned(ma, mb);
	DynamicMod x(pair.first, mb / gcd);			// mod に注意

	/* 答えは va + ma * (diff / gcd) * x となる */
	x *= DynamicMod(diff / gcd, x.mod_m);
	return DynamicMod(va + ma * x.value_m, lcm);
}

/* 中国剰余定理により複数の条件に共通して当てはまるものを求める */
DynamicMod DynamicMod::conjunction(const std::vector<DynamicMod>& values) {
	DynamicMod product(0, 1);
	for (const DynamicMod& value : values) {
		product = conjunction(product, value);
	}
	return product;
}

/* mod の初期値を設定 */
void DynamicMod::setDefaultMod(uint64_t mod) {
	assert(mod != 0);
	default_mod_m = mod;
}

DynamicMod operator+(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::add(left, right);
}

DynamicMod operator-(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::subtract(left, right);
}

DynamicMod operator*(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::multiply(left, right);
}

bool operator==(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::compare(left, right) == 0;
}

bool operator!=(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::compare(left, right) != 0;
}

bool operator<(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::compare(left, right) < 0;
}

bool operator<=(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::compare(left, right) <= 0;
}

bool operator>(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::compare(left, right) > 0;
}

bool operator>=(const DynamicMod& left, const DynamicMod& right) {
	return DynamicMod::compare(left, right) >= 0;
}

#define _CRT_SECURE_NO_WARNINGS
#define _SCL_SECURE_NO_WARNINGS
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>
#include <string>
#include <vector>
#include <list>
#include <utility>
#include <algorithm>
#include <functional>
#include <cmath>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <iomanip>
#include <sstream>
#include <bitset>
#include <limits>
#include <numeric>
#include <valarray>
#include <fstream>
#include <array>
#include <random>
#include <unordered_set>
#include <unordered_map>

using namespace std;
using uint = uint32_t;
using LL = int64_t;
using ULL = uint64_t;
using PP = pair<LL, LL>;
template <typename T> using PriorityQ = priority_queue<T, vector<T>, greater<T> >;
#define REP(i, a, n) for(LL i = (a), i##_max_ = (n); i < i##_max_; ++i)
#define REM(i, a, n) for(LL i = (LL)(n) - 1, i##_min_ = (a); i >= i##_min_; --i)
#define FLOAT fixed << setprecision(16)
#define SPEEDUP { cin.tie(NULL); ios::sync_with_stdio(false); }
const int INF = 0x3FFFFFFF;
const LL INFLL = 0x3FFFFFFF3FFFFFFF;
const double INFD = 1.0e+308;
const double EPS = 1.0e-9;

void YesNo(bool b) { cout << (b ? "Yes" : "No") << endl; }
void YESNO(bool b) { cout << (b ? "YES" : "NO") << endl; }
template <class T, class U> istream& operator>>(istream& ist, pair<T, U>& right) { return ist >> right.first >> right.second; }
template <class T, class U> ostream& operator<<(ostream& ost, const pair<T, U>& right) { return ost << right.first << ' ' << right.second; }
template <class T, class TCompatible, size_t N> void Fill(T(&dest)[N], const TCompatible& val) { fill(dest, dest + N, val); }
template <class T, class TCompatible, size_t M, size_t N> void Fill(T(&dest)[M][N], const TCompatible& val) { for (int i = 0; i < M; ++i) Fill(dest[i], val); }
template <class T> T Next() { T buf; cin >> buf; return buf; }
istream& Ignore(istream& ist) { string s; ist >> s; return ist; }
bool Inside(int i, int j, int h, int w) { return i >= 0 && i < h&& j >= 0 && j < w; }

#ifdef ONLY_MY_ENVIR
#include "Accumulator.h"
#include "Algebraic.h"
#include "BinaryMatrix.h"
#include "BinaryTree.h"
#include "Bipartite.h"
#include "BIT.h"
#include "Compressor.h"
#include "Decompositions.h"
#include "DiscreteLog.h"
#include "DynamicMod.h"
#include "Exponential.h"
#include "Factorization.h"
#include "FFT.h"
#include "FlowSolver.h"
#include "Geometric2D.h"
#include "Geometric2DFloat.h"
#include "Geometric3D.h"
#include "Geometric3DFloat.h"
#include "Graph.h"
#include "GraphUtil.h"
#include "Interpolation.h"
#include "IntMod.h"
#include "KDTree.h"
#include "LazySegmentTree.h"
#include "LIS.h"
#include "List.h"
#include "Math.h"
#include "MathUtil.h"
#include "Matrix.h"
#include "MinCostFlowSolver.h"
#include "MinMax.h"
#include "Numbers.h"
#include "Optimize.h"
#include "Permutation.h"
#include "Polynomial.h"
#include "Position.h"
#include "Random.h"
#include "Range.h"
#include "Rational.h"
#include "RollingHash.h"
#include "RuntimeMod.h"
#include "SegmentTree.h"
#include "SegmentTree2D.h"
#include "Sets.h"
#include "Shortest.h"
#include "SlidingWindow.h"
#include "SpanningTree.h"
#include "StringSearching.h"
#include "SuffixArray.h"
#include "SwitchList.h"
#include "Timer.h"
#include "Tree.h"
#include "TreeUtil.h"
#include "UnionFind.h"
#include "Util.h"
#include "VectorUtil.h"
#endif

#ifdef __GNUC__
typedef __int128 LLL;
istream& operator>>(istream& ist, __int128& val) { LL tmp; ist >> tmp; val = tmp; return ist; }
ostream& operator<<(ostream& ost, __int128 val) { LL tmp = val; ost << tmp; return ost; }
#endif

#include <cctype>
#include <bitset>
#include <iostream>
#include <string>
#include <vector>
using namespace std;


LL n;
int main() {
	cin >> n;
	vector<DynamicMod> ds;
	ds.push_back(DynamicMod(0, 1));
	for (auto p : factorize(n)) {
		LL pk = power(p.first, p.second);

		vector<DynamicMod> es;

		if (p.first == 2) {
			es.push_back(DynamicMod(0, 2 * pk));
			es.push_back(DynamicMod(2 * pk - 1, 2 * pk));
		} else {
			es.push_back(DynamicMod(0, pk));
			es.push_back(DynamicMod(pk - 1, pk));
		}

		vector<DynamicMod> nds;
		for (auto d : ds) {
			for (auto e : es) {
				nds.push_back(DynamicMod::conjunction(d, e));
			}
		}
		ds = move(nds);
	}

	LL ans = INFLL;
	for (auto d : ds) {
		LL k = d.getValue() == 0 ? d.getMod() : d.getValue();
		ans = min(ans, k);
	}
	cout << ans << endl;
	return 0;
}


Compiler output:
a.cc: In function 'int main()':
a.cc:977:25: error: 'power' was not declared in this scope
  977 |                 LL pk = power(p.first, p.second);
      |                         ^~~~~