摸了昨天长安杯的series,研究了一下模的Fibonacci,算了一下Fibonacci矩阵的阶,搞了个模n的Fibonacci的快速计算方法,记了个笔记.
题目代码:https://paste.ubuntu.com/p/ZsYfBZw6gh/
关于Fibonacci·
就是那个斐波那契数 (/数列),即:
F 1 = F 2 = 1 F n = F n − 1 + F n − 2 ( n > 2 ) F_1=F_2=1 \\
F_n = F_{n-1}+F_{n-2}\ (n > 2)
F 1 = F 2 = 1 F n = F n − 1 + F n − 2 ( n > 2 )
直接看递归式的话可能没这么方便,用矩阵形式的话会好用一点:
[ F k F k − 1 ] = [ 1 1 1 0 ] ∗ [ F k − 1 F k − 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}
[ F k F k − 1 ] = [ 1 1 1 0 ] ∗ [ F k − 1 F k − 2 ]
然后一直递归下去的话其实就是矩阵M f = [ 1 1 1 0 ] M_f=\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} M f = [ 1 1 1 0 ] 在做次方,也就得到通项:
[ F k F k − 1 ] = [ 1 1 1 0 ] k − 2 ∗ [ F 1 F 0 ] = [ 1 1 1 0 ] k − 2 ∗ [ 1 1 ] \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}
[ F k F k − 1 ] = [ 1 1 1 0 ] k − 2 ∗ [ F 1 F 0 ] = [ 1 1 1 0 ] k − 2 ∗ [ 1 1 ]
关于mod p·
先假设p是素数,mod p的Fibonacci其实就是简单的F i ( m o d p ) F_i\ (mod\ p) F i ( m o d p ) ,同样用矩阵的写法:
[ F k F k − 1 ] ≡ [ 1 1 1 0 ] k − 2 ∗ [ 1 1 ] ( m o d 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)
[ F k F k − 1 ] ≡ [ 1 1 1 0 ] k − 2 ∗ [ 1 1 ] ( m o d p )
如果可以实现M f x ( m o d p ) M_f^x\ (mod\ p) M f x ( m o d p ) 的快速计算的话,F i ( m o d p ) F_i\ (mod\ p) F i ( m o d p ) 就没什么问题了,即问题转化成了快速计算:
[ 1 1 1 0 ] x ( m o d p ) \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^x\ (mod\ p)
[ 1 1 1 0 ] x ( m o d p )
矩阵快速幂·
如果还记得整数的快速幂的话,类比一下其实矩阵可是可以做快速幂的,即可把暴力的x次矩阵乘法降到Θ ( x ) \Theta(x) Θ ( x ) 次,参考代码:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 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 x 是非常大的数(比如x = y y y . . . x=y^{y^{y^{...}}} x = y y y ... )的话,只用快速幂的优化其实是不够的。
先看一下在整数下的情况:
g x ≡ g x ( m o d ϕ ( p ) ) ≡ g x ( m o d p − 1 ) ( m o d p ) g^x \equiv g^{x\ (mod\ \phi(p))} \equiv g^{x\ (mod\ p-1)}\ (mod\ p)
g x ≡ g x ( m o d ϕ ( p )) ≡ g x ( m o d p − 1 ) ( m o d p )
是可以把x x x 降到x ( m o d p − 1 ) x\ (mod\ p-1) x ( m o d p − 1 ) 的,这里的ϕ \phi ϕ 即欧拉函数 ,也即群Z p Z_p Z p 的阶(order)。类比一下,如果如果能找到mod p矩阵的群的阶~~~~的话,就可以用类似的方法加快计算,结合快速幂,最终优化为(这里借用一下符号,假设阶是ϕ ( p ) \phi(p) ϕ ( p ) ):
Θ ( x m o d ϕ ( p ) ) = O ( ϕ ( p ) ) \Theta(x\ mod\ \phi(p))=O(\phi(p))
Θ ( x m o d ϕ ( p )) = O ( ϕ ( p ))
关于这个阶的计算,其实有挺多东西写的,所以打算以后再开一篇文章(挖个坑,逃)
mod n矩阵群的阶·
假设n = p q n=pq n = pq ,p 、 q p、 q p 、 q 是素数的话(这里是series的情况,其实n = p 1 a 1 p 1 a 1 . . . n=p_1^{a_1}p_1^{a_1}... n = p 1 a 1 p 1 a 1 ... 的情况也是通用的),有一些快速的计算方法,可以参考一下这里 的6.3节,里面说的Pisano periods其实等同Fibonacci的矩阵M f ( m o d p ) M_f\ (mod\ p) M f ( m o d p ) 的阶,大概如下:
ϕ ( p 1 a 1 p 1 a 1 . . . ) = L C M ( ϕ ( p 1 a 1 ) ϕ ( p 1 a 1 ) . . . ) \phi(p_1^{a_1}p_1^{a_1}...) = LCM(\phi(p_1^{a_1})\phi(p_1^{a_1})...) ϕ ( p 1 a 1 p 1 a 1 ... ) = L CM ( ϕ ( p 1 a 1 ) ϕ ( p 1 a 1 ) ... )
ϕ ( p a ) = p a − 1 ϕ ( p ) \phi(p^{a}) = p^{a-1}\phi(p) ϕ ( p a ) = p a − 1 ϕ ( p )
对于n来说的话,就是:
ϕ ( n ) = ϕ ( p q ) = ϕ ( p ) ϕ ( q ) \phi(n)=\phi(pq)=\phi(p)\phi(q)
ϕ ( n ) = ϕ ( pq ) = ϕ ( p ) ϕ ( q )
例题——长安杯的series·
题目实现的其实是一个类似“Triple-Fibonacci”的东西 (名字我乱取的) ,大概就是:
F 1 = F 2 = F 3 = 1 F n ≡ F n − 1 + F n − 2 ( m o d p ) , ( n > 3 ) F_1=F_2=F_3=1 \\
F_n \equiv F_{n-1}+F_{n-2}\ (mod\ p)\ ,\ (n > 3)
F 1 = F 2 = F 3 = 1 F n ≡ F n − 1 + F n − 2 ( m o d p ) , ( n > 3 )
同样可以用矩阵形式:
[ F k F k − 1 F k − 2 ] ≡ [ 1 1 1 1 0 0 0 1 0 ] k − 3 ∗ [ 1 1 1 ] ( m o d 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)
⎣ ⎡ F k F k − 1 F k − 2 ⎦ ⎤ ≡ ⎣ ⎡ 1 1 0 1 0 1 1 0 0 ⎦ ⎤ k − 3 ∗ ⎣ ⎡ 1 1 1 ⎦ ⎤ ( m o d p )
目标是要算出(e e e 是250位的随机数,实测是两个素数p e 、 p q p_e、p_q p e 、 p q 的积):
k = e e s ≡ ( ∏ i = 1 4 F k + i ) + F 202 1 k 4 ( m o d e ) k = e^e \\
s \equiv (\prod_{i=1}^4 F_{k+i}) + F_{2021^k}^4 \ (mod\ e)
k = e e s ≡ ( i = 1 ∏ 4 F k + i ) + F 202 1 k 4 ( m o d e )
(插句话:这里得暴打一下出题人,在抄代码的时候式子后面那部分忘记mod e了)
实际上如果可以mod e的Fibonacci的快速计算的话,上面问题是很容易解决的,按照上面说的一堆东西,使用“快速幂+模ϕ ( e ) \phi(e) ϕ ( e ) ”的方法,可以矩阵乘法的次数可以优化到:
O ( l o g ( ϕ ( e ) ) ) O(log(\phi(e)))
O ( l o g ( ϕ ( e )))
但是 ,怎么算这个阶呢?这里说一个结论算了:经过一些神奇的验算,矩阵M f 3 ≡ [ 1 1 1 1 0 0 0 1 0 ] ( m o d e ) M_{f3} \equiv \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}\ (mod\ e) M f 3 ≡ ⎣ ⎡ 1 1 0 1 0 1 1 0 0 ⎦ ⎤ ( m o d e ) 所在的子群的阶是e 2 − 1 e^2-1 e 2 − 1 (另:更好的阶应该是e ∗ ( e + 1 ) ∗ ( e − 1 ) ∗ ( e 2 + e + 1 ) , ( e ≠ 2 ) e*(e+1)*(e-1)*(e^2+e+1), \ (e \ne 2) e ∗ ( e + 1 ) ∗ ( e − 1 ) ∗ ( e 2 + e + 1 ) , ( e = 2 ) ,e 2 − 1 e^2-1 e 2 − 1 是这里其中一个因子,不过说了到时再补坑就到时再说吧,逃)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 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)
所以可以算出:
ϕ ( p e ) = p e 2 − 1 , ϕ ( q e ) = q e 2 − 1 ϕ ( e ) = ϕ ( p e ) ϕ ( q e ) = ( p e 2 − 1 ) ( q e 2 − 1 ) \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)
ϕ ( p e ) = p e 2 − 1 , ϕ ( q e ) = q e 2 − 1 ϕ ( e ) = ϕ ( p e ) ϕ ( q e ) = ( p e 2 − 1 ) ( q e 2 − 1 )
最终可以把模e e e 的三阶(rank)矩阵的乘法次数优化为:
O ( l o g ( ϕ ( e ) ) ) = O ( ( p e 2 − 1 ) ( q e 2 − 1 ) ) O(log(\phi(e))) = O((p_e^2-1)(q_e^2-1))
O ( l o g ( ϕ ( 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 from Crypto.Util.number import *from gmpy2 import next_primedef 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 ep = 33285073849485750791903437807279991921 eq = 38845988283830087557982578789883120419 assert ep*eq == ephi = 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 )) 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)))
要不是昨天早上去见大一小朋友了+中午出去吃了顿饭就可以拿一血了(逃