Project Euler 29（3） - inamori’s diary

例えば、16は16^1/4=2と16^1/2=4と16^3/4=8に底を移行できます。20までとして16が移行できる指数はそれぞれ、



    2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

1/4 o o o o x x x x  x  x  x  x  x  x  x  x  x  x  x

1/2 o o o o o o o o  o  x  x  x  x  x  x  x  x  x  x

3/4 x o x x o x x o  x  x  o  x  x  o  x  x  x  x  x

oがついているところが移行できる指数です。これを総合すると、



    2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

3/4 o o o o o o o o  o  x  o  x  x  o  x  x  x  x  x

12個が移行できることになります。しかし、このように数えるのは20までならいいですが、100万となると時間がかかりそうです。
もう少し複雑になる例、ここでは10乗を使います。2¹⁰は、2¹⁰の

1/10, 1/5, 3/10, 2/5, 1/2, 3/5, 7/10, 4/5, 9/10

乗の底に移行できます。分子が同じ分数は最も大きいものしか意味がないので、

2/5, 1/2, 3/5, 7/10, 4/5, 9/10

また、分子が1のものより小さい分数は意味がないので、

1/2, 3/5, 7/10, 4/5, 9/10

となります。1000までとして、領域を区切ります。2〜500, 501〜600, 601〜700, 701〜800, 801〜900とします。2〜500はすべてOKで、それを除いた各領域は、

501〜600 3または7または4または9の倍数
601〜700 7または4または9の倍数
701〜800 4または9の倍数
801〜900 9の倍数

となります。ただし、最初の領域で9は3で割り切れるので、これは除きます。
ここで、Problem 1で用いた手法を再び使います。最初の領域で、まず3の倍数、7の倍数、4の倍数の個数を数えて足し合わせます。その次に、21の倍数、28の倍数、12の倍数の個数を引きます。ここで、21は3と7の最小公倍数です。最後に、84の倍数を足します。
このようにすると、10⁷まででも十分に速く求められます。
この問題ももう少しで終わります。



from itertools import count, imap

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():

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

        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():

    a = [ num_overlap(e) * pows[e] for e in range(2, max_e + 1) ]

    return (N - 1) ** 2 - sum(a)
N = 10 ** 7

ary = pow_numbers()

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

pows = convert(ary)

print solve()