摸了昨天长安杯的series，研究了一下模的Fibonacci，算了一下Fibonacci矩阵的阶，搞了个模n的Fibonacci的快速计算方法，记了个笔记.

关于Fibonacci

就是那个斐波那契数（/数列），即：
F_1=F_2=1 \\ F_n = F_{n-1}+F_{n-2}\ (n > 2)
直接看递归式的话可能没这么方便，用矩阵形式的话会好用一点：
\begin{bmatrix} F_k \\ F_{k-1} \end{bmatrix}= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}* \begin{bmatrix} F_{k-1} \\ F_{k-2} \end{bmatrix}
然后一直递归下去的话其实就是矩阵M_f=\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}在做次方，也就得到通项：

\begin{bmatrix} F_{k} \\ F_{k-1} \end{bmatrix}= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{k-2} * \begin{bmatrix} F_1 \\ F_0 \end{bmatrix}= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{k-2} * \begin{bmatrix} 1 \\ 1 \end{bmatrix}

关于mod p

先假设p是素数，mod p的Fibonacci其实就是简单的F_i\ (mod\ p)，同样用矩阵的写法：
\begin{bmatrix} F_k \\ F_{k-1} \end{bmatrix} \equiv \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{k-2} * \begin{bmatrix} 1 \\ 1 \end{bmatrix} \ (mod\ p)
如果可以实现M_f^x\ (mod\ p)的快速计算的话，F_i\ (mod\ p)就没什么问题了，即问题转化成了快速计算：
\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^x\ (mod\ p)

矩阵快速幂

如果还记得整数的快速幂的话，类比一下其实矩阵可是可以做快速幂的，即可把暴力的x次矩阵乘法降到\Theta(x)次，参考代码：

# Sage
def fastPow(M, x, n=None):
    if x==0:
        return identity_matrix(len(M.columns()))
    if n!=None:
        res = fastPow(M, x//2, n)
    else:
        res = fastPow(M, x//2)
    if x%2 == 1:
        if n!=None:
            return res*res*M % n
        else:
            return res*res*M
    else:
        if n!=None:
            return res*res % n
        else:
            return res*res

PS：在sage中如果把矩阵定义在mod p环上的话，其实就不用再在后面做"% n"了

mod p矩阵群的阶

如果x是非常大的数（比如x=y^{y^{y^{...}}}）的话，只用快速幂的优化其实是不够的。

先看一下在整数下的情况：
g^x \equiv g^{x\ (mod\ \phi(p))} \equiv g^{x\ (mod\ p-1)}\ (mod\ p)
是可以把x降到x\ (mod\ p-1)的，这里的\phi即欧拉函数，也即群Z_p的阶（order）。类比一下，如果如果能找到mod p矩阵的群的阶的话，就可以用类似的方法加快计算，结合快速幂，最终优化为（这里借用一下符号，假设阶是\phi(p)）：
\Theta(x\ mod\ \phi(p))=O(\phi(p))

关于这个阶的计算，其实有挺多东西写的，所以打算以后再开一篇文章（挖个坑，逃）

mod n矩阵群的阶

假设n=pq，p、 q是素数的话（这里是series的情况，其实n=p_1^{a_1}p_1^{a_1}...的情况也是通用的），有一些快速的计算方法，可以参考一下这里的6.3节，里面说的Pisano periods其实等同Fibonacci的矩阵M_f\ (mod\ p)的阶，大概如下：

\phi(p_1^{a_1}p_1^{a_1}...) = LCM(\phi(p_1^{a_1})\phi(p_1^{a_1})...)
\phi(p^{a}) = p^{a-1}\phi(p)

对于n来说的话，就是：
\phi(n)=\phi(pq)=\phi(p)\phi(q)

例题——长安杯的series

题目实现的其实是一个类似“Triple-Fibonacci”的东西~~（名字我乱取的）~~，大概就是：
F_1=F_2=F_3=1 \\ F_n \equiv F_{n-1}+F_{n-2}\ (mod\ p)\ ,\ (n > 3)
同样可以用矩阵形式：
\begin{bmatrix} F_k \\ F_{k-1} \\ F_{k-2} \end{bmatrix} \equiv \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}^{k-3} * \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \ (mod\ p)
目标是要算出（e是250位的随机数，实测是两个素数p_e、p_q的积）：
k = e^e \\ s \equiv (\prod_{i=1}^4 F_{k+i}) + F_{2021^k}^4 \ (mod\ e)

（插句话：这里得暴打一下出题人，在抄代码的时候式子后面那部分忘记mod e了）

实际上如果可以mod e的Fibonacci的快速计算的话，上面问题是很容易解决的，按照上面说的一堆东西，使用“快速幂+模\phi(e)”的方法，可以矩阵乘法的次数可以优化到：
O(log(\phi(e)))
但是，怎么算这个阶呢？这里说一个结论算了：经过一些神奇的验算，矩阵M_{f3}\equiv\begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}\ (mod\ e) 所在的子群的阶是e^2-1（另：更好的阶应该是e*(e+1)*(e-1)*(e^2+e+1), \ (e \ne 2)，e^2-1是这里其中一个因子，不过说了到时再补坑就到时再说吧，逃）

# Sage
except_p = []
p = 2
for i in range(10000):
  p = next_prime(p)
  I = identity_matrix(IntegerModRing(p), 3)
  m = matrix(IntegerModRing(p), [[1, 1, 1], [1, 0, 0], [0, 1, 0]])
  order = p*(p+1)*(p-1)*(p^2+p+1)
  try:
    assert pow(m, order)==I
  except:
    print(p)
    print(pow(m, order))
    except_p.append(p)
print('done!')
print(except_p)  # []

所以可以算出：
\phi(p_e) = p_e^2-1 ,\ \phi(q_e) = q_e^2-1 \\ \phi(e) = \phi(p_e)\phi(q_e) = (p_e^2-1)(q_e^2-1)
最终可以把模e的三阶（rank）矩阵的乘法次数优化为：
O(log(\phi(e))) = O((p_e^2-1)(q_e^2-1))
参考代码：

# Sage
from Crypto.Util.number import *
from gmpy2 import next_prime

def fastPow(M, x, n=None):
    if x==0:
        return identity_matrix(len(M.columns()))
    if n!=None:
        res = fastPow(M, x//2, n)
    else:
        res = fastPow(M, x//2)
    if x%2 == 1:
        if n!=None:
            return res*res*M % n
        else:
            return res*res*M
    else:
        if n!=None:
            return res*res % n
        else:
            return res*res

e = 1292991588783542706506728336494377723983115217051171962646571511384590134899
c = 229797522574801936576076488492034448896863980731763047709941641260180597290800402814069755381965565755866855389082787759443816945304000719176334587540293777658369250939545994974691382505993209963323032684771922094686136104097942892330051349688373437571196103392801691879287264056022383484359551333197
# yafu!
ep = 33285073849485750791903437807279991921
eq = 38845988283830087557982578789883120419
assert ep*eq == e
phi = int((ep^2-1)*(eq^2-1))
phi2 = int(euler_phi(ep+1)*euler_phi(ep-1)*euler_phi(eq+1)*euler_phi(eq-1))  # euler_phi(phi)

one = matrix(ZZ, [1, 1, 1]).transpose()
Fme = matrix(IntegerModRing(e), [[1, 1, 1], [1, 0, 0], [0, 1, 0]])

def Fe(n, i=0):
  return int((Fme^n*Fme^(i-3)*one)[0][0])%e

def Fe2(n, i=0):
  n = int(n)
  return int((Fme^pow(n, n, phi)*Fme^(i-3)*one)[0][0])%e

def Fe3(m, n, i=0):
  m = int(m)
  n = int(n)
  return int((Fme^pow(m, pow(n, n, phi2), phi)*Fme^(i-3)*one)[0][0])%e

s1 = reduce(lambda a, b: a*b, [Fe2(e, i) for i in range(1, 5)])
s2 = Fe3(2021, e)^4
s = s1+s2
p = next_prime(s)
print('c = %s' % c)
print('p = %s' % p)

d = e.inverse_mod(p-1)
m = pow(c, d, p)
print(long_to_bytes(int(m)))

~~要不是昨天早上去见大一小朋友了+中午出去吃了顿饭就可以拿一血了（逃~~

原文链接：https://tover.xyz/2021/09/26/Fibonacci-series/
PS：怎么论坛公式的tag是[math]不是$$的啊！（