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

关于Fibonacci

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

[math] \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} [/math]

关于mod p

先假设p是素数，mod p的Fibonacci其实就是简单的 [imath]F_i\ (mod\ p)[/imath] ，同样用矩阵的写法：
[math] \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) [/math]
如果可以实现 [imath]M_f^x\ (mod\ p)[/imath] 的快速计算的话， [imath]F_i\ (mod\ p)[/imath] 就没什么问题了，即问题转化成了快速计算：
[math] \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^x\ (mod\ p) [/math]

矩阵快速幂

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

# 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矩阵环的阶

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

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

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

mod n矩阵环的阶

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

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

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

例题——长安杯的series

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

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

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

# 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)  # []

所以可以算出：
[math] \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) [/math]
最终可以把模 [imath]e[/imath] 的三阶（rank）矩阵的乘法次数优化为：
[math] O(log(\phi(e))) = O((p_e^2-1)(q_e^2-1)) [/math]
参考代码：

# 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]不是$$的啊！（