Project Euler #23 - Non-abundant sums

Non-abundant sums: 두 개의 abundant numbers(자기 자신을 뺀 약수의 합이 자기 자신 보다 커지는 숫자)를 더해서 만들 수 없는 숫자들의 합은?

문제에서 28123보다 큰 수는 모두 두 개의 a.n 합으로 표현할 수 있다고 했으니 그보다 작은 수만 체크하면 된다. 두 개의 a.n 로 만들 수 "없는" 것보다는 만들 수 "있는" 쪽을 찾는 게 쉬워 보인다.

sympy를 쓰면 어렵지 않다.

from sympy import divisors

abun = []
for i in range(12,28123):
    if sum(divisors(i)[:-1]) > i: abun.append(i)

allabuncomb = [True] * 28123
for i,a in enumerate(abun):
    for b in abun[i:]:
        if a+b >= 28123: break
        allabuncomb[a+b] = False

print sum([n for n in range(28123) if allabuncomb[n]])

역시나 다른 문제들처럼 소수-소인수분해-약수합 부분을 직접 구현하면 약간 빨라진다. 하지만 abundant number 를 찾는 부분보다, 그 숫자들로 조합을 만드는 부분이 오래 걸리는 거라 큰 개선은 없다.

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 factor(n,pn):
    rs = {}
    for p in pn:
        if n % p == 0:
            cnt = 0
            while n % p == 0:
                n /= p
                cnt += 1
            rs[p] = cnt
            if n == 1: return rs
        if p*p > n:
            rs[n] = 1
            return rs

def isAbun(n,pn):
    fDic = factor(n,pn)
    dn = 1
    for k in fDic:
        dn *= sum([k**i for i in range(fDic[k]+1)])
    return (dn - n > n)

pn = rwh_primes1(28123)
abun = []
for n in range(12,28123):
    if isAbun(n,pn): abun.append(n)

allabuncomb = [True] * 28123
for i,a in enumerate(abun):
    for b in abun[i:]:
        if a+b >= 28123: break
        allabuncomb[a+b] = False

print sum([n for n in range(28123) if allabuncomb[n]])

Project Euler #21 - amicable numbers

10000이하의 친화수를 찾는 문제

sympy의 divisors함수를 쓰면 쉽다.

from sympy import divisors

d_nDic = {}
for i in range(3,10000):
    d_nDic[i] = sum(divisors(i))-i

ami = []
for i in d_nDic:
    if d_nDic[i] > 3 and d_nDic[i] < 10000:
        if i == d_nDic[d_nDic[i]] and i != d_nDic[i]:
            ami.append(i)

print sum(ami)

10000까지 반복해서 소인수분해를 해야 하는데, 좀 더 빠르게 하려면 소수를 직접 찾고, 소인수 분해를 직접 한 다음에, 약수의 합을 계산하면 된다.
마지막 스텝은 중학교 때 배운 공식 24 = 2^3 * 3 => 24의 약수의 합=(1+2+4+8)*(1+3)=1+2+4+8+3+6+12+24 를 쓰면 된다.
코드는 길어지지만 실행시간은 1/4로 줄어든다.

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 factor(n,pn):
    rs = {}
    for p in pn:
        if n % p == 0:
            cnt = 0
            while n % p == 0:
                n /= p
                cnt += 1
            rs[p] = cnt
            if n == 1: return rs
        if p*p > n:
            rs[n] = 1
            return rs

def d_n(n,pn):
    fDic = factor(n,pn)
    dn = 1
    for k in fDic:
        dn *= sum([k**i for i in range(fDic[k]+1)])
    return dn - n

pn = rwh_primes1(10000)
d_nDic = {}
for i in range(3,10000):
    d_nDic[i] = d_n(i,pn)

ami = []
for i in d_nDic:
    if d_nDic[i] > 3 and d_nDic[i] < 10000:
        if i == d_nDic[d_nDic[i]] and i != d_nDic[i]:
            ami.append(i)

print sum(ami)

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

Project Euler #5 - 20 이하의 모든 자연수로 나누어지는 최소의 숫자

Smallest multiple - 20 이하의 모든 자연수로 나누어지는 최소의 숫자

크게 어렵지 않으니 그냥 계산하면 된다
2*2*2*2*3*3*5*7*11*13*17*19 = 232792560

괜히 어렵게 문제를 풀어보면..

20까지의 숫자를 소인수 분해하고, 그 결과 지수들의 최대값을 찾고, 최소공배수를 찾아보자.

def chk_pn(n,pn):
    for p in pn:
        if n % p == 0: return False
        if p**2 > n: return True

def inc_pn(pn):
    n = pn[-1] + 2
    while True:
        if chk_pn(n,pn): pn.append(n); return pn
        n += 2

def factor(n,pn):
    # n > 2
    i, fDic = 0, {}
    while True:
        p = pn[i]
        if p == pn[-1]: pn = inc_pn(pn)
        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:
            if n in fDic: fDic[n] += 1
            else: fDic[n] = 1
            return fDic

pn, SCM = [2,3], {}
for i in range(4,21):
    fDic = factor(i,pn)
    for k in fDic:
        if k in SCM: SCM[k] = max(SCM[k],fDic[k])
        else: SCM[k] = fDic[k]

op = 1
for x in SCM:
    op *= x**SCM[x]

print op

좀 부끄러운 지저분한 코드지만 답은 잘 나온다. 

마찬가지로 from sympy import factorint 를 이용하면 삽질을 줄일 수 있다.

Project Euler #3

600851475143를 소인수분해했을 때 가장 큰 소수는?

소인수분해하는 함수는 누군가 만들어 두었을테니 가져다 쓸 것인지, 직접 만들 것인지를 정해야 하네요. 일단, 처음이니 만든다고 하고, 다음은 소인수 분해를 위해 '소수'의 리스트를 직접 만들 것인지, 가져다 쓸 것인지, 소수의 리스트 없이 소인수 분해 할 것인지를 정해야겠어요. 이왕 하는 거, 소수를 먼저 찾고, 소인수 분해하는 것까지 해 보겠습니다.

def chk_pn(n,pn):  
    for p in pn:  
        if n % p == 0: return False  
        if p**2 > n: return True  
   
def inc_pn(pn):  
    n = pn[-1] + 2  
    while True:  
        if chk_pn(n,pn): pn.append(n); return pn  
        n += 2  
   
def factor(n,pn):  
    # n > 2  
    i, fDic = 0, {}  
    while True:  
        p = pn[i]  
        if p == pn[-1]: pn = inc_pn(pn)  
        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:  
            if n in fDic: fDic[n] += 1  
            else: fDic[n] = 1  
            return fDic  
   
print max(factor(600851475143,[2,3]))  

소수의 리스트가 어디까지 필요한지 몰라서 필요한 만큼 소수를 만들어 내도록 했는데, 별로 효율적인 방법은 아닙니다. 소수를 빨리 찾는 것은 "에라스토테네스의 체"를 이용해야죠. 아무튼, 1.5ms가 걸렸습니다.


sympy의 factorint를 쓰면 0.5ms가 걸립니다.

from sympy import factorint  
print max(factorint(600851475143))  

소수를 효율적으로 만들고 소인수 분해 결과를 {71: 1, 839: 1, 1471: 1, 6857: 1} 이런 식으로가 아니라 {71, 839, 1471, 6857} 요런 식으로 저장해 봤습니다(계산 시간은 차이가 거의 없지만, 코드가 간단해져서..)
factorint 가져다 쓰는 것보다 아주 약간 느리네요. 소수를 10000까지만 만드는 것은 찜찜하죠. 더 많이 만들어서 쓰려면 시간이 더 오래 걸리고.. --;

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 factor(n,pn):  
  rs = set()  
  for p in pn:  
    if n % p == 0:  
      while n % p == 0:  
        n /= p  
      rs.add(p)  
      if n == 1: return rs  
    if p*p > n:  
      rs.add(n)  
      return rs  
   
print max(factor(600851475143,rwh_primes1(10000)))  

소수인지 신경 쓰지 않고 그냥 다 나눠보면 이렇게도 할 수 있네요.

n, p = 600851475143, 2  
while p*p < n:  
  if n % p == 0:  
    n /= p  
  else:  
    p += 1  
   
print n  

7ms. 시간도 괜찮은데.. 소수니 소인수분해니 왜 이렇게 삽질을 했을까요..ㅋ
그래도 이렇게 공부하니 좋죠~

마지막으로, 짝수로 나눠볼 필요는 없으니..

n, p = 600851475143, 3  
while p*p < n:  
  if n % p == 0:  
    n /= p  
  else:  
    p += 2  
   
print n  

이렇게 하면 4ms!