【Project Euler】Problem 72: 分数の数

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題 72. 分数の数

原文 Problem 72: Counting fractions

問題の要約：分母$d \le 10^6$の既約真分数の個数を求めよ

例として$d=8$では以下のように21個あります。

問題の見た目は問題 71.と似ていますが内容は全く異なります。以下のように分母$d$の既約真分数の分子の数は「dと互いに素である1以上d以下の自然数の個数 」となり「オイラーのφ関数(Wikipedia)」の定義そのものとなります。

分母	分子	分子の数
2	1	1
3	1,2	2
4	1,3	2
5	1,2,3,4	4
6	1,5	2
7	1,2,3,4,5,6	6
8	1,3,5,7	4

従って答えはsympy.ntheoryのtotient関数を使えば簡単に出ますが、制限時間1分を超えてしました。

from sympy.ntheory import totient
print(f"Answer: {sum([totient(den) for den in range(2,10**6+1)])}")

このtotient関数を自分でfactorintを使ってインプリして何とか1分を切れました。

from sympy import factorint
def phi(n):
  ret = n
  for i in factorint(n).keys():
    ret = ret*(i-1)//i
  return ret

print(f"Answer: {sum([phi(den) for den in range(2,10**6+1)])}")

もしこれ以上のスピードアップを望むなら、ちょっと長いですが「ユークリッドの果樹園 とオイラーのファイ関数の和の高速化（その１）」という記事（その3まであります）を読んでいただけたらと思います。ちなみに今回の答えも瞬時に出てきます。

# totsum V3.1 totsum
from functools import lru_cache

def totsum(n):
  @lru_cache(maxsize=None)
  def G(n):
    if n == 0: return 0
    ret, m1 = 0, 2
    while m1 < n:
      n1 = n//m1
      m2 = n//n1 + 1
      ret += (m2-m1)*G(n1)
      m1 = m2
    return (n*(n-1))//2-ret
  return G(n)+1

n = 10**6
print(f"Answer: {totsum(n)-1}")

（開発環境：Google Colab）