[程式] 排容原理實作 (inclusion-exclusion principle implementation)

由於先前這一篇會用到排容原理，因此就花了點時間研究排容原理要怎麼實作比較好。不過這玩意兒因為只是個原則，實作方法會根據要處理的問題有些許差別，這邊的實作針對的是以下的問題：

給定兩個正整數 $n$ 跟 $x$，要找出有多少正整數 $y$ 使得 $y < n \ \& \ gcd(y, x) = 1$，其中 $gcd$ 代表最大公因數


詳細的說明有放在實作中了，這邊就不多做解釋，實作如下：

	/*
	* =====================================================================================
	*
	* Filename: inclusion_exclusion.cpp
	*
	* Description: Implementation of inclusion-exclusion principle
	*
	* Version: 1.0
	* Created: 2018/02/07 (yyyy/mm/dd)
	* Revision: none
	* Compiler: g++
	*
	* Author: lionking
	* Organization: None
	*
	* =====================================================================================
	*/
	#include <cstdio>
	#include <vector>
	/*****************************************************************************************
	* 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();
	}
	/****************************************************************************************************
	* 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 = 14 and x = 6, the result is 3 ({1, 3, 5})
	* 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;
	}
	int main(int argc, char *argv[])
	{
	#define MAX_VALUE 100
	std::vector<int> prime;
	GeneratePrime(prime, MAX_VALUE);
	constexpr int p = 14, x = 6;
	const auto ex_prime = factorization(prime, p);
	printf("Number of values less than %d and coprime with %d: %d\n", x, p, coprimeCount(ex_prime, x));
	// Sample output:
	// "Number of values less than 6 and coprime with 14: 3"
	return 0;
	}

view raw inclusion_exclusion.cpp hosted with ❤ by GitHub