零基础完全理解ECC加密算法（软件方法）

mlsily

1.算法简介
ECC加密算法是在椭圆曲线形成的加法群上的EIG算法（EIG原本是在乘法群上的），是我国二代身份证使用的加密算法。

2.前置知识
椭圆曲线：在实数域上，满足方程y^2 = x^3 + ax + b 的所有点构成的集合。记为椭圆曲线E(a,b)。系数4a^3 + 27b^2 ≠ 0 的椭圆曲线，被称为有效的椭圆曲线。

无穷远点：平行线相交于无穷远点。
E(a,b)只有一个无穷远点，即平行于Y轴的直线相交的点，记为｛O｝。
无穷远点是加法运算的单位元。

椭圆曲线上的加法运算：任取E(a,b)上两点P和Q，过P、Q的直线交E(a,b)于点R（当P、Q重合则在该点做E(a,b)的切线），过点R做Y轴的平行线交E(a,b)于R'，定义P+Q = R'。
由定义进行推导可得出公式：若P ≠ Q，斜率k = (y1-y2)/(x1-x2)，若P=Q，斜率k = (3x1^2 +a)/2y1。 P+Q = (k^2 -x1-x2,k(x1-x3)-y1 )。


3. 算法原理

1.生成椭圆曲线Ep(a,b)。
2.取曲线上一点G。
3.选取整数k，1<k<|Ep(a,b)|（|Ep(a,b)|表示曲线|Ep(a,b)|上的点的数量，即曲线Ep(a,b)构成的加法群的阶。）
4. 计算K = kG  (即k个G相加)
公钥：（Ep(a,b)，G，K)
私钥：k
加密过程：将明文m编码为Ep(a,b)上的一点M，生成随机整数r（1<r<|Ep(a,b)|）,
                计算C1 = rG, C2 = M+rK
                密文为C1,C2
解密过程：M = C2-rK = C2 - rkG = C2 - krG = C2 - kC1
                对M解码即得到明文。

4.椭圆曲线Ep(a,b)的生成

由于我们要使用的是有限元的加法循环群，所以我们要求定义的椭圆曲线上只有整数点，故使用mod p运算定义曲线。
生成一个大于3的素数p
7定义椭圆曲线Ep(a,b)：选取a,b属于mod p 的剩余系，有 (4a^3+27b^2)mod p ≠ 0, 则有Ep(a,b) = ｛（x，y）|y^2=x^3+ax+b，x，y∈Zp｝∪{O}
与上面的实数域的定义方式相同，只是挪到了mod p 的整数域上。这样的曲线被成为素曲线。

5.ECC算法的安全性
与EIG相同，ECC算法同样基于离散对数难题。这要求我们使用足够大的素数p，以保证我们选择的K 构成的群有足够大的阶，以保证离散对数表的复杂度。
所以，在生成椭圆曲线之后，我们需要对曲线的阶进行验证，以保证加密算法的安全性。一般我们使用SEA（ Schoof-Elkies-Atkin，不是图像处理的SEA算法）算法对椭圆曲线的阶进行求解。这是一直复杂度很高的算法，需要进行大量的计算。同时我们要求生成的椭圆曲线的阶不能被多个小质数分解，同时随机选择的基点G具有大素数阶，即保证我们算法使用的在该椭圆曲线上的子群应该有足够大的阶，这同样需要对这个子群的阶进行计算。
由于曲线的生成需要大量的验证运算，所以我们平时使用时并不自己生成椭圆曲线，而是选择一条已经经过验证的曲线使用。

6.明文嵌入
在对明文使用ECC加密前，需要将明文m编码为曲线Ep(a,b)上的一个点M。
1.计算点Q = rK = (qx，qy)
2.设将发送的两个明文分组分别为m1，m2。
c1 = rG
c2 = M+rK = (m1，m2)+(qx，qy) =  (m1*qx) mod p，(m2*qy)mod p （即点的加法转换为向量的乘法，实际上加密后的点并不在椭圆曲线上，但解密后依然具有等价性。）
解密过程：使用私钥k计算Q = rkG 得到qx，qy。
m1 =  （(m1*qx) mod p）* qx^-1 mod p，
m2 = （(m2*qy)mod p）*qy^-1 mod p


7.算法实现

def ecadd(P, Q):                    #椭圆加法
    if P == (0, 0):                   #其中一个点是单位元的情况
        return Q
    if Q == (0, 0):
        return P

    (x1, y1) = P
    (x2, y2) = Q

    if x1 == x2 and y1 == y2:                              
        m = (3 * x1 * x1 + a) * pow(2 * y1, -1, p) % p         # P、Q重合的情况  切线斜率  
    else:
        m = (y2 - y1) * pow(x2 - x1, -1, p) % p                    # P ≠ Q

    x3 = (m * m - x1 - x2) % p
    y3 = (m * (x1 - x3) - y1) % p

    return (x3, y3)

def ec_mul(k, P):                                          #椭圆上的标量乘法，即一个整数乘以一个点，即K = kG ,倍点的计算
    R = (0, 0)
    for i in range(k.bit_length()):                     #使用二进制长度计算数值
        if (k >> i) & 1:
            R = ec_add(R, P)                               #k 个P 相加，即多次调用点的加法
        P = ec_add(P, P)
    return R

def generate_key_pair():                                 #生成密钥
    k = int.from_bytes(os.urandom(32), byteorder='big') % n                        #直接调用python库生成随机的私钥k,这里的n是所选曲线基点G的阶数。
    K = ec_mul(k, G)                                       #公钥K = kG
    return k, K

def encrypt(public_key, plaintext):                    #加密
    r = int.from_bytes(os.urandom(32), byteorder='big') % n                   #生成随机数 r
    C1 = ec_mul(r, G)                                                                           # 计算 C1 = rG    C1 = ec_mul(r, G)
    Q=ec_mul(r, K)                                                                                   # 获取 Q 的 x 坐标
    q_x = Q[0]
    q_y = Q[1]
    m1, m2 = plaintext
    m1_prime = (m1 * q_x) % p                                                            #对明文编码加密
    m2_prime = (m2 * q_y) % p
    return C1, m1_prime, m2_prime

def decrypt(private_key, C1, m1_prime, m2_prime):                               #解密过程，输入私钥，C1，密文1，密文2

    Q = ec_mul(private_key, C1)                                                            # 计算 Q = rK
    q_x = Q[0]                                                                                        # 获取 Q 的 x 坐标
    q_y = Q[1]                                                                                    
    q_x_inv = pow(q_x, -1, p)                                                                # 计算 q_x 的逆元
    q_y_inv = pow(q_y, -1, p)
    m1 = (m1_prime * q_x_inv) % p                                                      # 恢复 m1 = m1' * q_x^-1 % p
    m2 = (m2_prime * q_y_inv) % p                                                      # 恢复 m2 = m2' * q_x^-1 % p

    return m1, m2