フィボナッチ数列とピサノ周期

ピサノ周期

フィボナッチ数列に関する投稿はかなりしてきましたが、その剰余を取った数列は周期的になりそれをピサノ周期と呼ぶということを初めて知りました。

まずフィボナッチ数$F_{i} \pmod{n}$を表にしてみます。赤い部分が周期的に現れているように見えます。

def fibmodgen(M): # Fibonacci number generator w/ modulo
  f1, f2 = 0, 1
  while True:
    yield f1
    f2, f1, = (f1+f2)%M, f2

for n in range(2,8+1):
  fib = fibmodgen(n)
  print(n, [next(fib) for i in range(20)])

n	$F_{0}$	$F_{1}$	$F_{2}$	$F_{3}$	$F_{4}$	$F_{5}$	$F_{6}$	$F_{7}$	$F_{8}$	$F_{9}$	$F_{10}$	$F_{11}$	$F_{12}$	$F_{13}$	$F_{14}$	$F_{15}$	$F_{16}$	$F_{17}$	$F_{18}$	$F_{19}$
2	0	1	1	0	1	1	0	1	1	0	1	1	0	1	1	0	1	1	0	1
3	0	1	1	2	0	2	2	1	0	1	1	2	0	2	2	1	0	1	1	2
4	0	1	1	2	3	1	0	1	1	2	3	1	0	1	1	2	3	1	0	1
5	0	1	1	2	3	0	3	3	1	4	0	4	4	3	2	0	2	2	4	1
6	0	1	1	2	3	5	2	1	3	4	1	5	0	5	5	4	3	1	4	5
7	0	1	1	2	3	5	1	6	0	6	6	5	4	2	6	1	0	1	1	2
8	0	1	1	2	3	5	0	5	5	2	7	1	0	1	1	2	3	5	0	5

このピサノ周期を以下のステップで求めていきたいを思います。

1. nが素数のときに0が現れる周期$\alpha(p)$を求める
2. nが$p^k$ののときに0が現れる周期$\alpha(p^k)$を求める
3. nが正の整数のときの0が現れる周期$\alpha(n)$を求める
4. ピサノ周期を0が現れる周期の倍数と仮定して、nのピサノ周期を求める

１．nが素数のときに0が現れる周期を求める

$\pmod{p}の$フィボナッチ数列の0が現れる周期を$\alpha(p)$と表すとき$p<40$で$\alpha(p)$を求めます

p	2	3	5	7	11	13	17	19	23	29	31	37
$\alpha$	3	4	5	8	10	7	9	18	24	14	30	19

これをよく見ると$p=5$以外では$\alpha(p)$は$(p+1)$または$(p-1)$の約数になっていそうな気がします。この定理とかを探してみるとありました。

Wall–Sun–Sun prime(Wikipedia)

この内容はかなり専門的なので要約すると。

\begin{align}
&p>=7のとき　\\
&F_{p-1} \equiv 0 \pmod p  : p \equiv 1 or 4 \pmod 5\\
&F_{p+1} \equiv 0 \pmod p  : p \equiv 2 or 3 \pmod 5\\
&このとき(d|p \pm 1)でF_d \equiv 0 \pmod p となる\\
&dの最小値が\alpha(p)
\end{align}

プログラムは結構シンプルです。fibn(n,p) は$F_n \pmod p$をもとめる関数で「フィボナッチ数列を行列を使ってO(log(n))で計算する」で作ったものを使っています

import sympy as sp
def alpha_prime(p):
  fn = [p, p-1, p+1, p+1, p-1][p%5]  # # fn depends on p mod 5 
  for d in sp.divisors(fn)[1:]:  # find minimum d fib(d,p) == 0 
    if fibn(d, p) == 0: break
  return d

print([alpha_prime(p) for p in sp.primerange(2,100)])
#[3, 4, 5, 8, 10, 7, 9, 18, 24, 14, 30, 19, 20, 44, 16, 27, 58, 15, 68, 70, 37, 78, 84, 11, 49]

２．nが素数のべき乗のときに0が現れる周期を求める

Pisano period->Number of zeros in the cycle (Wikipedia)によれば（まだ仮説とのことですが）

\begin{align}
&\alpha(p^k) = p^{k-1} \alpha(p)
\end{align}

これを使って$\alpha(p^k)$を求めるalpha_prime_powを作ります。ただ検証してみると１つだけ例外があって$p=2$で$k>=3$のときは$k=k-1$とする必要がありました。

def alpha_prime_pow(p,k):
  if p==2 and k >= 3: k -= 1  # p=2 and k>=3 is an exception case
  return p**(k-1)*alpha_prime(p)
  
print([alpha_prime_pow(2,k) for k in range(1,10+1)])
# [3, 6, 6, 12, 24, 48, 96, 192, 384, 768]

３．nが正の整数のときの0が現れる周期を求める

これも同じWikipediaからの式で、

\begin{align}
&n = p_1^{k_1} p_2^{k_2} \dots p_m^{k_m}と因数分解したとき\\
&\alpha(n) = lcm(\alpha( p_1^{k_1}),\alpha(p_2^{k_2}),\dots,\alpha( p_m^{k_m}))
\end{align}

これをコードにすると、

def alpha(n):
  fct = sp.factorint(n)
  al = 1
  for p in fct.keys():
    al = sp.lcm(al,alpha_prime_pow(p, fct[p]))
  return al 
print([alpha(n) for n in range(1,20+1)])
# [1, 3, 4, 6, 5, 12, 8, 6, 12, 15, 10, 12, 7, 24, 20, 12, 9, 12, 18, 30]

ピサノ周期を0が現れる周期の倍数と仮定して、nのピサノ周期を求める

さてやっとピサノ周期を求める準備ができました。やはりWikipediaによると0が現れる周期はピサノ周期に1か2か4回現れるということなので、いずれかの場合に0の次に1が来るかをチェックすればピサノ周期が求まると仮定してコードを書きました。

def pisano(n): 
  al = alpha(n)
  for m in [1,2,4]:
    if fibn(al*m+1, n) == 1: 
      return al*m
  return al

このコードを使ってWikipediaにあるのと同じ表を作って$n<=144$までは正しいことを確認しました。

$\pi(n)$	1	2	3	4	5	6	7	8	9	10	11	12
0+	1	3	8	6	20	24	16	12	24	60	10	24
12+	28	48	40	24	36	24	18	60	16	30	48	24
24+	100	84	72	48	14	120	30	48	40	36	80	24
36+	76	18	56	60	40	48	88	30	120	48	32	24
48+	112	300	72	84	108	72	20	48	72	42	58	120
60+	60	30	48	96	140	120	136	36	48	240	70	24
72+	148	228	200	18	80	168	78	120	216	120	168	48
84+	180	264	56	60	44	120	112	48	120	96	180	48
96+	196	336	120	300	50	72	208	84	80	108	72	72
108+	108	60	152	48	76	72	240	42	168	174	144	120
120+	110	60	40	30	500	48	256	192	88	420	130	120
132+	144	408	360	36	276	48	46	240	32	210	140	24

（開発環境：Google Colab）

この考え方はProject Euler、Problem399: Squarefree Fibonacci Numbersを解くのに役に立ちます。