[筆記] Codeforces 920G: List Of Integers

[題意] 給定任意正整數 p, 則 L 是一串與 p 互質的序列 L(p) = {y | gcd(y, p) = 1}，且 L 是由小到大排序的 (gcd 就是最大公因數)。再給任一正整數 x，求 L 中比 x 大的第 k 個數值。

這題基本上是數論 + 二分搜尋的整合題，只要能求得以下函數的值再搭配二分搜尋就好：
$$G(p, x) = |L(p, x)|, L(p, x) = \{y\ |\ y \leq x\ \&\ y \in L(p)\}$$
白話一點來說，$L(p, x)$ 就是 $L(p)$ 中比 $x$ 小的數值所成的序列，因此 $G(p, x)$ 就是比 $x$ 小的數值中有多少個跟 $p$ 互質。

接著，我們就可以利用二分搜尋的方式在 [1, inf] 中去尋找 $M$ 使得 $G(p, M) = k + G(p, x)$，而 $M$ 就是我們要的答案。因為對任意正整數 $x_1, x_2, G(p, x_1) \leq G(p, x_2)\ if\ x_1\leq x_2$

因此接下來我們就專注在怎麼找出 $G(p, x)$ 即可。而這個值可以簡單的利用排容原理 (inclusion-exclusion principle) 算出來：

令 $x = p_1 ^ {o_1} \times p_2 ^ {o_2} \times \ \dots \ \times p_m ^ {o_m}$，則
$$G(p, x) = p - \sum_{i=1}^{m} \lfloor{\frac{p}{x_i}}\rfloor \quad + \quad \sum_{i=1}^{m} \sum_{j=j+1}^{m} \lfloor{\frac{p}{x_i \ x_j}}\rfloor \quad - \quad \dots$$

上面那個式子可以利用莫比烏斯方程式 (Möbius function) 改寫：
$$G(p, x) = \sum_{d|p} \mu(d) \times \lfloor{\frac{x}{d}}\rfloor$$

實作在此：

	/*
	* =====================================================================================
	*
	* Filename: 920G.cpp
	*
	* Description: codeforces 920G: List Of Integers
	*
	* Version: 1.0
	* Created: 2018年02月07日 12時48分22秒
	* Revision: none
	* Compiler: gcc
	*
	* Author: lionking
	* Organization: None
	*
	* =====================================================================================
	*/
	#include <cstdio>
	#include <vector>
	using namespace std;
	/*****************************************************************************************
	* Generate prime table
	* Given a threshold, this function will generate list all primes less than the threshold.
	*
	* Parameter:
	* 1. prime: the resultant prime table
	* 2. limit: the given threshold
	*
	* Return value:
	* Number of primes in the generated prime table
	*****************************************************************************************/
	template <typename T>
	size_t GeneratePrime(std::vector<T> &prime, const size_t limit)
	{
	prime.clear();
	std::vector<bool> check(limit, true);
	/* check[i] = false: prime else not prime *
	* mask all even number to be not prime */
	check[0] = check[1] = false;
	for (int i = 4; i < limit; i += 2) check[i] = false;
	for (int i = 9; i < limit; i += 3) check[i] = false;
	prime.push_back(2);
	prime.push_back(3);
	for (int i = 5; i<limit; i += 2) {
	if (check[i]) {
	/* is prime */
	prime.push_back(i);
	for (int j = i*i; j<limit; j += i) { check[j] = false; }
	}
	}
	return prime.size();
	}
	/**********************************************************************************************
	* Euler Totient Function : phi(n)
	* This function will generate a table containing phi(i)
	*
	* Parameter list:
	* 1. phi: resultant table containing phi(i)
	* 2. limit: threshold. Only positive integers less than this threshold will generate phi(i)
	* 3. prime: prime tables with all primes less than the given threshold (namely, 2nd parameter)
	*
	* Return value:
	* None
	**********************************************************************************************/
	template <typename T>
	void euler_totient(std::vector<T>& phi, const size_t limit, const std::vector<T>& prime)
	{
	const auto primecount = prime.size();
	phi.clear();
	phi.resize(limit, 0);
	phi[1] = phi[2] = 1;
	for(int i=3; i<=limit; ++i){
	int temp = phi[i] = i;
	for(int j=0; (j<primecount) && (prime[j]*prime[j]<=temp); ++j){
	if(!(temp%prime[j])){
	phi[i] = (phi[i]/prime[j])*(prime[j]-1);
	while(!(temp%prime[j]))
	temp /= prime[j];
	}
	}
	if(temp!=1)
	phi[i] = (phi[i]/temp)*(temp-1);
	}
	}
	/****************************************************************************************************
	* Prime Factorization
	* Given a positive integer n = p1 ^ k1 * p2 ^ k2 * ... * pm ^ km, where p1, p2, ... , pm are primes,
	* this function will return the list {p1, p2, ... , pm}
	*
	* Parameters:
	* 1. prime: prime table
	* 2. value: the given value for prime factorization (namely, n described above)
	*
	* Return value:
	* Required prime list (namely, {k1, k2, ... , km} described above)
	****************************************************************************************************/
	template <typename T>
	std::vector<T> factorization(const std::vector<T>& prime, const T& value)
	{
	std::vector<T> result;
	T curr = value;
	for (size_t i = 0; i < prime.size(); ++i) {
	if (prime[i] * prime[i] > curr) { break; }
	if ((curr % prime[i]) == 0) {
	result.push_back(prime[i]);
	while ((curr % prime[i]) == 0) { curr /= prime[i]; }
	}
	}
	if (curr > static_cast<T>(1)) { result.push_back(curr); }
	return result;
	}
	/******************************************************************************************************************
	* Given two positive integers n and x, where x = p1 ^ k1 * p2 ^ k2 * ... * pm ^ km, and p1, p2, ... pm are primes,
	* this function will return the number of values y that y <= n and gcd(y, x) = 1, where
	* gcd(y, x) is the greatest common divisor of y and x.
	* Namely, it will find how many values are less than n are coprime with x.
	*
	* Example: Let n = 12 and x = 6, the result is 4 ({1, 5, 7, 11})
	* Note: When n = x, the result is equal to phi(n), where phi is Euler Totient Function.
	*
	* Parameters:
	* 1. prime: prime list {p1, p2, ... , pm}
	* 2. value: n
	*
	* Return value:
	* Number of values that are less than n and coprime with x.
	******************************************************************************************************************/
	template <typename T>
	T coprimeCount(const std::vector<T>& prime, const T value)
	{
	T ans = 0;
	// Use Inclusion-Exclusion principle
	const T bitmask = (static_cast<T>(1) << prime.size());
	for (T bit = 0; bit < bitmask; ++bit) {
	T curr = 1, flag = 1;
	// Let bit = b1 b2 ... bm, where bit is in [0, 2^ps), ps = prime.size()
	// bi = 1, if currently we need to process prime[i]; otherwise skip prime[i]
	for (size_t i = 0; i < prime.size(); ++i) {
	if ((bit >> i) & 1) { // bi = 1
	curr *= prime[i];
	flag *= -1;
	}
	}
	// flag = 1 when there is even number of bi = 1
	// otherwise, flag = -1
	ans += flag * (value / curr);
	}
	// Example: x = 6 = 2 ^1 * 3 ^1
	// prime = {2, 3}, and bit will be {00, 01, 10, 11}
	// bit = 00: ans += 1 * (value / 1)
	// bit = 01: ans += (-1) * (value / 3)
	// bit = 10: ans += (-1) * (value / 2)
	// bit = 11: ans += 1 * (value / (2 * 3))
	return ans;
	}
	/******************************************************************************************************************
	* Given two positive integers n and x, where x = p1 ^ k1 * p2 ^ k2 * ... * pm ^ km, and p1, p2, ... pm are primes,
	* this function will return the number of values y that y <= n and gcd(y, x) = 1, where
	* gcd(y, x) is the greatest common divisor of y and x.
	* Namely, it will find how many values are less than n are coprime with x.
	*
	* Example: Let n = 12 and x = 6, the result is 4 ({1, 5, 7, 11})
	* Note: When n = x, the result is equal to phi(n), where phi is Euler Totient Function.
	*
	* Parameters:
	* 1. prime: prime list {p1, p2, ... , pm}
	* 2. value: n
	* 3. x: as the x described above
	* 4. phi_x: phi(x), where phi is Euler Totient function.
	*
	* Return value:
	* Number of values that are less than n and coprime with x.
	******************************************************************************************************************/
	template <typename T>
	T coprimeCount(const std::vector<T>& prime, const T value, const T x, const T phi_x)
	{
	return coprimeCount(prime, value % x) + (value / x) * phi_x;
	}
	// return the kth value that is relative prime to p and larger than x
	int kthValue(const std::vector<int>& prime, const int x, const int p, int k)
	{
	#define INF 0x3fffffff
	const auto ex_prime = factorization(prime, p);
	k += coprimeCount(ex_prime, x);
	int left = 1, right = INF;
	// binary search
	while (left < right) {
	const auto middle = (left + right) / 2;
	if (coprimeCount(ex_prime, middle) < k) {
	left = middle + 1;
	}
	else {
	right = middle;
	}
	}
	return right;
	}
	// return the kth value that is relative prime to p and larger than x
	int kthValue(const std::vector<int>& prime, const std::vector<int>& phi, const int x, const int p, int k)
	{
	#define INF 0x3fffffff
	const auto ex_prime = factorization(prime, p);
	k += coprimeCount(ex_prime, x, p, phi[p]);
	int left = 1, right = INF;
	// binary search
	while (left < right) {
	const auto middle = (left + right) / 2;
	if (coprimeCount(ex_prime, middle, p, phi[p]) < k) {
	left = middle + 1;
	}
	else {
	right = middle;
	}
	}
	return right;
	}
	int main(int argc, char *argv[])
	{
	#define MAX_VALUE 1000000
	std::vector<int> prime, phi;
	const auto prime_count = GeneratePrime(prime, MAX_VALUE);
	euler_totient(phi, MAX_VALUE, prime);
	int t;
	while (scanf("%d", &t) == 1) {
	for (int i = 0; i < t; ++i) {
	int x, p, k;
	if (scanf("%d %d %d", &x, &p, &k) == 3) {
	printf("%d\n", kthValue(prime, x, p, k));
	//printf("%d\n", kthValue(prime, phi, x, p, k)); // much slower...
	}
	}
	}
	return 0;
	}

view raw codeforces-920G.cpp hosted with ❤ by GitHub

[後記]

當初本來想要用 Euler Totient Function 來解，雖然有找到下面這個 Lemma：

Lemma: 令 $x = d \times p$，其中 $d$ 為正整數，則 $G(p, x) = d \times \phi(p)$，其中 $\phi$ 為 。

不過問題會出在當 $x = p \times d + r,\ 0 < r < p$，此時 $G(p, x) = d \times \phi(p) + \phi(p, r)$，其中 $\phi(p, r) = |\{n\ |\ n \in \phi(p)\ \&\ n < r\}|$，也就是 $\phi(p)$ 中跟 $p$ 互質且比 $r$ 小的個數。

後來 Google 找了一陣子找不到數論中有什麼方法可以算 $\phi(p, r)$，因此最後還是得回歸排容原理來算 = =...上面的實作雖然也有把這方法做出來，但說實在的還比不用更慢...