Project Euler 1（2） - inamori’s diary

元の問題を一般的に、

nと互いに素でないm未満の自然数で総和を求めよ。

としましょう。例えば、nを60とすると、

60 = 2² • 3 • 5

なので、2または3または5の倍数の総和を求めればよいことになります。これを求めるには、まず2の倍数の総和と3の倍数の総和と5の倍数の総和を求めます。ここからダブっている6の倍数の総和と10の倍数の総和と15の倍数の総和を引きます。さらに、30の倍数を引きすぎているので、足します。一般に、包除原理より、奇数個の素数からなる倍数の場合は総和を足して、偶数個なら引けばよいことになります。

これをコーディングするには、どうすればいいでしょう。15で考えると、

3 5 -15

と出すことを考えます。3なら3の倍数の総和を足して、-15なら15の倍数の総和を引くようにします。
では、こう出すにはどうすればいいでしょう。一旦、符合は無視して、

3 5 15

と出すことを考えましょう。

3¹•5⁰ 3⁰•5¹ 3¹•5¹

と書けます。1をつけ加えると、

3⁰•5⁰ -> 00
3¹•5⁰ -> 01
3⁰•5¹ -> 10
3¹•5⁰ -> 11

0から3まで2進数に対応させることができます。コードにすると、



def value(code, a):

    n = 1

    for p in a:

        if code == 0:

            break

        if code & 1:

            n *= p

        code >>= 1

    return n
a = [ 3, 5 ]

for c in range(1, 1 << len(a)):

    print value(c, a),

再帰を使ったほうがきれいに書けます。



def value(code, a, k = 0):

    if code == 0:

        return 1

    

    n = 1

    if code & 1:

        n = a[k]

    code >>= 1

    return n * value(code, a, k + 1)

これに符号をつければ期待した結果が得られます。



def value(code, a, k = 0):

    if code == 0:

        return -1

    

    n = 1

    if code & 1:

        n = -a[k]

    code >>= 1

    return n * value(code, a, k + 1)

Pythonではこんなきれいな書き方もできます。



def gen_value(a, k = 0):

    if k == len(a):

        yield -1

    else:

        for n in gen_value(a, k + 1):

            yield n

            yield -a[k] * n
a = [ 3, 5, 7 ]

for n in filter(lambda n: n != -1, gen_value(a)):

    print n,

以下の全体のコードでは、別のもう少し分かりにくい書き方をしています。
これは次のように使います。

〜.py 15 1000



import sys
def sum_sequence(n, d):

    m = (n - 1) / d

    return d * m * (m + 1) / 2
def sum_seq(n, d):

    if d > 0:

        return sum_sequence(n, d)

    else:

        return -sum_sequence(n, -d)
def gen_divisors(ary_p, k0 = 0):

    if k0:

        yield -1

    

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

        for d in gen_divisors(ary_p, k + 1):

            yield -ary_p[k] * d
def get_exp(n, p):

    e = 0

    while n % p == 0:

        e += 1

        n /= p

    return e, n
def gen_prime_candidates(n):

    yield 2

    for p in range(3, n, 2):

        yield p
def factorize(n):

    f = [ ]

    for p in gen_prime_candidates(n):

        if p * p > n:

            break

        e, n = get_exp(n, p)

        if e:

            f.append(p)

    if n != 1:

        f.append(n)

    return f
def solve(n, m):

    f = factorize(n)

    return sum(map(lambda d: sum_seq(m, d), gen_divisors(f)))
if len(sys.argv) == 3:

    print solve(int(sys.argv[1]), int(sys.argv[2]))

else:

    print "usage : %s n upper_limit." % sys.argv[0]