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

Beatty数列とは

あまりなじみのない名前ですがビーティ数列 (Wikipedia)とは、

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

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

\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数列の面白い性質

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

\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)」と呼ぶそうです。

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

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

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

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

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

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

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

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

という漸化式が出来ました。これをコードにします。

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）