ECPP——利用有限域上的椭圆曲线，精确判定素数的算法

要说明的是，现在这篇东西是研究用椭圆曲线精确判定一个正整数是否为素数，而不是研究用椭圆曲线加解密的！

当然，后面提供的一些关于有限域上的椭圆曲线上的点的加法，乘法，除法的实装函数都是可以直接用于加解密的。

这次要研究的ECPP，它是接我以前的一篇关于素数的blog（http://blog.csdn.net/chen09/archive/2008/03/24/2214219.aspx ），“试除法”以外的素数精确判定的方法之一，它的算法复杂度为：O((log n)^6)。

 

“有限域上的椭圆曲线”，这个概念的子概念，顾名思义有两个：1，有限域，2，椭圆曲线。这两个子概念都很容易理解，但是合在一起，定义在有限域上的椭圆曲线，而不是定义在普通的实数上的椭圆曲线，是怎么一回事儿？我大概搞了整整一个月才明白。网上的中文和日文资料都太少了，而且错误百出，让我走了不少弯路。在这里想把我所知道的写清楚，让后来人可以多一份可以参考的资料。

 

什么是椭圆曲线？什么是有限域？什么是定义在有限域上的椭圆曲线？它有什么特征？和判断素数有什么关系？

这里有一篇《ECC加密算法入门介绍》（http://www.docin.com/p-853208.html ），写得非常的好，错误几乎没有。通读一边有益无损。

初中时，我们就学过，方程x^2+y^2=1是个圆，以(0,0)为中心，半径为1。
稍作变形，方程(x-a)^2+(y-b)^2=c^2也是个圆，以(a,b)为中心，半径为c。
然后说说椭圆方程，(x/a)^2+(y/b)^2(a>b>0)是个标准方程。
稍作变形，((x-h)/a)^2+((y-k)/b)^2=1。h,k自然是x轴和y轴上的位移了。
有没有都回想起来呢？如果回想起来了，那么把它们都忘了吧！因为椭圆曲线和椭圆方程没有什么关系。8-)

 

《ECC加密算法入门介绍》（http://www.docin.com/p-853208.html ）里面讲椭圆曲线时先介绍了射影平面，这个有点超纲（在我看来，所有计算机的题目都是用初中数学可以解决的，凡是用到高中数学以上的东西，都是超纲的）。所以还是直接考察笛卡尔坐标系上的，椭圆曲线：y^2+a1*x*y+a3*y=x^3+a2*x^2+a4*x+a6。

看上去有点复杂？其实已经很简单化了！要知道：两个元（x和y）且是3次的一般方程要复杂的多。而现在要求每个点都是光滑的才是椭圆曲线，所以方程的形式真的已经简单很多。

还是觉得有点复杂？没关系，实际运用中的椭圆曲线的形式将更为简单。

 

现在，这个椭圆曲线是定义在实数上的，对于曲线上的点都可以求导（因为我们要求它是光滑的），求导的方法很简单，两边分别对y和x求偏导就可以。求导自然是为了就切线，为什么要算切线？稍微过一会儿说。

 

现在说一下，曲线上的2个点的加法。

什么是加法？

大家都知道1+1=2(不是哥德巴赫猜想，是小学或者说幼儿园的数学)，但是1+1为什么等于2？因为这个“+”是人为定义的一个算法，在首先自然数上定义，然后扩展到整数，有理数，实数，复数。

罗素（Bertrand Russell, 1872-1970）和他的导师怀特黑(Alfred Whitehead, 1861-1947)用了几十页的篇幅，最终成功地证明1+1=2（参见，若奇的《孤独的破译者和他的计算机器》 http://xys4.dxiong.com/xys/netters/psi1b/ruoqi.txt ）。

 

现在，天才的数学家们也定义了椭圆曲线上的两个点的加法！
令椭圆曲线两个点P，Q，做一条过点P和点Q的直线（如果P，Q重合，则做P点上的切线（这就是刚才说到为什么要求导求切线的原因）），交椭圆曲线上的另一点R'，过R'做y轴的平行线交椭圆曲线与点R。规定：P+Q=R!

 

为什么要定义这么个不着边际的运算法则？
的确，它离我们的目标，无论是精确素数判断还是费马大定理都离得很远。
这里，我想可以运用人择原理：如果它是没有用的，那么就不会被我们注目；虽然它不是那么直观，但是它是有用的，所以我们现在在学习它。

 

就好象，1+1=2，原来也许没有什么用，但是要计算1个苹果加上1个生梨，共有几个水果时，起到了作用，所以1+1=2被保留了下来。

1+1威力如此之大，不但用于幼儿的数学启蒙，也创造了诺贝尔奖获得者罗素！大家都知道诺贝尔是没有数学奖的，罗素得到的是文学奖。罗素和怀特黑的两千多页的《数学原理》（和牛顿的伟大的著作同名，说明他们有如此的自信）耗尽了他的心血，哥德尔(Kurt G?del, 1906-1978)和他的“哥德尔不完备定理”的出现，使罗素则彻底丧失了对数理逻辑的兴趣，他改而倾心哲学和文学，关注社会事务。终于，在1950年他获得了诺贝尔文学奖。

 

1+1威力还不止于此，同样因为哥德尔的出现，也使冯·诺依曼(John von Neumann, 1903-1957)放弃了纯粹数学的研究。1944年，他与奥斯卡·摩根斯特恩合著出版《博弈论与经济行为》，标志着现代系统博弈理论的的初步形成，从而被称为“博弈论之父”；1945年6月，他戈德斯坦、勃克斯等人联名发表了著名的“101页报告”，而成了现代计算机之父。

 

哥德尔是有颠覆性的，他的不完备定理最初可以描述为“如果将算术运算公理引入一阶谓词系统，则由此而成的算术运算系统是不完备的，即存在这样的命题，它们的真伪是不可证明的。”而一阶谓词显然与50年后（20世纪70，80年代）的计算机届有直接的关系，以日本为首的推崇的第5代计算机就是以Prolog为主要语言，而Prolog就是以一阶谓词逻辑为基础的。虽然第5代计算机以失败而告终，但是不得不佩服日本敢为天下先的勇气，如果成功，那么我们现在的世界也许完全不是现在这个样子的了。

 

扯了半天，想说明的就是，P+Q=R!

 

定义在实数上的椭圆加法，还有图像可以看，但是定义在有限域上的椭圆曲线及其加法，就很难理解了。杨振宁曾说过，现在的数学论文，有的是看了第一页就看不下去的；有的是看了第一句就看不下去的。幸亏，有限域上的椭圆曲线还是可以看下去的，虽然它比“1个苹果加上1个生梨”要复杂些。但是只要有了抽象的头脑，它并不比“2头牛加上3头羊等于5头动物”要难多少。

刚才，我们把椭圆曲线定义在实数上，它是连续的。现在我们把它定义在有限域上。

有限域的定义很简单：包含有限个元素的域被称为有限域。

 

如果从有限域扯开去，那么和我们熟知的一些数学和计算机都有很大关系。比方说，每个方程对应于一个域，即含有方程全部根的域，称为这方程的伽罗瓦域，这个域对应一个群，即这个方程根的置换群，称为这方程的伽罗瓦群。伽罗瓦域的子域和伽罗瓦群的子群有一一对应关系；当且仅当一个方程的伽罗瓦群是可解群时，这方程是根式可解的。

这个理论是不是很难理解？没有关系，它的推论就容易理解的多：1，可以得出五次以上一般代数方程沒有公式解（华罗庚曾经指出过某人的5次方程公式解的错误）；2，尺规作图三等分任意角和作倍立方体不可能；等等。

 

不过我们不需要知道那么多。我们只要知道定义在有限域上的椭圆曲线就是：
那些点P或者Q的坐标(x,y)都不是实数范围那么大，而是有限的。比方说，（0，1，2，3，4）或者（0，1，2，3，4，5，6）。注意到，这连个范围的个数5或者7都是素数！

 

刚才说到，椭圆曲线函数太复杂了，实际运用中，我们用比较简单的形式：y^2 = x^3 + a*x +b 。（含有x*y项的椭圆曲线方程在这篇文章里面被俺完全忽略了，有兴趣的可以参阅网上的文章或者自己推导一下，公式比不含x*y项的要复杂一些，但是原理是一样的）
然后，很容易写出个程序，把这个方程的解都算出来。

 

但是，说到这里要稍微等等。我们虽然把椭圆函数定义到了有限域上，但是这时候有限域的威力还没有发挥。我们把有限域的个数（素数）也放到椭圆曲线方程里面去！

y^2 = x^3 + a*x +b (mod p)，p为素数，x,y 属于 (0，1，...，p-1)

好了！现在我们把有限域，椭圆曲线合成起来了。

举个例子，我们来最简单的y^2 = x^3 + x +1 (mod 5)
下面就是所有的解：
[x=3 y=1]
[x=3 y=4]
[x=4 y=3]
[x=4 y=2]
[x=2 y=1]
[x=2 y=4]
[x=0 y=4]
[x=0 y=1]
要验算很简单吧？把上面的x和y代入上面的方程，看看左右是否相等。别忘记，方程左右要对5取余。

 

那么上面的解，也就是椭圆曲线的各个点，P+Q=R这个加法怎么办呢？形式上很像的：
加法： a+b ≡ c (mod p)
乘法： a*b ≡ c (mod p)
除法： a/b ≡ c (mod p) 即 a * b^(-1) ≡ c (mod p)

 

乘法？第一次听说吧？椭圆曲线上的点和点是没法相乘的，这里的乘法就是连加，和初等数学类似。a*b 也就是 a + a + ... + a ，一共b个a相加。

乘法理解后，除法也就容易了，a/b也就是 a * b^(-1)，而b^(-1)也就是要求出一个数x（0<= x <p）让x*b ≡ 1 (mod p)

同样是y^2 = x^3 + x +1 (mod 5)

加法的例子：
[x=3 y=1] + [x=4 y=3] = [x=2 y=1]

 

乘法的例子：
[x=3 y=1] * 2 = [x=3 y=1] + [x=3 y=1] = [x=0 y=1]
[x=3 y=1] * 3 = [x=3 y=1] * 2 + [x=3 y=1] = [x=0 y=1] + [x=3 y=1] = [x=2 y=4]

 

这里的加法就是前面定义过的加法，那个画通过被加的2个点的直线，然后相交，然后再画平行y轴的直线，再相交的那个点，就是加法的结果。算法虽然比较复杂，一般公式推导却可以自己来。这里就不推导了，直接给出一般公式：

椭圆曲线： y^2 = x^3 + a*x +b (mod p)

P(x1,y1) + Q(x2,y2) = R(x3,y3)
x3 ≡ λ^2 - x1 - x2  (mod p)
y3 ≡ λ * (x1 - x3) - y1 (mod p)
其中
如果，PQ重合，则 λ ≡ (3 * x^2 + a) / (2 * y) (mod p)
PQ相异，则 λ ≡ (y2-y1)/(x2-x1) (mod p)

要注意，上面的除法，是我们刚刚定义的除法（就是除数乘以结果的对p取余等于被除数对p取余），而不是实数上的除法。

有了乘法，我们再引入一个概念，阶（或者秩（不是线性代数行列式里面的秩）或者级数）。

我们还是拿最简单的 y^2 = x^3 + x +1 (mod 5) 来做例子。
最简单的一个解P[x=0 y=1]，我们来计算，P * 2，P * 3，...
P * 1 : [x=0 y=1]
P * 2 : [x=4 y=2]
P * 3 : [x=2 y=1]
P * 4 : [x=3 y=4]
P * 5 : [x=3 y=1]
P * 6 : [x=2 y=4]
P * 7 : [x=4 y=3]
P * 8 : [x=0 y=4]
那么P * 9是多少？也许有人说，岂不是很简单？套用上面的公式就可以。但是......
P * 9 = P * 1 + P * 8，P * 1 和 P * 8的x都是0，也就是说通过它们的直线是垂直于x轴，平行于y轴的，就不能和椭圆曲线再相交了！用公式套用时，λ算不出了。
同样，P * 9 = P * 2 + P * 7 = P * 3 + P * 6 = ......可惜，P * 2 和 P * 7的x是一样的，P * 3 和 P * 6的x也是一样的。也就是说P * 9是不存在的，或者说交点在无限远。

而，这个数字9就是就是点P[x=0 y=1]的阶。

重新定义一下，椭圆曲线上的某个点的阶n就是，使P*n等于无限远。

上面关于 y^2 = x^3 + x +1 (mod 5)的点[x=0 y=1]的阶，可以参考下面的链接《椭圆曲线密码系统》的第11页。http://www.docin.com/p-17972799.html 但是很可惜，这个链接里面P * 3算错了，应该是[x=2 y=1]，而不是[x=4 y=2]。

 

大家是不是有点头晕了？至少我自学到这里，已经非常头晕了。但是，再忍一忍，我们的素数判断算法呼之欲出了。那就是：

定义在有限域p上的椭圆曲线上的每个点（解），必然有阶！
具体的证明就不写了，我们也不感兴趣，我们只要结论，只要有可操作的方法。
现在，we got it。

 

从头开始整理一遍：
1，简单的椭圆曲线：y^2 = x^3 + a*x +b (mod p)
2，我们可以取a = 1, b = 1 ，但是必须要保证 4*a^3 + 27*b^2 mod p 不等于零。比方说，我们要判定31是不是素数时，就不能用y^2 = x^3 + x + 1 (mod 31)，这个时候，我们可以取a = 1, b = 2。
3，我们计算出所有的解P(x,y)，其中 0 <= x < p，0 <= y < p。
4，计算每个解P(x,y)的乘法P * n , 其中n从2开始往上递增。每次计算出来的P * n的x值都保存起来，如果后来计算出来的P * n的x值和前面重复了，那就说明找到了这个解的阶，这个时候处理下一个解，如果所有的解都处理了，都找到了阶，那就这个p就是素数；如果计算P * n是发生了错误，一般来说，就是计算λ时，那个除法没有解，那么这个p就是合数（不是素数）。

 

上面就是计算的流程了。

 

然后，就让我们把它们变成Java的代码吧。

 

首先，先定义一下简单的Exception，当计算出错时可以被抛出。

 

package elliptic; public class EllipticException extends Exception { /** * */ private static final long serialVersionUID = 3528609824516918284L; public EllipticException() { } public EllipticException(String message) { super(message); } public EllipticException(Throwable cause) { super(cause); } public EllipticException(String message, Throwable cause) { super(message, cause); } }

 

这个没有什么好解释的吧？

 

然后我们来定义椭圆曲线：

 

package elliptic; /** * * y^2 = x^3 + a*x +b (mod p) * @author z1j2110 * * */ public class Line { int a; int b; int p; static public Line newInstance(int a, int b, int p) throws EllipticException { return new Line(a, b, p); } private Line(int a, int b, int p) throws EllipticException { if ((4 * a * a * a + 27 * b * b) % p == 0) throw new EllipticException("4*a^3 + 27*b^2 == 0 mod p : a=" + a + "b=" + b + "p=" + p); this.a = a; this.b = b; this.p = p; } @Override public String toString() { return "[a=" + a + " b=" + b + " p=" + b + "]"; } }

 

这个Line类就是我们的椭圆曲线了，里面的a,b,p都解释过了，还有印象吧？

 

然后就是解（点）了：

package elliptic; public class Point { int x; int y; static Point newInstance(int x, int y) { return new Point(x, y); } Point(int x, int y) { this.x = x; this.y = y; } @Override public boolean equals(Object obj) { if (!(obj instanceof Point)) return false; Point other = (Point) obj; if (this.x != other.x) return false; if (this.y != other.y) return false; return true; } @Override public String toString() { return "[x=" + x + " y=" + y+"]"; } }

 

这个也没有什么好说的，就是点(x,y)而已。

 

关键的东西来了：

 

package elliptic; import java.util.Collections; import java.util.HashMap; import java.util.HashSet; import java.util.Map; import java.util.Set; public class Elliptic { static private Map<Integer, Point> mapResult = Collections.synchronizedMap(new HashMap<Integer, Point>()); static public boolean isPrime(int x) { Line l = null; int i = 1; while (l == null) { try { l = Line.newInstance(1, i, x); } catch (EllipticException e) { i++; } } Set<Point> set = getSolution(l); for (Point p : set) { if (p.y == 0) continue; int n = 0; Set<Integer> setX = Collections.synchronizedSet(new HashSet<Integer>()); try { while (true) { Point point; point = multiply(p, ++n, l); if (setX.contains(point.x)) { break; } else { setX.add(point.x); } } } catch (EllipticException e) { return false; } } return true; } static public Set<Point> getSolution(Line l) { Set<Point> set = Collections.synchronizedSet(new HashSet<Point>()); for (int x = 0; x < l.p; x++) { for (int y = 0; y < l.p; y++) { Point p = Point.newInstance(x, y); if (isSolution(p, l)) set.add(p); } } return set; } static public boolean isSolution(Point p, Line l) { int left = (p.y * p.y) % l.p; int right = (p.x * p.x * p.x + l.a * p.x + l.b) % l.p; if (right < 0) right += l.p; return left == right; } static public Point plus(Point p, Point q, Line l) throws EllipticException { // if (!isSolution(p, l)) // throw new EllipticException(p + " is not the solution of line " + l); // // if (!isSolution(q, l)) // throw new EllipticException(q + " is not the solution of line " + l); int k = lambda(p, q, l); int x3 = (k * k - p.x - q.x) % l.p; if (x3 < 0) x3 += l.p; int y3 = (k * (p.x - x3) - p.y) % l.p; if (y3 < 0) y3 += l.p; return new Point(x3, y3); } static public Point multiply(Point p, int x, Line l) throws EllipticException { mapResult.clear(); // mapResult.put(1, p); return multiply(p, x, l, mapResult); } static private Point multiply(Point p, int x, Line l, Map<Integer, Point> mapResult) throws EllipticException { if (x == 1) return p; if (mapResult.containsKey(x)) return mapResult.get(x); int half = x / 2; while (half > 0) { try { Point point = plus(multiply(p, half, l, mapResult), multiply(p, x - half, l, mapResult), l); mapResult.put(x, point); return point; } catch (EllipticException ee) { half--; } } throw new EllipticException("can't multiply " + p + "*" + x + " mod " + l.p); } static public int lambda(Point p, Point q, Line l) throws EllipticException { if (p.x == q.x) { if (p.y == q.y) { return (3 * p.x * p.x + l.a) * inverse(2 * p.y, l.p); } else throw new EllipticException("p+q=∞ p=(" + p + ")q=(" + q + ")"); } else return (q.y - p.y) * inverse(q.x - p.x, l.p); } static public int inverse(int x, int p) throws EllipticException { while (x < 0) x += p; for (int i = 0; i < p; i++) { if (i * x % p == 1) return i; } throw new EllipticException(x + " * x == 1 mod " + p + " x is not exist."); } }

 

这个就是我们的主要的功能类了。稍微解释一下。

 

static public boolean isSolution(Point p, Line l)

这是用来判断，点p是否是椭圆曲线l的解？

 

static public Set<Point> getSolution(Line l)

这是用来计算椭圆曲线的所有解得函数，返回值是个集合，里面的元素是Point点

 

static public int lambda(Point p, Point q, Line l) throws EllipticException

顾名思义，这是计算λ值的函数。

 

static public int inverse(int x, int p) throws EllipticException

顾名思义，这是计算除法时要用到，求(^-1)的函数。

 

static public Point plus(Point p, Point q, Line l) throws EllipticException

顾名思义，这是加法。

 

static public Point multiply(Point p, int x, Line l) throws EllipticException { mapResult.clear(); // mapResult.put(1, p); return multiply(p, x, l, mapResult); } static private Point multiply(Point p, int x, Line l, Map<Integer, Point> mapResult) throws EllipticException { if (x == 1) return p; if (mapResult.containsKey(x)) return mapResult.get(x); int half = x / 2; while (half > 0) { try { Point point = plus(multiply(p, half, l, mapResult), multiply(p, x - half, l, mapResult), l); mapResult.put(x, point); return point; } catch (EllipticException ee) { half--; } } throw new EllipticException("can't multiply " + p + "*" + x + " mod " + l.p); }

顾名思义，这是乘法。

稍微说明一下，这里用到了递归。因为显然当计算P*n时，n非常大的话，n个P连加也是很费时间的，所以要用点乘法，学过CPU内核设计的都知道。不记得的话，回去看看本科计算机的《硬件原理》。

比方说要计算P*9，就是计算P*4+P*5；计算P*4，就是计算P*2+P*2；P*2 = P+P；P*5 = P*2+P*3；P*3 = P*1+P*2；把计算出的中间结果保存在mapResult中，可以反复利用，以加快速度。

 

最后，我们的判定函数就是：

static public boolean isPrime(int x) { Line l = null; int i = 1; while (l == null) { try { l = Line.newInstance(1, i, x); } catch (EllipticException e) { i++; } } Set<Point> set = getSolution(l); for (Point p : set) { if (p.y == 0) continue; int n = 0; Set<Integer> setX = Collections.synchronizedSet(new HashSet<Integer>()); try { while (true) { Point point; point = multiply(p, ++n, l); if (setX.contains(point.x)) { break; } else { setX.add(point.x); } } } catch (EllipticException e) { return false; } } return true; }

如果你耐心的看到这里，上面的函数一点不难理解。

就是上面的4步的流程了。

 

然后让我们来写一段测试用的代码：

public static void main(String[] args) throws EllipticException { for (int i = 10; i < 1000; i++) { boolean isPrime = Elliptic.isPrime(i); if (isPrime) System.out.print(i + " "); } System.out.println(); } }

 

这是用来，输出所有的10到1000之间的素数。结果为：

11 13 17 19 23 29 31 37 41 43 47 53 59 61...... 947 953 967 971 977 983 991 997

大家可以验证一下。

 

最后，说说算法复杂度。

首先要算出椭圆曲线的解，那就是n^2了，因为x从0到p-1，y也从0到p-1，逐一找出解。

然后，对于每个解都要看有没有阶，为了得到阶，要计算P*n，所以又多了个n^1，当然我们用了点乘法，所以是log2的n^1

然后就是，加法时要计算除法，除法时要从1到p检测，所以又有个n^1。

我的算法的复杂度是O((log n)*n^3)。说实话，我也不知道O((log n)^6)（应该是正解）这个值是怎么来的。

呵呵。

题外的一些话：
1，在相同的安全强度下，ECC的密钥长度和处理速度比RSA要短且快，随着RSA被破解，以后或许会是ECC的世界，它更适合放在IC卡或者手机里。在日本，上网的手机用户已经超过PC用户，以后网上商店或者网上拍卖的网站，将越来越倾向于移动用户，所以ECC应该大有可为。
2，费马大定理，也称费马最后定理（当整数 n>2 时，关于x,y,z的不定方程：x^n+y^n=z^n的自然数解是不存在的）的被证明过程中，诞生了代数几何上的椭圆曲线。在计算机上最大的应用就是前面讲的非对称加密。在数学上，谷山·志村猜想（椭圆曲线和模形式一一对应）的最后被证明直接可以反证出费马大定理。（因为可以证明，如果x^n+y^n=z^n n>2有自然数解，谷山·志村猜想将不成立）。