Tuesday, April 14, 2015

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!

No comments: