Wednesday, June 17, 2015

Project Euler #67 - Maximum path sum II

프로젝트 오일러 67번 - 최대합 경로 두번째 문제
첨부한 파일에서 경로를 따라 가며 숫자를 더했을 때 최대합은?

18번 문제와 완전히 같은 문제. 18번은 숫자가 몇 개 안 되어서 '모든 경로'를 검토할 수 있었지만 여기서는 그게 불가능하다는 정도의 차이.
풀이방법도 18번에서 설명했던 것과 같다.

tri = [map(int, line.strip().split(' ')) for line in open('p067_triangle.txt').readlines()][::-1]

for i in range(1,100):
    for j in range(100-i):
        tri[i][j] += max(tri[i-1][j],tri[i-1][j+1])

print tri[99][0]

0.006초.

Tuesday, June 16, 2015

Project Euler #66 - 디오판틴 방정식

프로젝트 오일러 66번 - Diophantine equation
1000보다 작은 D 중에서(단, D는 제곱수가 아니다)
x^2 - D y^2 = 1 의 최소정수해 중 x를 가장 크게 만드는 것은?

모든 D에 대해서 x의 최소해를 구해 봐야 한다.
D를 고정하고 x와 y를 늘려가며 찾아보는 건 쉬운데.. 답이 안 나온다. 한참 고민하다가 구글링..
D = 61일 때의 최소해는 x=1766319049, y = 226153980 이라고 한다. brute force는 포기. 인터넷 글을 읽다 보니, 페르마도 풀려다 못 풀었다고 한다. 페르마가 못 푼 문제를 내가 풀 수 있을 것 같지는 않고.. 풀이 방법까지 찾아보게 됐다.

디오판틴 방정식 중에서 특별히 이런 모양의 방정식을 펠방정식 - Pell's equation이라고 부른다는 것도 알게 됐다.
http://en.wikipedia.org/wiki/Chakravala_method#Examples
여기에 있는 풀이 방법을 따랐다. 우변이 ±1, ±2, or ±4 일 때는 좀 더 쉬운 특별한 방법이 있다고 하지만, 그 '특별한 방법'을 쓰지 않아도 iteration을 조금 더 진행하면 답이 나온다고 위키피디아에 친절하게 나와 있다.

#http://en.wikipedia.org/wiki/Chakravala_method#Examples
from math import sqrt

def dio_init(D):
    Di = int(sqrt(D))
    tmpdic = {Di:abs(Di**2-D),Di+1:abs((Di+1)**2-D)}
    a = min(tmpdic,key=tmpdic.get)
    return a,1,a*a-D

def dio_update(a,b,k,D):
    tmpdic = {}
    for m in range(int(sqrt(D))-abs(k),int(sqrt(D))+abs(k)+1):
        if (a + b*m) % abs(k) == 0:
            tmpdic[m] = abs(m*m - D)
    m = min(tmpdic,key=tmpdic.get)
    return (a*m+D*b)/abs(k),(a+b*m)/abs(k),(m*m-D)/k

def dio(D):
    a,b,k = dio_init(D)
    while True:
        a,b,k = dio_update(a,b,k,D)
        if k == 1: return a, b

mx, argmx = 0,0
for D in [x for x in range(2,1001) if x not in [y*y for y in range(1,32)]]:
    if dio(D)[0] > mx: mx = dio(D)[0]; argmx = D

print argmx

ps. x**.5 와 sqrt(x)는 같다. 후자는 from math import sqrt가 필요해서 보통은 0.5제곱하는 쪽을 선호한다. 하지만 가끔, sqrt로 쓰는 게 더 잘 읽히는 날은 sqrt로 쓴다.

psps. 수학이나 알고리즘 공부가 아니라 python 코딩공부를 목적으로 프로젝트 오일러 문제를 풀고 있었다. 이 문제를 접하고 나서 동기부여가 안 된다. 66번에서 그만 두기는 거시기하니까 100번까지 하고 그만 두는 걸로...



Monday, June 15, 2015

Project Euler #65 - e의 수렴

프로젝트 오일러 65번 - convergent of e

e 를 64번 문제와 같은 continued fraction 형태로 표현하면 1,2,1, 1,4,1, 1,6,1, ... 형태로 전개된다. 100까지 전개한 다음 분수로 표현했을 때 분자의 자리수의 합은?

가장 작은 부분부터 차례로 분수를 풀어보면 된다. 이 때 1+1/n 형태를 (n+1)/n 으로 풀거나 뒤집어도 약분되지 않는다. 마음놓고 100번째까지 풀어 보면,

nume,deno=0,1
for k in range(33,0,-1):
    nume,deno=deno,deno+nume
    nume,deno=deno,deno*2*k+nume
    nume,deno=deno,deno+nume
nume = 2*deno+nume
print sum([int(x) for x in str(nume)])


Sunday, June 14, 2015

Project Euler #64 - Odd period square roots

프로젝트 오일러 64번 - 제곱근을 분수로 풀었을 때 주기가 홀수인 경우는?

모든 수는 무한분수(continued fraction을 뭐라고 번역하는지 모르겠다)로 표현할 수 있다. 자연수의 제곱근은 continued fraction에서 일정한 주기를 갖는데, 이 주기가 짝수일 수도 있고 홀수일 수도 있다. 10000 이하 자연수 중 주기가 홀수인 경우는 몇 번?

continued fraction은 57번에서 처음 나왔다. 이 때는 대충 넘어갈 수 있었는데, 64번과 그 이후 문제들을 풀려면 continued fraction만드는 법을 제대로 알아야 한다.
문제에서 예로 든 sqrt(23)의 전개가 바로 이해가 되면 good, 그렇지 않으면 아래의 과정을 차근차근 따라가 보는 걸로..
(열심히 타이핑하다가 너무 번거로워서 연습장에 쓰는 걸로.. )


분자, 분모를 뒤집고, 자기보다 작은 가장 큰 수를 찾고, 켤레수를 곱하고, 약분하는 과정을 반복하다가 앞에서 나왔던 것과 같은 패턴을 찾으면 끝.
한가지.. 언제나 약분이 되니까 소인수 분해 과정 거치지 않고 나누어 버렸다.(왜 그런지는 모르겠다)

이걸 코드로 구현하면 이렇게 된다.

from math import sqrt

# (sqrt(a)-b)/c
def findnext(a,b,c):
    h = int(c/(sqrt(a)-b))
    return h, a, h*(a-b*b)/c - b, (a-b*b)/c

def period(n):
    if sqrt(n) == int(sqrt(n)): return 0
    rs = 0
    first = int(sqrt(n))
    a,b,c = (n,first,1)
    while (a,b,c)!=(n,first,1) or rs == 0:
        h,a,b,c = findnext(a,b,c)
        rs+= 1
    return rs%2

result = 0
for i in range(2,10001):
    result += period(i)

print result

실행 시간은 0.35초

Saturday, June 13, 2015

Project Euler #63 - powerful digit counts

프로젝트 오일러 63번 - 제곱 자리수 세기
7^5은 5자리, 8^9은 9자리처럼 어떤 수의 n제곱이 n자리가 되는 경우는 총 몇 가지인가?

1^1, 2^1, ... 9^1
1^2 안 되고 ... 4^2 되고.. 9^2까지 모두 되고
4^3 안 되고, 5^3 되고.. 9^3까지 모두 되고
5^4 안 되고...

두 가지 성질..
1. k^n이 되면 k+1, k+2, ... 9까지 다 된다
2. k^n이 되었으면, k^(n+1)부터 체크하면 된다. (n+1)을 위해 1, 2, ... (k-1)은 체크할 필요 없다.

k, n, rs = 1, 1, 0
while k < 10:
    if len(str(k**n)) == n: rs += 10-k; n += 1
    else: k += 1
print rs

실행시간은.. 순식간.

Friday, June 12, 2015

Project Euler #62 - 세제곱 순열

프로젝트 오일러 62번 - cubic permutations

345^3 = 41063625 의 자리수를 바꾸면, 56623104 = 384^3 도 되고, 66430125 = 405^3도 된다. 실제로, 345^3은 이런 성질을 갖는 가장 작은 수이다.
5개의 permutation을 갖는 가장 작은 세제곱수는?


계속 문제가 길어지고 코드는 지저분해져서 재미가 떨어지고 있었는데 오랜만에 쉽게 풀리는 문제를 만났다. 60번, 61번에서 시간을 너무 많이 써서 더 그런 느낌..

그냥 n=1,2,3,4,5.. 에 대해 세제곱수들의 개수를 세다가 5가 되면 return. 단, permutation에 invariant하게 저장해야 하므로 각 자리수의 숫자들을 정렬하여 숫자를 세었다.

n = 1
rs = {}
while True:
    tmp = ''.join(sorted(`n**3`))
    if tmp in rs: rs[tmp].add(n)
    else: rs[tmp] = {n}
    if len(rs[tmp]) == 5:
        print min(rs[tmp])**3
        break
    n += 1

0.05초.
좀 걸리는 부분이 있다. 문제에서 "which exactly five"라고 했는데 위의 코드로는 6개가 되어도 5개 찾은 순간 return하게 된다. 엄밀하게 하려면, 자리수가 바뀌기 전까지 break하지 말고 계속 체크해야 하지만.. five permutations도 매우 드문데, 갑자기 six permutations는 되지 않을 것 같아서 굳이 코드에 넣지 않았다.

Thursday, June 11, 2015

Project Euler #61 - cyclical figurate numbers

프로젝트 오일러 61번 - 순환 다항 수
세 숫자 8128, 2882, 8281 는 두가지 독특한 성질이 있다.
1. 두자리씩 끊었을 때 28, 82, 81 이렇게 연결이 되는 성질(cyclical)
2. 8128은 3각수  - n(n+1)/2의 형태, 8281은 4각수 - n^2의 형태, 2882는 5각수 - n(3n-1)/2의 형태태

네자리 숫자 6개로 아래의 성질을 갖는 경우가 단 하나 있다. 이 6개 숫자의 합은?
1. 두자리씩 끊었을 때 연결된다
2. 숫자 6개가 각각 삼각, 사각, 오각, 육각, 칠각, 팔각수


예시에도 있지만 cycle의 순서는 아직 모른다. 3각->4각->5각->... 일 수도 있고, 3각->8각->4각->...일 수도 있다. 순환의 속성상 시작을 어디로 하든 마음대로라는 점은 다행이다. 8각수의 개수가 제일 적으니 거기에서 시작해 보자.
모든 8각수에 대해, 뒤 2자리(꼬리)를 끊고,
이 2자리로 시작하는(머리가 되는) 모든 케이스를 3각~7각에서 찾고
이 때의 꼬리를 찾고,
이 꼬리가 머리가 되는 모든 케이스를 3각~7각 중에서 위에서 선택하지 않았던 곳에서 찾고,
이 때의 꼬리를 찾고,
이 꼬리가 머리가 되는.....
이렇게 3각~7각을 한번씩 돌고 나면 순환이 완성된다.

늘 그렇듯 긴 코드가 실행시간을 줄여준다.

p, ph, pt = [],[],[]

p.append([`n*(n+1)/2`   for n in range(1,200) if 1000 <= n*(n+1)/2   < 10000])
p.append([`n*n`         for n in range(1,200) if 1000 <= n*n         < 10000])
p.append([`n*(3*n-1)/2` for n in range(1,200) if 1000 <= n*(3*n-1)/2 < 10000])
p.append([`n*(2*n-1)`   for n in range(1,200) if 1000 <= n*(2*n-1)   < 10000])
p.append([`n*(5*n-3)/2` for n in range(1,200) if 1000 <= n*(5*n-3)/2 < 10000])
p.append([`n*(3*n-2)`   for n in range(1,200) if 1000 <= n*(3*n-2)   < 10000])

big = {}
#{'10': {0: ['35', '81'],  1: ['24', '89']}, '11': {0: ['28', '76'],  4: ['60']},...}

for i in range(5):
    for e in p[i]:
        head, tail = e[:2], e[-2:]
        if tail[0] != '0':
            if head in big:
                if i in big[head]: big[head][i].append(tail)
                else: big[head][i] = [tail]
            else:
                big[head] = {}
                big[head][i] = [tail]


def prob61(p,big):
    for e5 in p[5]:
        e5h, e5t = e5[:2], e5[-2:]
        if e5t in big:
            for s1 in big[e5t]:
                for h1 in big[e5t][s1]:
                    if h1 in big:
                        for s2 in big[h1]:
                            if s2 not in [s1]:
                                for h2 in big[h1][s2]:
                                    if h2 in big:
                                        for s3 in big[h2]:
                                            if s3 not in [s1,s2]:
                                                for h3 in big[h2][s3]:
                                                    if h3 in big:
                                                        for s4 in big[h3]:
                                                            if s4 not in [s1,s2,s3]:
                                                                for h4 in big[h3][s4]:
                                                                    if h4 in big:
                                                                        for s5 in big[h4]:
                                                                            if s5 not in [s1,s2,s3,s4]:
                                                                                for h5 in big[h4][s5]:
                                                                                    if h5 == e5h:
                                                                                        return int(e5)+int(e5t+h1)+int(h1+h2)+int(h2+h3)+int(h3+h4)+int(h4+h5)
                                                                                    
print prob61(p,big)

0.0015초