素性检验算法系列（三）：椭圆曲线素性检验

系列目录：

• 【素性检验算法系列（一）：试除、Fermat小定理、Lucas序列】
• 【素性检验算法系列（二）：基于因数分解的群阶数检验】
• 【素性检验算法系列（三）：椭圆曲线素性检验】
• 【素性检验算法系列（四）：高次扩域检验、APR算法】
• 【素性检验算法系列（五）：Frobenius检验、AKS算法】

在前两篇文章的介绍中，我们已经发现了除试除法之外的所有素性检验方法的共同本质：验证与 $N$ 有关的某个群的阶数 。例如，Fermat、Miller-Rabin、Solovay-Strassen、Lucas、Polington等方法都是在验证整数模群 $(\mathbb{Z}/N\mathbb{Z}{)}^{\ast }$ 的阶数等于 $N-1$ ，而Lucas拟素数检验及其因数分解版本、包括其变体与推广方法则是在验证二次域整数模群的商群 $((\mathbb{Z}/N\mathbb{Z})[\sqrt{D}]{)}^{\ast }/(\mathbb{Z}/N\mathbb{Z}{)}^{\ast }$ 的阶数等于 $N+1$ . 本篇及后面要介绍的素性检验算法将依然遵循这一本质。

在通过验证群阶数进行素性检验时，我们一共要搞清楚三件事情： （1）什么群？阶数是多少？（2）群元素是什么？怎么取群元素？怎么做群的运算？（3）对群阶数进行因数分解的界限是什么？ 在前两篇文章中，我们针对整数模群和二次域商群这两种情形，都比较清晰地回答了这三个问题。在本篇文章及后面的介绍中，我们将重点回答前两个问题；第三个问题通常涉及繁琐的推导过程，我们就只描述结论，将具体证明留给文献和参考资料。

本篇文章中，我们将介绍基于椭圆曲线群的素性检验方法。

本文为所有算法提供了对应的python 3代码，放在了github上：【GitHub - pyprimality】.

四、基于椭圆曲线的方法

我们介绍一种新的代数结构——有限域上的椭圆曲线群（elliptic curve group）。通过检验椭圆曲线群的阶数来进行素性检验的方法叫做椭圆曲线素性检验（elliptic curve primality proving，以下简称ECPP）。

先来看实数域 $\mathbb{R}$ 上的椭圆曲线群。首先，定义如下集合为实数域上的椭圆曲线：

$$E(\mathbb{R})=\{(x,y)|x,y,a,b\in \mathbb{R};y^2=x^3+ax+b;4a^3+27b^2\ne 0\}$$

它包含 $a,b$ 两个参数，且要求 $4a^3+27b^2\ne 0$ ，这是为了保证曲线上没有不可导的奇异点。在不同的 $a,b$ 下，曲线在平面直角坐标系中的图像有所差别，但永远关于 $x$ 轴对称。

在特定的加法下，曲线 $E(\mathbb{R})$ 上的所有点与无穷远点可以构成一个阿贝尔群，称为椭圆曲线群。群中的加法运算法则如下：

• 规定无穷远点为加法单位元 $0$ ；
• 非无穷远点 $P(x,y)\in E$ 的逆元是其关于 $x$ 轴对称的点，即 $-P(x,-y)$ ；
• 二元加法运算(1)：对于两个不同的非无穷远点 $P,Q\in E$ ，直线 $PQ$ 与椭圆曲线交于第三个非无穷远点 $R$ ，则有 $P+Q+R=0$ ，即 $P+Q=-R$ ；
• 二元加法运算(2)：对于非无穷远点 $P\in E$ ，椭圆曲线在 $P$ 处的切线与曲线交于另一个非无穷远点 $Q$ ，则有 $P+P=-Q$ ，也可记为 $2P=-Q$ ；
• 当 $P+Q+R=0$ 时三点共线，由韦达定理可得 $x_P+x_Q+x_R=k^2$ ，其中 $k$ 为直线斜率；由这一等式，可在已知 $P,Q$ 坐标时计算出 $R$ 的坐标；其中，斜率 $k$ 可根据 $P,Q$ 坐标算出，若 $P,Q$ 重合，则斜率为切线斜率 $k=(3x_P^2+a)/(2y_P)$ .

以上所定义的加法满足结合律和交换律。

椭圆曲线还可以被限定在一个比实数更小的域内，如有限域。令 $p$ 为素数，则模 $p$ 整数可以构成一个有限域 ${\mathbb{F}}_p=\mathbb{Z}/p\mathbb{Z}$ ，在该域上也可以定义椭圆曲线：

$$E(\mathbb{Z}/p\mathbb{Z})=\{(x,y)|x,y,a,b\in \mathbb{Z}/p\mathbb{Z};y^2=x^3+ax+b(\mathrm{mod}p);4a^3+27b^2\ne 0(\mathrm{mod}p)\}$$

此时的椭圆曲线不再是一条连续的线，而是一系列离散的有限数量的点。这些点与无穷远点依然可以构成一个阿贝尔群，其加法运算法则与实数域椭圆曲线的情形基本相同，但具体细节有如下改变：

• 非无穷远点 $P(x,y)$ 的逆元的坐标须对 $p$ 取模，即 $-P(x,-y(\mathrm{mod}p))$ ，因此实际效果为 $P$ 与 $-P$ 关于 $y=p/2$ 对称；
• 同一直线上的三个点累加之和为 $0$ ，直线的定义须改为 $\alpha x+\beta y+\gamma =0(\mathrm{mod}p)$ ，直线斜率仍可根据 $P,Q$ 坐标计算，计算过程中的除法定义为模 $p$ 乘法的逆运算；
  计算 $P+P=2P$ 需要寻找椭圆曲线在 $P$ 处的切线，而离散点集无法定义切线；但数学证明表明，当 $p\ne 2,3$ 时，离散点集的“切线”斜率表达式与连续曲线切线斜率一致，对 $p$ 取模即可；
• 当 $P+Q+R=0$ 时三点共线，韦达定理表达式为 $x_P+x_Q+x_R=k^2(\mathrm{mod}p)$ .

给定待检验的正整数 $N$ ，若 $N$ 是素数，则 $E(\mathbb{Z}/N\mathbb{Z})$ 是有限阶数的群。通过对该群的阶数进行检验，就可以对 $N$ 的素性进行判断。这就是ECPP的基本思路：

定理4.1（ECPP）：对于正整数 $N$ ，令 $E$ 为定义在 $\mathbb{Z}/N\mathbb{Z}$ 上的椭圆曲线 $y^2=x^3+ax+b(\mathrm{mod}N)$ ，曲线上的点的个数（包含无穷远点）为 $m$ . 若 $m$ 存在因数 $s>{\left(N^{1/4}+1\right)}^2$ ，且存在 $E$ 上的一点 $P$ ，使得以下条件成立：
（1） $mP=0$ ；
（2）对于 $s$ 的任意素因数 $q$ 都有 $(m/q)P\ne 0$ ；
则 $N$ 是素数。

上述定理的具体证明可参考英文维基百科词条【Elliptic curve primality】。如果 $N$ 是素数，则 $E$ 可构成群，$m$ 是群的阶数，所以一定有 $mP=0$ （ $mP$ 表示 $m$ 个 $P$ 按照群加法累加的结果）；另一方面，$(m/q)P\ne 0$ 则说明群元素 $P$ 的阶数不是小于 $m$ 的因数。这两个条件一起保证了椭圆曲线群的阶数一定是 $m$ . 从这一角度理解，定理4.1与Pocklington素性检验（定理3.4）最为接近；而椭圆曲线群的优点则在于其阶数会随着曲线而变化，并非固定的 $N-1$ 或 $N+1$ ，因此，如果符合条件的因数分解不存在或难以找到，则可以重新寻找一条椭圆曲线，这样就避免了因数分解的难度受限于 $N$ .

11. Goldwasser-Kilian算法

在定理4.1的理论基础上，第一个可执行的ECPP算法被提出，称为Goldwasser-Kilian算法。其主要思路为：先随机生成椭圆曲线，再计算其阶数，然后根据 $mP$ 和 $(m/q)P$ 的计算结果检验群的阶数。如果在计算过程中出现了无意义/无定义的表达式，或出现其他不符合椭圆曲线群性质的情况，则说明 $E$ 无法构成群，则 $N$ 是合数。

算法4.2（Goldwasser-Kilian算法）：输入正整数 $N>3$ ；
【步骤1】随机选择整数 $a,x,y\in [0,N-1]$ ，令 $b=y^2-x^3-ax\mathrm{mod}N$ ，并要保证 $4a^3+27b^2\ne 0(\mathrm{mod}N)$ ，由此得到椭圆曲线 $E$ 和曲线上的点 $P(x,y)$ ；
【步骤2】利用特定算法计算 $\mathbb{Z}/N\mathbb{Z}$ 上的椭圆曲线 $E$ 的点的个数（推荐使用Schoof算法，此处不介绍）；若计算过程中出现了无意义/无定义表达式或其他不符合椭圆曲线群性质的情况，则 $N$ 是合数；
【步骤3】将椭圆曲线上点的个数（包含无穷远点）记为 $m$ ，寻找 $m$ 的一个因数 $s>{\left(N^{1/4}+1\right)}^2$ ，并对 $s$ 进行完全的素因数分解（如果无法确定每个因数都是素数，则也可以分解为拟素数，即通过某些方式认为可能是素数的数）；
【步骤4】计算 $mP$ 和所有 $(m/q_i)P$ ，其中 $q_i$ 是 $s$ 的所有（拟）素因数；若计算过程中出现了无意义/无定义表达式或其他不符合椭圆曲线群性质的情况，则 $N$ 是合数；
【步骤5】若 $mP\ne 0$ ，则 $N$ 是合数；
【步骤6】若存在 $(m/q_i)P=0$ ，则从步骤1重新开始；
【步骤7】若 $mP=0$ 且所有 $(m/q_i)P\ne 0$ ，则在假定所有 $q_i$ 都是素数的情况下，$N$ 是素数；如果存在某个/某些 $q_i$ 无法准确判定为素数，则需要对其进行进一步的素性检验，可递归使用本算法，也可使用其他方法；若所有 $q_i$ 都确定是素数，则 $N$ 是素数；若某个 $q_i$ 确定是合数，则从步骤1重新开始；
【步骤8】重复足够多次步骤1后，以足够高的概率相信 $N$ 是合数。

该算法复杂度期望是 $\mathcal{O}({\log }^{5+ϵ}N)$（注意这是期望，实际上有一定随机性）。它是第一个多项式时间的（相对于 $\log N$ ）、适用于任意正整数的、不产生伪素数的素性检验算法（三个buff），也是目前最广泛使用的快速素性检验算法之一。它可以给被检验的素数生成对应的素性证书（certificate），通过这个证书可以更快速地验证这个数的素性，而不需要再经过一遍检验。ECPP不会产生伪素数，但由于其具有一定的随机性，如果输入的是素数，那么理论上会有一点点非常微小的概率判断不出素性（即步骤8），不过在实用中可以忽略这种情况。

Goldwasser-Kilian算法复杂度的最主要来源是步骤2中用来计算椭圆曲线群阶数的Schoof算法。尽管Schoof算法是一个多项式时间的算法，但其过程极其繁琐（也就是说多项式前面的常数很大）。在实际应用中，对于较大的 $N$ ，算法也会非常耗时。

补充：关于Schoof算法的具体原理和步骤，可以参考原始文献【Schoof, R. 1985】，也以可参考一篇硕士论文里的综述部分【高付强. 2010】，其中后者对算法的描述更为具体且清晰（也可能是因为中文比较好读）。但是需要注意，这两篇文章中的公式都有很多错误！即使是原始文献也有错误！错误种类包括但不限于幂次错误、缺少符号、奇偶情况混淆、符号打印错误、括号位置错误、定义冲突等。如果要参阅这两篇文献，请务必自己推一遍公式！！！如果懒得推又想知道正确结果，请直接参考笔者推导并总结好的笔记【有限域椭圆曲线求阶的Schoof算法】 ~ 另外，Schoof算法默认所输入的有限域特征是素数，如果要将Schoof算法用在素性检验中，输入的整数不一定是素数，则可能会在运算过程中出现错误或无意义表达式，此时就说明输入的整数其实是合数。

对算法进行测试后发现，当 $N$ 很大时，Schoof算法所需要的时间要远大于步骤3的因数分解，因此Schoof算法是限制Goldwasser-Kilian算法效率的最关键因素。除此之外还有另一个问题——原则上说这是一个确定性算法，即如果 $N$ 是素数则一定存在满足定理4.1的曲线以及曲线上的点，但没有理论保证算法能在可接受的重复次数内找到合适的曲线和点，因此在较坏的情况下，算法所需要的尝试次数是无法估计的。

再补充一句：Goldwasser-Kilian算法在做因数分解时使用了拟素因数分解+递归的策略，然而这种策略并不能用于上一篇文章中的Lucas、Pocklington等素性检验方法，因为这些方法所分解的 $N\pm 1$ 是固定的，如果 $N\pm 1$ 的某个拟素因数被判定为合数，我们除了继续分解它以外没有任何办法。相对而言，ECPP所分解的阶数 $m$ 则是可以变化的，因此如果某个拟素因数被判定为合数，只要重新生成一个 $m$ 就可以。

12. Atkin-Morain算法

为了解决Goldwasser-Kilian算法中的Schoof算法太慢的问题，另一种ECPP算法被提出，称为Atkin-Morain算法。该算法利用了一种特殊的椭圆曲线——具有复乘（complex multiplication）的椭圆曲线，这种椭圆曲线可以根据指定的阶数构造出来（构造方法详见文献【IEEE Std 1363-2000】附录A14）。因此Atkin-Morain算法是先给出一个阶数，再构造指定阶数的椭圆曲线，这个顺序与Goldwasser-Kilian算法是相反的，省去了繁琐的Schoof算法，从而更加高效。

算法4.3（Atkin-Morain算法）：输入正整数 $N>3$ ；
【步骤1】在一定范围内选择无平方因子的正整数 $D$ ，使得方程 $4N=u^2+D{v}^2$ 有整数解 $(u,v)$ ；
【步骤2】构造阶数为 $m=N+1\pm u$ 的椭圆曲线 $E$（方法详见文献【IEEE Std 1363-2000】附录A14），随机选择曲线上的一点 $P$ ；
【步骤3】检验 $m$ 取何值时 $mP=0$ ；若计算 $mP$ 的过程中出现了无意义/无定义表达式或其他不符合椭圆曲线群性质的情况，或者对于所有可能的 $m$ 都有 $mP\ne 0$ ，则 $N$ 是合数；
【步骤4】对于满足 $mP=0$ 的 $m$ ，寻找它的一个因数 $s>{\left(N^{1/4}+1\right)}^2$ ，并对 $s$ 进行完全的素因数分解（或分解为拟素数）；
【步骤5】计算所有 $(m/q_i)P$ ，其中 $q_i$ 是 $s$ 的所有（拟）素因数；若计算过程中出现了无意义/无定义表达式或其他不符合椭圆曲线群性质的情况，则 $N$ 是合数；
【步骤6】若存在 $(m/q_i)P=0$ ，则从步骤1重新开始（或者也可以回到步骤2，在同一条曲线上选择另外一个点 $P$ ）；
【步骤7】若 $mP=0$ 且所有 $(m/q_i)P\ne 0$ ，则在假定所有 $q_i$ 都是素数的情况下，$N$ 是素数；如果存在某个/某些 $q_i$ 无法准确判定为素数，则需要对其进行进一步的素性检验，可递归使用本算法，也可使用其他方法；若所有 $q_i$ 都确定是素数，则 $N$ 是素数；若某个 $q_i$ 确定是合数，则从步骤1重新开始；
【步骤8】重复足够多次检验后，以足够高的概率相信 $N$ 是合数。

由于Atkin-Morain算法不需要利用Schoof算法求阶，因此总体运行时间要远少于Goldwasser-Kilian算法。Atkin-Morain算法运行过程中的绝大部分时间都消耗在了步骤4的因数分解上。

针对Goldwasser-Kilian和Atkin-Morain两种算法的代码进行了一些测试，统计了不同输入数据下两种算法各自进行单次检验所需要的时间，得到的结果如下表（“单次检验所需要的时间”指的是算法从步骤1开始到步骤6结束所用的时间，不包括重复检验和递归的时间，也不代表算法给出最终结果所需要的时间）：

通过对比能看出Atkin-Morain算法具有相对而言更高的效率。对于同一个 $N$ 重复进行测试，两种算法的运行时间都有波动，这是由于算法产生的椭圆曲线是随机的，对阶数进行因数分解的难度也不同。此外在测试中还发现一个规律：如果输入的是合数，那么Goldwasser-Kilian算法基本上在步骤2（Schoof算法）就能判断出来，而Atkin-Morain算法基本上在步骤1（求解二元二次方程部分）就能判断出来，且后者做出判断所需要的时间更短，通常可忽略不计。

13. Tsumura算法（针对特定数）

对于具有特定形式的数字，存在更高效的ECPP算法。这里我们研究 $N=2^{k}n-1$ 形式的数字，针对该形式数字的ECPP算法由Yu Tsumura所提出（此处称之为Tsumura算法，原始文献为【Yu Tsumura, 2009】）。

假设有一个素数 $p=2^{k}n-1$ ，其中 $n$ 是奇数且 $k\ge 2$ . 在有限域 ${\mathbb{F}}_p=\mathbb{Z}/p\mathbb{Z}$ 上定义一条特殊的椭圆曲线：$E:y^2=x^3-mx(\mathrm{mod}p)$ ，其中 $m\ne 0(\mathrm{mod}p)$ . 通过数学证明，可以得到如下结论（证明过程见原始文献）：

• 椭圆曲线群 $E({\mathbb{F}}_p)$ 的阶数等于 $p+1=2^{k}n$ ；
• 若 $m$ 是 $p$ 的二次非剩余，则 $E({\mathbb{F}}_p)$ 同构于 ${\mathbb{Z}}_{2^{k}n}$ ，否则 $E({\mathbb{F}}_p)$ 同构于 ${\mathbb{Z}}_2\oplus {\mathbb{Z}}_{2^{k-1}n}$ （符号 ${\mathbb{Z}}_d$ 表示 $d$ 阶循环群）；
• 若 $m$ 是 $p$ 的二次非剩余，$P(x,y)$ 是曲线 $E$ 上的一点，且 $x$ 也是 $p$ 的二次非剩余，则 $P$ 在群中的阶数一定是 $2^k$ 的倍数。

Tsumura算法是以这些结论为理论基础而构建的。算法分为两种情形：$n$ 较小情形（ $n<N/{\left(N^{1/4}+1\right)}^2$ ）和 $n$ 较大情形（ $2^k<N/{\left(N^{1/4}+1\right)}^2$ ）。对于 $n$ 较小的情形，有如下定理：

定理4.4： 正整数 $N=2^{k}n-1$ ，其中 $n$ 是奇数且 $k\ge 2$ ，且有 $n<N/{\left(N^{1/4}+1\right)}^2$ . 令 $E$ 为定义在 $\mathbb{Z}/N\mathbb{Z}$ 上的椭圆曲线 $y^2=x^3-mx(\mathrm{mod}N)$ ，其中 $m$ 是 $N$ 的二次非剩余。则 $N$ 是素数的充要条件为：
（1）曲线 $E$ 上存在一点 $P(x,y)$ ，使得 $S_k=0(\mathrm{mod}N)$ ，且
（2）对于所有 $1\le i\le k-1$ 都有 $\text{gcd}(S_i,N)=1$ .
其中，数列 $S_i(i>0)$ 定义为：令数列 $X_0=x$ 且 $i>0$ 时 $X_i=\frac{{\left(X_{i-1}^2+m\right)}^2}{4\left(X_{i-1}^3-m{X}_{i-1}\right)}\mathrm{mod}N$ ，则 $S_i=4\left(X_{i-1}^3-m{X}_{i-1}\right)$ .

定理的具体证明见原始文献。其中数列 $S_i$ 的含义可以这样理解：设点 $P$ 的坐标为 $(x,y)$ ，则根据椭圆曲线群的点加法法则，可以计算出 $2P$ 的横坐标为 $\frac{{\left(x^2+m\right)}^2}{4\left(x^3-mx\right)}\mathrm{mod}N$ ，因此数列 $X_i$ 就是点 $2^{i}P$ 的横坐标；而 $S_i$ 则是 $X_i$ 的分母，$S_k=0(\mathrm{mod}N)$ 意味着点 $2^{k}P$ 是无穷远点，即 $2^{k}P=0$ ；$\text{gcd}(S_i,N)=1$ 则意味着 $2^{i}P\ne 0$ . 因此定理4.4实际上表明曲线 $E$ 上存在阶数等于 $2^k$ 的点 $P$ . 此外，数列 $S_i$ 的定义也可以舍弃常系数 $4$ .

对于 $n$ 较大的情形，有 $n>N/{\left(N^{1/4}+1\right)}^2$ ，这时则又回到了最初的定理4.1，$n$ 就是定理中的 $s$ ，此时需要检验是否有 $2^{k}nP=0$ ，以及是否对于 $n$ 的所有素因数 $q$ 都有 $(2^{k}n/q)P\ne 0$ .

综合以上两种情形，给出Tsumura算法的具体描述：

算法4.5（Tsumura算法）：输入正整数 $N=2^{k}n-1$ ，其中 $n$ 是奇数且 $k>2$ ；
【步骤1】选择整数 $x$ 使得Jacobi符号 $(x|N)=-1$ ，再选择整数 $y$ 使得 $(x^3-y^2|N)=1$ . 令 $m=(x^3-y^2)/x\mathrm{mod}N$ ，定义椭圆曲线为 $E:y^2=x^3-mx(\mathrm{mod}N)$ ，则点 $P(x,y)$ 为曲线 $E$ 上的一点；

若 $n<N/{\left(N^{1/4}+1\right)}^2$ ，执行步骤2~4：

【步骤2】计算 $R=nP$ ，若 $R=0$ ，则 $N$ 是合数；
【步骤3】以点 $R$ 为初始条件，计算数列 $S_i(1\le i\le k-1)$ ，若存在 $\text{gcd}(S_i,N)>1$ ，则 $N$ 是合数；
【步骤4】若 $S_k=0(\mathrm{mod}N)$ ，则 $N$ 是素数，否则 $N$ 是合数；

若 $n>{\left(N^{1/4}+1\right)}^2$ ，执行步骤5~7：

【步骤5】若 $2^{k}P=0$ ，则回到步骤1，重新选择 $y$ ；
【步骤6】对于 $n$ 的所有素因数 $q$ ，若存在 $(2^{k}n/q)P=0$ ，则回到步骤1，重新选择 $y$ ；
【步骤7】若 $2^{k}nP=0$ ，则 $N$ 是素数，否则 $N$ 是合数。

下面笔者根据自己的理解做一些补充说明：(1) 在执行上述算法的过程中，如果在任何步骤出现了任何无意义/无定义表达式，或其他不符合椭圆曲线群性质的情况，都可以判定 $N$ 是合数；(2) 对于 $n>{\left(N^{1/4}+1\right)}^2$ 的情形，如果 $n$ 的所有素因数难以计算，并不能像前面的ECPP算法一样使用拟素因数分解+递归策略，因为这里的 $n$ 在 $N$ 确定的情况下是固定的，不会随着所选取的曲线而改变，除了分解它以外没有任何办法；这就要求我们所输入的 $n$ 的素因数分解不能太难，或者需要提前知道 $n$ 的素因数信息，否则算法会非常低效；(3) 事实上，任何大于 $1$ 的奇数都可以写成 $N=2^{k}n-1$ 的形式，这使得Tsumura算法看上去是普适的，但实际上算法并没有给出 $N/{\left(N^{1/4}+1\right)}^2\le n\le {\left(N^{1/4}+1\right)}^2$ 的情形，而实际情况中这种情形的发生并不能忽略不计；例如 $k=5,n=31$ 时有 $N=991$ ，此时 $22.677\approx N/{\left(N^{1/4}+1\right)}^2<n<{\left(N^{1/4}+1\right)}^2\approx 43.702$ ；按照我的理解，如果出现这种情况，Tsumura算法就不再适用。

在 $n$ 较小的情形下，Tsumura算法所需要的时间几乎可忽略不计。而在 $n$ 较大的情形下，Tsumura算法需要的时间几乎完全取决于 $n$ 的素因数分解的难度。如果想利用Tsumura算法对任意奇数 $N$ 进行检验，需要先反算出 $k$ 和 $n$ 再输入到算法中，这时大多数情况下算法无法高效给出结果，其中一部分情况是 $n$ 不满足算法的适用条件，更多的情况是 $k$ 很小而 $n$ 很大，且没有 $n$ 的素因数信息可以利用，算法的运行时间基本上不可接受。

下一篇：【素性检验算法系列（四）：高次扩域检验、APR算法】