結果

問題 No.2181 LRM Question 2
コンテスト
ユーザー MasKoaTS
提出日時 2022-08-28 23:51:47
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 793 ms / 2,000 ms
コード長 8,759 bytes
コンパイル時間 1,506 ms
コンパイル使用メモリ 97,832 KB
最終ジャッジ日時 2025-02-06 23:21:36
ジャッジサーバーID
（参考情報） 		judge4 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 23
権限があれば一括ダウンロードができます

ソースコード

diff # 

/*
* Copyright (c) 2019 Ryotaro Sato
* All functions and data structures prior to the main function are released under the MIT license.
* https://github.com/hitonanode/cplib-cpp/blob/master/LICENSE
*/

#include <iostream>
#include <algorithm>
#include <cassert>
#include <tuple>
#include <utility>
#include <vector>
#include <map>
using ll = long long;

// CUT begin
// Solve ax+by=gcd(a, b)
template <class Int> Int extgcd(Int a, Int b, Int& x, Int& y) {
	Int d = a;
	if (b != 0) {
		d = extgcd(b, a % b, y, x), y -= (a / b) * x;
	}
	else {
		x = 1, y = 0;
	}
	return d;
}
// Calculate a^(-1) (MOD m) s if gcd(a, m) == 1
// Calculate x s.t. ax == gcd(a, m) MOD m
template <class Int> Int mod_inverse(Int a, Int m) {
	Int x, y;
	extgcd<Int>(a, m, x, y);
	x %= m;
	return x + (x < 0) * m;
}

// Require: 1 <= b
// return: (g, x) s.t. g = gcd(a, b), xa = g MOD b, 0 <= x < b/g
template <class Int> /* constexpr */ std::pair<Int, Int> inv_gcd(Int a, Int b) {
	a %= b;
	if (a < 0) a += b;
	if (a == 0) return { b, 0 };
	Int s = b, t = a, m0 = 0, m1 = 1;
	while (t) {
		Int u = s / t;
		s -= t * u, m0 -= m1 * u;
		auto tmp = s;
		s = t, t = tmp, tmp = m0, m0 = m1, m1 = tmp;
	}
	if (m0 < 0) m0 += b / s;
	return { s, m0 };
}

template <class Int>
/* constexpr */ std::pair<Int, Int> crt(const std::vector<Int>& r, const std::vector<Int>& m) {
	assert(r.size() == m.size());
	int n = int(r.size());
	// Contracts: 0 <= r0 < m0
	Int r0 = 0, m0 = 1;
	for (int i = 0; i < n; i++) {
		assert(1 <= m[i]);
		Int r1 = r[i] % m[i], m1 = m[i];
		if (r1 < 0) r1 += m1;
		if (m0 < m1) {
			std::swap(r0, r1);
			std::swap(m0, m1);
		}
		if (m0 % m1 == 0) {
			if (r0 % m1 != r1) return { 0, 0 };
			continue;
		}
		Int g, im;
		std::tie(g, im) = inv_gcd<Int>(m0, m1);

		Int u1 = m1 / g;
		if ((r1 - r0) % g) return { 0, 0 };

		Int x = (r1 - r0) / g % u1 * im % u1;
		r0 += x * m0;
		m0 *= u1;
		if (r0 < 0) r0 += m0;
	}
	return { r0, m0 };
}

// 蟻本 P.262
// 中国剰余定理を利用して，色々な素数で割った余りから元の値を復元
// 連立線形合同式 A * x = B mod M の解
// Requirement: M[i] > 0
// Output: x = first MOD second (if solution exists), (0, 0) (otherwise)
template <class Int>
std::pair<Int, Int>
linear_congruence(const std::vector<Int>& A, const std::vector<Int>& B, const std::vector<Int>& M) {
	Int r = 0, m = 1;
	assert(A.size() == M.size());
	assert(B.size() == M.size());
	for (int i = 0; i < (int)A.size(); i++) {
		assert(M[i] > 0);
		const Int ai = A[i] % M[i];
		Int a = ai * m, b = B[i] - ai * r, d = std::__gcd(M[i], a);
		if (b % d != 0) {
			return std::make_pair(0, 0); // 解なし
		}
		Int t = b / d * mod_inverse<Int>(a / d, M[i] / d) % (M[i] / d);
		r += m * t;
		m *= M[i] / d;
	}
	return std::make_pair((r < 0 ? r + m : r), m);
}

int pow_mod(int x, long long n, int md) {
	if (md == 1) return 0;
	long long ans = 1;
	while (n > 0) {
		if (n & 1) ans = ans * x % md;
		x = (long long)x * x % md;
		n >>= 1;
	}
	return ans;
}
#line 5 "number/combination.hpp"

// nCr mod m = p^q (p: prime, q >= 1)
// Can be used for n, r <= 1e18, m <= 1e7
// Complexity: O(m) (construction), O(log(n)) (per query)
// https://ferin-tech.hatenablog.com/entry/2018/01/17/010829
struct combination_prime_pow {
	int p, q, m;
	std::vector<int> fac, invfac, ppow;

	long long _ej(long long n) const {
		long long ret = 0;
		while (n) ret += n, n /= p;
		return ret;
	}

	combination_prime_pow(int p_, int q_) : p(p_), q(q_), m(1), ppow{ 1 } {
		for (int t = 0; t < q; ++t) m *= p, ppow.push_back(m);
		fac.assign(m, 1);
		invfac.assign(m, 1);
		for (int i = 1; i < m; ++i) fac[i] = (long long)fac[i - 1] * (i % p ? i : 1) % m;
		invfac[m - 1] = fac[m - 1]; // Same as Wilson's theorem
		assert(1LL * fac.back() * invfac.back() % m == 1);
		for (int i = m - 1; i; --i) invfac[i - 1] = (long long)invfac[i] * (i % p ? i : 1) % m;
	}

	int nCr(long long n, long long r) const {
		if (r < 0 or n < r) return 0;
		if (p == 2 and q == 1) return !((~n) & r); // Lucas
		long long k = n - r;
		long long e0 = _ej(n / p) - _ej(r / p) - _ej(k / p);
		if (e0 >= q) return 0;

		long long ret = ppow[e0];
		if (q == 1) { // Lucas
			while (n) {
				ret = __int128(ret) * fac[n % p] * invfac[r % p] * invfac[k % p] % p;
				n /= p, r /= p, k /= p;
			}
			return (int)ret;
		}
		else {
			if ((p > 2 or q < 3) and (_ej(n / m) - _ej(r / m) - _ej(k / m)) & 1) ret = m - ret;
			while (n) {
				ret = __int128(ret) * fac[n % m] * invfac[r % m] * invfac[k % m] % m;
				n /= p, r /= p, k /= p;
			}
			return (int)ret;
		}
	}
};

// nCr mod m
// Complexity: O(m) space worst (construction), O(log(n) log(m)) (per query)
// Input: pairs of (prime, degree), such as vector<pair<int, int>> and map<int, int>
// https://judge.yosupo.jp/problem/binomial_coefficient
struct combination {
	std::vector<combination_prime_pow> cpps;
	std::vector<int> ms;

	template <class Map> combination(const Map& p2deg) {
		for (auto f : p2deg) {
			cpps.push_back(combination_prime_pow(f.first, f.second));
			ms.push_back(cpps.back().m);
		}
	}

	int operator()(long long n, long long r) const {
		if (r < 0 or n < r) return 0;
		std::vector<int> rs;
		for (const auto& cpp : cpps) rs.push_back(cpp.nCr(n, r));
		return crt(rs, ms).first;
	}
};


struct Sieve {
	std::vector<int> min_factor;
	std::vector<int> primes;
	Sieve(int MAXN) : min_factor(MAXN + 1) {
		for (int d = 2; d <= MAXN; d++) {
			if (!min_factor[d]) {
				min_factor[d] = d;
				primes.emplace_back(d);
			}
			for (const auto& p : primes) {
				if (p > min_factor[d] or d * p > MAXN) break;
				min_factor[d * p] = p;
			}
		}
	}
	// Prime factorization for 1 <= x <= MAXN^2
	// Complexity: O(log x)		   (x <= MAXN)
	//			 O(MAXN / log MAXN) (MAXN < x <= MAXN^2)
	template <class T> std::map<T, int> factorize(T x) const {
		std::map<T, int> ret;
		assert(x > 0 and
			x <= ((long long)min_factor.size() - 1) * ((long long)min_factor.size() - 1));
		for (const auto& p : primes) {
			if (x < T(min_factor.size())) break;
			while (!(x % p)) x /= p, ret[p]++;
		}
		if (x >= T(min_factor.size())) ret[x]++, x = 1;
		while (x > 1) ret[min_factor[x]]++, x /= min_factor[x];
		return ret;
	}
	// Enumerate divisors of 1 <= x <= MAXN^2
	// Be careful of highly composite numbers https://oeis.org/A002182/list
	// https://gist.github.com/dario2994/fb4713f252ca86c1254d#file-list-txt (n, (# of div. of n)):
	// 45360->100, 735134400(<1e9)->1344, 963761198400(<1e12)->6720
	template <class T> std::vector<T> divisors(T x) const {
		std::vector<T> ret{ 1 };
		for (const auto p : factorize(x)) {
			int n = ret.size();
			for (int i = 0; i < n; i++) {
				for (T a = 1, d = 1; d <= p.second; d++) {
					a *= p.first;
					ret.push_back(ret[i] * a);
				}
			}
		}
		return ret; // NOT sorted
	}
	// Euler phi functions of divisors of given x
	// Verified: ABC212 G https://atcoder.jp/contests/abc212/tasks/abc212_g
	// Complexity: O(sqrt(x) + d(x))
	template <class T> std::map<T, T> euler_of_divisors(T x) const {
		assert(x >= 1);
		std::map<T, T> ret;
		ret[1] = 1;
		std::vector<T> divs{ 1 };
		for (auto p : factorize(x)) {
			int n = ret.size();
			for (int i = 0; i < n; i++) {
				ret[divs[i] * p.first] = ret[divs[i]] * (p.first - 1);
				divs.push_back(divs[i] * p.first);
				for (T a = divs[i] * p.first, d = 1; d < p.second; a *= p.first, d++) {
					ret[a * p.first] = ret[a] * p.first;
					divs.push_back(a * p.first);
				}
			}
		}
		return ret;
	}
	// Moebius function Table, (-1)^{# of different prime factors} for square-free x
	// return: [0=>0, 1=>1, 2=>-1, 3=>-1, 4=>0, 5=>-1, 6=>1, 7=>-1, 8=>0, ...] https://oeis.org/A008683
	std::vector<int> GenerateMoebiusFunctionTable() const {
		std::vector<int> ret(min_factor.size());
		for (unsigned i = 1; i < min_factor.size(); i++) {
			if (i == 1) {
				ret[i] = 1;
			}
			else if ((i / min_factor[i]) % min_factor[i] == 0) {
				ret[i] = 0;
			}
			else {
				ret[i] = -ret[i / min_factor[i]];
			}
		}
		return ret;
	}
	// Calculate [0^K, 1^K, ..., nmax^K] in O(nmax)
	// Note: **0^0 == 1**
	template <class MODINT> std::vector<MODINT> enumerate_kth_pows(long long K, int nmax) const {
		assert(nmax < int(min_factor.size()));
		assert(K >= 0);
		if (K == 0) return std::vector<MODINT>(nmax + 1, 1);
		std::vector<MODINT> ret(nmax + 1);
		ret[0] = 0, ret[1] = 1;
		for (int n = 2; n <= nmax; n++) {
			if (min_factor[n] == n) {
				ret[n] = MODINT(n).pow(K);
			}
			else {
				ret[n] = ret[n / min_factor[n]] * ret[min_factor[n]];
			}
		}
		return ret;
	}
};


/*
* Main Code
*/

int main(void) {
	ll l, r, m;  std::cin >> l >> r >> m;
	combination nCr(Sieve(1 << 20).factorize(m));

	ll ans = 0;
	for (ll i = l; i <= r; i++) {
		ans += (ll)nCr(2 * i, i) + m - 2ll;
		ans %= m;
	}
	std::cout << ans << std::endl;

	return 0;
}
0