高效因子分解

本文主要包括了一些常见和不常见的高效因子分解算法，重点在于椭圆曲线法的讲解（虽然Dixon’s 的篇幅更多一些貌似）

主要参考为 Wikipedia 的 Integer Factorization 页面

Warning:多项内容需修正

前言

质因数分解一直是一个密码学领域相当热门的话题，因为他直接关系着 RSA 等使用大素数相乘来加密的算法的安全性

相比于

• 求素数（筛法或者特殊的生成法）

• 素性测验(能做到$O((\log n))^{6+\epsilon}$)

• 素数计数（求$[1,n]$的素数个数，可以使用Meissel-Lehmer在亚线性复杂度内求解）

因数分解一直是一个很困难的内容

现有的对于一般数字分解的最好算法是GNFS（General Number Field Sieve，一般数域筛）的复杂度为$L_n[\frac{1}{3},\sqrt[3]{\frac{64}{9}}]=\exp((\sqrt[3]{\frac{64}{9}}+o(1)(\ln n)^\frac{1}{3}(\ln\ln n)^{\frac{2}{3}})$

能够分解上百位的合数

而数域筛涉及到大量的代数数论的内容，故笔者在此介绍一些更加简单易懂的算法包括 ：

Dixon’s ，CFRAC 和 Lanstra 椭圆曲线 算法

大概只需要一些基本的数论知识（和一小点群论），同时给出了一个椭圆曲线法的Python Demo

基础的想法

试除

一个基础的想法是在$[2,\sqrt n]$的范围内试除，复杂度为$O(\sqrt n)$，虽然在 OI 里很多情况下都已经够用，但是对于分解数十甚至上百位整数的需求就难以接受

Fermat

而 Fermat 分解提供了一个新的思路

注意到若奇合数 $n=pq$ 那一定可以表示为两数平方差 $a^2-b^2$ 仅需取 $a=\frac{p+q}{2} \quad b=\frac{p-q}{2}$

变形一下得到 $a^2-n=b^2$ 那我们枚举 $a$ 测试 $a^2-n$ 是否为完全平方数即可

但是这样的复杂度为 $O(n)$ 甚至更慢，但是能够对更快的的质因数分解做出一定的启示

比如我们考虑如果放宽要求 $a^2\equiv b^2 \pmod n$ 那么$\gcd(a - b,n)$ 就很可能是一个非平凡因子

因为有$n|(a-b)(a+b)$，所以$\gcd(a-b,n) ,\gcd(a+b,n)$一定有一个因子 (只要不是 $a=\pm b$)

Pollard’s Rho算法

该算法是 OI 中不常用的算法之一，但是是一个极好的高效因数分解的例子

Pollard’s Rho算法的思想是在接近随机的数列$\{x_i\}:x_{i+1}= (x_i^2+c) \bmod n$

出现$x_a\equiv x_b\pmod n$的期望步数为$O(\sqrt n)$ （根据生日悖论）

而对于$n=pq$ 的 $\{x_i \bmod p\}$ 依然是接近随机行为的数列，同时期望为$O(\sqrt p)$

于是只要$x_a\ne x_b$计算$\gcd(x_a-x_b,n)$就得到了一个$n$的非平凡因子

于是我们得到了$O(\sqrt p)$或者说$O(n^{\frac{1}{4}})$的优秀算法

亚指数级

$O(n^{\frac{1}{4}})$已经很快的，在 OI 中大概也是极限了（再大建议物理解决出题人）

但是它仍然不！够！快！

下面是 Linux 下 factor 的速度（内部实现为 Pollard’s Rho）

time factor 60000000000000003100000000000000033
60000000000000003100000000000000033: 200000000000000003 300000000000000011
factor 60000000000000003100000000000000033  38.38s user 0.00s system 99% cpu 38.383 total

可以看到分解仅 $35$ 位的合数就用了 $38$ s 这显然是不能接受的 （实际上 Linux 的跑的相当快,我写的丢人 Pollard’s Rho 跑了两分钟还没跑出来

而我们完全可以用 Fermat 的想法做到更快的分解速度

Dixon’s 与 连分数 (CFRAC)

为了方便的叙述先定义两个概念

因子基

一个因子基意味前几个素数的集合$\{p_1,p_2,p_3,\dots\}$

L记号

L 记号同样是衡量时间复杂度的，$L_n[a,c]=\exp\{(c+o(1)(\ln n)^a(\ln \ln n)^{1-a})\}$

其中 $c > 0,a \in [0,1]$

主要用来数论中算法的复杂度

Dixon’s

Dixon’s 的主要思想是：

选定一个大小为 $b$ 的因子基 $B$，选取一些数字$x_i$，使得 $x_i^2\bmod n=\prod p_i^{k_i}$ 的所有 $p_i \in B$ （显然不会所有的$k_i$都为偶数，否则就找到了一个$a^2\equiv b^2$）

然后通过奥妙重重的方法选取一系列$x_i$使得 $\prod x_i \equiv (\prod p_i^{k_i})^2 \pmod n$

这样我们就得到了一个$a^2 \equiv b^2 \pmod n$

那我们怎么才能寻找到这个奥妙重重的方法呢？

考虑 $\prod x_i = \prod p_i^{a_{i,1}+a_{i,2}+a_{i,3}+\dots} \pmod n$

令 $\vec{a_i}=(a_{i,1} \bmod 2,a_{i,2} \bmod 2, \dots ,a_{i,k} \bmod 2) \in(\mathbb{Z}_{2})^k$

可以用高斯消元解出来向量 $\vec{a}=(0,0,0,\dots,0) \in (\mathbb{Z}_2)^k$

把对应的数字乘起来即可

而 $x_i$ 的选取则是完全随机的（

大概分析下复杂度

对于确定的 $b=n^{-\alpha}$ 而言 找到一个合适的 $x$ 的概率大概是是 $a^{-\alpha}$ （这实际上Dickman function $\rho$ 的一个近似，之后我们会在椭圆曲线那里讨论）

而对于这个 $b$ 我们仅需找到 $\pi(b)$ 个 $x$ 即可总共随机约 $\exp\{(\alpha+1)\ln \alpha +\frac{\ln n}{\alpha} -\ln\ln n\}$ 次数 每次需要 $\pi(b)$ 次试除

取 $\alpha=2\sqrt{\frac{\ln n}{\ln \ln n}}$ 得到复杂度约为$L_n[\frac{1}{2},2\sqrt2]$

连分数(CFRAC)

笔者并没有系统的学习过连分数相关的理论，故仅在此给出简单的介绍

连分数因数分解实际上在 Dixon’s 上的一个优化

我们发现如果 $x^2 \bmod n$ 比较小的话就容易分解

那么用周期连分数表示 $\sqrt {kn}$ 来得到比较好的精度 这样可以做到 $L_n[\frac{1}{2}, \sqrt2]$

(抄的Wikipedia，如果有懂的大佬希望可以给出相关文献QAQ)

二次筛

二次筛法同样是对 Dixon’s 的优化，可以做到 $L_n[\frac{1}{2},1]$ 在 GNFS 之前一直是最快的因子分解算法

具体原理在学了，学会了也许会单独写一篇

椭圆曲线法

虽然椭圆曲线法同样是亚指数级别的，但是因为要给出比较详细的讲解故将其分离了出来

一些定义

椭圆曲线

椭圆曲线是域上亏格为1的光滑射影曲线

我们可以理解为形如 $y^2=x^3+ax+b$ 的曲线

B-smooth

当我们说 $n$ 是 B-smooth 的指的是 $n$ 所有的质因子都小于等于 $B$

B-powersmooth

同样的，B-powersmooth 指的是 $n$ 的任意因子的次幂不超过 $B$ 例如 $2^4 \times 3^2$ 是 16-powersmooth的

椭圆曲线群

我们在椭圆曲线上的点 $P,Q$ 定义加法为 $P+Q=T$，$T$是椭圆曲线上 $P,Q$ 两点连线的第三个交点的关于 $x$ 轴的对称点（这样定义是为了保持群的结构和交换）

定义自己加自己是沿着自己的切线的交点

再定义无穷远点记作$O$，容易验证椭圆曲线上的点加上点 $O$ 和加法构成结构是一个阿贝尔群，点 $O$ 是其幺元

Pollard’s p-1

虽然 Pollard’s p-1 算法甚至不是正确的（不一定能得到答案）

但是它仍然有着借鉴意义

Pollard’s p-1 算法的思想是如果 $n=pq$

那么 $a^{k(p-1)} \equiv 1 \pmod p$ （费马小定理）

如果 $a$ 的 某个幂次 $k$ 是模 $n$ 的因子同余于 $1$ 那么 $\gcd(a^k-1,n)$ 大概率是一个 $n$ 的非平凡因子

我们选取一个 B-smooth 的数字 $M$ ，然后计算 $\gcd(x^M,n)$ 重复多次以寻找因子

而 $x$ 需要选取互质于 $n$ 的数字 对于奇合数来说选 $2$ 即可

椭圆曲线法

我们将椭圆曲线限制在有理点上，可以发现群的结构仍然保持

我们可以考虑有限域$\mathbb{F}_q$而非有理域，因为我们不需要有理数的结构

考虑椭圆曲线上的加法过程

令 $s=\frac{3x^2+a}{2y}$

则 $2(x,y)=(s^2-2x,y+s(u-x))$

如果我们在 $E(\mathbb{Z}/n\mathbb{Z})$ 下出现了无法除 $2y$ 的情况 ，那么我们在模 $n$ 意义下出现了不存在的逆元，也就是 $\gcd(n,t)\neq1$ 我们就大概率找到了一个因子

流程

1. 选定一条曲线 $y^2=x^3+ax+b$ 与非平凡点 $P$

2. 选定素数上界 $B$

3. 计算 $(\prod\limits_{p_i<B}p_i^{\ln B/\ln p_i})P$ （$\ln B/\ln p_i$ 是为了保证这个数 B-powersmooth）

4. 加法过程中寻找因子

算法原理

可以发现我们的基本思想是和 Pollard’s p-1 相同，只不过把乘法改成了在椭圆曲线群上的加法

所以它为什么比Pollard’s p-1 要好呢？

Pollard’s p-1 的本质在于在群 $\mathbb{Z}_q \times \mathbb{Z}_p$ 中找到了一个 $\mathbb{Z}_p$ 的幺元 $e_p$ 但是非 $\mathbb{Z}_q$ 的幺元

而这个群的阶 $p$ 不一定是一个 B-powersmooth 的数

而根据 Hasse 定理，椭圆曲线群 $|E(\mathbb{Z}/p\mathbb{Z})|\in[p+1-2\sqrt p,p+1+2\sqrt p]$

这样群的阶就会有变化，避免了不 powersmooth 的问题

进一步我们考虑曲线 $y^2=x^3+ax+b$ 分别在模 $p$ $q$ 意义下的群 $E(\mathbb{Z}/p\mathbb{Z})$ 与 $E(\mathbb{Z}/q\mathbb{Z})$ （假设 $n=pq$）

记他们的元素数量为 $N_q$ $N_p$

根据拉格朗日定理（$|G|=[G:H]\times |H|$）若 $kP=O \pmod p$

那么$N_pP=O \pmod p$

于是和 Pollard’s p-1 一样，如果我们找到了某个 $k$ 使得 $kP=O \pmod p$ 但不是$E(\mathbb{Z}/q\mathbb{Z})$ 中 $P$ 的阶，此时我们就找到一个因子

复杂度

找到一个可分解的椭圆曲线群的概率为 $B\ln^2n$

若设$B=p^{-a}$ 终止的概率为 $a^{-a}$

和 Dixon’s 一样的分析能得到复杂度大致为 $L_p[\frac{1}{2},\sqrt2]$ $p$ 为 $n$ 的最小素因子

可以看到椭圆曲线法的复杂度主要取决于 $n$ 的最小素因子大小而非 $n$

所以椭圆曲线法可以快速去除一些小的素因子（大约 $30 \sim 40$ 位的素因子）

具体实现

操作中选取曲线上的非平凡点这个行为可以变为随机一个点然后生成一个曲线

对于点 $(x,y)$ 选取一个 $a$ 计算得到 $b=y^3-x^3-ax$

Demo

from random import randint
from math import gcd

def EulerSieve(n):
    vis = [True] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if vis[i]:
            primes.append(i)
        for j in primes:
            if i * j > n:
                break
            vis[i * j] = False
            if i % j == 0:
                break
    return primes

def exgcd(a,b):
    if b == 0:
        return 1, 0, a
    q, r = divmod(a, b)
    x, y, d = exgcd(b, r)
    return y, x - q * y, d

def elliptic_add(p, q, a, b, m):
    if p[2] == 0: 
        return q
    if q[2] == 0: 
        return p
    if p[0] == q[0]:
        if (p[1] + q[1]) % m == 0:
            return 0, 1, 0
        num = (3 * p[0] * p[0] + a) % m
        denom = (2 * p[1]) % m
    else:
        num = (q[1] - p[1]) % m
        denom = (q[0] - p[0]) % m
    inv, _, g = modular_inv(denom, m)
    if g > 1:
        return 0, 0, denom  
    z = (num * inv * num * inv - p[0] - q[0]) % m
    return z, (num * inv * (p[0] - z) - p[1]) % m, 1

def elliptic_mul(k, p, a, b, m):
    res = (0, 1, 0)
    while k > 0:
        if p[2] > 1:
            return p
        if k % 2 == 1:
            res = elliptic_add(p, res, a, b, m)
        k = k // 2
        p = elliptic_add(p, p, a, b, m)
    return res

def lenstra(n, limit):
    g = n
    while g == n:
        q = randint(0, n - 1), randint(0, n - 1), 1
        a = randint(0, n - 1)
        b = (q[1] * q[1] - q[0] * q[0] * q[0] - a * q[0]) % n
        g = gcd(4 * a * a * a + 27 * b * b, n) 
    if g > 1:
        return g
    primes=EulerSieve(limit)
    for p in primes:
        pp = p
        while pp < limit:
            q = elliptic_mul(p, q, a, b, n)
            if q[2] > 1:
                return gcd(q[2], n)
            pp = p * pp
    return False

def main():
    n=998244354
    B=7000
    curve_count=200
    for i in range(curve_count):
        res=lenstra(n, B)
        if res != False:
            print(res)
            break
if __name__ == '__main__':
    main()

一些优化

B的选定

如果定义 $\Phi(x,y)$ 为 $\le x$ 的 y-smooth 的数的个数，那么 $\Phi(x,x^{-a}) \sim x\rho(a)$

其中 $\rho$ 是之前提到过的 Dickman Function ，定义为 $x\rho'(x)+\rho(x-1)=0$

我们希望最小化 $\Phi(p,B)^{-1}\pi(B)$ 的值

如果采用近似 $\Phi(n,n^{-a}) \sim n\rho(a) \quad \pi(B) \sim \mathrm{Li}(B)$

可以用牛顿迭代一类的方法找极值，当然通常使用经验对 $30$ 位左右的数字用 $7000$ 的 B 即可

多曲线

我们可以一次性生成多条曲线，一起在这多条曲线上做加法，这样只要有一个曲线上找到了因子就可以退出

求逆元的时候也可以先求一次然后然后用乘法来算其他逆元，你一个个求要快

蒙哥马利取模

Min_25 的博客有过介绍

核心步骤Reduce: （ $x \in [0,rn)$ 返回 $xn^{-1} \bmod r$）

def reduce(x):
    m=((x%n)*(inv(r)))%n
    x=(x-m*r)//n
    if x>=r:
        return x-r
    return x

究其本质大概是在 $n,r$ 互质的情况下映射 $f(a)=an$ 是 $\mathbb{F}_r$ 的自同态？

大量乘法取模的时候实际效率很高

C++实现

太菜了写不动这个东西，给一份 zball 大佬的实现，效率极高

Here

0.20788920s consumed
200000000000000003 300000000000000011

比 factor 快上百倍、

同时 wikipedia 上也列举了许多 ecm 的高效实现，可以在 Wikipedia中查看（需独特的上网技巧）

Asusetic eru quionours