제곱수의 합

자연수 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도 두개의 제곱수로 표현 가능하다 (디오판투스 항등식)
  - 2 또는 4k+1 형태의 소수는 2개의 제곱수의 합으로 표현 가능하다. 4k+3 형태의 소수는 불가능하다.
- 다양한 증명법들이 있다. 보통은 오일러의 증명으로 설명하지만, 그림을 이용한 직관적인 증명도 있다.
- 시간복잡도는 소인수분해의 시간복잡도에 바운드된다. 폴라드로 방법으로 소인수분해를 하면 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))

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개의 제곱수의 합으로 나타내는 것은 (디오판투스 항등식)을 통해서 쉽게 찾을수 있다.
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/ 이다. 해당되는 백준 문제 제곱수의 합 2 (More Huge)의 출제자가 직접 쓴 글이다.
- Hermite-Serret 알고리즘 (또는 Cornacchia 알고리즘)은위키피디아에 소개된 방법이다. 구현하기에는 더 간단한데, 증명은 좀 찾기가 어렵다

가우스 정수를 이용한 알고리즘

알고리즘 자체는 간단하게 요약된다. 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)가 걸린다. 총 시간복잡도는 소인수분해의 시간복잡도에 바운드된다.
이렇게 코드를 작성해서 제출하면 제곱수의 합 2 (More Huge)를 통과할수 있다.

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))

목차

제곱수의 합

제곱수의 합으로 표현할 때 필요한 제곱수의 최소 갯수를 찾기

관련 문제

n을 최소 갯수의 제곱수의 합으로 표현하기

2개의 제곱수의 합으로 표현

가우스 정수를 이용한 알고리즘

Hermite-Serret 알고리즘

관련 문제