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

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

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

$\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

次に$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}φ(N)$を求める最終的な関数totsumをインプリしました。この結果$10^{10}$まで約1分で計算できるようになりました。

# 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**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	10s	64s

このアルゴリズムはEuler ProjectのProblem 351, 202等に役にたつと思います。

（開発環境：Google Colab）