Project Euler 39

PE 数学 プログラミング

Problem 39
pを直角三角形の周囲の長さとすると、p = 120の解は3つある。
（中略）p ≤ 1000に対して、解の個数が最大のものは？
http://projecteuler.net/index.php?section=problems&id=39

ピタゴラス数は、

a = 2lmn b = l(m2 - n2) c = l(m2 + n2)
m > n
mとnは互いに素
m + nは奇数

と書けることを知っているのが前提の問題です。周囲の長さpは、

p = 2lm(m + n)
mとm + nは互いに素
m + nは奇数
m < m + n < 2m

となります。
例えば、p = 120で考えます。m(m + n) ≥ 6だから、l ≤ 10。lは60の約数だから、l = 1, 2, 3, 4, 5, 6, 10。l = 1とすると、m(m + n) = 60 = 22 • 3 • 5。これをmと(m + n)で分け合うことになります。22はmが取るので、残りを分け合います。この例は適切ではないですが、互いに素なのでどちらかが一つ素数を取るとその素数は全て取ることになります。15の分け合い方は、(1, 15), (3, 5), (5, 3), (15, 1)となり、22を掛けると、(4, 15), (12, 5), (20, 3), (60, 1)となりますが、条件に合いません。このように計算していけば答えは出ますが、遅いです（下のCode1）。


この問題は、やはり逆方向から考えましょう。lとmとm + nを回してpを計算し、0に初期化したリストに1を足していきます。こうすれば速くなります（Code2）。問題が大きくなるとメモリが足りなくなりますが、そういうときはうまく分割しましょう。




# Code1

from itertools import ifilter, imap
def make_primes(max_p):

    a = [ True ] * (max_p + 1)

    for p in ifilter(lambda n: a[n], xrange(2, int(max_p ** 0.5 + 0.5))):

        for n in xrange(p * 2, max_p + 1, p):

            a[n] = False

    return list(ifilter(lambda n: a[n], xrange(2, max_p + 1)))
def div_pow(n, p):

    while n % p == 0:

        n /= p

    return n
def factorize(max_n):

    a = range(max_n + 1)

    b = [ () for n in a ]

    max_p = int(max_n ** 0.5 + 0.5)

    for p in ifilter(lambda n: a[n] == n, xrange(2, max_p + 1)):

        for n in xrange(p, max_n + 1, p):

            a[n] = div_pow(a[n], p)

            b[n] = b[n] + (p,)

    

    for n in xrange(2, max_n + 1):

        if a[n] != 1:

            b[n] = b[n] + (a[n],)

    

    return b
def calc_exp(n, p):

    e = 0

    while n % p == 0:

        e += 1

        n /= p

    return e
def gen_divisors(f, k0 = 0):

    yield ()

    if k0 < len(f):

        for k in range(k0, len(f)):

            p, e_max = f[k]

            for e in range(1, e_max + 1):

                f1 = ( (p, e),)

                for f2 in gen_divisors(f, k + 1):

                    yield f1 + f2
def index(f, p):

    for k in range(len(f)):

        if f[k][0] == p:

            return k

    return -1
def sub(f1, f2):

    f = ()

    for p1, e1 in f1:

        k = index(f2, p1)

        if k != -1:

            e2 = f2[k][1]

            if e2 < e1:

                f = f + ( (p1, e1 - e2),)

        else:

            f = f + ( (p1, e1),)

    return f
def value(f):

    return reduce(lambda x, y: x * y[0] ** y[1], f, 1)
def gen_mn_core(f, k = 0):  # (m, n) = 1

    if k == len(f):

        yield (), ()

    else:

        for f1, f2 in gen_mn_core(f, k + 1):

            yield (f[k],) + f1, f2

            yield f1, (f[k],) + f2
def gen_mn(f):

    if len(f) > 0 and f[0][0] == 2:

        f1 = (f[0],)

        f2 = f[1:]

    else:

        f1 = ()

        f2 = f

    

    for mf, mnf in gen_mn_core(f2):

        m = value(f1 + mf)

        mn = value(mnf)     # m + n

        n = mn - m

        if n > 0 and n < m:

            yield m, n
def num_right_triangles(n):

    nn = n

    counter = 0

    n /= 2

    f = reduce(lambda x, y: x + ( (y, calc_exp(n, y)), ), fac[n], ())

    for lf in gen_divisors(f):

        l = value(lf)

        for m, n in gen_mn(sub(f, lf)):

            counter += 1

    return counter
N = 10 ** 5

M = N / 2

fac = factorize(M)

g = imap(lambda n: (n, num_right_triangles(n)), xrange(2, N + 1, 2))

print reduce(lambda x, y: x if x[1] >= y[1] else y, g, (0, 0))[0]



# Code2

from fractions import gcd
N = 10 ** 6

M = N / 2

a = [ 0 ] * (M + 1)

for l in xrange(1, M / 6 + 1):

    for m in xrange(2, int( (M / l) ** 0.5) + 1):

        begin = (m + 1) / 2 * 2 + 1         # odd

        end = min(m * 2, M / (l * m) + 1)

        for mn in xrange(begin, end, 2):

            if gcd(m, mn) == 1:

                k = l * m * mn

                a[k] += 1

print reduce(lambda x, y: x if a[x] >= a[y] else y, xrange(M + 1)) * 2