====== 제곱수의 합 ======
* 자연수 n을 최소갯수의 제곱수의 합으로 표현하는 것에 관련된 이론들이다.
* 크게 두가지 주제가 있다.
* n을 제곱수의 합으로 표현할 때 필요한 제곱수의 최소 갯수를 찾기
* n을 최소 갯수의 제곱수의 합으로 표현할 수 있는 제곱수를 찾기
* 'n을 k개의 제곱수의 합으로 표현하는 방법의 갯수를 찾기' 에 관해서도 이론들이 있는데.. 아직 이것을 필요로 하는 문제는 본적이 없으니 패스
* 필요하면 [[https://en.wikipedia.org/wiki/Sum_of_squares_function]] 를 참고하자.
* 참고자료
* [[https://00ad-8e71-00ff-055d.tistory.com/63]] - 제곱수의 최소 갯수를 찾는 것에 대한 설명.
* [[https://infossm.github.io/blog/2019/10/18/sum-of-squares/]] - 제곱수의 최소 갯수를 찾는 것 + 제곱수들을 실제로 찾는 것에 대한 상세한 설명 (파이썬으로 구현된 코드 포함)
===== 제곱수의 합으로 표현할 때 필요한 제곱수의 최소 갯수를 찾기 =====
* 우선 간단히 알아둘 것은, 모든 자연수는 최대 4개의 제곱수의 합으로 표현 가능하다
* 그러면 최소 갯수를 찾기 위해서는 1개,2개,3개로 표현 가능한지 여부를 각각 확인하면 된다.
* 1개의 제곱수로 표현 가능? : 당연히 n이 자체로 제곱수여야 한다. O(1)에 확인 가능
* 2개의 제곱수로 표현 가능?
* 단순한 방법은 n이하의 모든 제곱수 x^2에 대해서 n-x^2 이 제곱수인지 확인해보면 된다. O(n^(1/2))에 구할수 있다.
* 큰 n에 대해서는, 페르마의 두 제곱수 정리 (=페르마의 크리스마스 정리)를 이용한다
* n을 소인수분해했을때, 4k+3 꼴의 소인수의 지수가 2의 배수로 존재하지 않으면 2개의 제곱수의 합으로 표현할수 있다. 역도 성립.
* 관련된 소정리들
* a와 b가 둘다 2개의 제곱수의 합으로 표현되는 수라면 a*b도 두개의 제곱수로 표현 가능하다 ([[https://en.wikipedia.org/wiki/Brahmagupta%E2%80%93Fibonacci_identity|디오판투스 항등식]])
* 2 또는 4k+1 형태의 소수는 2개의 제곱수의 합으로 표현 가능하다. 4k+3 형태의 소수는 불가능하다.
* 다양한 증명법들이 있다. 보통은 오일러의 증명으로 설명하지만, [[https://mathoverflow.net/questions/31113/zagiers-one-sentence-proof-of-a-theorem-of-fermat/299696#299696|그림을 이용한 직관적인 증명]]도 있다.
* 시간복잡도는 소인수분해의 시간복잡도에 바운드된다. 폴라드로 방법으로 소인수분해를 하면 O(n^(1/4))에 처리할수 있다.
* n이 크지 않아서 O(n^(1/2)) 소인수분해 알고리즘을 돌리는게 가능할정도라면, 그냥 위에서 말한 O(n^(1/2)) 방법을 쓰면 된다.
* 3개의 제곱수로 표현 가능?
* n=4^a(8b+7) 꼴이 아니면 3개의 제곱수로 표현 가능하다. (르장드르의 세 제곱수 정리).
* 증명은 꽤 복잡하다고 한다
* 4개의 제곱수로 표현 가능?
* 모든 자연수는 4개의 제곱수로 표현 가능하다 (라그랑주의 네 제곱수 정리)
* 따라서 3개 이하로 표현 불가능하다면, 4개가 최소이다
* 결국 최소 갯수를 찾을때, 시간이 걸리는 부분은 2개의 제곱수로 표현 가능한지 여부를 확인하는 부분이고 (O(n^(1/4))), 나머지는 O(1)에 구할수 있다. 총 O(n^(1/4))
==== 관련 문제 ====
* [[ps:problems:boj:17626]] (n<50,000)
* [[ps:problems:boj:1699]] (n<100,000)
* [[ps:problems:boj:17633]] (n<1,000,000,000,000,000,000) - n범위가 크므로, O(n^(1/4)) 소인수분해 알고리즘을 써서 찾아야 한다
* [[ps:problems:boj:4913]] 페르마의 두 제곱수 정리에 따라, 4k+1꼴의 소수들의 갯수를 세면 된다.
===== n을 최소 갯수의 제곱수의 합으로 표현하기 =====
* 1개의 제곱수의 합으로 표현하는 것은 그냥 sqrt(n)을 구하면 된다.
* 2개의 제곱수의 합으로 표현하는 것은, 조금 복잡하므로 [[#2개의 제곱수의 합으로 표현]]에서 따로 설명. 시간복잡도는 O(n^(1/4))이다
* 3개의 제곱수의 합으로 표현하는 것
* 한개의 제곱수를 고정시키고, 나머지 두개의 제곱수를 [[#2개의 제곱수의 합으로 표현]] 방법을 사용해서 찾는다.
* 즉 (n-1^2), (n-2^2), (n-3^2), ... 에 대해서 2개의 제곱수의 합으로 표현가능할때까지 시도해본다는 말이다
* 무식해보이고 많은 시간이 걸릴것 같지만, 그렇지도 않다.
* N 이하의 자연수중 두 제곱수의 합으로 표현될 수 있는 수의 개수는 대략 $ K \cdot \frac{n}{\sqrt{\ln n}} $ 이라는 것이 알려져있다.
* 따라서 n-x^2이 2개의 제곱수의 합으로 표현될 확률은 O(1/sqrt(logn)) 이므로, 평균적으로 O(sqrt(logn))번의 시도만에 찾을수 있다.
* 4개의 제곱수의 합으로 표현하는 것
* 4개의 제곱수의 합으로 표현하는 것은 4^a(8b+7) 형태일때이다. 제곱수 하나를 (2^a)^2 로 고정시키고 4^a(8b+6)에 대해서 3개의 합으로 표현하는 방법을 찾아주면 된다.
* 좀더 효율적으로는 4^a(8b+6)이 아니라, (8b+6)을 세제곱수의 합으로 표현한 뒤에 각각에 4^a을 곱해주면 된다.
* 결국 총 시간복잡도는 O(n^(1/4) * (logn)^(1/2)) 이 된다
==== 2개의 제곱수의 합으로 표현 ====
* n을 소인수분해해서, 4k+1꼴의 p를 2개의 제곱수의 합으로 표현하는 방법을 찾는다면, n을 2개의 제곱수의 합으로 나타내는 것은 ([[https://en.wikipedia.org/wiki/Brahmagupta%E2%80%93Fibonacci_identity|디오판투스 항등식]])을 통해서 쉽게 찾을수 있다.
* 4k+1꼴의 p를 2개의 제곱수의 합으로 표현하는 방법을 찾는 알고리즘은 몇가지가 있다.
* 우선 p가 작을 경우에는 그냥 브루트포스 방식으로 구해도 된다. 2개의 제곱수로 표현 가능여부를 찾을때도 언급했지만, n이하의 모든 제곱수 x^2에 대해서 n-x^2 이 제곱수인지 확인하는 방식으로 처리해도 O(n^(1/2))에 계산 가능하다. n의 범위가 작은 경우 (n의 소인수분해를 O(n^(1/2))으로 처리해도 돌아가는 경우)에는 굳이 소인수분해해서 각 소인수들에 대해 일일히 제곱수의 합으로 표현하는 방법을 찾지 말고, 그냥 이 방법을 바로 n에 대해서 돌리자
* p가 큰 경우에는, 가우스 정수를 사용하는 방법과 Hermite-Serret 알고리즘을 사용하는 방법이 있다.
* 가우스 정수를 사용하는 방법으로, 증명 + 알고리즘 + 코드 까지 가장 친절하게 설명된 자료는 [[https://infossm.github.io/blog/2019/10/18/sum-of-squares/]] 이다. 해당되는 백준 문제 [[ps:problems:boj:17646]]의 출제자가 직접 쓴 글이다.
* Hermite-Serret 알고리즘 (또는 Cornacchia 알고리즘)은[[https://en.wikipedia.org/wiki/Fermat's_theorem_on_sums_of_two_squares#Algorithm|위키피디아]]에 소개된 방법이다. 구현하기에는 더 간단한데, 증명은 좀 찾기가 어렵다
=== 가우스 정수를 이용한 알고리즘 ===
* 알고리즘 자체는 간단하게 요약된다. x^2 % p = -1 이 되는 quadratic residue x를 찾은 다음에, p와 x+i의 최소공약수를 구한다. 이렇게 최소공약수 a+bi를 찾는다면 a^2+b^2=p가 성립한다.
* %% quadratic residue를 찾는 것은.. 1≤a≤p-1인 a에 대해서 a^((p-1)/2)%p 이 절반은 1, 절반은 -1이라는 사실을 이용해서 쉽게 찾을수 있다. a를 1부터 증가시키면서 저 값을 구해보면 평균 2번만에 mod p가 -1이 되는 a를 찾을수 있다. x= a^((p-1)/4) 로 잡으면 된다 %%
*
def quadratic_residue(p):
k = (p - 1) // 4
return next(x for a in range(2, p - 1) if (x := pow(a, k, p)) * x % p == p - 1)
* p와 x+i의 최소공약수라는 개념이 생소한데.. 이것을 이해하려면 가우스정수에 대한 배경지식이 필요하다. a와 b가 모두 정수일때 a+bi를 가우스 정수라고 하는데.. 여기에서 나머지연산, gcd연산을 정의해 줄 수가 있다. 자세한 설명은 링크를 참고하고, 코드만 적어보면 이런식이다
*
def gcd(a: GaussianInt, b: GaussianInt) -> GaussianInt:
while b.real != 0 or b.imag != 0:
a, b = b, a % b
return a
class GaussianInt:
__slots__ = ('real', 'imag')
def __init__(self, real, imag=0):
self.real = real
self.imag = imag
def __sub__(self, other):
return GaussianInt(self.real - other.real, self.imag - other.imag)
def __mul__(self, other):
return GaussianInt(
self.real * other.real - self.imag * other.imag,
self.real * other.imag + self.imag * other.real,
)
def __floordiv__(self, other):
norm = other.real * other.real + other.imag * other.imag
r = self.real * other.real + self.imag * other.imag
i = self.imag * other.real - self.real * other.imag
return GaussianInt(round(r / norm), round(i / norm))
def __mod__(self, other):
return self - self // other * other
* 빌트인 complex 클래스를 상속받아서 구현하면, 덧셈,뺄셈,곱셈 등은 이미 정의되어있으므로, 정수나눗셈과 나머지연산만 만들어주면 되기 때문에 편리하다. 정수나눗셈도 일반 나눗셈의 계산 결과를 round만 해주면 되므로 더 쉽게 구현 가능하다. 하지만 문제는 complex 클래스의 실수부와 허수부의 값은 부동소숫점으로 저장되기 때문에, 값이 커지면 실수 오차가 생긴다는 것이다.. 그래서 아쉽지만 complex 클래스를 상속받아서 쓰는것은 불가능하고, 결국은 바닥부터 다 구현해야 했다..
* 위에서 정의한 가우스정수 클래스를 이용하면, n을 2개을 제곱수의 합으로 나타내는 함수는 아래처럼 구현 가능하다
*
def sum_of_two_squares_if_possible(n) -> tuple[int,int] | None:
factoriaztion = numtheory.prime_factorization(n)
primes = [p for p, e in factoriaztion.items() if e % 2 == 1]
if any(p % 4 == 3 for p in primes):
return None
ans = GaussianInt(math.isqrt(n // math.prod(primes)))
for p in primes:
if p == 2:
ans *= GaussianInt(1, 1)
else:
x = quadratic_residue(p)
ans *= gcd(GaussianInt(p), GaussianInt(x, 1))
return (abs(ans.real), abs(ans.imag))
* 2개의 제곱수의 합로 표현 불가능 하면 None를 리턴한다
* quadratic residue를 구하는 것이나 각 연산들은 모두 O(1)이고, gcd를 구하는 것에만 O(logp)가 걸린다. 총 시간복잡도는 소인수분해의 시간복잡도에 바운드된다.
* 이렇게 코드를 작성해서 제출하면 [[ps:problems:boj:17646]]를 통과할수 있다.
=== Hermite-Serret 알고리즘 ===
* Hermite-Serret 알고리즘 이라는 이름 없이 그냥 알고리즘 설명만 되어있는 자료도 많다. 방법은 다음을 참고
* [[https://en.wikipedia.org/wiki/Fermat's_theorem_on_sums_of_two_squares#Algorithm]]
* [[https://projecteuler.chat/viewtopic.php?p=30262&sid=65c363ba34e03e9b3ad2be21c1abe1ea#p30262]]
* 대충 결과적으로는 Gaussian GCD를 이용하는것과 동일한 원리라고 하는것 같다..
* Cornacchia's 알고리즘은 x^2+dy^2=p 을 만족하는 x와 y를 찾는 알고리즘이다. Hermite-Serret의 일반화된 버전이라고 생각해도 될거 같은데.. Hermite-Serret보다는 이쪽이 더 유명한것 같다.. (따로 위키페이지도 있고)
* [[https://en.wikipedia.org/wiki/Cornacchia%27s_algorithm]]
* [[https://math.stackexchange.com/questions/5877/efficiently-finding-two-squares-which-sum-to-a-prime/5895#5895]]
* 여기에서는 (x^2)%p=-d 인 x를 찾기 위해서, Tonelli-Shanks Algorithm 을 사용해야 하는 것으로 설명되어있다.
* 하지만 우리가 풀려는 문제는 d=1로 고정되어있으므로, 앞에서 설명했던 방식으로 quadratic residue x를 간단하게 찾을수 있다.
* 이 방법을 사용하면, 복잡한 가우스정수 클래스 없이도 코드를 좀더 간단하게 짤수 있다.
*
def hermite_serret(p):
if p == 2:
return (1, 1)
q = quadratic_residue(p)
sqrt_p = math.sqrt(p)
c = None
while q > 0:
if q < sqrt_p:
if c is None:
c = q
else:
return (c, q)
p, q = q, p % q
def sum_of_two_squares_if_possible(n):
factoriaztion = prime_factorization(n)
primes = [p for p, e in factoriaztion.items() if e % 2 == 1]
if any(p % 4 == 3 for p in primes):
return None
a, b = math.isqrt(n // math.prod(primes)), 0
for p in primes:
c, d = hermite_serret(p)
a, b = a * c + b * d, a * d - b * c
return (abs(a), abs(b))
==== 관련 문제 ====
* [[ps:problems:boj:23287]] (n<1,000,000,000,000,000,000) 이지만, 소인수의 크기가 100,000 이하라는 제한이 있어서, trial division을 이용한 소인수분해와, brute force를 이용한 두 제곱수의 합 표현 찾기 방법을 써서도 풀수 있다.
* [[ps:problems:boj:17646]] (n<1,000,000,000,000,000,000)