Project Euler #12 - 500개 이상의 약수를 갖는 최소의 삼각수

"500개 이상의 약수를 갖는 최소의 삼각수" - 우리말로 써 놓기는 했지만 이걸 보고 문제를 이해하기는 어렵겠다.


요렇게 하면 0.6초 걸려서 답이 나온다.

from sympy import factorint
from operator import mul

i = 10
while True:
    f = factorint(i*(i+1)/2)
    if reduce(mul,[f[k]+1 for k in f]) > 500:
        print i*(i+1)/2
        break
    i += 1

i를 굳이 10에서 시작할 필요는 없다. 10000 정도에서 시작해도 무방.
i*(i+1)/2 로 삼각수를 만들었는데.. 아래가 더 낫다.

i, t = 10, 55
while True:
    f = factorint(t)
    if reduce(mul,[f[k]+1 for k in f]) > 500:
        print i*(i+1)/2
        break
    i += 1
    t += i

하지만 병목은 '삼각수'를 만드는 게 아니라 '소인수 분해'하는 쪽이어서 이렇게 한다고 뭐가 개선되거나 하지는 않는다.

reduce - mul에 대해서는 8번 문제에서 간단히 썼다. reduce를 쓰지 않고 더 짧은 코드로 문제를 풀 수도 있다. sympy에는 약수를 모두 찾아주는 함수(divisors)도 있다. 대신 이렇게 하면 실행 시간이 두배 정도 걸린다. 약수의 개수가 필요한데 약수를 모두 계산하려니 시간이 더 걸리나보다.

from sympy import divisors

i = 0
while True:
    if len(divisors(i*(i+1)/2)) > 500:
        print i*(i+1)/2
        break
    i += 1

더 나은 방법은 sympy를 쓰지 않고 5번에서 만들어 두었던 소인수 분해 함수를 응용하는 것이다. 소인수분해를 하려면 소수(prime numbers)의 리스트가 필요하다. 소인수분해를 할 때마다 이 리스트를 반복해서 만들면 효율이 떨어질 수 밖에 없다. 필요한만큼 넉넉하게 소수를 준비해 두고 소인수분해에서 활용하면 좋은데, sympy의 factorint를 써서는 이게 불가능한 듯 하다. 그럼 직접 만들어 써야지. 아래에서는 10만까지의 소수만 만들었으니, 만약 답이 100억이 넘는다면 소수를 더 만들어서 다시 돌려야 한다. 다행히 답은 1억 아래에서 나왔고 실행 시간은 0.3초 정도로 sympy.factorint를 썼을 때의 절반으로 줄었다.

from operator import mul

def rwh_primes1(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]

def ndiv(n,pn):
    # n > 2
    i, fDic = 0, {}
    while True:
        p = pn[i]
        if n % p == 0:
            if p in fDic: fDic[p] += 1
            else: fDic[p] = 1
            n /= p
            i -= 1
        i += 1
        if pn[i]**2 > n and n > 1:
            if n in fDic: fDic[n] += 1
            else: fDic[n] = 1
            return reduce(mul,[x+1 for x in list(fDic.values())])

pn = rwh_primes1(100000)
        
tn, lasttn = 1, 2
while True:
    tn += lasttn
    if ndiv(tn,pn) > 500:
        print tn
        break
    lasttn += 1