Project Euler 241 - inamori’s diary

Q241.
nの全ての約数の和をσ(n)、p(n) = σ(n) / nとする。
p(n) = k + 1 / 2(kは正の整数)となるn ≤ 10¹⁸の総和を求めよ。

n = p₁^e₁...p_m^e_mなら、σ(n) = (p₁^e₁+1-1)/p₁^e₁(p-1)...(p_m^e_m+1-1)/p_m^e_m(p_m-1)だから、あとは再帰的に素因数を選んでいく。素因数は小さいほうから選ぶ。2がなければならない。意外に大きな素数は使えなくて、例えばnが1000以下で、今の素数が11で、13/10を作ろうとしても、2つしか使えないから、13/9 > (14/13)²となるからダメ。それから、分母に3の素因数があるから、もう3を超えているのでダメ。これでも10¹⁶まではそんなに時間がかからなかったが、10¹⁸にしたらなかなか進まなくなってしまった。それでもほっておいたら、正しい答えが出た。



from itertools import count

from math import sqrt

from fractions import Fraction

import time
def is_prime(n):

    for p in primes:

        if p * p > n:

            return True

        elif n % p == 0:

            return False

    return True
def make_primes(n):

    m = int(sqrt(n + 0.5))

    for p in xrange(3, m + 1, 2):

        if is_prime(p):

            primes.append(p)

    

    a = [ 0 ] * ((n + 15) / 16)

    for p in primes:

        if p * p > n:

            break

        if p == 2:

            continue

        for k in xrange(3 * p, n, 2 * p):

            a[k>>4] |= 1 << (k & 15)

    

    for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):

        if (a[k>>4] & (1 << (k & 15))) == 0:

            primes.append(k)
def calc_quotient(p, e):

    q = p ** e

    return Fraction( (p * q - 1) / (p - 1), q)
def gen_quotient(n0, k0 = 0):

    if n0 > 1:

        for k in xrange(k0, len(primes)):

            p = primes[k]

            n = n0 / p

            if n == 0:

                break

            for e in count(1):

                f = calc_quotient(p, e)

                yield (p ** e, f)

                for m, q in gen_quotient(n, k + 1):

                    yield (p ** e * m, f * q)

                n /= p

                if n == 0:

                    break
t = time.clock()

N = 10 ** 7

primes = [ 2 ]

make_primes(N)

for n, q in gen_quotient(N):

    if q.denominator == 2:

        print n, q

print time.clock() - t