引入
给定一正整数
N
∈
N
∗
N\in \mathbb{N} ^{*}
N∈N∗,试快速找到它的一个因数。
很久很久以前,我们曾学过试除法来解决这个问题。很容易想到因数是成对称分布的:即
N
N
N的所有因数可以被分成两块:
[
1
,
N
]
[1,\sqrt{N} ]
[1,N]和
[
N
,
N
]
[\sqrt{N},N ]
[N,N]。这个很容易想清楚,我们只需要把区间
[
1
,
N
]
[1,\sqrt{N} ]
[1,N]扫一遍就可以了,此时试除法的时间复杂度为
O
(
N
)
O(\sqrt{N} )
O(N)。
对于
N
≥
1
0
18
N\ge 10^{18}
N≥1018的数据,这个算法运行起来无疑是非常糟糕的。我们希望有更快的算法。对于这样的算法,一个比较好的想法是去设计一个随机程序,随便猜一个因数。如果你运气好,这个程序的时间复杂度下限为
O
(
1
)
O(1)
O(1)。但对于一个
N
≥
1
0
18
N\ge 10^{18}
N≥1018的数据,这个算法给出答案的概率是
1
1000000000000000000
\frac{1}{1000000000000000000}
10000000000000000001。当然,如果在
[
1
,
N
]
[1,\sqrt{N} ]
[1,N]里面猜,成功的可能性会更大。那么,有没有更好的改进算法来提高我们猜的准确率呢?那就是Pollard’s rho算法。
算法流程
Pollard’s rho算法的重要思想就是:最大公约数一定是某个数的约数。也就是说, ∀ k ∈ N ∗ \forall k\in \mathbb{N}^{*} ∀k∈N∗, g c d ( k , n ) ∣ n gcd(k,n) |n gcd(k,n)∣n。只要选适当的k使得 g c d ( k , n ) > 1 gcd(k,n)>1 gcd(k,n)>1就可以求得一个约数 g c d ( k , n ) gcd(k,n) gcd(k,n)。这样的k数量还是蛮多的:k有若干个质因子,而每个质因子的倍数都是可行的。
如果这样的k我们随机取的话,就体现不出Pollard’s rho算法的精妙之处了。我们不妨考虑构造一个伪随机数序列,然后取相邻的两项的差来求gcd。为了生成一串优秀的随机数,Pollard设计了这样一个函数: f ( x ) = ( x 2 + c ) m o d N f(x)=(x^{2}+c ) \mod N f(x)=(x2+c)modN,其中 c c c是一个随机的常数。
我们随便取一个 x 1 x_{1} x1,令 x 2 = f ( x 1 ) x_{2} =f(x_{1} ) x2=f(x1), x 3 = f ( x 2 ) x_{3} =f(x_{2} ) x3=f(x2),……, x i = f ( x i − 1 ) x_{i} =f(x_{i-1} ) xi=f(xi−1)。在一定的范围内,这个数列是基本随机的,可以取相邻两项作差求gcd。
这样构造的数列会有一些性质:
(1)相邻两数的差会包含前面所有相邻两数差的乘积,这些乘积中可能会包含n的因子。考虑下面这个式子:
x
i
−
x
i
−
1
=
f
(
x
i
−
1
)
−
f
(
x
i
−
2
)
=
x
i
−
1
2
−
x
i
−
2
2
=
(
x
i
−
1
+
x
i
−
2
)
(
x
i
−
1
−
x
i
−
2
)
(
m
o
d
n
)
x_{i}-x_{i-1}= f(x_{i-1} )-f(x_{i-2})={x_{i-1}}^{2} -{x_{i-2}}^{2}=(x_{i-1}+x_{i-2})(x_{i-1}-x_{i-2}) \pmod n
xi−xi−1=f(xi−1)−f(xi−2)=xi−12−xi−22=(xi−1+xi−2)(xi−1−xi−2)(modn)
在等式最右边
x
i
−
1
x_{i-1}
xi−1和
x
i
−
2
x_{i-2}
xi−2又可以继续分解,直到
x
1
x_{1}
x1。所以会有上面那个性质。
(2)序列中距离为k的两项的差,一定为前面任意距离为k的两项的差的倍数。这个性质算是上面那个性质的推广,当
k
=
1
k=1
k=1时,其实和性质(1)是一样的。在序列中取
x
i
,
x
j
x_{i}, x_{j}
xi,xj,假设
i
>
j
i>j
i>j,则有:
x
i
−
x
j
=
f
(
x
i
−
1
)
−
f
(
x
j
−
1
)
=
x
i
−
1
2
−
x
j
−
1
2
=
(
x
i
−
1
+
x
j
−
1
)
(
x
i
−
1
−
x
j
−
1
)
(
m
o
d
N
)
x_{i}-x_{j}=f(x_{i-1})-f(x_{j-1})={x_{i-1}}^{2}-{x_{j-1}}^{2}=(x_{i-1}+x_{j-1})(x_{i-1}-x_{j-1}) \pmod N
xi−xj=f(xi−1)−f(xj−1)=xi−12−xj−12=(xi−1+xj−1)(xi−1−xj−1)(modN)
其中等式最右边
x
i
−
1
x_{i-1}
xi−1和
x
j
−
1
x_{j-1}
xj−1又可以继续分解。
(3)序列中任意两数的差,也一定可以转换为相邻两个数的差的倍数。在序列中取
x
i
,
x
j
x_{i}, x_{j}
xi,xj,假设
i
>
j
i>j
i>j,则有:
x
i
−
x
j
=
x
i
−
x
i
−
1
+
x
i
−
1
−
⋯
−
x
j
+
x
j
−
x
j
−
1
=
b
∗
(
x
j
−
x
j
−
1
)
(
m
o
d
N
)
x_{i}-x_{j}=x_{i}-x_{i-1}+x_{i-1}-\dots -x_{j}+x_{j}-x_{j-1}=b*(x_{j}-x_{j-1}) \pmod N
xi−xj=xi−xi−1+xi−1−⋯−xj+xj−xj−1=b∗(xj−xj−1)(modN)
根据性质(1),相邻两数的差会包含前面所有相邻两数差的乘积。所以,
x
i
−
x
i
−
1
x_{i}-x_{i-1}
xi−xi−1一直到
x
j
+
1
−
x
j
x_{j+1}-x_{j}
xj+1−xj会包含
x
j
−
x
j
−
1
x_{j}-x_{j-1}
xj−xj−1,将所有的无关项提取出来形成一个b就得到了上面等式的最右边。
根据前面的分析,序列有一些特殊的性质,即下标距离为k的两个数之差,是前面的下标距离为k的某两个数之差的倍数。所以,对于距离为k的两个数的差,我们只需要检测最后那一对即可。这样,每一个距离都只需要检测一次。
为了判断环的存在,可以用一个简单的Floyd判圈算法,也就是"龟兔赛跑"。假设乌龟为 t t t,兔子为 r r r,初始时 t = r = 1 t=r=1 t=r=1。假设兔子的速度是乌龟的二倍。过了时间 i i i后, t = i , r = 2 i t=i,r=2i t=i,r=2i。此时两者得到的数列值 x t = x i , x r = x 2 i x_{t}=x_{i},x_{r}=x_{2i} xt=xi,xr=x2i。假设环的长度为 c c c,在环内恒有: x i = x i + c x_{i}=x_{i+c} xi=xi+c。 如果龟兔"相遇",此时有: x r = x t x_{r}=x_{t} xr=xt,也就是 x i = x 2 i = x i + k c x_{i}=x_{2i}=x_{i+kc} xi=x2i=xi+kc。此时两者路径之差正好是环长度的整数倍。
出现这种乌龟被套圈的情况,再继续下去只会一直重复,所以需要通过调整函数 f ( x ) f(x) f(x)里常数项的值来继续“碰运气”。
这样一来,我们得到了一套基于Floyd判圈算法的Pollard Rho 算法。Python实现如下
import random
from math import *
from Crypto.Util.number import isPrime
def f(x, c, n):
return (x * x + c) % n
def Pollard_rho(n):
if isPrime(n):
return n
while True:
c = random.randint(1, n - 1)
t = f(0, c, n)
r = f(t, c, n)
while t != r:
d = gcd(abs(t - r), n)
if d > 1:
return d
t = f(t, c, n)
r = f(r, c, n)
r = f(r, c, n)
n = int(input())
fac = Pollard_rho(n)
print(fac)
椭圆曲线上的Pollard’s rho算法
设椭圆曲线为 G G G, P , Q P,Q P,Q为椭圆曲线上两点,欲求出满足 d P = Q dP=Q dP=Q的 d d d。
选择哈希函数将
G
G
G分成尺寸大致相同的三部分
S
1
,
S
2
,
S
3
S_{1},S_{2} ,S_{3}
S1,S2,S3,其中
O
∉
S
2
O\notin S_{2}
O∈/S2;定义一个随机步数的迭代函数
f
f
f:
R
i
+
1
=
f
(
R
i
)
=
{
Q
+
R
i
,
R
i
∈
S
1
2
R
i
,
R
i
∈
S
2
P
+
R
i
,
R
i
∈
S
3
R_{i+1}=f(R_{i})=\left\{\begin{matrix} Q+R_{i}, & R_{i} \in S_{1}\\ 2R_{i}, & R_{i} \in S_{2} \\ P+R_{i}, & R_{i} \in S_{3} \end{matrix}\right.
Ri+1=f(Ri)=⎩
⎨
⎧Q+Ri,2Ri,P+Ri,Ri∈S1Ri∈S2Ri∈S3
令
R
i
=
a
i
∗
P
+
b
i
∗
Q
R_{i}=a_{i}*P+b_{i}*Q
Ri=ai∗P+bi∗Q,则:
a
i
+
1
=
{
a
i
,
R
i
∈
S
1
2
a
i
m
o
d
N
,
R
i
∈
S
2
a
i
+
1
,
R
i
∈
S
3
a_{i+1}=\left\{\begin{matrix} a_{i}, & R_{i} \in S_{1}\\ 2a_{i} \mod N, & R_{i} \in S_{2} \\ a_{i}+1, & R_{i} \in S_{3} \end{matrix}\right.
ai+1=⎩
⎨
⎧ai,2aimodN,ai+1,Ri∈S1Ri∈S2Ri∈S3
b
i
+
1
=
{
b
i
+
1
,
R
i
∈
S
1
2
b
i
m
o
d
N
,
R
i
∈
S
2
b
i
,
R
i
∈
S
3
b_{i+1}=\left\{\begin{matrix} b_{i}+1, & R_{i} \in S_{1}\\ 2b_{i} \mod N, & R_{i} \in S_{2} \\ b_{i}, & R_{i} \in S_{3} \end{matrix}\right.
bi+1=⎩
⎨
⎧bi+1,2bimodN,bi,Ri∈S1Ri∈S2Ri∈S3
初始化参数
R
0
=
P
,
a
0
=
1
,
b
0
=
0
R_{0}=P,a_{0}=1,b_{0}=0
R0=P,a0=1,b0=0,产生配对
(
R
i
,
R
2
i
)
(R_{i},R_{2i})
(Ri,R2i) ,直到对某个
m
m
m有
R
m
=
R
2
m
R_{m}=R_{2m}
Rm=R2m,此时有:
R
m
=
a
m
P
+
b
m
Q
R_{m}=a_{m}P+b_{m}Q
Rm=amP+bmQ
R
2
m
=
a
2
m
P
+
b
2
m
Q
R_{2m}=a_{2m}P+b_{2m}Q
R2m=a2mP+b2mQ
所以有
d
=
a
2
m
−
a
m
b
m
−
b
2
m
m
o
d
N
d=\frac{a_{2m}-a_{m}}{b_{m}-b_{2m}} \mod N
d=bm−b2ma2m−ammodN