Project Euler #25 - 피보나치 수열은 몇번째 항에서 1000자리를 넘어가는가?

피보나치 수열에서 1000자리를 넘어가는 첫번째 항은?

포럼에는 역시 기발한 방법들이 있지만, 여기서는 나이브하게...

i, t, f = 1, 0, 1

while len(str(f)) < 1000:
    i += 1
    t, f = f, f+t

print i

Project Euler #24 - 100만번째 permutation

0~9의 permutation 으로 만들 수 있는 숫자들을 차례로 늘어놓았을 때, 100만번째 수는?

1번째는 0123456789

9!+1번째는 1023456789

9!*2+1번째는 2013456789

그럼 제일 앞에 2를 고정하고, 9개의 숫자로 permutation을 만들었을 때  100만 - 9!*2 번째 숫자를 찾는 문제로 바뀐다.

코드로 표현하면,

factorial = {1:1}
for i in range(2,10):
    factorial[i] = i*factorial[i-1]

answer = ''
d = range(10)
N = 1000000-1
for i in range(9,0,-1):
    answer += str(d.pop(N / factorial[i]))
    N %= factorial[i]

print answer + str(d[0])

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 #22 - 이름 숫자

알파벳을 숫자로 변환하고 주어진 공식으로 순위와의 가중합을 구하는 문제

포럼에 올라온 글에 자극받아서 아예 한 줄 코드로 만들어 봤다.

print sum([(i+1)*sum([ord(c)-64 for c in name]) for i,name in enumerate(sorted([name.strip('"') for name in open('p022_names.txt').read().split(',')]))])

ord는 문자에 해당하는 아스키번호를 알려주는 함수라고 한다.

아래와 같이 정상적인(?) 코드를 작성할 수도 있다.

ABC = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
alphanum = {}
for i,c in enumerate(ABC):
    alphanum[c] = i + 1

def calcname(name):
    return sum([alphanum[s] for s in name])

with open('/Users/Dongug/Downloads/p022_names.txt') as f:
    names = [name.strip('"') for name in f.readline().split(',')]

rs = 0
for i,name in enumerate(sorted(names)):
    rs += (i + 1) * calcname(name)

print rs

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 #20 - 100!의 숫자 합

100!의 각 자리수의 합은?

파이썬이 아니라면 고민스럽겠지만...

from operator import mul
print sum([int(s) for s in str(reduce(mul,range(1,100)))])

마지막에 range(1,101)보다는 range(1,100)이 조금 더 이뻐 보여서..

Project Euler #19 - 1일 일요일이 몇 번?

1901년부터 2000년까지 매달 1일이 일요일인 경우는 몇 번?

1. 1200번 중 몇 번인지를 물어보는 문제인데, 100년이면 충분히 긴 시간이고 중심극한정리를 생각하면 정답은 1200/7 = 171.xx 근처에 있을 거다.

2. 1번으로 풀면 cheating인 듯 하고.. 이런 코드를 짜 봤다. 매달 1일이 무슨 요일인지 다 찾아보는 방법이다.

d = 2 #Monday - 1 Jan 1901
c = 0
for year in range(1901,2001):
    for mon in range(1,13):
        if d == 0: c += 1
        if mon in (1,3,5,7,8,10,12): d = (d + 3) % 7
        elif mon in (4,6,9,11): d = (d + 2) % 7
        elif year % 4 == 0: d = (d + 1) % 7

print c

3. 포럼에서 본 글 중 인상적인 것.
(일단, 윤년이 아닌 경우)
1월 1일이 일요일인 경우 그 해에 1일이 일요일은 몇 번?
1월 1일이 월요일인 경우 그 해에 1일이 일요일은 몇 번?

1월 1일이 화요일인 경우 그 해에 1일이 일요일은 몇 번?

...

1월 1일이 토요일인 경우 그 해에 1일이 일요일은 몇 번?

이걸 미리 찾아두고, 100년 동안 1월 1일이 일요일이었던 횟수, 월요일이었던 횟수, ... 를 세어 본다. 2번에서는 y*m번의 요일 계산이 필요했지만 3번은 y+m번의 요일 계산이면 충분하다.

Project Euler #18 - Maximum path sum I

최대합 경로 문제

처음엔 이렇게 풀었다.

st = '''75
95 64
17 47 82
18 35 87 10
20 04 82 47 65
19 01 23 75 03 34
88 02 77 73 07 63 67
99 65 04 28 06 16 70 92
41 41 26 56 83 40 80 70 33
41 48 72 33 47 32 37 16 94 29
53 71 44 65 25 43 91 52 97 51 14
70 11 33 28 77 73 17 78 39 68 17 57
91 71 52 38 17 14 91 43 58 50 27 29 48
63 66 04 68 89 53 67 30 73 16 69 87 40 31
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23'''

og = []
for line in st.split('\n'):
    og.append([int(e) for e in line.split()])

for i in range(len(og)-1):
    og[i+1][0] = og[i][0] + og[i+1][0]
    for j in range(len(og[i])-1):
        og[i+1][j+1] = max(og[i][j],og[i][j+1]) + og[i+1][j+1]
    og[i+1][-1] = og[i][-1] + og[i+1][-1]
    
print max(og[-1])

포럼 글을 읽다보니.. 삼각형을 뒤집으면 문제가 훨씬 쉽게 풀린다.

st = '''75
95 64
17 47 82
18 35 87 10
20 04 82 47 65
19 01 23 75 03 34
88 02 77 73 07 63 67
99 65 04 28 06 16 70 92
41 41 26 56 83 40 80 70 33
41 48 72 33 47 32 37 16 94 29
53 71 44 65 25 43 91 52 97 51 14
70 11 33 28 77 73 17 78 39 68 17 57
91 71 52 38 17 14 91 43 58 50 27 29 48
63 66 04 68 89 53 67 30 73 16 69 87 40 31
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23'''

og = [[int(e) for e in line.split()] for line in st.split('\n')][::-1]

l = len(og)
for i in range(l):
    for j in range(l - i-1):
        og[i+1][j] += max(og[i][j],og[i][j+1])

print og[-1][-1]

처음 이런 코드를 작성할 때는 어떤 리스트의 몇번째 위치에 어떤 값이 들어 있는지 헷갈려서 머리 싸매고 했었는데.. 100번까지 풀고 블로그 정리용으로 다시 작성하니까 훨씬 쉽다. 조금이나마 발전이 있는 것 같아서 기분 up.

Project Euler #17 - 숫자 읽기

1000까지의 숫자를 영어로 쓸 때 필요한 알파벳의 수는?

성실과 근면을 요구하는 문제.. forty인지 fourty인지 헷갈리는 ESL 입장에선 약간 짜증나는 문제.

d1 = ['one','two','three','four','five','six','seven','eight','nine']
d2 = ['ten','eleven','twelve','thirteen','fourteen','fifteen','sixteen','seventeen','eighteen','nineteen']
d3 = ['twenty','thirty','forty','fifty','sixty','seventy','eighty','ninety']
d4 = ['hundred','and','thousand']

l1 = [len(e) for e in d1]
l2 = [len(e) for e in d2]
l3 = [len(e) for e in d3]
l4 = [len(e) for e in d4]
# 1~99
s1 = sum(l1)*9 + sum(l2) + sum(l3)*10

# 1~99 101~199 ... 901~999
s2 = s1*10 + sum(l1)*99 + sum(l4[:2])*99*9

# 1~1000
s3 = s2 + sum(l1) +  l4[0]*9 + l1[0] + l4[2]
print s3

글타래를 읽다 보니.. 요런 스타일의 코드가 있다.

def read(n):
    ...
    return ...

sm = 0
for i in range(1,1001):
    sm += len(read(n))

print sm

귀찮은 정도는 똑같지만, 답이 틀렸을 때 틀린 부분을 찾기가 좀 더 쉬울 것 같다.

Project Euler #16 - Power digit sum

2^1000 의 숫자들을 다 더하면?

300자리 정도 되겠지만 큰 정수에 강한 파이썬을 믿어보자~

print sum([int(s) for s in str(2**1000)])

2**1000000 계산은 위의 코드보다 빠른 방법이 있는데.. 1000제곱 정도는 이것보다 빠른 방법을 못 찾겠다.

Project Euler #15 - Lattice paths

격자 이동 - 20*20인 격자의 왼쪽 위에서 오른쪽 아래로 가는 경로는 몇가지인가? 단, 오른쪽, 아래로만 이동할 수 있다.

오른쪽으로 20칸, 아래로 20칸 이동하는 문제이다. 40번 이동/20번 오른쪽/20번 아래쪽은 변함이 없다. 40칸의 자리가 있고, 그 중 오른쪽으로 이동할 20칸을 어떻게 할당하느냐의 문제이고 간단히 40C20을 계산하면 된다.

40!/20!/20! 계산 문제인데.. 40!는 큰 숫자이고 다른 언어라면 큰 숫자의 메모리/정확도 문제 때문에.. 40/20*39/13/3*38/19/2*.... 이런 걸 만들어야 할 수도 있다. 하지만 파이썬은 큰 정수를 다루는 데 전혀~ 문제가 없다.

from operator import mul
print reduce(mul,range(21,41))/reduce(mul,range(1,21))

파이썬 만세!

Project Euler #14 - Longest Collatz sequence

Longest Collatz sequence

13->40->20->10->5->16 ... 이렇게 해서 13:10을 찾았다면, 40:9, 20:8 도 메모리에 저장해 두고, 나중에 20이나 40이 나왔을 때 다시 계산하지 않도록 한다. 이런 식이라면, 이미 10:7이 메모리에 있었을테니, 13에 대해 10번 계산이 아니라 세 번 계산만에 13:10을 찾을 수 있다.

def collatz(n,cDic):
    m = n
    i = 0
    tl = [n]
    while True:
        if m in cDic:
            for s, t in enumerate(tl[::-1]):
                cDic[t] = s + 1 + cDic[m]
            return
        if m % 2 == 0: m = m / 2
        else: m = 3 * m + 1
        tl.append(m)
        i += 1

cDic = {1:1}
for i in range(2,1000000):
    collatz(i,cDic)

print max(cDic, key=cDic.get)

이렇게 하면 4.3초...
그런데, 13:10만 저장하고 40:9, 20:8 을 저장 안 하면 속도가 더 빨라진다. 아래 코드는 2.7초..
tl 이라는 리스트에 값을 하나씩 더하고 다시 꺼내는 데 드는 비용이 cache 비용보다 더 많이 드나보다.

def collatz(n,cDic):
    m = n
    i = 0
    while True:
        if m in cDic:
            cDic[n] = i + cDic[m]
            return
        if m % 2 == 0: m = m / 2
        else: m = 3 * m + 1
        i += 1

cDic = {1:1}
for i in range(2,1000000):
    collatz(i,cDic)

print max(cDic, key=cDic.get)

Project Euler #13 - 큰 수의 합

50자리 숫자 100개를 더했을 때 제일 앞자리 숫자 10개는?

파이썬이 큰 정수 다룰 때 얼마나 편한지 보여주는 좋은 예!

st = '''37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538

...

53503534226472524250874054075591789781264330331690'''

print str(sum([int(s) for s in st.split()]))[:10]

---
파이썬에서는 큰 정수를 편하고 안전하게 쓸 수 있다는 점을 다른 곳에 활용할 수도 있다. 예를 들어 sqrt(2)를 소수점 이하 100자리까지 구할 때, float로 다루는 것보다 2*10^200 의 sqrt를 찾는 게 더 쉽다.

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