yafu为什么事神

这篇文章主要是想总结一下四个分解因子算法 向大家推荐一下yafu 。

下面所有算法都是基于分解 $n=pq$ （$p,q$ 是素数）这个问题。

$Pollard\ p-1$ 算法

原理

​ 很简单啊，就事在自然数里找一个 $ s $ ，满足 $ s=1(mod p) $ ，这不就说明 $ p $ 是 $ （s-1） $ 的一个因子吗，然后计算 $ gcd(s-1,n) $ 。好了，讲完了

实现

​ 虽然原理看上去是个原理，但是问题在于如何寻找这个 $ s $ 。在不知道 $ p $ 的情况下，我们很难找到满足上式的 $ s $ ，不过伟大的费马给出了他的公式：
$$
a^{p-1}\equiv1\ mod\ p
$$
所以我们要计算p-1（混乱 虽然直接计算p-1还是很困难，但是我们可以直接修改费马公式啊（
$$
a^{k}\equiv1\ mod\ p\ \ \ \ if\ \ \ (p-1)|k
$$
根据这个原理，Pollard给出了他的算法。我们假设 $ （p-1） $ 有一个最大因子 $ B $ ，那么就一定满足$(p-1)|B!$ ，同时我们再令 $ a=2 $ 方便计算。最终我们得出结论，最多计算 $ B $ 次模幂运算就可以得到我们想要的 $ s $ （也就是$ a^{B!} $）。

代码

import math
def Pollard_p_1(n,B):
    print('Pollard_(p-1) start !')
    a=2
    for i in range(1,B):
        a=pow(a,i,n)
        p=math.gcd(a-1,n)
        if p!=1 and p!=n:
            print('p=',p,'\nPollard_(p-1) finish !')
            return p
    print('Pollard_(p-1) failure !')
    return

评价

​ 该算法实现的要求是B不能太大，否则和试除法无异，所以一般只适合 $ （p-1） $ 光滑的情况。

$Williams\ p+1$ 算法

​ 这是一个比较抽象的算法，适用于 $p+1$ 光滑的情况，需要用到数论的 $Lucas$ 序列。先介绍一下这个序列叭。
$$
U_n(P,Q)=PU_{n-1}(P,Q)-QU_{n-2}(P,Q)\\
V_n(P,Q)=PV_{n-1}(P,Q)-QV_{n-2}(P,Q)\\
U_0(P,Q)=0, \ U_1(P,Q)=1\\
V_0(P,Q)=2, \ V_1(P,Q)=P
$$
​ 看上去很简单（雾

​ 这个数列如果在有限域 $p$ 内进行，会得到一些很奇妙的结论（我不会证明反正）。这里拿出其中一个来，即
$$
若\ p\ 是奇素数，且\ p\nmid PQD,则\ p|U_{l}(P,Q)。其中\ D=P^2-4Q,\ l=p-(\frac{D}{p})\\ (\frac{D}{p})是勒让德符号。
$$
​ 我们于是利用这个性质分解 $n$ 。

​ 第一步，我们保证 $\frac{D}{p}=-1$ ，这样 $l=p+1$ 。根据前面的算法，我们是可以找到 $(p+1)|B!$ 这样的 $B$ 的。现在的问题是 $p|U_l$ ，那么 $p$ 是否整除 $U_{kl}$ ?

​ 我们烧为用特征方程证一下。
$$
数列的特征方程为 x^2-Px+Q=0\\
解得 a=\frac{P+\sqrt{D}}{2}\ b=\frac{P-\sqrt{D}}{2} 。\\通解为U(n)=\frac{a^n-b^n}{\sqrt{D}} 。\\
\frac{U_kl}{U_l}=\frac{a^{kl}-b^{kl}}{a^l-b^l}=\frac{(a^{l}-b^{l})(\sum_{i=0}^{i=n-1} a^ib^{n-i-1})}{a^l-b^l}
$$
​ 这个通解是不是和斐波那契数列很像？

​ 现在我们已经可以通过在一定范围内寻找一个 $B$ ，然后计算 $U_{B!}$ 的方法来分解 $n$ 了。最后一个问题是如何计算 $U_{B!}$ ?，显然一个一个计算太浪费时间了。于是我们有了矩阵快速幂，利用等式
$$
\left[\begin{matrix}
U_{i+1}\\
U_{i}
\end{matrix}\right]=
\left[\begin{matrix}
P&-Q\\
1&0
\end{matrix}\right]^i*
\left[\begin{matrix}
U_1\\
U_{0}
\end{matrix}\right]
$$

from Crypto.Util.number import *
import math
n=
def Williams(n,A_list):
    v=vector(Zmod(n),[0,1])
    for A in A_list:
        L=Matrix(Zmod(n),[[A,2],[1,0]])
        for j in range(3000):
            L=L^(j+1)
            V=int((L*v)[0])
            p=math.gcd(V,n)
            if(p!=1):
                q=n//p
                return p,q
A_list = [5, 7, 9, 11, 13]
p,q=Williams(n,A_list)
print(f"p={p}\nq={q}")

$Pollard\ rho$ 算法

原理

​ 这个就更简单了，咱就是说找两个数 $ x $ ， $ y $ ,也没什么限制，就随便找。如果运气好找出的两个数满足 $ x= y （mod \ p） $ ，那么恭喜，你可以计算 $ gcd(x-y,n) $ 来获取 $ p $ 。

实现

​ 根据上面的原理。很容易看出难点在于 如何寻找满足条件的数。 $ Pollard $ 给出了一种”随机”算法。首先确定一个函数 $ f $ ，一般取 $ f=x^21 $ ，这个函数在这里相当于随机数生成器。我们先通过一个实例来感受一下为什么这个算法叫做“$\rho$”算法。

实例

​ 令 $ n=71*101 $ ，起始值为 $ x_1=1 $ ， $ f=x^2+1 $ ，每次记录 $ x_i=f(x_{i-1})\ mod \ n $ ，得到前30个数的序列：

​ 1,2,5,26,677,6557,4105,6347,4903,2218,219,4936,4210,4560,4872,375,4377,4389,2016,5471…

​ 对它们模 $ 71 $ ，

​ 1,2,5,26,38,25| 58,28,4,17,6,37,21,16,44,20,46 $ | $ 58,28,4…

​ 我们发现第二条序列是一些无规则数字和一群循环数字的组合，很像希腊字母中的$\rho$，该算法由此得名。

​ 那么是不是任何n都可以这么做呢？

​ 首先我们可以保证，在有限域 $ n $ 中肯定可以找到两个数 $ x_i $ , $ x_j $ （比如 $ 1 $ 和 $ p+1 $ ），满足 $ x_i=x_j（mod\ p） $ , 我们设 $ j-i=l $ 。由于我们选择的函数是一个多项式函数，满足了 $ x_{i+1}=x_{j+1}（mod\ p） $ （不会证但是应该很显然），进而有 $ x_{i+s}=x_{j+s}（mod\ p） $ , $ 0<s<l $ （归纳法），也就证实了从 $ x_i $ 开始，所有的数构成一个长度为 $ l $ 且不断循环的圈。

​ 这个结论的好处是，只要我取的 $ x $ 足够多，最终它们基本都会落在环上，而环的长度是有限的（就是 $ j-x=l $ ），一旦落在环上的数大于 $ l $ ，那么就一定存在两个数满足 $ x_i=x_j（mod\ p） $ 。这时就只要计算 $ gcd（x_j-x_i，n） $ 即可完成分解。

​ 接下来的问题是该如何确定这样一个数对。如果我们假定其中一个不变，去寻找它的匹配数，由于我们是不知道数环的具体情况，也就不能判断我们所假设的数是不是在环上，如果我们假定的数不在数环上，那么就一定找不到它的匹配数，所以“一定一动”的方法不太合适。而如果认为两个数都是变化的，就需要对它们的轨迹进行一定的约束，否则完全随机的取数很可能会浪费太多时间。

​ 这时候我们需要学习一下小学数学中非常常见的追及问题。现在假设小明和小红（太对了）要从教学楼出发去操场跑步，其中小明的速度是 $ 1 $ ，小红是 $ 2 $ 。然后问题是小明和小红在一直能不能在操场的某处相遇。（这个模型其实就是把该算法实体化，从教学楼到操场的距离相当于$\rho$的尾巴，两人在操场跑步就是$\rho$的圈）这个问题的答案显然是肯定的，而且和教学楼到操场的距离无关，甚至和他们俩开始所在的位置无关。

​ pollard给出的算法也类似，两个动点一开始的分别是 $ x=x_0 $ ， $ y=x_1 $ 。其中第一个点每次映射一次而第二个点映射两次，这么做保证了它们首先不会在$\rho$的尾部相遇。而进入循环后，不管它们相对位置如何，最多再进行 $ l $ 次操作必定可以相遇。那么我们要做的就是在每一次映射后计算 $ gcd（x-y，n） $ 。

​ 这个算法有个问题，如果你运气很好导致公因子算出来恰好是 $ n $ ，那么算法就不能继续，要更换初始值来算。

代码

import math
def f(x):
    return x**2+1
#可能会陷入死循环哦
def Pollard_rho(s,n,f):
    print('Pollard_rho start !')
    x=s
    y=f(s)
    g=math.gcd(abs(y-x),n)
    while g==1:
        x=f(x)
        y=f(f(y))
        g=math.gcd(abs(x-y),n)
        if g==n:
            print("failure\nPollard_rho finish !")
            return
        elif g!=1 :
            print('p=',g,'\nPollard_rho finish !')
            return g

评论

​ 不太清楚为什么感觉实际操作起来很慢（？），不过思路确实好玩。

$Dixon$ 随机平方算法

原理

​ 这个原理就更更简单了，平方差公式嘛，小学生都会（润主要是用到了公式
$$
a^2\equiv b^2 \ mod\ n\ <=>\ (a-b)(a+b)\equiv\ 0\ mod\ n
$$
​ 在此基础上如果 $ （a-b） $ ， $ （a+b） $ 不是 $ 0 $ 或者 $ n $ 的倍数，那么 $ p $ ， $ q $ 就一定包含在 $ （a+b） $ 或者 $ （a-b） $ 当中，而且是每一个因子对应一个式子。

​ 这个很好证明，看等式 $ a^2-b^2=kpq=(a+b)*(a-b) $ 就可以了。 $ k $ 放哪无所谓， $ p $ ， $ q $ 已经是素数不能拆了，所以肯定是整个放在 $ （a+b） $ 或者 $ （a-b） $ 里面。

实现

​ 那么问题就只剩下怎么去找两个同余的平方数了。 $ Dixon $ 算法给出了一种利用素数基，也就是若干个小素数的集合（比如前5~15个素数）的方法。该算法说我们需要找到一些整数，使得它们的平方模n后是素数基中的元素的乘积。

​ 那么压力来到了如何去寻找符合条件的数。首先我们认识到，素数基中的数的乘积往往不是很大（至少相对于n来说比较小），这就给了我们一个思路，去寻找 $ kn+\varepsilon$ 的平方根，更确切地说，寻找 $ kn $ 的平方根附近的整数。于是我们可以依次求出n,2n…的平方根的近似值，然后去尝试能否分解为素数基中的数，如果成功，我们得到等式
$$
a^2_j\equiv\prod_{i=1}^{|B|} \ x^{k_i}_i\ mod\ n\ \ \ \ x_i是素数基B中的元素,a_j表示这是我们找到的第j个符合的整数,k_i是该元素的次数
$$
​ 接下来我们发现，现在等式左边已经是平方的形式了，但是右边还不一定，这就是我们为什么要找多个满足条件的数的原因。如果我们将得到的某几个等式相乘，可以使得右边也变成平方数的形式。

​ 为了使处理更加简洁，我们利用矩阵的方法，在GF(2)上定义一个 $1*|B|$ 的矩阵（ $|B|$ 是素数基中元素的个数），称之为幂矩阵。我们根据上式右边的 完善矩阵 $A_j$ ，即矩阵中第 $ i $ 列的值是 $ k_i mod 2 $ 。之后我们对一系列幂矩阵进行加减运算（相当于对对等式进行乘法运算）直到其中不含有 ​ 1 ，现在右边所有元素的幂都是偶数（模2为0），换句话说右边是一个平方数。

​ 根据线性代数的知识我们可以一步到位。如果我们已经得到了 $ s $ 个这样的等式，我们把 $ A_j $ 合并成一个$s*|B|$的矩阵 $ A $ ，求解方程 $ xA=0 $ 就是我们需要的组合系数（这样算出来一般会有很多解，但有些不是我们想要的）。

代码

#sage

B = []#素数基
n =
k=  #kn中的k值，影响最终获取的等式数目,此处取sqrt（kn）附近的两个整数
def dixon(B, n, k):
    print("Dixon start !")
    if -1 not in B:
        B.append(-1)
    A = []
    kn = []
    for i in range(1, k + 1):
        temp = [floor((i * n).sqrt()), ceil((i * n).sqrt())]
        for a in temp:
            a_ = a
            a = a**2 % n
            if a > (n / 2):
                a -= n
            mat = [0 for _ in range(len(B))] #mat相当于前文提到的幂矩阵
            f = a.factor()
            q = 1
            for x in f:
                if x[0] not in B:
                    q = 0
                else:
                    if a < 0:
                        mat[0] = 1
                    mat[B.index(x[0])] = x[1] % 2
            if q:
                A.append(mat)
                kn.append((a_, f))
    A = matrix(GF(2), A)
    x = A.left_kernel()
    #下面开始检验哪些x是我们需要的
    for xi in x:
        a2 = 1
        b2 = 1
        xi = list(xi)
        if 1 not in xi:
            continue
        b=[0 for _ in range(len(B))]
        for k in range(len(xi)):
            if xi[k]==1:
                a2=(a2*kn[k][0])
                for t in range(len(kn[k][1])):
                    b[B.index(kn[k][1][t][0])]+=kn[k][1][t][1]
        b2=1
        for i in range(len(b)):
            assert b[i]%2==0
            b2*=B[i]**(b[i]/2)
        p=gcd(a2%n+b2%n,n) 
        if p!=1 and p!=n:
            print("Dixon finish !")
            return p
    print("Dixon failure !")
    return
p=Dixon(B,n,k)

评价

​ 这个算法看起来很有趣但是代码实现有点困难，首先是没有合理的公式告诉我们素数基B应该取多大，其次是我们不知道需要多少个这样的等式。上面代码是自己理解的，实际操作后最多能快速分解50位左右n。（我的评价是不如yafu）

好了，那我们现在正式开始这篇文章，这篇文章的内容是 大数分解网站（别老惦记着你那yafu了，在线网址不好吗