Project Euler 200

プロジェクトオイラー
http://projecteuler.net/index.php

Q200.
p²q³(pとqは異なる素数）をsqubeと呼ぶことにする。また、200のようにどの一つの桁の数字を変えても素数にならない数をprime-proofと呼ぶことにする。prime-proofであるsqubeで1992008のように"200"を含む数で200番目のものを求めよ。

ほとんど問題文をそのままコードにするだけ、のはずだったが、squbeを生成するコードを一文字間違えていて、それをなかなか見つけられなかった。

def is_prime(n):

    if n < 2:

        return False

    for p in primes:

        if p * p > n:

            return True

        elif n % p == 0:

            return False

    return True
def make_primes(n):

    m = int((n + 0.5) ** 0.5)

    for p in xrange(3, m + 1, 2):

        if is_prime(p):

            primes.append(p)

    

    a = [ 0 ] * ((n + 15) / 16)

    for p in primes:

        if p * p > n:

            break

        if p == 2:

            continue

        for k in xrange(3 * p, n, 2 * p):

            a[k>>4] |= 1 << (k & 15)

    

    for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):

        if (a[k>>4] & (1 << (k & 15))) == 0:

            primes.append(k)
def calc_sqube(i, j):

    return primes[i] ** 2 * primes[j] ** 3
def next_sqube(a, i):

    j = a[i][0] + 1

    if j == i:

        j += 1

    a[i] = (j, calc_sqube(j, i))
def new_row(a, i):

    a.append( (0, calc_sqube(0, i)) )
def gen_sqube():

    a = [(1, calc_sqube(1, 0)) ]

    yield a[0][1]

    next_sqube(a, 0)

    new_row(a, 1)

    

    while True:

        i = reduce(lambda x, y: min(x, y,

                        key = lambda k: a[k][1]), xrange(len(a)))

        yield a[i][1]

        next_sqube(a, i)

        if i == len(a) - 1:

            new_row(a, i + 1)
def encode(n):

    a = [ ]

    while n:

        a.append(n % 10)

        n /= 10

    a.reverse()

    return a
def decode(a):

    n = 0

    for d in a:

        n = n * 10 + d

    return n
def is_prime_proof(n):

    a = encode(n)

    for k in range(len(a)):

        orig_d = a[k]

        for d in range(10):

            a[k] = d

            if is_prime(decode(a)):

                return False

        a[k] = orig_d

    return True
N = 200

subs = str(N)

primes = [ 2 ]

make_primes(10 ** 6)

counter = 0

for sq in gen_sqube():

    if str(sq).find(subs) != -1 and is_prime_proof(sq):

        counter += 1

        print counter, sq

        if counter == N:

            break