사용자 도구

사이트 도구


ps:제곱수의_합

제곱수의 합

  • 자연수 n을 최소갯수의 제곱수의 합으로 표현하는 것에 관련된 이론들이다.
  • 크게 두가지 주제가 있다.
    • n을 제곱수의 합으로 표현할 때 필요한 제곱수의 최소 갯수를 찾기
    • n을 최소 갯수의 제곱수의 합으로 표현할 수 있는 제곱수를 찾기
  • 'n을 k개의 제곱수의 합으로 표현하는 방법의 갯수를 찾기' 에 관해서도 이론들이 있는데.. 아직 이것을 필요로 하는 문제는 본적이 없으니 패스
  • 참고자료

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

  • 우선 간단히 알아둘 것은, 모든 자연수는 최대 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 알고리즘을 사용하는 방법이 있다.

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

  • 알고리즘 자체는 간단하게 요약된다. 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 알고리즘 이라는 이름 없이 그냥 알고리즘 설명만 되어있는 자료도 많다. 방법은 다음을 참고
  • Cornacchia's 알고리즘은 x^2+dy^2=p 을 만족하는 x와 y를 찾는 알고리즘이다. Hermite-Serret의 일반화된 버전이라고 생각해도 될거 같은데.. Hermite-Serret보다는 이쪽이 더 유명한것 같다.. (따로 위키페이지도 있고)
  • 이 방법을 사용하면, 복잡한 가우스정수 클래스 없이도 코드를 좀더 간단하게 짤수 있다.
    • 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))

관련 문제

  • 23287 (n<1,000,000,000,000,000,000) 이지만, 소인수의 크기가 100,000 이하라는 제한이 있어서, trial division을 이용한 소인수분해와, brute force를 이용한 두 제곱수의 합 표현 찾기 방법을 써서도 풀수 있다.
  • 제곱수의 합 2 (More Huge) (n<1,000,000,000,000,000,000)

토론

댓글을 입력하세요:
I O I G X
 
ps/제곱수의_합.txt · 마지막으로 수정됨: 2023/06/14 04:53 저자 teferi