오일러 프로젝트 72

오일러 프로젝트 72 번 문제는 여태껏 나왔던 문제에서의 최고 난이도를 또 한 번 갱신했다. 오일러 피 함수(\phi )의 1에서 100만까지의 자연수에 대한 피함수 값의 합을 구해야하는 문제이며, 피 함수를 빠르게 작성하는 것이 얼마나 고된(?)일인지 알고 있다면 이 문제를 brute force로 푸는 것은 정말 답이 없다는 점에서 마음을 단단히 먹어야 한다.

n과 d가 양의 정수이고 n < d 인 분수 n/d을 GCD(n, d) = 1 일 때 기약 진분수라고 부르기로 합니다. d <= 8인 기약 진분수들을 커지는 순으로 늘어놓으면 아래와 같고 개수는 모두 21개입니다.

1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8

d가 백만 이하의 경우에 기약 진분수는 모두 몇 개나 있습니까?

접근

분모가 a인 기약진분수의 개수는 1…a 사이에서 a와 서로 소인 자연수의 개수와 같다. 이는 곧 오일러의 피 함수(phi)와 같은 의미이다. 따라서 이 문제는 1 \le n \le 1,000,000 범위일 때 \phi(n) 의 합계를 구하라는 문제와 같다. 이를 위해서 제법 준수한 수준의 피함수 구현이 필요하다.

  1. 소인수의 지수는 필요없고 소인수 목록만을 리턴하는 소인수분해함수 작성
  2. 1.로부터 생성된 소인수를 이용해서 각 소인수의 배수를 제거해나가는 체를 작성
  3. 체에서 남은 숫자를 카운트

이러한 구현을 파이썬으로 해보면 다음과 같다.

def factorize(n):
  def divisor(a):
    if a in (1, 2):
      if n % (a + 1) is 0: return a + 1
      return divisor(a+1)
    if a is 3:
      if n % 5 is 0: return 5
      return divisor(5)
    k, l = a+2, n ** 0.5
    ds, i = [4, 2] if k % 3 == 2 else [2, 4], 0
    while k <= l:
      if n % k is 0: return k
      k, i = k + ds[i], (i + 1) % 2
    return n
    result, a = [], 1
    while n > 1:
      a = divisor(a)
      result.append(a)
      while n % a is 0:
        n //= a
    return result

def phi(n):
  fs = factorize(n)
  s = [False] + [True] * n
  for i in fs:
    if s[i]:
      s[i::i] = [False] * (n // i)
  return sum(1 for x in s if s)

은근 복잡해 보이지만, 피 함수 자체가 복잡한 함수라는 점을 감안하면 제법 준수한 성능이 나온다. 그런데 문제는 우리가 찾아야 하는 값의 범위가 1백만이라는 점이다. 이 함수로는 10000까지의 피함수의 합을 계산하는데만도 약 5초 정도가 소요된다. (테스트 환경의 CPU – E6700, 3.2Ghz) 아마 백만까지 계산하려 한다면 몇 시간.. 가지고는 어림도 없을 수준의 시간이 걸릴 것 같다.

그래서 좀 다른 접근법이 필요하다. 예를 들면 동적 계획법이라든지, 동적 계획법이나 동적 계획법 따위의 것들 말이다.

동적 계획법을 적용하여 계산 시간을 단축하기

이 알고리듬의 기본적인 아이디어는 자연수 N이 분모일 때, 1이하의 (분자가 자연수인) 분수는 N개 만들 수 있다는 점이다.

배열 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] 이 있다고 가정해보자. 이 중에 한 원소 3을 보자. 3은 3인데 (?) 이렇게 해석할 수 있다. “3이하의 자연수의 개수는 3개이다.”

임의의 자연수 N과 그 소인수인 A를 생각해보자. 1부터 N까지의 자연수 중에서 A의 배수들은 N을 분모로 하는 분수의 분자가 될 때, A에 의해서 분모와 약분된다. 즉 A의 배수들은 모두 우리가 구하고자하는 대상에 포함되지 않는다.

이러한 약분되어 없어질 분자들은 N이하에서 (N / A)개 만큼 있다. 이 때, N의 범위가 정해져 있으므로,  1…N 범위의 자연수 배열을 만들고 각 소인수의 배수의 개수를 제거해 나간다. 그러면 최종적으로 남는 것들은 모두 N과 서로 소 인 수만 남게 될 것이다.

시간을 빠르게 하기 위해서 소인수 A를 어떻게 찾으면 좋을까? 이건 소수체를 만들 때와 비슷한데, 인덱스 번호 2부터 출발해서 한 번도 나눠진 적이 없는 칸을 처음만나면,  그 칸의 인덱스가 소수라는 의미가 될 것이다.

  1. 분모가 N인 경우, 1이하의 분수는 분자 1~N까지 총 N개가 있을 수 있다.
  2. 만약 N의 소인수 a (단 a > 1)이 있다면, 1~N 사이에는 a의 배수가 N / a 개 만큼 있을 것이다. 그런데, 이 수들은 모두 a로 나눠질 것이므로, 분모 N과 함께 약분될 것이다. 이는 다른 소인수 b, c, d…에 의해서도 마찬가지 이다.
  3. 따라서 약분가능한 분자들 (각 소인수이 배수)의 개수를 N에서 점점 빼 나가면 최종적으로 기약진분수를 만드는 분자만 남게 된다.
  4. 이 과정은 동적 계획법으로 풀어낼 수 있다.

풀이는 다음과 같다.

def number_of_P(N):
  phi = list(range(N+1)) ## 0~N 까지 분모가 n 인 분수의 개수를 담은 리스트
  for a in range(2, N+1):
    if phi[a] == a : ## 지금까지 손댄적이 없는, 즉 n이 소수이면
      for k in range(a, N+1, a):
        phi[k] -= phi[k] // a ## 모든 a의 배수에 대해서 a으로 약분되는 분자 개수를 뺀다.
   return sum(phi) - 1 ## 분모 1을 제외한 나머지 분수들의 개수를 더해서 리턴.

%time print(number_of_P(100_0000))
## wall time: 1.73s