콘텐츠로 건너뛰기
Home » 오일러 프로젝트 53

오일러 프로젝트 53

1, 2, 3, 4, 5 다섯 숫자 중에서 세 개를 고르는 것에는 다음과 같은 10가지 경우가 있습니다.

>> 123, 124, 125, 134, 135, 145, 234, 235, 345

조합론이라는 분야에서는 이것을 5C3 = 10 이라고 표현하며, 일반적인 식은 아래와 같습니다.

\binom{n}{r} =\frac{n!}{r!(n-r)!}

이 값은 n=23에 이르러 23C10 = 1144066으로 처음으로 백만을 넘기게 됩니다. 그렇다면 1 ≦ n ≦ 100 일 때 nCr의 값이 1백만을 넘는 경우는 모두 몇 번입니까?

http://euler.synap.co.kr/prob_detail.php?id=53

접근

조합의 공식은 간단하기 때문에 팩토리얼 함수만 있다면 간단하게 만들 수 있다. n = 100이면 팩토리얼을 계산하는 중간에 매우 큰 값이 등장하지만 파이썬에서는 그냥 계산하기만 하면 된다.

from functools import reduce

def factorial(n: int) -> int:
  if n < 2:
    return 1
  return reduce(lambda x, y: x * y, range(1, n+1))

def ncr(n, r):
  return factorial(n) // factorial(r) // factorial(n - r)

%%time
print(sum(1 for n in range(1, 101) for r in range(1, n+1) if ncr(n, r) > 100_0000))

늘 말하지만 큰 수를 사용하여 문제를 푸는 것은 너무 반칙같아서, 큰 정수의 마법을 사용하지 않고 피해보자. 이전의 덧셈이나 곱셈은 문자열을 사용해서 손으로 계산하는 방법을 그대로 적용했지만, 나눗셈은 구현을 고민하기가 싫어서, 계산과정에서 값이 너무 커지지 않도록 하는 nCr계산 함수를 작성해서 미안함을 보상해보도록 하겠다. nCr의 계산 결과는 언제나 자연수이기 때문에, nCr을 구하는 이 분수식의 분모의 모든 인수는 약분되어 1만 남게 된다. 분모의 인수들을 차례로 1이 될때까지 분자의 인수들과 약분하는 것을 반복한 후 분자에 남은 수들만 모두 곱하면 된다.

def gcd(a: int, b: int) -> int:
  if a < b:
    return gcd(b, a)
  r = a % b
  return b if r == 0 else gcd(b, r)

def nCr(n: int, r:int) -> int:
  ns = [*range(1, n + 1)]
  rs = [*range(1, r + 1), *range(1, n - r + 1)]
  for k in range(len(rs)):
    for j in range(len(ns)):
      if rs[k] == 1:
        break
      g = gcd(ns[j], rs[k])
      if g > 1:
        ns[j], rs[k] = ns[j] // g, rn[k] // g
  return reduce(lambda x, y: x * y, ms)

이 코드를 사용하면 nCr를 계산하는 동안 큰 수가 발생하지 않는다. 같은 알고리듬을 사용하면 큰 수를 지원하지 않는 프로그래밍 언어에서도 빠르게 계산할 수 있을 것이다.

파이썬의 큰 수 계산 능력은 깡패수준이기 때문에 이 계산은 순식간에 처리되고 답이 나온다. 큰 수 연산을 사용하지 않는다면 어떻게 계산할 수 있을까? 문제에서 제시한 23C10을 계산하기 위해 23!을 계산하면 ‘25,852,016,738,884,976,640,000’이 나오는데, 이 함수를 사용하면 17 x 19 x 7 x 22 x 23 만 계산하면 된다.

sum(nCr(n, r) > 100_0000 for n in range(1, 101) for r in range(1, n+1))

다만 nCr을 계산할 때 루프를 제법 돌기 때문에 수행 시간은 다소 오래 걸린다.

이항정리 – 큰 수 없이 동적 계획법으로 풀기

조합의 가지수를 구할 때 위 공식에서 r ≦ n 이므로, 어떤 n에 대해서 가능한 모든 nCr의 값을 계산해 볼 수 있고, 이 때의 각각의 값은 (x + y)n을 전개한 식의 각각의 항의 계수와 일치한다. nCr은 n개 중에서 r 개를 고르는 가짓수를 계산하는 방법이다. n 개 중에서 r개를 고르는 상황을 좀 자세히 들여다보자.

n 개 중에서 한 개 요소를 고를 때 우리가 선택할 수 있는 것은 그 한 개를 선택하는 경우와 선택하지 않는 경우가 있다. 만약 한 개를 선택했다면, 남은 n-1 개 중에서 r-1개를 선택하면 되고, 선택하지 않았다면 남은 n-1 개 중에서 r 개를 선택하면 된다. 결국 하나의 원소를 선택한 경우와 선택하지 않은 경우의 남은 원소로 조합을 만드는 경우의 수를 결합하여 nCr을 구할 수 있는 셈이다. 이 규칙을 식으로 나타낸 것이 아래와 같고, 이를 ‘이항정리’라고 부른다.

_nC_r = \space_{n-1}C_{r-1} + _{n-1}C_r 

그 외에 r = 1 인 경우에는 nCr = n, r = n 인 경우에는 nCr = 1 이다. 이 기준값들과 이항정리 공식을 사용하면 nCr 은 간단한 동적 계획법으로 풀이할 수 있다. 대각선을 효과적으로 다루기 위해서 그리드를 2차원 리스트 대신 1개의 1차원 리스트를 사용해서 표현했다. 또한 큰 수를 계산하지 않기 위해서 100만을 넘어가는 시점에 그저 1백만 + 1 값이 되도록 상한을 정해버린다.

%%time

s = [0] * 10000
s[::100] = [i + 1 for i in range(100)]
s[::101] = [1] * 100
for n in range(1, 100):
    for r in range(1, n):
        s[n*100+r] = min(100_0001, s[n*100+r-100] + s[n*100+r-101])
print(sum(x > 100_0000 for x in s ))

1차원 배열을 2차원 그리드처럼 사용하는 기술은 numpy 배열에도 동일하게 적용된다. 역시 numpy를 사용하여 위 로직과 비슷하게 구현할 수 있다. 기본적으로 워낙 간단하고, 루프 수는 생각보다 많지 않기 때문에 단순한 코드를 쓰는 위쪽 코드가 더 나은 성능을 보이는 경우가 종종 있다는 점이 놀랍다.

%%time

s = np.zeros((100, 100), dtype='int')
s[:, 0] = np.arange(1, 101)
s[np.diag_indices(100)] = np.ones(100)
for n in range(1, 100):
    for r in range(1, n):
        s[n, r] = min(100_0001, s[n - 1, r - 1] + s[n - 1, r])
(s > 100_0000).sum()