Project Euler 73

PE 数学 プログラミング

Problem 73
n/dという分数を考える。ここで、nとdは正の整数。n<dでHCF(n,d)=1のとき、既約真分数と呼ぶ。
d≤8の既約真分数の集合を昇順に並べると次のようになる：
1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
1/3と1/2の間に3つの分数があるのが分かる。
d ≤ 12000に対し、1/3と1/2の間に既約真分数がいくつあるか。
http://projecteuler.net/index.php?section=problems&id=73

dに対し、互いに素なnがいくつあるか数えていけば答えが出ます。

from fractions import gcd

N = 10 ** 3 * 12
counter = 0
for d in xrange(1, N + 1):
    for n in xrange(d / 3 + 1, (d + 1) / 2):
        if gcd(n, d) == 1:
            counter += 1
print counter

速くなるように工夫しましょう。
0から1ならdを分母とする既約分数の個数はφ(d)ですが、1/3から1/2の間にある既約分数の個数はφ(d)/6にならないから難しいのです。しかし、dを素因数分解したときに、因子のいずれかのφが6の倍数なら、それが成り立ちます。例えば104 = 23 * 13を考えると、φ(23) = 4,φ(13) = 12だから、後者が6で割り切れるので、φ(104) / 6 = 4 * 12 / 6 = 8個の既約分数が1/3と1/2の間にあることになります。

35/104 37/104 41/104 43/104 45/104 47/104 49/104 51/104

これだけでもかなり速くなります。
実は、dを分母とする既約分数は、dの因子のφの約数分はきれいに分かれます。例えば、d = 22 = 2 * 11とすると、φ(11) = 10の約数5に対し、1/5ずつできれいに分かれます。

1/22 3/22 | 5/22 7/22 | 9/22 13/22 | 15/22 17/22 | 19/22 21/22

|で区切って、最初の領域は1/5まで、次の領域は2/5まで、などとなっています。
別の例でd = 100を考えましょう。素因数分解すると100 = 22 * 52、それぞれの因子について、φ(22) = 2、φ(52) = 20なので、大きいほうの20を取ります。1/3より小さくて最も大きい20を分母とする分数を考えます。それは6/20です。(6/20,1/2) = (30/100,50/100)の区間に分母が100の既約分数はφ(100)*(50-30)/100 = 8個あることになります。余分に数えた分数は一つずつ調べて31/100と33/100の2つで、これを差し引いて6個となります。
これで最初のコードに比べて2桁近く速くなりました。

from itertools import imap, ifilter
from fractions import gcd

def count_if(f, a):
    counter = 0
    for e in ifilter(f, a):
        counter += 1
    return counter

def factorize(max_n):
    def sieve(max_n):
        a = range(max_n + 1)
        for p in ifilter(lambda n: a[n] == n, xrange(2, max_n / 2 + 1)):
            for k in ifilter(lambda k: a[k] == k,
                            xrange(p * 2, max_n + 1, p)):
                a[k] = p
        return a
    
    def div_pow(n, p):
        e = 0
        while n % p == 0:
            e += 1
            n /= p
        return n, e
    
    a = sieve(max_n)
    a[1] = []
    for n in xrange(2, len(a)):
        if a[n] == n:   # prime
            a[n] = [(n, 1)]
        else:
            p = a[n]
            m, e = div_pow(n, p)
            a[n] = [(p, e)] + a[m]
    return a

def phi_f(f):
    p, e = f
    return (p - 1) * p ** (e - 1)

def phi(f):
    return reduce(lambda x, y: x * phi_f(y), f, 1)

def is_div(f, d):
    return any(imap(lambda g: g[0] % d == 1, f))

def count_reduced(d):
    f = facts[d]
    if is_div(f, 6):
        return phi(f) / 6
    else:
        max_f = reduce(max, imap(phi_f, f))
        n_min = max_f / 3
        s = phi(f) * (max_f / 2 - n_min) / max_f
        return s - count_if(lambda n: gcd(n, d) == 1,
                        xrange(n_min * d / max_f + 1, d / 3 + 1))

N = 10 ** 3 * 12
facts = factorize(N)
print sum(imap(count_reduced, xrange(2, N + 1)))