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

题目代码:https://paste.ubuntu.com/p/ZsYfBZw6gh/

关于Fibonacci·

就是那个斐波那契数(/数列),即:

F1=F2=1Fn=Fn1+Fn2 (n>2)F_1=F_2=1 \\ F_n = F_{n-1}+F_{n-2}\ (n > 2)

直接看递归式的话可能没这么方便,用矩阵形式的话会好用一点:

[FkFk1]=[1110][Fk1Fk2]\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}

然后一直递归下去的话其实就是矩阵Mf=[1110]M_f=\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}在做次方,也就得到通项:

[FkFk1]=[1110]k2[F1F0]=[1110]k2[11]\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其实就是简单的Fi (mod p)F_i\ (mod\ p),同样用矩阵的写法:

[FkFk1][1110]k2[11] (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)

如果可以实现Mfx (mod p)M_f^x\ (mod\ p)的快速计算的话,Fi (mod p)F_i\ (mod\ p)就没什么问题了,即问题转化成了快速计算:

[1110]x (mod p)\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^x\ (mod\ p)

矩阵快速幂·

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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矩阵群的阶·

如果xx是非常大的数(比如x=yyy...x=y^{y^{y^{...}}})的话,只用快速幂的优化其实是不够的。

先看一下在整数下的情况:

gxgx (mod ϕ(p))gx (mod p1) (mod p)g^x \equiv g^{x\ (mod\ \phi(p))} \equiv g^{x\ (mod\ p-1)}\ (mod\ p)

是可以把xx降到x (mod p1)x\ (mod\ p-1)的,这里的ϕ\phi欧拉函数,也即群ZpZ_p的阶(order)。类比一下,如果如果能找到mod p矩阵的群的阶~~~~的话,就可以用类似的方法加快计算,结合快速幂,最终优化为(这里借用一下符号,假设阶是ϕ(p)\phi(p)):

Θ(x mod ϕ(p))=O(ϕ(p))\Theta(x\ mod\ \phi(p))=O(\phi(p))

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

mod n矩阵群的阶·

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

  • ϕ(p1a1p1a1...)=LCM(ϕ(p1a1)ϕ(p1a1)...)\phi(p_1^{a_1}p_1^{a_1}...) = LCM(\phi(p_1^{a_1})\phi(p_1^{a_1})...)
  • ϕ(pa)=pa1ϕ(p)\phi(p^{a}) = p^{a-1}\phi(p)

对于n来说的话,就是:

ϕ(n)=ϕ(pq)=ϕ(p)ϕ(q)\phi(n)=\phi(pq)=\phi(p)\phi(q)

例题——长安杯的series·

题目实现的其实是一个类似“Triple-Fibonacci”的东西 (名字我乱取的) ,大概就是:

F1=F2=F3=1FnFn1+Fn2 (mod p) , (n>3)F_1=F_2=F_3=1 \\ F_n \equiv F_{n-1}+F_{n-2}\ (mod\ p)\ ,\ (n > 3)

同样可以用矩阵形式:

[FkFk1Fk2][111100010]k3[111] (mod p)\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)

目标是要算出(ee是250位的随机数,实测是两个素数pepqp_e、p_q的积):

k=ees(i=14Fk+i)+F2021k4 (mod e)k = e^e \\ s \equiv (\prod_{i=1}^4 F_{k+i}) + F_{2021^k}^4 \ (mod\ e)

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

实际上如果可以mod e的Fibonacci的快速计算的话,上面问题是很容易解决的,按照上面说的一堆东西,使用“快速幂+模ϕ(e)\phi(e)”的方法,可以矩阵乘法的次数可以优化到:

O(log(ϕ(e)))O(log(\phi(e)))

但是,怎么算这个阶呢?这里说一个结论算了:经过一些神奇的验算,矩阵Mf3[111100010] (mod e)M_{f3} \equiv \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}\ (mod\ e) 所在的子群的阶是e21e^2-1(另:更好的阶应该是e(e+1)(e1)(e2+e+1), (e2)e*(e+1)*(e-1)*(e^2+e+1), \ (e \ne 2)e21e^2-1是这里其中一个因子,不过说了到时再补坑就到时再说吧,逃)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 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) # []

所以可以算出:

ϕ(pe)=pe21, ϕ(qe)=qe21ϕ(e)=ϕ(pe)ϕ(qe)=(pe21)(qe21)\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)

最终可以把模ee的三阶(rank)矩阵的乘法次数优化为:

O(log(ϕ(e)))=O((pe21)(qe21))O(log(\phi(e))) = O((p_e^2-1)(q_e^2-1))

参考代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
# 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)))

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