[백준 17646] 제곱수의 합 2 (More Huge)

문제 링크

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

 

사용 알고리즘

• 폴라드 로(Pollard's rho)
• 밀러 라빈 소수 판별법(Miller Rabin Primality Test)
• Cornacchia's algorithm
• Tonelli-Shanks algorithm

 

풀이 과정

이 문제에서 제곱수의 개수를 구하는 방법은 '제곱수의 합 (More Huge)'와 동일합니다.

따라서 이 부분은 생략하고 설명하겠습니다.

 

먼저 제곱수의 개수가 4인 경우부터 생각해 보겠습니다.

제곱수의 개수가 4인 경우는 \(n = 4^a(8b + 7)\) 입니다.

그리고 \(4^0(8b + 6)\)은 세 제곱수로 나타낼 수 있기 때문에 이것을 \(x^2 + y^2 + z^2\)이라 하면 다음과 같이 나타낼 수 있습니다.

$$n = {(2^a)}^2(x^2 + y^2 + z^2 + 1)$$

$${(2^a)}^2(x^2 + y^2 + z^2 + 1) = (x * 2^a)^2 + (y * 2^a)^2 + (z * 2^a)^2 + {(2^a)}^2$$

따라서 4개인 경우를 해결하기 위해서는 3개인 경우인 \({n\over4^a} - 1\)의 결과 값이 필요합니다.

 

다음으로 제곱수의 개수가 3인 경우를 살펴보겠습니다.

먼저 \(x\)이하의 자연수중 두 제곱수의 합으로 표현될 수 있는 수는 대략 \(K * {x \over {\sqrt{ln x}}}\)라는 것을 란다우가 밝혀냈으며, 여기서 \(K\)는 란다우 라마누잔 상수입니다.

이것을 활용해서 임의의 \(t\)를 잡고, \(n - t^2\)을 했을 때, 이 수가 두 제곱수로 나올 확률이 그렇게 나쁘지 않다는 것을 알 수 있습니다.

따라서 3개인 경우는 \(t\)를 증가시키면서 \(n - t^2\)이 두 개로 표현될 수 있는 수로 나올 때까지 반복하면,

\(x^2 + y^2 + t^2\)으로  3개인 경우도 해결할 수 있습니다.

 

마지막으로 제곱수의 개수가 2인 경우입니다.

2인 경우에는 \(x^2 + dy^2 = m\)에서 \(x\)와 \(y\)를 찾아주는 알고리즘인 Cornacchia's algorithm을 사용해서 해결할 수 있고,

중간에 계산하는 과정에서 \(r^2 \equiv -1 \ (mod \ p)\)를 구하는 부분이 있는데, 이 부분을 계산할 때 시간초과가 나기 때문에 Tonelli-Shanks algorithm을 사용해서 시간을 단축할 수 있습니다.

 

이렇게 해서 문제를 해결할 수 있습니다.

 

사담

이 문제도 상당히 어려웠던 문제같습니다.

Cornacchia's algorithm과 Tonelli-Shanks algorithm을 모른다면 풀 수 없었을 것 같은데, 

다른 분의 블로그에서 두 알고리즘의 존재를 알게 되어서 겨우 풀 수 있었던 것 같습니다.

 

Cornacchia's algorithm과 Tonelli-Shanks algorithm은 잘 모르는 사람이 많을 것 같아서 코드를 한 번 남겨볼까합니다.

다음 글에서 뵙겠습니다.

읽어주셔서 감사합니다.

 

코드

• Cornacchia's algorithm

Two_sqrts cornacchia(ll n) {
	ll r1 = n, r2;

	r2 = tonelli_shanks(r1 - 1, r1);

	while ((__int128)(r1) * r1 >= n) {
		ll t = r1 % r2;
		r1 = r2;
		r2 = t;
	}

	ll x = r1, y = sqrt(n - r1 * r1);
	return { x, y };
}

 

• Tonelli-Shanks algorithm

ll tonelli_shanks(ll n, ll p) {
	ll s = 0;
	ll q = p - 1;

	while ((q & 1) == 0) { q /= 2; ++s; }
	if (s == 1) {
		ll r = poww(n, (p + 1) / 4, p);
		if (mul(r, r, p) == n) return r;
		return 1;
	}
	
	ll z = 1;
	while (poww(++z, (p - 1) / 2, p) != p - 1);

	ll c = poww(z, q, p);
	ll r = poww(n, (q + 1) / 2, p);
	ll t = poww(n, q, p);
	ll m = s;
	while (t != 1) {
		ll tt = t;
		ll i = 0;
		while (tt != 1) {
			tt = mul(tt, tt, p);
			++i;
			if (i == m) {
			    return 1;
			}
		}
		ll b = poww(c, poww(2, m - i - 1, p - 1), p);
		ll b2 = mul(b, b, p);
		r = mul(r, b, p);
		t = mul(t, b2, p);
		c = b2;
		m = i;
	}
	if (mul(r, r, p) == n) return r;
	return 1;
}