Project Euler 29(4)

この問題は、べき乗数を列挙する、重複する数を数える、という2段階に分かれています。今まで前段は後段にくらべて時間がかからなかったのですが、後段が速くなって前段のほうが時間がかかるようになったので前段を改良します。
例えば100までのべき乗数を考えます。ほしいのは指数ごとの個数です。具体的にべき乗数を列挙する必要はありません。べき乗の指数は2〜6です。2乗を考えると、22〜102となります。ただし、4,8,9はべき乗数ですので、これらは取り除かなければなりません。取り除く個数は自分自身を呼び出して得られます。このように簡単に再帰的に求められます。
これで、107まででも十分に速く求められます。
後段はなかなか速くならないようです。



from itertools import count, imap, takewhile
from fractions import gcd, Fraction

def index(a, n):
for k in xrange(len(a)):
if a[k][0] == n:
return k
return -1

def pow_numbers(max_n):
if max_n <= 3:
return [ 0 ]

ary = [ 0, 0]
max_e = max(takewhile(lambda e: N >> e, count(1)))
for e in range(2, max_e + 1):
max_m = int( (max_n + 0.5) ** (1. / e) )
ary.append(max_m - sum(pow_numbers(max_m)) - 1)

return ary

def num_overlap(e):
a = [ 0 ] * (N + 1)
for k in range(1, e):
d = gcd(e, k)
e1 = e / d
k1 = k / d
for l in xrange(max(2, k1), N * k1 / e1 + 1, k1):
a[l] = 1
return sum(a)

def max_denominator(e):
a = { }
for k in range(1, e):
f = Fraction(k, e)
k1 = f.numerator
a[k1] = f

b = [ f for f in a.itervalues() ]
b.sort()

while True:
if b[0].numerator != 1:
b.pop(0)
else:
break

return b

def degenerate_nums(nums):
for k in range(len(nums) - 1, 1, -1):
for n in nums[:k]:
if nums[k] % n == 0:
nums.pop(k)
break

def lcm(a, b):
a2 = abs(a)
m = a2 * b / gcd(a2, b)
return m if a > 0 else -m

def gen_divisors(ary_d, k0 = 0):
if k0:
yield -1

for k in range(k0, len(ary_d)):
for d in gen_divisors(ary_d, k + 1):
yield -lcm(d, ary_d[k])

def count_range(lower, upper, d):
if d > 0:
plus = True
else:
d = -d
plus = False
s = upper / d - (lower + d - 1) / d + 1
return s if plus else -s

def num_overlap_range(b, k):
upper = int(N * b[k])
if k == 0:
return upper - 1
lower = int(N * b[k-1]) + 1
nums = [ f.numerator for f in b[k:] ]
degenerate_nums(nums)
return sum(imap(lambda d: count_range(lower, upper, d), gen_divisors(nums)))

def num_overlap(e):
b = max_denominator(e)
return sum([ num_overlap_range(b, k) for k in range(len(b)) ])

def solve():
max_e = max(takewhile(lambda e: N >> e, count(1)))
pows = pow_numbers(N)
a = [ num_overlap(e) * pows[e] for e in range(2, max_e + 1) ]
return (N - 1) ** 2 - sum(a)

N = 10 ** 9
print solve()