有趣的椭圆曲线加密

一、概述

椭圆曲线加密算法依赖于椭圆曲线理论，后者理论涵盖的知识比较深广，而且涉及数论中比较深奥的问题。经过数学家几百年的研究积累，已经有很多重要的成果，一些很棘手的数学难题依赖椭圆曲线理论得以解决（比如费马大定理）。

本文涉及的椭圆曲线知识只是抽取与密码学相关的很小的一个角落，涉及到很浅的理论的知识，同时也是一点比较肤浅的总结和认识，重点是利用椭圆曲线结合数学技巧阐述加密算法的过程和原理。

本文特意构造有比较多的实例方便理解其过程和原理。

二、椭圆曲线

椭圆曲线方程来源于椭圆积分，后者来最初来源于计算椭圆周长的问题，有一段时间的历史了，在欧拉时期就开始研究。椭圆周长没有精确的初等函数的公式表示，只有近似的公式表示，精确的椭圆周长可以用不定积分表示。

现在一般将形如如下形式的积分定义为椭圆积分：

其中R是其两个参数的有理函数，P是一个无重根的3或4阶多项式，而c是一个常数。椭圆曲线方程与P(t)表现形式比较相像。

数学上的椭圆曲线一般由如下形式给出：

E:y^{2}=x^{3}+ax^{2}+bx+c，其中判别式\Delta (E)=-4a^{3}c+a^{2}b^{2}-4b^{3}-27c^{2}+18abc\neq 0

椭圆曲线都是关于X轴对称的曲线。

典型的椭圆曲线如：y^{2}=x^{3}-4x^{2}+16，其图像为：

更多的椭圆曲线图像：

限定\Delta不为零有特殊的意义。如果判别式\Delta (E)等于零，由三次方程判别式判定理可知，方程x^{3}+ax^{2}+bx+c=0存在二重根或者三重根，曲线表现为”自相交”或者有“尖点”。

两个典型的例子是:

C1有三重根，表现为有”尖点”；C2有二重根，表现为“自相交”，它们都不是椭圆曲线，其图像分别如下：

在密码学中用到的椭圆曲线方程一般限定为：

E:y^{2}=x^{3}+ax+b，其中4a^{3}+27b^{2}\neq 0。

也即是这里的二次项系数为0。

三、椭圆曲线算术

椭圆曲线上可以定义一些很有意思的特殊运算规则。一般来说会定义两种运算：加法和数乘运算。加法运算是点与点之间的运算；数乘运算基于加法运算，重复的加法运算就是数乘。

1、实数域的加法运算

1.1、加法运算的几何解释

已知椭圆曲线上两个不同的点P和Q，则这两个点之和R=P+Q可以通过如下操作得到：过P、Q两点做直线L，与椭圆曲线相交于第三点，该点关于X轴的对称点即是所求的R点。椭圆曲线的这种加法运算有比较明确的几何含义。如下所示：

以这种比较奇特的规则来定义加法运算会让人觉得比较怪异，其思想很可能是借鉴于求椭圆曲线有理解的方法（没有去严格考据）。

求椭圆曲线有理解考虑的问题是寻找有理点(x,y)使其满足椭圆曲线方程。其求解过程是在有限的已知有理点的集合中，选两个点P1,P2，作直线与椭圆曲线相交于第三点个有理点Q3。此时如果再利用P1,P2,Q3三个点中的任意两点作直线不能再产生新的有理解（因为他们本身是已经在一条直线上，不会产生新的交点），但是考虑Q3关于X轴对称的点Q'3必定也是有理点，于是可以利用Q3'与P1或者Q3'与P2继续做直线与椭圆曲线相交得到新的有理解，对新的交点再取对称点，以此迭代下去。由此利用交点的对称点作直线来生成新的交点，进而可逐步求解满足椭圆曲线的有理解。

椭圆曲线加法运算的规则中“取交点的对称点”正是与上述求解过程极其相似。

对于加法运算也有另外一种描述：若椭圆曲线上三个点在同一直线上，则他们的和为O，也即是P+Q+R'=O，其中的O是无穷远点或者零点。

更完整的椭圆曲线加法运算规则如下：

1、O+O=O，对任意的P，有P+O=P；O看做零点，对加法运算没有实际贡献（类似于四则运算加法运算中的0）。

2、P=(x,y)的负元是关于X轴对称的点-P=(x,-y)（而不是关于原点对称），P+(-P)=O。过P和-P的直线与X轴垂直，实际上可以看做它与椭圆曲线相交于无穷远点（射影平面，也即是在欧式平面上添加了无穷远点和无穷远直线的平面），因此将也将O视作无穷远点。

3、计算P和Q的和是通过做过P和Q两点的直线，与椭圆曲线相交于第三点，再取该点关于X轴的对称点以此作为P,Q之和，正如上面的几何图形展示的那样。

4、计算P点（P \neq O）的两倍时，是做该点的切线，再取交点S关于X轴的对称点-S，也即是2P=P+P=-S

容易验证，对于椭圆曲线上的点和O点组成的集合，以及集合上定义的二元加法运算，构成一个Abel群。单位元是O点，P(x,y)的逆元是P(x,-y)，封闭性，结合性以及交换性也是显然满足的。

1.2、加法运算的代数解释

几何解释更直观，代数解释更有利于数值计算。

过曲线上P(x_{p},y_{p})和Q(x_{Q},y_{Q})两点（P和Q不互为负元）做直线，求与曲线的第三个交点的问题是很容易用代数的方法来描述的。

也即是求：

其中斜率k=\frac{y_{Q}-y_{P}}{x_{Q}-x_{P}}。

将（2）代入（1）再利用次数对齐的方式容易求得第三个交点的对称点也即P,Q之和R(x_{R},y_{R})为：

如果需要计算倍乘，可以让多个点自身重复相加得到。例如P+P=2P=R，当y_{P}\neq 0时，代数描述为：

x_{R}=\left ( \frac{3x^{2}_{P}+a}{2y_{P}} \right )^{2}-2x_{P}

2、模素数P的加法运算

密码学中普遍采用的是有限域上的椭圆曲线，也即是变量和系数均在有限域中取值的椭圆曲线。使用模素数p的有限域Z_{p}，将模运算引入到椭圆曲线算术中，变量和系数从集合{0,1,2,...,p-1}中取值而非是在实数上取值。

此时讨论椭圆曲线形式如下：

其中(4a^{3}+27b^{2})\mod p\neq 0\mod p，变量和系数均在Z_{p}中取值。

将满足上式的所有非负整数对和O点记为集合E_{p}(a,b)，这是一个有限的离散点集。由此可知集合中的点分布在(0,0)到(p-1,p-1)的象限中，集合中的点有可能刚好也在椭圆曲线上，更多的可能是在椭圆曲线外。例如点(13,7)是满足y^{2} \mod 23=(x^{3}+x+1)\mod 23的点，但是(13,7)并不在椭圆曲线上。

实际上，集合E_{p}(a,b)与模p的加法运算构成循环阿贝尔群，其生成元，阶和生成子群问题在本节后面会讨论。

对于较小的素数p，完全可以暴力穷举找出集合E_{p}(a,b)中的点。比如参数a=1,b=3,p=23，E_{23}(1,3)有27个点（包含O点），暴力穷举这些点分别为（第八节给出了一些分析椭圆曲线问题的demo实现）：

(0,7)        (6,15)        (15,9)
(0,16)       (7,10)        (15,14)
(2,6)        (7,13)        (19,2)
(2,17)       (10,1)        (19,21)
(4,5)        (10,22)       (21,4)
(4,18)       (12,8)        (21,19)
(5,8)        (12,15)       (22,1)
(5,15)       (14,1)        (22,22)
(6,8)        (14,22)        O

E_{p}(a,b)上的加法规则和实数域上的加法基本一致，只是多加了模运算。但是模p的加法没有显而易见的几何解释，只有代数描述。

求解(x_{R},y_{R})的代数表达式为：

其中

例如a=1,b=1,p=23,P(3,10),Q(13,16)，求R=P+Q.

此时P\neq Q，计算\lambda=(\frac{y_{Q}-y_{P}}{x_{Q}-x_{P}})\mod p=(\frac{16-10}{13-3})\mod 23=6\times 10^{-1}\mod 23.

要计算上式首先要计算10^{-1}\mod 23.

令x\equiv 10^{-1}(\mod 23)，由于10 \equiv 10(\mod 23)，所以10x\equiv 1(\mod 23)，利用扩展欧几里德算法求得x=7.

所以

所以R=(0,1).

还可以按照以上规则计算2P,3P等等倍乘点。

实际上E_{23}(1,1)中共有28个点（包含无穷远点O），以P(3,10)开始的所有倍乘点：P,2P,3P...27P,28P可以暴力计算得出:

P=(3,10)
2P=(7,12)
3P=(19,5)
4P=(17,3)
5P=(9,16)
6P=(12,4)
7P=(11,3)
8P=(13,16)
9P=(0,1)
10P=(6,4)
11P=(18,20)
12P=(5,4)
13P=(1,7)
14P=(4,0)
15P=(1,16)
16P=(5,19)
17P=(18,3)
18P=(6,19)
19P=(0,22)
20P=(13,7)
21P=(11,20)
22P=(12,19)
23P=(9,7)
24P=(17,20)
25P=(19,18)
26P=(7,11)
27P=(3,13)
28P=O

容易验证，上述计算过程中Q(13,16)点就是8P，P+Q=P+8P=9P=(0,1)，与上述计算结果是吻合的，读者也可以验证更多的结果。

现在提出一个问题：E_{p}(a,b)中有多少个点呢？这个问题的准确答案并不好回答，但是有一些粗略的规律。

设E_{p}(a,b)中点的个数为N_{p}，并且a_{p}=p-N_{p}，实际上N_{p}与p比较接近。

Hasse定理表明：

这实际上解释了a_{p}的误差，从而有-2\sqrt{p}< a_{p}=p-N_{p}<2\sqrt{p}

也就是p-2\sqrt{p}<N_{p}<p+2\sqrt{p}，等价于p+1-2\sqrt{p} \leqslant N_{p}\leqslant p+1+2\sqrt{p}

接下来的一个问题是子群、生成元和阶：

不难得知，有限域上的椭圆曲线的点和加法运算构成一个有限交换群S。

以一个点G作为生成元，进行重复的加法运算，能够生成一个子群S'。群S'有可能与S相同，更有可能是S的一个真子群。

比如椭圆曲线为y^{2}=x^{3}+x+1，p=23

以G=(3,10)生成的子群刚好等于S，因为28G=O，说明此时S还是一个有限循环群，群中每个元素都能由(3,10)通过重复的加法运算得到，群中元素的个数为28，其阶为28(最小的使得nG=O成立的n)，生成元为(3,10)。

但是以G=(6,19)生成的子群是S的真子群，因为14G=O，它也是一个有限循环群，因为群中每个元素都能由(6,19)通过重复的加法运算得到，该子群的元素个数为14，阶为14（最小的使得nG=O的成立n），生成元为(6,19)。

拉格朗日定理还告诉我们，生成的子群S'的阶与原来的群S的阶存在约数关系，也就是说子群S'的阶只能是1,2,4,7,14,28。群S的阶与子群S'的阶的比值一般称为协因子(cofactor)：h=\frac{\left | S \right |}{\left | S' \right |}。一般会取h=1，以保证生成子群的阶比较大。

以G=(6,19)生成子群过程（为方便计算，取O=(-1,-1)）：

1P=(6,19)        15P=(6,19)
2P=(13,16)       16P=(13,16)
3P=(7,11)        17P=(7,11)
4P=(5,19)        18P=(5,19)
5P=(12,4)        19P=(12,4)
6P=(17,20)       20P=(17,20)
7P=(4,0)         21P=(4,0)
8P=(17,3)        22P=(17,3)
9P=(12,19)       23P=(12,19)
10P=(5,4)        24P=(5,4)
11P=(7,12)       25P=(7,12)
12P=(13,7)       26P=(13,7)
13P=(6,4)        27P=(6,4)
14P=(-1,-1)      28P=(-1,-1)

在后面的ECC加密算法过程中会有一个给定的基点G（也就是生成元）生成一个子群，然后秘钥空间在此子群取得，一般会要求保证子群的阶会尽量大，基点及其子群的阶n都是公开的信息。

四、椭圆曲线密码学中的离散对数问题

构造一个数学难题来保证加密的安全性是现代密码学中加密算法的主要思想。类似RSA算法中大数的质因子分解难题一样，椭圆曲线也有类似的数学难题。

考虑Q=kP，其中Q,P\in E_{p}(a,b),k<p。

对于给定的k,p计算Q是很容易的；反过来给定Q,P，计算k是相当困难的，这就是椭圆曲线的离散对数问题（这里之所以称之为离散对数问题大概是为了与其他加密算法的说法保持一致，便于理解）。

正因为如此，可以将Q作为公钥，公开出去；k作为私钥，秘密保管，通过公钥来破解私钥十分困难。

目前由椭圆曲线公钥求解私钥的最有效算法复杂度为O(\sqrt{p})，其中p是阶数n的最大素因子。

五、椭圆曲线实现秘钥协商（ECDHE）

1、原理和流程

利用模p有限域上椭圆曲线算术规则可以用来实现秘钥协商，其流程如下：

（1）Alice和Bob会共享一些椭圆曲线的参数信息：(p, a, b, G, n, h)，其中a,b确定了椭圆曲线方程，p确定了模p的有限域，G是中的生成元，用于生成子群，要求G的阶n应该尽量大，G的阶是使得nG=O成立的最小的整数，也即是(n-1)G所代表的点与G点的横坐标刚好相同，协因子h一般等于1。

（2）Alice选择小于n整数n_{A}，然后计算P_{A}=n_{A}\times G，将P_{A}发送给Bob。

（3）同理，Bob也选择小于n整数n_{B}，然后计算P_{B}=n_{B}\times G，将P_{B}发送给Alice。

（4）Alice通过如下计算得出协商的秘钥：K_{A}=n_{A}\times P_{B}；Bob通过如下计算得出协商的秘钥：K_{B}=n_{B}\times P_{B}，容易证明，K_{A}=K_{B}。

因为K_{A}=n_{A}\times P_{B}=n_{A}\times (n_{B}\times G)=n_{B}\times (n_{A}\times G)= K_{B}。

由此Alice和Bob计算得到相同的秘钥，达到秘钥协商的目的。

2、实例

例如，参数a=0,b=-4，椭圆曲线方程为y^{2}=x^{3}-4，p=211,G=(2,2)，因为240G=(2,209)所以n=241。

Alice选择n_{A}=151，P_{A}=151G=(62,59)；

Bob选择n_{B}=171，P_{B}=171G=(209,153)；

他们协商的秘钥K=151P_{B}=171P_{A}=151(209,153)=171(62,59)=(95,194)=(151*171\%n)G=34G。

六、椭圆曲线实现加解密

1、原理和流程

利用椭圆曲线来进行加密和解密也是与RSA一定程度的类似，每一个用户都有属于自己的公钥和私钥。私钥就是用户选定的数字n，私钥自己保存；公钥就是由P=nG，计算出来的点，公钥公开。

假设Alice与Bob进行加密通信，其加密的流程如下：

（1）Alice首先将明文消息转换（编码）为E_{p}(a,b)中的P_{m}(x,y)，然后随机选定一个正整数k，并且利用Bob的公钥P_{B}通过如下计算出密文：

因此，密文实际上是有两个点组成。

（2）Bob收到密文C_{m}，利用自己的私钥n_{B}进行如下计算，可以解密得到明文：

也就是用第二个点P_{m}+kP_{B}减去第一个点kG与自己的私钥n之积。

2、实例

考虑参数a=0,b=-4,p=199,G=(2,2)，椭圆曲线方程为y^{2}=x^{3}-4，

Bob选定的私钥为n_{B}=119，其公钥P_{B}=119G=(183,173)。

Alice希望将消息P_{m}=(76,66)加密后发送给Bob，于是Alice随机选定正整数k=133，并通过Bob的公钥加密得到密文：

Bob收到密文消息利用自己的私钥n_{B}=119进行解密：

P_{m}+kP_{B}-n_{B}(kG)=(180,163)-119(40,147)=(180,163)-(98,52)=(76,66)=P_{m}.

由此，Bob顺利解密得到明文消息，Alice与Bob之间完成加密通信。

七、杂项

公钥加密算法中除了ECC，还有另外一个广泛使用的加密算法–RSA公钥加密算法。

ECC与RSA相比，主要的优点是在相同的安全级别下ECC使用的秘钥长度要短很多，由此带来处理速度、带宽和存储空间上的额外优势。下表展示了不同加密算法秘钥的位数对比情况：

比特币中使用椭圆曲线密码算法，以保证比特币网络中信息安全性和签名认证问题。每个用户都有属于自己的公钥和私钥。利用私钥可以对交易信息进行签名（ECDSA），其他人可以利用其公钥进行认证，公钥也用来构造钱包地址。所使用的椭圆曲线是采用了Certicom推荐的椭圆曲线secp256k1 ，其参数由六元组构成：D=(p,a,b,G,n,h)，其中p是很大的质数，h就是协因子，表征由生成元G生成子群的阶与母群的阶的关系，前文有解释；参数信息都是公开的，有兴趣可以到网上查到。

八、数值计算Demo实现

以下实现一些比较常用的方法，主要涉及暴力穷举E_{p}(a,b)中的点（无穷远点O=(-1,-1)），求有限交换群的阶，椭圆曲线加法运算和乘法运算，以及相关的求模反元素和扩展欧几里德算法等。主要在学习椭圆曲线加密过程中，方便快速的构造用例和实验、验证。主要实现核心的计算原理，细节部分的处理难免有疏漏。

/**
 * Date:2018.05.16
 * Author:Qcer
 * Function:ECC
 * */
package com.demo;
import java.util.ArrayList;
public class ECCImpl {
    private long a = 0;
    private long b = 0;
    private long p = 0;
    private ECCImpl() {}
    public ECCImpl(long a, long b, long p) {// 参数初始化，椭圆曲线方程:y^2=x^3+ax+b的形式，
        this.a = a;
        this.b = b;
        this.p = p;
    }
    // 暴力搜索(0,0)到(p-1,p-1)象限且满足模p方程的所有点
    public ArrayList search() {
        ArrayList<Long[]> list =new ArrayList<>();
        long left,right;
        for (long x = 0; x < p; x++) {//i=x   from 0 to p-1
            for (long y = 0; y < p; y++) {//j=y form 0 to p-1
                left = y * y % p;
                right = (x * x * x + a * x + b) % p; 
                if (left == right) {
                    list.add(new Long[] {x,y});
                    System.out.println("("+x+","+y+")");
                }
            }
        }
        return list;
    }
    //获取椭圆曲线上有限交换群的阶，包含无穷云巅O
    public long getOrder() {
        return search().size()+1;
    }
    //自定义求模运算，主要考虑a小于0的情况
    public long mod(long a, long p) {
        long r = Math.abs(a) % p;
        return a >= 0 ? r : p - r;
    }
    //椭圆曲线上的加法运算,R=P+Q
    public long[] add(long x1, long y1, long x2, long y2) {
        if (x1 == -1) {// O+P=P
            return new long[]{x2,y2};
        }
        if (x2 == -1) {// P+O=P
            return new long[]{x1,y1};
        }
        long[] point = new long[2];
        long ratio = 0;//系数    0<=ratio<=p-1
        if (x1 == x2) {
            if (y1 == y2) {// P==Q
                ratio = mod((3  *x1 * x1 + a) * (getMMI(2 * y1,p)),p);
            }else { // P+P=O
                return new long[] {-1,-1};
            }            
        }else {// P!=Q
//            long flag = (y2-y1)*(x2-x1)>0?1:-1;//锚定正负号
            long mask = 1<<63;
            long flag = (((y2 - y1) & mask) ^ ((x2 - x1) & mask)) == 0 ? 1 : -1;
            ratio =  mod(flag * Math.abs(y2 - y1) * (getMMI(Math.abs(x2 - x1),p)),p);
        }
        point[0] = mod(ratio * ratio - x1 - x2, p);
        point[1] = mod(ratio * (x1 - point[0]) - y1,p);
//        System.out.println("("+point[0]+","+point[1]+")");
        return point;
    }
    //椭圆曲线上的乘法运算,乘法运算定义为重复的加法运算:kP=P+P+P+...+P
    public long[] multiply(int x, int y, int times) {
        long[] point = {x,y};
        for (long i = 2; i <= times; i++) {
            point = add(point[0],point[1],x,y);
            System.out.println(i+"P="+"("+point[0]+","+point[1]+")");
            if (point[0] == -1) {
                return point;
            }
        }
        System.out.println(times+"P="+"("+point[0]+","+point[1]+")");
        return point;
    }
    //利用扩展欧几里算法得求模范元素,ax=1(mod b) => ax+b(-y)=1=gdc(a,b)
    public long getMMI(long a, long b) {//modular multiplicative inverse
        long mmi = gcdExt(a,b)[1];
        while (mmi < 0) {
            mmi += b;
        }
        return mmi;
    }
    //扩展欧几里得算法:ax+by=g=gcd(a,b)  =>  tuple[1]x+tuple[2]y=gcd(a,b)
    public long[] gcdExt(long a, long b) {
        long ans;
        long[] tuple = new long[3];
        if (b == 0) {
            tuple[0] = a;
            tuple[1] = 1;
            tuple[2] = 0;
            return tuple;
        }
        long[] temp = gcdExt(b,a%b);
        ans = temp[0];
        tuple[0] = ans;
        tuple[1] = temp[2];
        tuple[2] = temp[1] - (a / b) * temp[2];
        return tuple;
    }
    public static void main(String[] args) {
        // TODO Auto-generated method stub
//        ECCImpl ecc = new ECCImpl(0, -4, 199);//0, -4, 211
//        ecc.search();
//        ecc.add(3,10,13,16);
        ECCImpl ecc = new ECCImpl(1, 1, 23);
        System.out.println(ecc.getOrder());
//        ecc.search();
//        ECCImpl ecc = new ECCImpl(4, 20, 29);
//        ecc.search();
        
//        ecc.multiply(12,19, 28);
    }
}

九、References

1、密码编码学与网络安全原理与实践

2、http://andrea.corbellini.name/

3、数论概述.第四十三章-第四十八章

4、https://mp.weixin.qq.com/s/jOcVk7olBDgBgoy56m5cxQ

5、https://blog.csdn.net/mrpre/article/details/72850598