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

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

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

```python 
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}$|
|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
|<b>2</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|0|1|1|0|1|1|0|1|1|0|1|1|0|1|1|0|1|
|<b>3</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|<font color='Red'>2</font>|<font color='Red'>0</font>|<font color='Red'>2</font>|<font color='Red'>2</font>|<font color='Red'>1</font>|0|1|1|2|0|2|2|1|0|1|1|2|
|<b>4</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|<font color='Red'>2</font>|<font color='Red'>3</font>|<font color='Red'>1</font>|0|1|1|2|3|1|0|1|1|2|3|1|0|1|
|<b>5</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|<font color='Red'>2</font>|<font color='Red'>3</font>|<font color='Red'>0</font>|<font color='Red'>3</font>|<font color='Red'>3</font>|<font color='Red'>1</font>|<font color='Red'>4</font>|<font color='Red'>0</font>|<font color='Red'>4</font>|<font color='Red'>4</font>|<font color='Red'>3</font>|<font color='Red'>2</font>|<font color='Red'>0</font>|<font color='Red'>2</font>|<font color='Red'>2</font>|<font color='Red'>4</font>|<font color='Red'>1</font>|
|<b>6</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|<font color='Red'>2</font>|<font color='Red'>3</font>|<font color='Red'>5</font>|<font color='Red'>2</font>|<font color='Red'>1</font>|<font color='Red'>3</font>|<font color='Red'>4</font>|<font color='Red'>1</font>|<font color='Red'>5</font>|<font color='Red'>0</font>|<font color='Red'>5</font>|<font color='Red'>5</font>|<font color='Red'>4</font>|<font color='Red'>3</font>|<font color='Red'>1</font>|<font color='Red'>4</font>|<font color='Red'>5</font>|
|<b>7</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|<font color='Red'>2</font>|<font color='Red'>3</font>|<font color='Red'>5</font>|<font color='Red'>1</font>|<font color='Red'>6</font>|<font color='Red'>0</font>|<font color='Red'>6</font>|<font color='Red'>6</font>|<font color='Red'>5</font>|<font color='Red'>4</font>|<font color='Red'>2</font>|<font color='Red'>6</font>|<font color='Red'>1</font>|0|1|1|2|
|<b>8</b>|<font color='Red'>0</font>|<font color='Red'>1</font>|<font color='Red'>1</font>|<font color='Red'>2</font>|<font color='Red'>3</font>|<font color='Red'>5</font>|<font color='Red'>0</font>|<font color='Red'>5</font>|<font color='Red'>5</font>|<font color='Red'>2</font>|<font color='Red'>7</font>|<font color='Red'>1</font>|0|1|1|2|3|5|0|5|

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

1. nが素数のときに0が現れる周期$\alpha(p)$を求める
1. nが$p^k$ののときに0が現れる周期$\alpha(p^k)$を求める
1. nが正の整数のときの0が現れる周期$\alpha(n)$を求める
1. ピサノ周期を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)](https://en.wikipedia.org/wiki/Wall%E2%80%93Sun%E2%80%93Sun_prime)

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

```math 
\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))で計算する](https://qiita.com/masa0599/items/5d412daa96695dea74bf)」で作ったものを使っています

```python 
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)](https://en.wikipedia.org/wiki/Pisano_period#Number_of_zeros_in_the_cycle)によれば（まだ仮説とのことですが）

```math 
\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$とする必要がありました。

```python
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からの式で、

```math 
\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}
```
これをコードにすると、

```python 
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が来るかをチェックすればピサノ周期が求まると仮定してコードを書きました。

```python 
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を解くのに役に立ちます。