Showing posts with label 소수만들기. Show all posts
Showing posts with label 소수만들기. Show all posts

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


Monday, April 20, 2015

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


Sunday, April 19, 2015

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)


Saturday, April 18, 2015

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


Tuesday, April 14, 2015

파이썬으로 소수 만들기

소수를 만드는 가장 좋은 방법은 "에라스토테네스의 체"입니다. 한글 위키피디아 페이지에서 오른쪽 상단의 움직이는 이미지를 보면 충분히 이해할 수 있을 겁니다.

이런 내용을 파이썬 코드로 구현하면 이렇게 됩니다. n 이하의 모든 소수를 리스트로 만들어 주는 거죠.

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

로버트 윌리엄 행스(rwh)라는 사람이 stackoverflow에 올린 걸 가져온 겁니다. n이하 홀수들 중에서 체(seive)로 걸러지지 않는 것을 매우 효율적으로 찾아줍니다. 위의 url에 들어가 보면, 소수를 찾는 다른 방법들도 나와 있습니다. 각 알고리즘을 비교하는 표도 찾을 수 있고요. 저는 이 함수가 길이도 짧고, 의미하는 바도 간결하고, 실행 속도도 매우 빨라서 요즘에는 이것만 씁니다.


=====

이전에는 http://www.python-course.eu/list_comprehension.php 요기 페이지 마지막 예제로 나온 함수를 사용했었습니다. 위의 rwh_primes1보다 속도도 느리고, 무엇보다 버그가 있습니다. primes(50)을 찾아보면 49를 소수라고 판단하죠. 버그 고치라고 사이트 관리자에게 메일을 보낸 지 한참 되었는데 안 고치네요.