Project Euler 29（2）

a^b（2 ≤ a ≤ 100, 2 ≤ b ≤ 100）は99×99個あるのですが、同じ数になるものがあるので問題になっているのです。例えば、

4² = 2⁴

これらを99×99から取り除けばよいです。
ここでは、底が小さくなれるものを数えます。上の例なら、

4² = (2²)² = 2⁴

すなわち、底がべき乗数ならこういうことが可能なわけです。
まず、10以下で考えましょう。べき乗数は、

4 8 9

です。4については、

4² = 2⁴
4³ = 2⁶
4⁴ = 2⁸
4⁵ = 2¹⁰

の4つが底を小さくできます。9についても、

9² = 3⁴
9³ = 3⁶
9⁴ = 3⁸
9⁵ = 3¹⁰

と同様になります。元の底の指数が最大2のものはこうなります。すなわち、指数だけで取り除ける個数は決まるわけです。さて、8については、

8² = 2⁶
8³ = 2⁹

のほかに、

8⁴ = 4⁶
8⁶ = 4⁹

があります。答えは、81 - (4 + 4 + 4) = 69となります。
もう少し例を見てみましょう。2⁶ = 64を取り上げます。代わりの底としては、2,4,8,16,32が考えられます。100までとして、

64² = 2¹² 64³ = 2¹⁸ … 64¹⁶ = 2⁹⁶
64² = 4⁶ 64³ = 4⁹ … 64³³ = 2⁹⁹
64² = 8⁴ 64³ = 8⁶ … 64⁵⁰ = 2¹⁰⁰
64² = 16³ 64⁴ = 8⁶ … 64⁶⁶ = 2⁹⁹
64⁵ = 32⁶ 64¹⁰ = 8¹² … 64⁸⁰ = 2⁹⁶

が重複していることになります。しかしこの重複もまた重複していて、結局、64の指数だけ示すと、

2 〜 50 52 54 55 56 58 60 62 64 66 70 75 80

となります。
まとめると、64は64^1/6,64^2/6,64^3/6,64^4/6,64^5/6に底を変えることができて、指数の上限は、

64ⁿ ≤ (64^k/6)¹⁰⁰

だから、

n ≤ 100k / 6

となります。また、k/6を約分してp / qとなったとき、64の指数はpの倍数となります。k/6を約分して並べると、

1/6 1/3 1/2 2/3 5/6

上の考察より、分子が同じなら分母が最も小さい指数以外は意味がなくなります。結局、

1/2 2/3 5/6

だけ考えればいいことになります。100までの配列を作り、これらに対して出てくる64の指数を記録していきます。
下のコードでは100万まででも十分に速く答えが出ます。

from itertools import count

from fractions import gcd
def index(a, n):

    for k in xrange(len(a)):

        if a[k][0] == n:

            return k

    return -1
def pow_numbers():

    ary = [ ]

    for a in count(2):

        n = a

        for b in count(2):

            n *= a

            if n > N:

                break

            if index(ary, n) == -1:

                ary.append( (n, b) )

        if b == 2:

            break

    

    return ary
# [ (4, 2), (9, 2), (8, 3) ] -> [ 0, 0, 2, 1 ]

def convert(a):

    b = [ 0 ] * (max_e + 1)

    for e in a:

        b[e[1]] += 1

    return b
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):

        d = gcd(k, e)

        k1 = k / d

        e1 = e / d

        a[k1] = e1

    return a
def num_overlap(e):

    b = max_denominator(e)

    a = [ 0 ] * (N + 1)

    for k, e in b.iteritems():

        for l in xrange(max(2, k), N * k / e + 1, k):

            a[l] = 1

    return sum(a)
def solve():

    return (N - 1) ** 2 \

        - sum(map(lambda e: num_overlap(e) * pows[e], range(2, max_e + 1)))
N = 10 ** 6

ary = pow_numbers()

max_e = max(map(lambda e: e[1], ary))

pows = convert(ary)

print solve()