ファレイ数列の応用（ユークリッドの果樹園、フォードの円）

### ファレイ数列とは

[ファレイ数列の４つの性質とその証明(高校数学の美しい物語)](https://manabitimes.jp/math/1167)によると、

> 0 以上 1 以下であり，分母が$n$以下であるような**既約分数**を小さい順に並べた数列を（$n$に対応する）ファレイ数列という。

まずファレイ数列を生成するプログラムを作ります。これには上記のリンクの性質４を使います。

> $F_n$の隣り合う二項で分母の和が$n+1$のもの $\Large \frac{q_i}{p_i}, \frac{q_{i+1}}{p_{i+1}}$について、この間に$\Large \frac{q_i+q_{i+1}}{p_i+p_{i+1}}$を入れる（該当箇所全部にこの操作をする）とが$F_{n+1}$できる

文章にするとやや分かりずらいですが、このアルゴリズムは[Stern-Brocot Tree (Wikipedia)](https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree)という名前がついていてサンプルプログラムも多く存在します。

このファレイ数列を生成する関数**sternBrocot**を作りました。f1, f2 で両端の分数の値、dmaxで分母の最大値を指定します。

```python 
import numpy as np
from collections import deque
def sternBrocot(f1, f2, dmax):
  f1, f2 = np.array(f1), np.array(f2)  # convert to np.array to make addtion easier
  q = deque([(1, f1, f2)])  # queue (pos, f1, f2), pos: position to insert next level
  ret = [f1,f2]             # intitial list: both edge values
  while len(q)>0:
    pos, f1,f2 = q.pop()
    fm = f1+f2              # next leve: add nom/denom respectively
    if fm[1] <= dmax:       # check denominator value
      ret.insert(pos,fm)    # insert next level in the middle 
      q.append((pos, f1, fm)) 
      q.append((pos+1, fm, f2))
  return list(map(list,ret))

def fform(ll):
  return ", ".join(list(map(lambda l2: str(l2[0])+"/"+str(l2[1]),ll)))

for n in range(2,5+1):
  print(f"n = {n}: {fform(sternBrocot((0,1),(1,1),n))}")
```
| N| ファレイ数列 |
|:---:|:---|
| 2 | 0/1, 1/2, 1/1 |
| 3 | 0/1, 1/3, 1/2, 2/3, 1/1 |
| 4 | 0/1, 1/4, 1/3, 1/2, 2/3, 3/4, 1/1 |
| 5 | 0/1, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 1/1 |
| 6 | 0/1, 1/6, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 5/6, 1/1 |
| 7 | 0/1, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 2/5, 3/7, 1/2, 4/7, 3/5, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 1/1 |

これだけ見るとそれが何の役に立つのかと思うかもしれませんが、意外なところに役に立つという例を示して行きたいと思います。

重要な点は　**「約分をしていないのに既に既約分数になっている」**　ということです。したがって見方を変えると以下のことが言えます。

* 分母がN以下の正の有理数をすべて列挙している
* N以下の「互いに素(coprime)」な自然数のペアをすべて列挙している

### 範囲を指定して既約分数を列挙する

[Project Euler, Problem 73: Counting fractions in a range](https://projecteuler.net/problem=73)がこの例となります。

ここで例題として挙げられている$1/3$と$1/2$の間の分母が最大$8$のものを上記の **sternBrocot**を使って（両端を含めて）リストしてみます。

```python 
print(f"{fform(sternBrocot((1,3),(1,2),8))} ")
# 1/3, 3/8, 2/5, 3/7, 1/2 
```

元の問題の解法は[「Problem 73: 範囲内の分数の数」](https://qiita.com/masa0599/items/4f72f86a66814f704c7c)を参考にしてください。

### ユークリッドの果樹園

詳細は以下の記事を参考にしてください。

> [ユークリッドの果樹園 とオイラーのファイ関数の和の高速化（その１）](https://qiita.com/masa0599/items/6968c99e0e928763cb42)

ここでユークリッドの果樹園の条件をもう一度書くと以下のようになり、その座標はまさに「互いに素(coprime)」な自然数のペア」ということになります。

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

従って5以下の原点から見える木の座標（x,y）は$dmax=5$として求められます。分子分母を逆にすると、下の図の青の点（原点から見える木）を表しているのが分かります。

```python 
ut = sorted(list(map(lambda x: x[::-1], sternBrocot((0,1),(1,1),5))))
print(f"{ut[1:]} ")
# [[1, 1], [2, 1], [3, 1], [3, 2], [4, 1], [4, 3], [5, 1], [5, 2], [5, 3], [5, 4]] 
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1922159/ef528f01-84ec-47c5-ea3c-572f0d8e06af.png" width=40%>

ここまで来るとファレイ数列の分母が$n$の分数の数はオイラーのファイ関数$φ(n)$で表されるということも推測できると思います。

### フォードの円

[フォードの円(Wikipedia)](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A9%E3%83%BC%E3%83%89%E3%81%AE%E5%86%86) は以下のように説明されます。

> ファレイ数列で得られた既約分数$p/q$に対して中心が$(p/q, 1/2q^2)$半径が$1/2q^2$の円を描いたもの

具体的にN=30の時のファレイ数列の各分数に対してこのルールで円を描くと以下にように**隣り合う円が互いに外接している**のが分かります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1922159/ac8daed3-ac02-9970-9c73-1c47320f71ff.png)

### フォードの円の応用

フォードの円の原理から以下のことが言えます。

* 互いに素な正の整数$a,b$が与えられた時$c=a+b$ と置くと、各々$1/a^2$、$1/b^2$、$1/c^2$ の半径の円を直線の上に互いに外接するように置くことが出来る。
* さらに半径をすべて$(abc)^2$倍することより、各々$(bc)^2$、$(ca)^2$、$(ab)^2$ の半径の円になるので、すべて半径が整数のフォードの円を書くことが可能になります。

例として$a,b,c = 11,14,25$のときに互いに外接する円を書いてみました。3つの円の半径を整数で表すことができます。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1922159/b8f5d488-9986-96d2-7191-ebf3a2598a6c.png)


（開発環境：Google Colab）

この考え方はProject Euler、Problem 510: Tangent Circlesを解くのに役に立ちます。