ユークリッドの果樹園

ユークリッドの果樹園(wiki)について色々調べてみました。

座標の格子点$(x,y)\space (x, yは自然数)$に細い木を立てたとき、原点から見える木の条件は$gcd(x,y)=1$である

というものでWikiにもある画像をN=25まで描いてみました。

ここからオイラーのファイ関数につながって行きます。まず次のような問題を考えます。

【問題】ユークリッドの果樹園で$0<x,y<N$のとき原点から見える木の数を求めよ。

すなわち上の図で青の点の数を数えよという問題になります。
まず素直にN=25で$gcd(x,y)=1$となる点$(x,y)$を数えてみます。

# V1.0
import math
N, count = 25. 0
for i in range(1,N+1):
  for j in range(1,N+1):
    if math.gcd(i, j) == 1:
      count += 1

print(f"answer = {count}")
# answer = 339

計算量は$O(n^2)$となるのでかなり効率が悪いですね。そこでオイラーのファイ関数$φ(n)$の登場となります。Wikiによれば、

$φ(n)$は n と互いに素である1以上n以下の自然数の個数

となります。これは下の図のようにユークリッドの果樹園の右下の三角形の縦の列の青の数と同じになります。したがってこれを2倍して1引けば答が出るはずです。（$(x,y)=(1,1)$が重複）

今回は組込みのsympy.ntheoryのtotient関数を使って計算すると同じ答えになりました。

# V2.0 use sympy totient function
import sympy as sp
N = 25  
print(f"answer = {sum(sp.totient(x) for x in range(1, N+1))*2 - 1}")
# answer = 339

計算時間は以下のようになります。V1.0(gcd)は$10^4$まで、V2.0(totient)では$10^6$まで100秒以内で計算できました。それ以上になるとさらなる高速化が必要になります。続編でオイラーのファイ関数の和の高速化について書きたいと思います。

	$10^3$	$10^4$	$10^5$	$10^6$	$10^7$
V1.0(gcd)	1s	38s	x	x	x
V2.0(totient)	1s	1s	5s	70s	x

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

（開発環境：Google Colab）