Project Euler 70 #Python

代数的に解きたかった

Project Eulerの問70について、$Φ$関数の実装を自分なりにしてみた所、ある意味中途半端で面白いものになったので載せてみます。

<$Φ$関数>

import sympy, itertools, collections

def prime_factorization(value):
    factors = set()
    for i in range(2, value + 1):
        while value % i == 0:
            factors.add(i)
            value //= i
        if value == 1:
            break
    return factors

def calculate_baisu(n):
    prime_factors = list(prime_factorization(n))
    baisu = 0
    for r in range(1, len(prime_factors) + 1):
        combinations = itertools.combinations(prime_factors, r)
        for combination in combinations:
            product = 1
            for divisor in combination:
                product *= divisor
            if size % 2 != 0:
                baisu += int(n / product)
            else:
                baisu -= int(n / product)
    return baisu

def calculate_phi(n):
    phi_value = n - calculate_baisu(n)
    return phi_value

if __name__ == "__main__":
    n = 87109 # 問題文中の例
    result = calculate_phi(n)
    print(f"phi({n}) = {result}")

# phi(n) = 79180

基本的なアイデアとしては
・自然数n以下のnと互いに素な自然数の個数Φ(n)を求める
ために、その背反である
・n以下の自然数でnと共通する約数を持つものの数を求め
これは
・nを割り切る素数の倍数がそれぞれn以下に何個あるか
からわかり、この計算は
・倍数の重複を考慮し、組み合わせ$nCr$を使えば求まる

ということなのですが、それを上手く計算するために

nを素因数分解しその素因数達の倍数を数えていく
nの約数の倍数達には重複があるので、これらの数について合併集合を作る。

例えば３つの集合A, B, Cの合併集合は
$A\cup B\cup C:=$ {${x\mid x\in A{\mbox{ or }}x\in B{\mbox{ or }} x\in C}$}
で定義され、その個数について$∣A∪B∪C∣=∣A∣+∣B∣+∣C∣−∣A∩B∣−∣A∩C∣−∣B∩C∣+∣A∩B∩C∣$
が成り立つ。

任意数個の集合についてもその個数は包除原理という以下の積和型の式で表されます。(*1)

イメージとしては重複している積集合の個数を引いて、そうすると引きすぎてしまうので足してあげて、そうすると今度は足しすぎてしまうのでまた引いてあげて...という感じです

という事を考えました。1. は関数prime_factorization(n) で定義されているように、nを2から割れるだけ割り切っていきながらその分だけ因数としてリストに入れていくことで実現しています。 2. ではnの素因数達についての組み合わせを得るため、引数に素因数のリスㇳをとるitertools.combinations(divisors, r)を使っています

具体例

n=60 の場合について見てみます。prime_factorization(60)では2で2回、3と5で一回ずつ割られた60が{2, 3, 5}という集合になって返ってきます。calculate_baisu(60)ではこの素因数の集合を使い、それらの倍数の個数を求めます。int(60/2) = 30, int(60/3) = 20, int(60/5) = 12なので$1, 2, 3, 4,..., 60$ の中には2, 3, 5 の倍数がそれぞれ上の数ずつ存在します。しかしその中には$6 =2×3$や$30 =2×3×5 $のように、複数の素因数で割り切れる数が含まれます。その数を考慮して、合併集合
$S=$ {${x\mid x<60, \ x\in A{\mbox{ or }}x\in B{\mbox{ or }} x\in C} \ $}

$where \ A= ${$x\mid x=0 \mod2$}$, \ B= ${$x\mid x=0 \mod3$}$,
\ C= ${$x\mid x=0 \mod5$}
の要素数$|S|$を求めると、これは次の組み合わせの集合の要素数の加減算で求めることができるのがわかります。
2, 3, 5 で割り切れる → $ \ _3C_1 $
6, 10, 15 で割り切れる → $ \ _3C_2 $
30 で割り切れる → $ \ _3C_3 $
これらの組み合わせを素因数のリストprime_factorsから計算しているのがcalculate_baisu(n)内のtertools.combinations(prime_factors, r)です。

時間計算量

このアルゴリズムでは数$n$の素因数分解で$O(\sqrt{n})$ かかります。これは$n$の素因数は最悪でも$\sqrt{n}$ 以下を走査すれば全て求まるからです。これが求まってしまえば、素因数の数等によっても左右されるもののあとの倍数の計算には大きな計算量はかからないです。しかしこの見積もりはとても大雑把で、場合によっては大きなオーダーの計算量になるかもしれません。

一方、n以下の数mについてmを割り切る数がnを割り切るか見ることで(n, m)の組が互いに素かを判定するといったとてもnaiveな実装では時間計算量が$O(n^2)$となり非現実的です。

今回の問題はここからさらに$\phi(n)$の permunation-性の判定と$\phi(n)/n$の計算などがいるので、もっと高速なアルゴリズムの発案が求められました。

もう少し

実は今回オイラーが求めていたであろう正解はもう少し踏み込んだところにありました。ここでは概要だけを説明します。基本的な考えは同じで、素数$p$で割り切れる自然数$m$以下の数を数えます。そのような数は $p, 2p, 3p,..., (m/p)p$ なので $p$で割り切れない$m$以下の数の個数は$m-m/p$とかけます。ここで、別の素数$q$でも割ることを考えると、この中で更に$q$でも割り切れないものの数は$m(1-\frac1p)-\bigl(\frac mq - \frac m{pq} \big) = m(1-\frac1q)(1-\frac1p)$と書けることが、簡単な議論からわかります。これを一般化すると、今回の問題は高々$n$の素因数個回だけの乗算で計算できるようになります。

詳しくはWolframの下記記事を御覧ください。
Toient Function

とても面白い記事になっています。自分の回答では直感的に進め過ぎたことで煩雑な計算を経由しなくてはいけませんでした。もう少し落ち着いて何が必要かを考えられるようになりたいです。

Wikipedia 含除原理