프로젝트 오일러 039

문제

39번 문제

가장 많은 직각삼각형이 만들어지는 둘레(≤ 1000)의 길이는?

사이냅소프트

풀이

이 문제는 피타고라스 트리플(피타고라스 삼조)를 이용하면 쉽게 풀 수 있는 문제이긴 합니다만, 일반적인 삼각형의 변의 길이에 대한 성질과 피타고라스 정리만 이용해도 범위를 크게 줄일 수 있기에 이 방법을 먼저 살펴보도록 하겠습니다.

• 삼각형의 가장 긴 변의 길이는 다른 두 변의 길이의 합보다 작습니다.
• 세 변의 길이가 같으면 직각 삼각형이 될 수 없습니다.
• 두 변의 길이가 같아도 직각 삼각형이 될 수 없습니다.

만약 세 변의 길이를 짧은 순으로 a, b, c 라 하고, 그 합을 p라 하겠습니다. 그러면 일정한 p에 대해서 a, b, c는 다음과 같은 관계를 가질 수 있습니다.

1. a는 가장 짧은 변이므로 p의 1/3 보다는 작습니다.
2. b는 a보다는 큰 자연수이며, ( p - a)의 1/2 보다는 작습니다.
3. a, b가 어떤 값을 갖는다면 c는 (p - a - b)로 결정되며, 이 세 수는 피타고라스 정리를 만족해야 합니다.

따라서 어떤 p에 대해서 세 변의 길이가 자연수인 직각 삼각형을 만들 수 있는 함수를 아래와 같이 작성해보겠습니다.

def f(p: int) -> tuple(int, int):
  res = 0
  for a in range(1, p // 3):
    for b in range(a + 1, (p - a) // 2):
      c = (p - a - b)
      if a * a + b * b == c * c:
        res += 1
  return (res, p)

그러면 6 - 1,000 범위에서 위 함수를 실행한 후 가장 큰 값을 고르면 됩니다.

개선

위 코드에서 a와 b가 결정되면 c 값은 p값에 따라 고정됩니다. 하지만 직각 삼각형의 빗변이 될 수 있는 c는 p에 무관하게 a, b만으로 결정됩니다. 그래서 아래와 같이 코드를 개선합니다. 결과적으로는 3중 루프인 것을 2중 루프로 변경했고, 함수 호출에 대한 오버헤드도 줄이게 돼서 20배 가량 성능이 올라가는 것을 확인할 수 있었습니다.

from math import hypot

def g() -> list[int]:
  limit = 1000
  res = [0] * limit
  for b in range(2, limit // 2):
    for a in range(1, b):
      c = hypot(a, b)
      d = a + b + int(c)
      if int(c) == c and d < limit:
        ps[d] += 1
  return res

  print(max(enumerate(g(1000)), key=lambda x: x[1]))

피타고라스 삼조

세 자연수 (a, b, c)에 대해 이 값들이 피타고라스 정리를 만족하면 이를 피타고라스 삼조라고 합니다. 이 때 a, b, c 가 모두 서로 소인 경우를 '원시 피타고라스 삼조'라고 합니다. 모든 피타고라스 삼조들은 원시 피타고라스 삼조의 배수입니다.

이 문제는 피타고라스 삼조를 쓰라는 의도가 보이는 문제라, 참고 차원에서 언급하고 넘어가겠습니다.

피타고라스 삼조는 서로 소인 두 정수 m, n (m > n)이 있을 때, (m² + n², 2mn, m² - n²)으로 구할 수 있습니다. 이 때 둘 중 하나가 짝수라면 원시 피타고라스 삼조가 됩니다. 원시 피타고라스 삼조를 만들었을 때, 세 요소 값의 합이 직각삼각형의 합이 되고, 한계값인 1,000보다 작은 그의 배수들 모두 직각 삼각형을 이룹니다. 따라서 배수와 그 합에서의 직각 삼각형의 개수를 기록하면 됩니다.

피타고라스 삼조의 합은 2m² + 2mn 으로, 이 값이 1000을 넘기지 않을 때 m의 최대값은 n을 고려하지 않는다면 한계값의 절반인 500의 제곱근보다 작거나 같습니다. 이 상의 내용으로 아래와 같은 코드를 작성하여 풀이할 수 있습니다. 물론 반복 횟수가 가장 작아서 이 코드가 가장 빠를 겁니다.

def gcd(a, b):
    if a < b:
        return gcd(b, a)
    r = a % b
    if r == 0:
        return b
    return gcd(b, r)


def main(limit=1000):
    ps = [0] * limit
    for m in range(2, int((limit / 2) ** 0.5 + 1)):
        for n in range(1, m):
            if gcd(m, n) == 1 and m * n % 2 == 0:
                p = 2 * m * (m + n)
                for i in range(p, limit, p):
                    ps[i] += 1
    return max(enumerate(ps), key=lambda x: x[1])


print(main())