Saturday, May 30, 2015

Project Euler #50 - 연속된 소수의 합

프로젝트 오일러 50번 - consecutive prime sum
41 = 2+3+5+7+11+13 처럼 연속된 소수의 합으로 표현되는 소수가 있다.
100만 이하 소수 중, 가장 긴 소수 합으로 표현되는 것은?

두가지 방법
첫번째.
A. 소수 n에 대해 n이 최대 몇 개의 연속된 소수의 합으로 표현되는지를 찾는 함수를 정의한다
 A-1. 소수들을 더해가면서 합이 n보다 크면 앞쪽 소수를 하나씩 빼고, 합이 n보다 작으면 뒤쪽에 하나씩 더한다.
 A-2. 소수의 개수가 너무 적어지면 그만 둔다.
B. A에서 정의한 함수를 모든 n에 대해 적용한다.

두번째.
A. 2부터 시작해서 합이 100만 넘기 직전까지 소수들을 더한다. '합'들을 저장해 둔다.
 A-1. 만약, 100만 넘기 직전의 합이 소수이면 그것이 정답
 A-2. 그게 아니라면, 2부터 어디까지 더했을 때 가장 긴 sequence가 되면서 합이 소수가 되는지 찾는다. 그 길이를 maxlen에 저장해 둔다.
B. 3부터 시작해서 합이 100만 넘기 직전까지 소수들을 더한다. '합'들을 저장해 둔다.
 B-1. 만약, 100만 넘기 직전의 합이 소수이고 그것이 maxlen보다 길면 정답
 B-2. 그게 아니라면, 3부터 어디까지 더했을 때 가장 긴 sequence가 되면서 합이 소수가 되는지 찾는다. 그 길이를 maxlen에 저장해 둔다.
C. 5부터 시작해서 합이 100만 넘기 직전까지 소수들을 더한다. '합'들을 저장해 둔다.
 C-1. 만약, 100만 넘기 직전의 합이 소수이고 그것이 maxlen보다 길면 정답
 C-2. 그게 아니라면, 5부터 어디까지 더했을 때 가장 긴 sequence가 되면서 합이 소수가 되는지 찾는다. 그 길이를 maxlen에 저장해 둔다.
D. E. F. ... - 7, 11, 13부터 시작해서...
끝까지 갔으면 maxlen이 정답

첫번째 방법을 해 봤더니 예상보다 긴 sequence에서 예상보다 큰 소수가 나왔다. 그래서 두번째 방법으로 다시 코딩. 두번째 방법은 우연히 D-1에서 끝나서 빨리 답을 찾은 거라.. 정석은 아닌 듯 하다.

첫번째 방법의 코드(6초 남짓)
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 maxcon(n,pn,mx):
    st = ed = 0
    sm = 0
    while True:
        if sm < n:
            sm += pn[ed]
            ed += 1
        elif sm > n:
            sm -= pn[st]
            st += 1
        elif sm == n:
            return ed - st
        if ed > mx and ed - st < mx:
            return 0

pn = rwh_primes1(1000000)
mx, mxp = 1, 1

for p in pn:
    if maxcon(p,pn,mx) > mx:
        mx = maxcon(p,pn,mx)
        mxp = p

print mxp


두번째 방법의 코드(0.05초)
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]]

N = 1000000
pnl = rwh_primes1(N)
pns = set(pnl)

l = len(pnl)

def prob():
    mxlen, mx = 1, 1
    for i in range(l):
        sm, sml = 0, []
        for j in range(i,l):
            sm += pnl[j]
            sml.append(sm)
            if sm > N:
                if sm - pnl[j] in pns and len(sml) >= mxlen:
                    return sm - pnl[j]
                break
        for k,q in enumerate(sml[::-1]):
            if q in pns:
                if len(sml) - k > mxlen:
                    mxlen, mx = len(sml)-k, q
                break
    return mx

print prob()





저장하다 날아가서 다시 쓰려니 의욕이 안 생긴다 --;

Friday, May 29, 2015

Project Euler #49 - 자리바꿈 소수

프로젝트 오일러 49번 - prime permutations
4자리 소수 중 3330, 6660을 더해도 소수이고 이 세 소수가 각 자리수의 자리바꿈인 경우가 딱 두 가지 있다. 하나는 1487, 다른 것은?

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]]

pn = rwh_primes1(10000)

for p in [n for n in pn if n > 1000 and n < 3330]:
    if p + 3330 in pn and p + 6660 in pn:
        p2, p3 = p + 3330, p + 6660
        if set(`p`) == set(`p2`) == set(`p3`) and p != 1487:
            print `p`+`p2`+`p3`
            break


Thursday, May 28, 2015

Project Euler #48 - 제곱수의 합

프로젝트 오일러 48번 - self powers
x가 1~1000일 때, x^x 를 모두 더하면?

1000^1000 이면 3천자리 숫자가 된다. But, python 으로 숫자 다룰 때 overflow를 걱정하는 건 기우..

print str(sum([x**x for x in range(1,1000)]))[-10:]

이런게 뚝딱 답이 나오는 게 아직도 신기하다. 솔직히 이게 최선인가 싶기도 하고.. 그래서 좀 더 조심하는 코드를 짜 봤다.

def dkp(n):
    m = 1
    for i in range(n):
        m = m*n % 10000000000
    return m

sm = 0
for i in range(1000):
    sm += dkp(i+1)

print `sm`[-10:]

똑같은 답이 나오고 시간은 6배...
믿고 쓰자 파이썬!

Wednesday, May 27, 2015

Project Euler #47 - 네 개의 소인수를 갖는 연속된 네 수

프로젝트 오일러 47번 - distinct primes factors
연속되는 네 수이면서 각각의 소인수가 4개인 것 중 첫번째 숫자는?

sympy의 factorint를 이용하면 간단히 코딩할 수 있다.

from sympy import factorint

chk4 = 0
n = 10
while chk4 < 4:
    if len(factorint(n)) == 4: chk4 += 1
    else: chk4 = 0
    n += 1

print n -4

약 2초 동안 13만번의 소인수분해를 했다.
12번에서 썼던 소인수분해 함수를 다시 가져와 보면, 1.2초 정도로 시간이 줄어든다.

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 fDic

pn = rwh_primes1(1000)

chk4 = 0
n = 10
while chk4 < 4:
    if len(ndiv(n,pn)) == 4: chk4 += 1
    else: chk4 = 0
    n += 1

print n -4


Tuesday, May 26, 2015

Project Euler #46 - 골드바흐의 다른 추론

프로젝트 오일러 46번
Goldbach's other conjecture

유명한 "골드바흐의 추론" 말고 다른 추론 - 소수가 아닌 모든 홀수는 (소수+제곱수)로 표현된다-은 틀렸다. 소수와 제곱수의 합으로 표현되지 않는 소수가 아닌 홀수 중 가장 작은 것은?

일단, 소수를 모두 찾아두고, 임의의 (소수가 아닌) 홀수에 대해 그보다 작은 제곱수들에 대한 loop을 돌면서 차이가 소수가 되는지 체크한다.

그냥.. 하면 된다 -.-;

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 Gold(n,pn,sq):
    for s in sq:
        if n - 2*s in pn: return True
        if n - 2*s < 1: return False
    return False

def findNGold(MX):
    pn = set(rwh_primes1(MX))
    sq = [x*x for x in range(1,int((MX/2)**.5)+1)]
    for n in [n for n in range(9,MX,2) if n not in pn]:
        if Gold(n,pn,sq) == False:
            return n

print findNGold(10000)

이 추론이 어디서 깨질 지 몰라서 10만까지 체크해 봤다. 그러면 0.02초
10000까지만 체크하면 0.005초

Monday, May 25, 2015

Project Euler #45 - Triangular, pentagonal, and hexagonal

프로젝트 오일러 45번
삼각, 오각, 육각수가 모두 되는 수는?
40755는 세 성질을 모두 만족한다. 세 성질을 모두 만족하는 40755 다음의 수는?

예전에 작성했던 코드 - 약 0.1초

i,j,k=285,165,143
a, b, c = 40755, 40755, 40755
while True:
    mn = min(a,b,c)
    if a == mn: a += i+1; i += 1
    elif b == mn: b += 3*j+1; j += 1
    elif c == mn: c += 4*k+1; k += 1
    if a == b == c and a > 40755:
        print i,j,k,a
        break


그런데, "모든 육각수는 삼각수"라는 성질(포럼에서 봤어요~)을 이용하면 a, i와 관련된 코드를 모두 없앨 수 있다. 위에서 가장 많이 업데이트 되는 게 a라는 점 때문에.. 이렇게 하면 반 이하로 실행시간이 줄어든다.

j,k=165,143
b, c = 40755, 40755
while True:
    mn = min(b,c)
    if b == mn: b += 3*j+1; j += 1
    elif c == mn: c += 4*k+1; k += 1
    if b == c and b > 40755:
        print j,k,b
        break


Sunday, May 24, 2015

Project Euler #44 - Pentagonal numbers

프로젝트 오일러 44번
pentagonal 수 두 개를 더해도 pentagonal, 차이도 pentagonal이 되는 것 중, 최소의 "차이"는?

p_k - p_j = p_l
p_k = p_l + p_j

p_k + p_j = p_m
라고 해 보자.

l < j 조건을 만족하는 (j,l) 쌍에 대해서( (2,1),(3,1),(3,2),...) 적당한 p_k, p_m 이 존재하는지 찾아본다.
p_k - p_j = p_l 이면 p_k - p_l = p_j 도 따라서 성립한다. 혹시 p_k + p_l = p_m이 되는 걸 찾아도 좋다.

from math import sqrt
j, cont = 1, 1
while cont == 1:
    for l in range(1,j):
        c = 3*(l*l + j*j) - l - j
        k = (1 + sqrt(1 + 12*c))/6
        if k == int(k):
            c2 = 3*(k*k + j*j) - k - j
            m = (1 + sqrt(1 + 12*c2))/6
            if m == int(m):
                print l*(3*l-1)/2
                cont = 0
            c2 = 3*(k*k + l*l) - k - l
            m = (1 + sqrt(1 + 12*c2))/6
            if m == int(m):
                print j*(3*j-1)/2
                cont = 0
    j += 1

1.6초. 예상보다 큰 숫자에서 답이 나왔으니 이 정도 시간이면 받아들이기로 하자.

Saturday, May 23, 2015

Project Euler #43 - sub-string divisibility

프로젝트 오일러 43번
2,3,4번째 자리는 2로 나누어 떨어지고, 3,4,5번째 자리는 3으로, 4,5,6번째 자리는 7로... 나누어 떨어지는 10자리 pandigital 수를 모두 더하면?

답을 아는 입장에서.. 그런 조건을 만족하는 경우는 6개 밖에 안 되고, 엑셀을 써서 답을 찾을 수도 있다.
문제의 의도가 이게 맞는지 모르겠지만, 좀 무식한 코딩으로...
(보기에 안 좋아서 그렇지.. 이런 스타일이 속도는 괜찮다. 앞의 동전 문제에서도 마찬가지)

sm = 0
for d8 in range(17,1000,17):
    if d8 < 100 or len(set(`d8`)) == 3:
        r7 = set('1234567890') - set(`d8`)
        for d7 in r7:
            if (int(d7)*100 + d8/10) % 13 == 0:
                r6 = r7 - {d7}
                for d6 in r6:
                    if (int(d6)*100 + int(d7)*10 + d8/100) %11 == 0:
                        r5 = r6 - {d6}
                        for d5 in r5:
                            if (int(d5)*100 + int(d6)*10 + int(d7)) % 7 == 0:
                                r4 = r5 - {d5}
                                for d4 in r4:
                                    if d6 in {'5','0'}:
                                        r3 = r4 - {d4}
                                        for d3 in r3:
                                            if (int(d3)+int(d4)+int(d5)) % 3 == 0:
                                                r2 = r3 - {d3}
                                                for d2 in r2:
                                                    if d4 in set('24680'):
                                                        d1 = list(r2-{d2})[0]
                                                        sm += int(d1+d2+d3+d4+d5+d6+d7+`d8`)

print sm


41번 문제에서 썼던 permutation을 활용하면 코드가 깔끔해진다.
from itertools import permutations

sm = 0
for p in permutations('1234567890'):
    n = ''.join(p)    
    if int(n[7:]) % 17 == 0 and \
        int(n[6:9]) % 13 == 0 and \
        int(n[5:8]) % 11 == 0 and \
        int(n[4:7]) % 7 == 0 and \
        int(n[3:6]) % 5 == 0 and \
        int(n[2:5]) % 3 == 0 and \
        int(n[1:4]) % 2 == 0:
        sm += int(n)

print sm

그런데, 처음 코드는  0.0014초, 두번째 코드는 4.2초. 이 정도 차이면 지저분한 코드쪽이 나아 보인다.

Friday, May 22, 2015

Project Euler #42 - Coded triangle numbers

프로젝트 오일러 42번
알파벳을 숫자로 변환했을 때, 단어의 알파벳 숫자합이 삼각수가 되는 경우는 몇 번인가?
A => 1, B => 2, ... S => 19.. 
이렇게 할당하면 알파벳 합은 ABC => 6, SKY => 55가 된다. 첨부한 파일의 대문자 단어들 중 알파벳 합이 삼각수(1,2,3,6,10,...)이 되는 경우는 몇 번?

ord('A') = 65, ord('S') = 83, ... 을 이용하면 쉽다.

words = open('Euler/p042_words.txt').read().strip('"').split('","')
#print max([len(w) for w in words])

tri, n = {1}, 3
for i in range(3,30):
    n += i; tri.add(n)

cnt = 0
for word in words:
    nm = 0
    for s in word:
        nm += ord(s)-64
    if nm in tri: cnt += 1

print cnt




Thursday, May 21, 2015

Project Euler #41 - largest pandigital prime

프로젝트 오일러 41번
1,2,..,n 이 한번씩만 나오는 n자리수 소수 중 가장 큰 것은?

1,2,3,...,9를 한번씩 사용하는 9자리 소수가 존재할까?
=> 각 자리수를 더하면 45, 3의 배수이므로 소수가 될 수 없다.
그럼 8자리 소수는?
=> 각 자리수를 더하면 36, 역시 소수가 될 수 없다.
7자리 소수는? => 찾아보자!
6자리, 5자리 소수는? => 역시 3의 배수

7자리에서 찾아보고, 안 되면 4자리수를 찾아봐야지~

코딩하기 쉬운 방법: 7자리 모든 소수를 찾고, 8이나 9가 들어간 걸 빼고, pandigital인지 체크하고, 이들 중 가장 큰 것을 찾는다. => 천만까지의 소수를 찾아야 한다. 약 7^7 만큼의 수(에서 소수인 것들 중에서) pandigital인지 체크해야 한다.

연산을 줄일 수 있는 방법: 1,2,..,7의 순열(permutation)으로 만들 수 있는 모든 수에 대해서 소수인지 체크한다. 이왕이면 큰 수부터 나오도록 순열을 만든다. => 소수인지 '체크'를 위해서는 sqrt(7654321) 이하의 소수를 미리 찾아두어야 한다.
약 0.023초 걸린다.

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 isprime(n, primes):
    for p in primes:
        if n % p == 0: return False
        if p*p > n: return True
    return True

def perm(lst):
    if len(lst) == 1:
        return [(lst[0],)]
    rst = []
    for i, e in enumerate(lst):
        rst += [(e,)+x for x in perm(lst[:i]+lst[i+1:])]
    return rst

pn = rwh_primes1(int(7654321**.5))
seven = '7654321'
for pm in perm(range(7)):
    n = int(''.join([seven[pm[i]] for i in range(7)]))
    if isprime(n,pn):
        print n
        break

permutation의 13번째만에 필요한 답이 나왔다. perm함수를 리스트로 return하는 대신 하나씩 yield하는 방법을 찾으면 더 좋을텐데(지금은 쓸데 없이 7!=5040길이의 리스트를 return), 재귀함수에 yield를 어떻게 쓰는지 모르겠다.

perm함수를 정의하는 대신, itertools의 permutations함수를 가져오면 이 문제가 해결되는 듯 하다. 실행시간이 1/50로 줄어든다.

from itertools import permutations as pm

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 isprime(n, primes):
    for p in primes:
        if n % p == 0: return False
        if p*p > n: return True
    return True

pn = rwh_primes1(int(7654321**.5))
for pm in pm('7654321'):
    if isprime(int(''.join(pm)),pn):
        print ''.join(pm)
        break