소인수분해하기

소인수분해와 약수의 개수

어떤 자연수 n을 소인수 분해하였을 때 나오는 소인수의 모든 지수에 1을 더하고, 그 값들을 모두 곱하면 n의 약수의 개수가 된다.

이는 n의 모든 약수는 각각의 소수인 인수 a, b, c.... 들을 0번에서 각 지수번까지 사용한 조합을 만드는 것과 동일하기 때문이다.

따라서 어떤 자연수 n을 소인수 분해 할 수 있다면 그 약수의 개수를 빠르게 계산할 수 있을 것이다.

문제

오일러 프로젝트 12번 문제는 약수의 개수를 빠르게 구하는 것에 중점을 두는 문제이다.

1부터 n까지의 자연수를 차례로 더하여 구해진 값을 삼각수라 합니다. 예를 들어 7번째 삼각수는 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 입니다. 이런 식으로 삼각수를 구해 나가면 다음과 같습니다.

1, 3, 6, 10, 15, 21, 28 ...

이 삼각수들의 약수를 구했을 때 5개 이상의 약수를 갖는 첫번째 삼각수는 28(6개)입니다. 그러면 500개 이상의 약수를 갖는 가장 작은 삼각수는 무엇입니까?

지난 번에 약수의 개수, 약수의 합을 구하는 루프를 도는 로직으로 이 문제는 풀 수 있다. 관건은 삼각수 수열의 일반항을 쉽게 구할 수 있는가인데, 삼각수의 일반항은 n 까지 자연수의 합이므로 파스칼의 공식을 통해서 다음과 같이 구할 수 있다.

$$ a_{n} = \frac {n \times (n + 1)} {2}$$

이렇게 구한 삼각수의 일반항의 값이 약수의 개수가 500 이상이 되는 최초의 케이스를 찾아서 보고하면 된다. 약수의 개수를 셀 때에는 곱셈식과 약수의 관계에 주목해볼 필요가 있다. 어떤 두 자연수 P, Q 의 곱이 N 이라하면 이 수들의 관계는 다음과 같이 나타낼 수 있다.

$$ p \times q = N $$

이 때, p, q는 그 자체가 소수인지와는 무관하게 N의 약수가 된다. 만약 p <= q 인 관계가 성립한다면, N의 제곱근 r 은 p와 q 사이에 오게 된다.

$$ p \le r \le q $$

따라서 N이 p 로 나눠진다면 우리는 약수 2개를 찾은 것이 된다. 다만 p가 N의 제곱근이라면 이 경우만 1개를 찾았다고 하면 된다. (즉 완전제곱수의 약수의 개수는 홀수이며, 그 외에는 항상 짝수 개의 약수를 가짐을 알 수 있다.)

정리하면 소수를 포함한 모든 정수는 기본적으로 2개의 약수(1과 그 자기자신)를 가진다. 그리고 2부터 N의 제곱근까지의 정수들을 후보로 나눠지는지 검사하고, 나눠진다면 약수의 개수에 2씩 더해준다. 예외적으로 완전제곱수의 제곱근은 1만 더해야 하며, 이 처리를 빼먹지 않아야 한다.

def number_of_divisors(n):
    # 경계값
    l = int(n**0.5 + 0.5)
    s = 1 if l * l == n else 2
    a = 2
    while a <= l:
        if n % a == 0:
            s += 2
        a += 1
    return s

def main():
    k = 1
    while True:
        # a는 k번째 삼각수
        a = k * (k + 1) // 2
        if number_of_divisors(n) >= 500:
            # 약수의 개수가 500개 이상이면 출력하고 중지
            print(a)
            break

# %time 명령은 ipython이나 jupyter notebook에서 사용가능
%time main()

이렇게 풀어도 답은 나온다만, 시간이 생각보다 매우 오래 걸린다. 이 문제에서는 500개 이상의 약수를 가질 수 있는 삼각수 자체가 워낙 큰 값이 될 것이기 때문에 약수의 개수를 계속 구하면서 조건을 만족하는 값을 찾기 까지는 시간이 제법 걸리게 될 것이다. (여러분의 좋은 컴퓨터에서는 모르겠지만, 여기선 8~9초 가량이 걸렸다.)

'너무 크지 않은 수'에 대해서는 약수의 개수를 구할 때 위 방법이 매우 효율적이긴하나, 대상 숫자가 너무 큰 값이라면 시간이 오래 걸린다. 여기서 정답은 7천 6백만이 조금 넘는 수이기 때문에 루프를 어마어마하게 많이 돌아야 한다. 그리고 이미 약수의 개수 구하는 함수는 n의 제곱근까지만 검사하기 때문에 어느 정도 최적화가 되어 있는 상태이다. 그래서 이렇게 큰 수의 약수의 개수를 구하려면 소인수 분해가 조금 더 나은 성능을 제공할 수 있다.

소인수분해

어떤 자연수 N을 소인수 분해하는 과정에 대해 생각해보자.

1. 가장 작은 소수인 2로 N을 나눠본다.
2. 한번 나눌 때마다 N은 2로 나눈 몫이되고, N을 나눌 수 있을 때까지 여러 번 나눈다. 이 때 나눈 횟수는 2의 지수가 된다.
3. 2보다 큰 가장 작은 소수인 3으로 2의 과정을 반복한다.
4. 선택하는 소수를 계속 늘려나가면서 2의 과정을 반복해 나간다.
5. 최종적으로 N이 작아지다 1이되면 소인수 분해가 끝난다.

보통 4의 과정에서 많은 선택지 중 하나를 택해야 한다. 소인수의 후보를 어떻게 찾을 것인가 하는 것이다. 반복의 횟수를 줄이는 것이 성능에서 가장 중요하게 고려해야 할 문제이다.

• 2부터 3, 4, 5 ... 로 계속 나눠본다.
• 2, 3, 5, 7, 11 .. 로 소수들로만 나눠본다.
• 2, 3, 7, 9, 11 ... 2와 홀수들로만 나눠본다.

어느쪽이 더 좋을지는 사실 본인의 선택이겠다. 소수들로만 나눠보는 것이 시행의 횟수를 가장 줄일 수 있겠지만, N 보다 작은 소수를 구하거나, 혹은 나누려는 수가 소수인지를 매번 체크하는 것이 추가적인 부담이 될 수 있다. 여기서는 세 번째 방법을 조금 개선한 것을 사용할 것이다.

튜플

예를 들어 240을 소인수 분해하면 다음과 같다.

$$240 = 2^{4} \times 3^{1} \times 5^{1} $$

이 때 각 소수와 그 지수를 어떻게 표현할까? 여기서는 밑과 지수를 하나의 튜플로 표현해서 튜플의 리스트로 표현해 볼 것이다.

[(2, 4), (3, 1), (5, 1)]

튜플은 순서가 정해진 요소의 묶음이다. 리스트와 달리 내부를 변경하지 않기 때문에, 리스트보다는 조금 가볍게 사용할 수 있다.

나눌 수 있는 값 찾기

문제는 N을 나눌 수 있는, 그리고 방금 처리를 거친 값보다 큰 소수를 어떻게 찾느냐 하는 점이다. 이를 찾는 헬퍼 함수를 하나 작성할 것이다.

def find_divisible(n, m):
    '''n을 나눌 수 있는 m 보다 큰 값
       m은 2이상이어야 한다.'''
    if m == 2:
        return 3
    if m == 3:
        return 5
    # 5부터 시작해서, 5, 7, (9), 11, 13, (15), ... 
    # 순으로 검사한다. 만약 n의 제곱근 이하에서 나눌 수 있는 수가 없다면
    # n은 소수이므로 n을 리턴하면 된다.
    k = 5
    while k * k <= n:
        if n % k == 0:
            return k
        if n % (k + 2) == 0:
            return k + 2
        k += 6
    return n

파이썬에서는 함수내에서 함수를 정의할 수 있는데, 이를 중첩된(nested) 함수라 한다. 아래와 같은 식으로 작성하는 경우 helper() 함수는 그 내부에서 n 값에 대한 참조를 유지할 수 있다.

def factorize(n):
    def helper(m):
        if m == 2:
            return 3
        if m == 3:
            return 5
        k = 5
        ...

헬퍼 함수를 사용해서

def factorize(n):
    # 헬퍼함수 정의
    def helper(m):
        if m == 2:
            return 3
        if m == 3:
            return 5
        k = 5
        while k * k <= n:
            if n % k is 0:
                return k
            if n % (k + 2) is 0:
                return k + 2
            k += 6
        return n

    # result는 결과를 담는 리스트, a는 나눠볼 소인수 후보, b는 소인수의 지수
    result, a, b = [], 2, 0
    while n > 1:
        # a로 나눌 수 있는 만큼 반복하면서 지수를 늘려나간다.
        while n % a is 0:
            n, b = n // a, b + 1
        # 한 번 이상 나눴다면 (a, b) 쌍을 결과에 추가한다.
        if b > 0:
            result.append((a, b))
            b = 0
        # 다음 후보 값을 찾는다.  
        a = helper(a)
    return result

이렇게 해서 소인수 분해 함수를 완성했다. 그럼 소인수 분해의 결과로 어떻게 약수의 개수를 구할 수 있을까? 어떤 자연수 N의 소인수 분해 결과가 a^3 * b^2 이라고 하자. 그러면 그 약수들은 항상 다음의 모양이 될 수 있다.

$$ a^{i} \times b^{j} (i=0, 1, 2, 3 j=0, 1, 2) $$

따라서 a를 03개 고르고 b를 02개 고르는 경우의 수가 되므로 4 * 3 = 12 개의 약수가 나오게 된다. 즉 소인수 분해한 결과에서 모든 지수에 1씩을 더해서 곱하면 되는 것이다.

소인수 분해 함수의 결과는 [(소수, 지수)]의 리스트이므로, 여기에서 지수 부분만 1을 더해서 모두 곱해주는 계산은 개인적인 취향으로 for 루프를 돌지 않고 reduce() 함수를 사용할 것이다.

def number_of_divisors(n):
    fs = factorize(n)
    return reduce(lambda x, y: x * y, (f[1] + 1 for f in fs))

reduce() 함수를 쓰려면 functools에서 해당 함수를 반입해주는 절차가 필요하다. 따라서 최종 결과는 다음과 같다.

def factorization(n):
    def helper(m):
        if m == 2:
            return 3
        if m == 3:
            return 5
        k = 5
        while k * k <= n:
            if n % k == 0:
                return k
            if n % (k + 2) == 0:
                return k + 2
            k += 6
        return n

    res, a, b = [], 2, 0
    while n > 1:
        while n % a == 0:
            n, b = n // a, b + 1
        if b > 0:
            res.append((a, b))
            b = 0
        a = helper(a)
    return res


def number_of_divisors(n):
    fs = factorization(n)
    return reduce(lambda x, y: x * y, (f[1] + 1 for f in fs), 1)



def main():
    i = 1
    while True:
        k = i * (i + 1) // 2
        if number_of_divisors(k) >= 500:
            print(k)
            break
        i += 1


if __name__ == "__main__":
    main()

# elapsed: 0.219s

코드가 다소 복잡해보이기는 하지만 결과는 매우 만족스럽다.

보너스 - 모든 약수의 합을 구하기

이전에 약수의 합을 구하는 문제에서 소인수분해로 모든 약수의 합을 구하는 방법이 있다고 하였다. 자연수 N을 소인수 분해한 결과가 [(a1, e1), (a2, e2), (a3, e3),...]일 때 N의 모든 약수의 합은 다음과 같이 계산한다.

1. 각 소수에 대해서 0부터 해당 지수까지의 지수를 갖는 수열의 합을 각각 구하고
2. 1의 각 값들을 곱한다.

즉 $$(1 + a_{1} + a_{1}^{2} + ... + a_{1}^{e_{1}}) \times (1 + a_{2} + a_{2}^{2} + .... + a_{2}^{e_{2}}) \times (a_{n}^{0} + a_{n}^{1} + ... + a_{n}^{e_{n}})$$이 된다.

따라서 약수의 합을 소인수분해를 통해서 구하려면 다음과 같이 코드를 작성하면 된다.

def sum_of_divisors(n):
    xs = factorize(n)
    ys = [sum((a ** i for i in range(e+1))) for (a, e) in xs]
    return reduce(lambda x, y: x * y, ys, 1)