Beatty数列とその和を求めるアルゴリズム

### Beatty数列とは

あまりなじみのない名前ですが[ビーティ数列 (Wikipedia)](https://ja.wikipedia.org/wiki/%E3%83%93%E3%83%BC%E3%83%86%E3%82%A3%E6%95%B0%E5%88%97)とは、

>  1より大きい無理数の整数倍の床関数をとることによって得られる整数列である。ビーティ列の名称は、1926年にそれらについて著したサミュエル・ビーティに因む。

と言われてもピンとこないので具体的に無理数$r$を黄金比$\phi$として例を示します。

```math
\begin{align}
& r = \phi = \frac{1+\sqrt{5}}{2} =  1.6180339... として\\
& \lfloor nr \rfloorをn=1,2,3,...,10で求めると以下のようになります \\
\end{align}
```

n |1 |2 |3 |4 |5 |6 |7 |8 |9 |10 |
:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
nr |1.62|3.24|4.85|6.47|8.09|9.71|11.33|12.94|14.56|16.18|
$\lfloor nr \rfloor$ |1|3|4|6|8|9|11|12|14|16|

### Beatty数列の面白い性質

これだけだとどこが面白いのか分かりませんが、

```math
\begin{align}
& \frac{1}{r} + \frac{1}{s} = 1となるs、すなわち s = \frac{r}{r-1}として \\
& s = \frac{\phi}{\phi-1} =   \frac{\phi^2-1}{\phi-1} = \phi+1 = 2.6180339...\\ 
& 同様に\lfloor sr \rfloorをn=1,2,3,...,10で求めると\\
\end{align}
```
n |1 |2 |3 |4 |5 |6 |7 |8 |9 |10 |
:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
$\lfloor ns \rfloor$ |2|5|7|10|13|15|18|20|23|26|

これでも分かりづらいと思うので$\lfloor nr \rfloor$と$\lfloor ns \rfloor$を並べてリストして見ます。そして数字の青と赤の部分を見てみると。**１から１６までのすべての整数がすべてしかも重複せずに**現れているのが分かります。


n |1 |2 |3 |4 |5 |6 |7 |8 |9 |10 |
:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
$\lfloor nr \rfloor$ |$\color{blue}{1}$|$\color{blue}{3}$|$\color{blue}{4}$|$\color{blue}{6}$|$\color{blue}{8}$|$\color{blue}{9}$|$\color{blue}{11}$|$\color{blue}{12}$|$\color{blue}{14}$|$\color{blue}{16}$|
$\lfloor ns \rfloor$ |$\color{red}{2}$|$\color{red}{5}$|$\color{red}{7}$|$\color{red}{10}$|$\color{red}{13}$|$\color{red}{15}$|18 |20 |23 |26 |

これがBeatty数列の面白い性質であり、これを[「レイリーの定理(Wikipedia)」](https://ja.wikipedia.org/wiki/%E3%83%AC%E3%82%A4%E3%83%AA%E3%83%BC%E3%81%AE%E5%AE%9A%E7%90%86)と呼ぶそうです。

> 数学におけるレイリーの定理とは、1より大きい無理数が、床関数によって自然数全体を互いに素な2つの集合に分ける方法を与える定理である

### 黄金比のBeatty数列の和をレイリーの定理を使って高速に求める

この性質を利用した問題を考えます。

**【例題】黄金比$\phi$に対して$S(\phi, N) = \sum_{k=1}^{N}\lfloor n\phi \rfloor$とするとき$S(\phi, 10^{16})$を求めよ**

まず$N=10$として$S(\phi, 10)$を考えます。これは上の表の<font color="blue">青の数字</font>の和です。<font color="blue">青の数字</font>と<font color="red">赤の数字</font>を合わせると１から１６までのすべての整数があるのでその和から<font color="red">赤の数字</font>を引けば求まります。

```math
\begin{align}
& 赤の数字の最大値m = n \times \frac{r}{s} = n \times \frac{\phi}{\phi+1}=6\\
& 1からnまでの和をT(n)とすると、\\
& S(\phi, N) = \sum_{k=1}^{N}\lfloor n\phi \rfloor = T(N\phi)-\color{red}{S(\phi+1,m)}...(1) \\
\end{align}
```

<font color="red">赤の数字$\lfloor ns \rfloor$</font>ですが同じ列の$\lfloor nr \rfloor$と比べるとその差が$n$になっています。これは$s=r+1$から$ns=nr+n$となるので当然ですね。

```math
\begin{align}
&したがって\color{red}{S(\phi+1,m)} = S(\phi, m) + T(m) ...(2) \\
& (2)を (1)に代入して \\
& S(\phi, N) =　T(N\phi)-T(m)-S(\phi,m) \\
\end{align}
```
という漸化式が出来ました。これをコードにします。

```python
from math import isqrt
def T(n):
  return n*(n+1)//2

def S1(N):  # 単純計算
  ret = 0
  for k in range(1,N+1):
    ret += (k+isqrt(5*k**2))//2
  return ret
  
def S2(n):  # 漸化式利用
  if n == 0: return 0
  s5n = isqrt(5*n**2)
  m = (-n+s5n)//2    # m = int(n*(phi-1))
  nr = (n+s5n)//2    # n = int(n*phi)
  return T(nr) - T(m) - S2(m)

k = 8
N = 10**k
print(f"N=10^{k}: S1(N)={S1(N)}, S2(N)={S2(N)}")
# N=10^8: S1(N)=8090169974651174, S2(N)=8090169974651174
```
これだと$N$の値が$\times 0.6188..$で減っていくので$O(log(n))$で計算できます。実際に単純計算$S1$と比べても差は歴然です。

N |$10^{9}$ |$10^{9}$ |$10^{10}$ |$10^{16}$ |
:---:|:---:|:---:|:---:|:---:|
S1| 28s| 240s| - | - | 
S2 | 0s| 0s| 0s | 0s |

（開発環境：Google Colab）