백준 26309번: Etched Emerald Orbs

https://www.acmicpc.net/problem/26309

26309번: Etched Emerald Orbs

---

문제의 상황은 주어지는 자연수 $k$ 에 대해 다음 조건을 만족하는 서로 다른 자연수 $x$, $y$ 를 찾으라는 것입니다.

$$ \frac{2}{k} = \frac{1}{x} + \frac{1}{y} $$

만약에 이런 $x$, $y$가 존재하지 않으면 $-1$ 을 출력해야 하는데, 다행히 이런 $x$, $y$ 는 항상 존재합니다.

---

존재성을 확인하기 위해 식을 변형합니다.

$$ 2xy = k(x+y) $$

$$ \rightarrow (2x-k)y = kx $$

$$ \rightarrow (2x-k)2y - 2kx + k^{2} = k^{2} $$

$$ \rightarrow (2x-k)(2y-k) = k^{2} $$

이제 $x$ 와 $y$ 는 서로 다른 자연수이므로 $k^{2}$ 은 최소한 $k$ 가 아닌 다른 약수를 가져야 합니다.

 

그런데 $k^{2}$ 꼴의 자연수 중 1을 제외하여 약수의 개수가 제일 적은 것은 소수의 제곱입니다.

 

그리고 소수 $p$에 대해 $p^{2}$ 은 $p$ 를 제외한 약수를 2개 가지고 있습니다. 따라서 $x$ 와 $y$ 는 적어도 한 가지 이상 존재합니다.

---

문제를 풀기 위해 $k^{2}$ 의 어떤 약수 $d$ 를 가정하면

$$ x = \frac{k+d}{2} $$

$$ y = \frac{k + \frac{k^{2}}{d}}{2} $$

라고 할 수 있습니다.

 

이때 $x$, $y$ 는 모두 자연수이며 $x < y$ 이므로 다음을 확인할 수 있습니다.

 

1. $d < k$ 이다.

2. $k$ 와 $d$, $\frac{k^{2}}{d}$ 는 모두 parity 가 일치한다.

 

이런 조건을 만족하는 최대의 $d$ 를 구해주면 문제가 해결됩니다. 참고로, $x < y$ 에 의해 $d$ 는 자연수입니다.

---

그런데 $d$는 $k^{2}$ 의 약수이며, $k$는 최대 $4 \times 10^{18}$ 까지 커질 수 있는 값입니다.

 

이 정도 범위의 정수를 소인수분해하기 위해서는 Pollard-rho 알고리즘이 필요하겠습니다.

 

또한, $k^{2}$ 의 약수를 순회하는 시간이 클 수 있습니다. 공식 해설에 따르면 최악의 경우는 약 4천만개 정도라고 언급됩니다. $k$ 보다 큰 약수를 탐색하지 않는 걸로 커팅하는 게 좋습니다.

 

시간 복잡도는 $k$ 의 소인수분해와, $k^{2}$ 의 약수의 개수의 합으로, 대략 $O(k^{\frac{1}{4}} + k^{O(1)})$ 이 되겠습니다. 약수의 개수의 asymptotic 한 분석은 후술하겠습니다.

 

답이 128bit 정수 범위에서 나오므로 출력 함수를 따로 만들어야 하겠습니다. 제 코드는 다음과 같습니다.

#include <iostream>
#include <cstdlib>
#include <numeric>
#include <algorithm>
#include <vector>
#include <chrono>
#include <random>
#include <cassert>

void print(__int128 x) {
	// reference : codeforces.com/blog/entry/75044
	if (x < 0) {
		putchar('-');
		x = -x;
	}
	if (x > 9) print(x / 10);
	putchar(x % 10 + '0');
}

template <typename Type>
class Random
{
	/* implementation : https://blog.naver.com/devmachine/221240867011
	*
	* Type : int, long long, double
	* Sample code :
	* (1) Random<int> gen(-1000, 1000); // integer between -1000 and 1000
	* 
	* (2) Random<double> gen(0.0, 1.0); // floating number between 0.0 and 1.0
	* 
	* (3) Random<__int64> gen(0, 100000000000); // large integer between 0 and 1e11
	* 
	* (4) Random<int> gen(-1000, 1000, 42); // Generate by user seed
	*/
public:
	using Distribution = typename std::conditional<std::is_floating_point_v<Type>,
	std::uniform_real_distribution<Type>,
	std::uniform_int_distribution<Type>>::type;
	
	using Engine = typename std::conditional<(std::is_integral_v<Type> && sizeof(Type) > 4),
	std::mt19937_64,
	std::mt19937>::type;
	
	Random(
		Type min = (std::numeric_limits<Type>::min)(), 
		Type max = (std::numeric_limits<Type>::max)(),
		typename Engine::result_type seed = std::random_device()()
		)
	: _engine(seed), _distribution((std::min)(min, max), (std::max)(min, max))
	{
	}
	
	Type operator()()
	{
		return static_cast<Type>(_distribution(_engine));
	}
	
private:
	Engine _engine;
	Distribution _distribution;
};
Random<long long> gen(0, 1e18, std::chrono::steady_clock::now().time_since_epoch().count());

namespace Primality {
	using u32 = unsigned int;
	using u64 = unsigned long long;
	using i128 = __int128;
	
	bool sqrt_test (u64 N) {
		/* input value : N is odd and small 
		* return value : {false := composite; true := prime}
		* time complexity : sqrt(N)
		* restriction : N < 10^9
		*/
		for (u32 i = 3; i * i <= N; i += 2) {
			if (N % i == 0) return false;
		}
		return true;
	}
	
	template <bool BigInt = false>
	u64 power_modular (u64 X, u64 Y, const u64 M) {
		u64 result = 1;
		while (Y) {
			if (Y % 2) {
				result = static_cast<typename std::conditional<BigInt, i128, u64>::type>(result) * X % M;
			}
			X = static_cast<typename std::conditional<BigInt, i128, u64>::type>(X) * X % M;
			Y /= 2;
		}
		return result;
	}
	
	template <bool BigInt = false>
	bool Miller_Rabin_Internal (u64 N, u64 K) {
		/* input value : N is odd, K < N
		* return value : {false := composite; true := probable-prime}
		*
		* explanation :
		* calculate k^d, k^(2d), k^(4d), ... , k^(2^s-1 * d)
		* where 2^s * d = N - 1
		* 
		* Then, N is probably prime if k^d = +1/-1 or k^(2^k * d) = -1 mod N
		* 
		* time complexity : O(log^2 N)
		*/
		u64 d = N - 1;
		for (; d % 2 == 0; d /= 2) {
			if (power_modular<BigInt>(K, d, N) == N - 1) {
				return true;
			}
		}
		if (u64 temp = power_modular<BigInt>(K, d, N); temp == 1 || temp == N - 1) {
			return true;
		}
		else {
			return false;
		}
	}
	
	template <bool BigInt = false>
	bool Miller_Rabin_test (u64 N) {
		/* input value : N is odd and larger than list[i]
		* return value : {false := composite; true := prime}
		*
		* explanation :
		* Originally, Miller_Rabin_Test can be wrong by 4^-k where k is number of rounds
		* But, for small N, there is a deterministic way which has 100% accuracy
		*
		* time complexity : O(k log^2 N)
		*/
		std::vector<u64> list;
		if (BigInt) list =  {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
		else list = {2, 7, 61};
		
		for (auto& i : list) {
			if (N % i == 0) return false;
			if (!Miller_Rabin_Internal<BigInt>(N, i)) return false;
		}
		return true;
	}
	
	bool is_Prime (u64 N) {
		if (N < 2) return false;
		if (N == 2 || N == 3) return true;
		if (N % 6 != 1 && N % 6 != 5) return false;
		if (N <= 1e4) return sqrt_test(N);
		if (N <= ((u64)1 << 32)) return Miller_Rabin_test(N);
		return Miller_Rabin_test<true>(N);
	}
}

namespace Factorize {
	void Pollard_Rho_Internal (long long N, std::vector<long long>& result) {
		/* result := <p1, p2, ... > where product of all elements is equal to N
		* (But there can be the same elements in result vector)
		*
		* Caution :
		* (1) include <algorithm>, <vector>, <numeric>, <cstdlib>, <random>, <chrono>
		* (2) for u64 or larger types, implement your own gcd and abs
		* (3) for i64 or larger types, miller-rabin should be safe "enough"
		*     (especially for multiplication)
		*
		* time complexity : approximately O(n^(1/4)), but actual proof is still open
		* 
		* implementation : https://blog.naver.com/PostView.naver?blogId=jinhan814&logNo=222141831551
		*/
		if (N == 1) return;
		if (Primality::is_Prime(N)) {
			result.push_back(N);
			return;
		}
		if (N % 2 == 0) {
			result.push_back(2);
			Pollard_Rho_Internal(N / 2, result);
			return;
		}
		
		// cycle-detection on functional graph (hare and tortoise)
		long long A = 0, B = 0, C = 0, G = N;
		do{
			auto f = [&](long long X) -> long long {
				return ((__int128)X * X + C) % N;
			};
			
			if (G == N) {
				A = B = gen() % (N - 2) + 2;
				C = gen() % 20 + 1;
			}
			
			A = f(A);
			B = f(f(B));
			G = std::gcd(std::abs(A - B), N);
		}while (G == 1);
		
		Pollard_Rho_Internal(G, result);
		Pollard_Rho_Internal(N / G, result);
	}
	
	std::vector<long long> Pollard_Rho (long long N) {
		std::vector<long long> result;
		Pollard_Rho_Internal(N, result);
		std::sort(result.begin(), result.end());
		return result;
	}
	
	std::vector<std::pair<long long, long long>> Multiplicity (long long n) {
		auto pF = Pollard_Rho(n);
		
		std::vector<long long> uPF;
		uPF.assign(pF.begin(), pF.end());
		uPF.erase(unique(uPF.begin(), uPF.end()), uPF.end());
		
		std::vector<std::pair<long long, long long>> result;
		for (auto& i : uPF) {
			result.push_back(std::make_pair(i,
				upper_bound(pF.begin(), pF.end(), i) -
				lower_bound(pF.begin(), pF.end(), i))); // prime number and its multiplicity
		}
		
		return result;
	}
}

using namespace std;
int main() {
	long long k;
	cin >> k;
	
	auto temp = Factorize::Multiplicity(k);
	for (auto& i : temp) {
		i.second *= 2;
	}
	
	auto exp = [](__int128 x, long long y) -> __int128 {
		__int128 result = 1;
		while(y--) result *= x;
		return result;
	};
	
	long long D = 1;
	std::function<void(int, int, long long, __int128)> Internal = [&](int pos, int last, long long power, __int128 cur) -> void {
		cur *= exp(temp[pos].first, power);
        if (cur >= k) return;
		if (pos == last) {
			if (D < cur) {
				// divisor of k^2, less than k, larger than current value D
				if (k % 2 == cur % 2) {
					if (((__int128)k * k / cur) % 2 == k % 2) {
						D = cur;
					}
				}
			}
		}
		else
			for (long long i = 0; i <= temp[pos + 1].second; i++) {
				Internal(pos + 1, last, i, cur);
			}
	};
	for (long long i = 0; i <= temp[0].second; i++) Internal(0, temp.size() - 1, i, 1);
	
	__int128 y = (__int128)k * k / D;
	y = (y + k) / 2;
	
	__int128 x = ((__int128)k + D) / 2;
	
	assert((__int128)k * k == (2 * x - k) * (2 * y - k));
	
	print(x);
	std::cout << ' ';
	print(y);
	
	return 0;
}

(C++20, 532ms, 2036KB, 제출번호 69669250)

---

약수의 개수에 관한 글들을 읽어보았는데, math.stackexchange.com 의 글이 좋은 것 같습니다.

 

OEIS 에서 highly composite numbers 를 찾아보니, $k$ 가 다음 값으로 주어질 때

$$ k = 3066842656354276800 = 2^{6} \times 3^{3} \times 5^{2} \times 7^{2} \times 11 \times 13 \times 17 \times 19 \times 23 \times 29 \times 31 \times 37 \times 41 $$

$k^{2}$의 약수의 개수가 44778825 개가 나옵니다. 아마 공식 해설은 이걸 언급한 것 같습니다.