Project Euler 250

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

Q250.
{1¹, 2², 3³, ..., 250250²⁵⁰²⁵⁰}の空でない部分集合で、その要素が250で割り切れるものの個数を、下16桁で。

これも前問とほとんど同じ。
nⁿの250の剰余は、(n%250)^kが250で割って1余る最小のkをあらかじめリストにしておくと速い。

import time
def pow_core(n, e):

    if e == 0:

        return 1

    

    m = pow_core(n, e / 2)

    if e % 2 == 1:

        return (m * m * n) % N

    else:

        return (m * m) % N
def pow(n, e):

    n %= N

    if min_e[n][0] < e:

        return pow_core(n, e) % N

    else:

        e = (e - min_e[n][0]) % (min_e[n][1] - min_e[n][0])

        return min_e[n][2] * pow_core(n, e) % N
def calc_min_e(n):

    m = n

    p = [ 0, n ]

    for e in range(2, N + 1):

        m = m * n % N

        if p.count(m):

            k = p.index(m)

            return (k, e, m)

        p.append(m)

    return -1
t = time.clock()

M = 250250

N = 250

L = 10 ** 16
min_e = [ calc_min_e(k) for k in range(N) ]

prev = [1] + [ 0 ] * (N - 1)

for n in xrange(1, M + 1):

    if n % 1000 == 0:

        print n

    S = prev[:]     # copy

    m = pow(n, n)

    for k in filter(lambda k: prev[k], range(N)):

        e = prev[k]

        S[(k+m)%N] += e

    for k in range(N):

        if S[k] >= L:

            S[k] %= L

    prev = S

print (S[0] - 1) % L

print time.clock() - t