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