ユークリッドの果樹園 とオイラーのファイ関数の和の高速化（その３）

オイラーのファイ関数の和の高速化

前回の投稿で以下の漸化式が出ましたので、これをプログラムにして行きたいと思います。

$\displaystyle \sum_{m=1}^{N} φ(N) = G(N,1)+1 $
$\displaystyle = \frac{N(N-1)}{2}-\sum_{m=2}^{N} G(\lfloor \frac{N}{m} \rfloor ,1)+1$

まず再帰関数にしやすいように＋１のところは後にして、以下の式を考えます

$\displaystyle G(N,1)= \frac{N(N-1)}{2}-\sum_{m=2}^{N} G(\lfloor \frac{N}{m} \rfloor ,1)$

これをプログラムにしたのが下のV3.0です、それなりに速くなっていると思いますが、まだ1分以内で計算できるのは$10^7$までです。さらなる高速化を考えます。

# V3.0 G(N,1) recursive with lru_cache
from functools import lru_cache
@lru_cache(maxsize=None)
def G(n):
  if n == 0: return 0
  ret = 0
  for m in range(2,n+1):
    ret += G(n//m)
  return n*(n-1)//2-ret 

n = 10**8
print(n, G(n))
# 100000000 3039635516365907

	$10^4$	$10^5$	$10^6$	$10^7$	$10^8$
V1.0(gcd)	8s	x	x	x	x
V2.0(totient)	0s	0s	70s	x	x
V3.0(G(N,1))	0s	0s	1s	16s	180s

さらなる高速化

これは前述のリンクHow to calculate these totient summation sums efficiently?で紹介されているSimplified hyperbola method 簡易双曲線法を使ったものです。

難しそうなので、なるべく分かりやすく説明します。

再帰的に呼ぶ引数の$N//m$に注目すると、その値は同じ値を取る場合が多いことに気づきます。例えば$N=10,m=2$の場合だと以下のようになります。すなわち$m\leqq x < m'$で同じ値になっています。従ってこの区間は再帰呼び出しを1回だけ行い$(m'-m)$倍してして、その後$m$を$m'$にスキップすればよいことが分かります。

N	m	N//m	m'
10	2	5	3
10	3	3	4
10	4	2	6
10	5	2	6
10	6	1	11
10	7	1	11
10	8	1	11
10	9	1	11
10	10	1	11

このように変更したコードが下のV3.1です。さらに+1するコードをかぶせて$ \sum_{m=1}^{N}φ(m)$を求める最終的な関数totsumを作りました。これを使って$10^{10}$まで約1分で計算できるようになりました。

# V3.1 Simplified hyperbola method #1
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**8
print(n, totsum(n))

	$10^4$	$10^5$	$10^6$	$10^7$	$10^8$	$10^9$	$10^{10}$
V1.0(gcd)	8s	x	x	x	x	x	x
V2.0(totient)	0s	0s	70s	x	x	x	x
V3.0(G(N,1))	0s	0s	1s	16s	180s	x	x
V3.1(totsum)	0s	0s	0s	0s	1s	3s	26s

（追記）簡易双曲線法の別の実装

上記のリンクでは以下の漸化式が出てくるので、これをそのまま実装したコードも添付します。コードはかなり違いますがやっていることは同等なのでパフォーマンスはV3.1とほぼ同じでした。

\begin{align}
\Phi(n) &= T(n) - \sum_{2\le i \le n/c}\Phi(\frac{n}{i}) - \sum_{j \le c-1}(\lfloor \frac{n}{j} \rfloor - \lfloor \frac{n}{j+1} \rfloor) \Phi(j) \\
& ここで　c = \lfloor \sqrt{n} \rfloor
\end{align}

# V3.2 Simplified hyperbola method #2
from functools import lru_cache
def totsum2(N):
  @lru_cache(maxsize=None) 
  def G(N):
    c = int(N**(0.5))
    return (N*(N-1))//2 \
          - sum(G(N//i) for i in range(2,c+1)) \
          - sum(((N//j) - (N//(j+1)))*G(j) for j in range(1,N//(c+1)+1))
  return 1 + G(N)

n = 10**8
print(n, totsum2(n))

（開発環境：Google Colab）

このアルゴリズムはEuler ProjectのProblem 351, 202, 787等を解くのに役に立ちます。