sm2椭圆曲线公钥密码算法

背景

国密算法介绍

国密即国家密码局认定的国产密码算法。主要有SM1，SM2，SM3，SM4。密钥长度和分组长度均为128位。

SM1 为对称加密。其加密强度与AES相当。该算法不公开，调用该算法时，需要通过加密芯片的接口进行调用。
SM2为非对称加密，基于ECC。该算法已公开。由于该算法基于ECC，故其签名速度与秘钥生成速度都快于RSA。ECC 256位（SM2采用的就是ECC 256位的一种）安全强度比RSA 2048位高，但运算速度快于RSA。
SM3 消息摘要。可以用MD5作为对比理解。该算法已公开。校验结果为256位。
SM4 无线局域网标准的分组数据算法。对称加密，密钥长度和分组长度均为128位。

由于SM1、SM4加解密的分组大小为128bit，故对消息进行加解密时，若消息长度过长，需要进行分组，要消息长度不足，则要进行填充。

SM2算法简介

SM2椭圆曲线公钥密码算法是我国自主设计的公钥密码算法，包括SM2-1椭圆曲线数字签名算法，SM2-2椭圆曲线密钥交换协议，SM2-3椭圆曲线公钥加密算法，分别用于实现数字签名密钥协商和数据加密等功能。SM2算法与RSA算法不同的是，SM2算法是基于椭圆曲线上点群离散对数难题，相对于RSA算法，256位的SM2密码强度已经比2048位的RSA密码强度要高。

sm2主要满足电子认证服务系统等应用需求。

椭圆曲线介绍

参考书籍

基本数学形式

椭圆曲线是下面方程的一组解： \[ Y^2 = X^3 + AX + B \] 这种类型的方程也被称为魏尔斯特拉斯方程，得名于19世纪对其进行广泛研究的数学家。下面给两个曲线的实例并画给出图片： \[ E1 : Y^2 = X^3 − 3X + 3 \;and\; E2 : Y ^2 = X^3 − 6X + 5 \]

实数域上椭圆曲线

加法运算

令E为如下椭圆曲线： \[ Y ^2 = X^3 − 15X + 18. \] 点P = (7, 16) 和 Q = (1, 2)为椭圆曲线上的两点，并构成直线L: \[ L : Y = \frac{7}{3}X − 1/3. \] 如下图所示，直线L和椭圆曲线E交于三个点P,Q,R

点R关于x轴对称得到R',我们定义椭圆曲线d上的加法运算, \[ P \oplus Q = R', \] 和我们通常理解的加法不同，这里的加法是在曲线上进行的，也就是说给出两个点，相加一定可以得到椭圆曲线上的第三个点R'，这个点就是加法的结果。

加法运算瓶颈——无穷远点O

还是按照上面的图片，我们尝试做这样的加法 \[ R \oplus R'=? \] 按照上一节的加法运算，直线RR‘应该和椭圆曲线交于三个点，但是这里没有第三个点，这时候该怎么办？数学家定义了一个无穷远点$O$,并假设点$O$也在椭圆曲线上，这样RR'就能交于椭圆曲线的点$O$了,于是有 \[ R \oplus R'=O. \] 基于无穷远点有这样的性质

无穷远点$O$和椭圆曲线上任意点P的连线一定是垂直于x轴的

结合之前的加法运算，于是有公式： \[ P \oplus O = P \] 这是不是很类似于一个零点，任何点加这个点都是其本身

除此之外，基于无穷远点还有如下性质 \[ P + O = O + P = P \; for \;all\; P ∈ E\\ P + (−P) = O \; for \;all\; P ∈ E\\ (P + Q) + R = P + (Q + R) \; for \;all\;P, Q, R ∈ E\\ P + Q = Q + P \; for \;all\; P, Q ∈ E \] 四条公式分别代表Albel群(交换群)的4条性质：存在零元，存在逆元，结合律，交换律

加法运算公式

这里给一张图片，很好的讲述了加法的运算公式

减法运算公式

定义负数元的概念：如果椭圆曲线上的一个点P=(a,b)，那么点P的负数元为$\ominus P = (a,-b)$

这样定义有很好的性质：

$P\ominus P = O$

相当于实数域加减法的相反数的概念。

有限域上椭圆曲线

有限域上椭圆曲线定义

从上一节的介绍上可以看到，椭圆曲线上的加法运算构成一个群，所有的点可以是小数，也可以整数。那我们能不能加以限制，令椭圆曲线上的点全部都集中在某个域中呢？

定义1.令$p$是一个奇素数，在有限域$\mathbb{F}_p$上的椭圆曲线有如下形式, \[ E : Y ^2 = X^3 + AX + B \;with\; A, B \in \mathbb{F}_p \;satisfying\; 4A^3 + 27B^2 \neq 0 \] 这里$4A^3 + 27B^2 \neq 0$是保证椭圆曲线上没有奇异点。

则定义在有限域上的椭圆曲线坐标是集合. \[ E(\mathbb{F}_p) = \{(x, y) | x, y \in \mathbb{F}_p \; satisfying \; y^2 = x^3 + Ax + B\}∪ \{O\} \]

一个例子,考虑下面的椭圆曲线： \[ E(\mathbb{F}_{13}) : Y ^2 = X^3 + 3X + 8 \;over \;the \;field \;\mathbb{F}_{13} \] 我们取X=0,得到$Y^2=8\pmod{13}$，但我们解不出这个Y，因为方程无解

我们再取X = 1,得到$Y^2=12\pmod{13}$，解出两个解Y=5或者Y=8，那么得到在$E(\mathbb{F}_{13})$有两个点(1,5)和(1,8)

X遍历有限域$\mathbb{F}_{13}$,我们可以用类似的方法得到下面所有点 \[ E(\mathbb{F}_{13}) = \{O,(1, 5),(1, 8),(2, 3),(2, 10),(9, 6),(9, 7),(12, 2),(12, 11)\} \]

有限域上椭圆曲线运算

类似与实数域上的加法运算，有限域上椭圆曲线的运算，除法用模逆运算代替，所有运算均在有限域上进行（算完之后要模一个p）。这里还是举个例子

椭圆曲线： \[ E(\mathbb{F}_p) = \{(x, y) | x, y \in \mathbb{F}_p \; satisfying \; y^2 = x^3 + Ax + B\}∪ \{O\} \] 计算P=(9,7)和Q=(1,8)的和 $$ λ = (y_2 − y_1)/(x_2 − x_1) = (8 − 7)/(1 − 9) = 1/(-8) = 1/5 = 5^{-1} =8\

x_3 = λ^2 − x_1 − x_2 = 64 − 9 − 1 = 54 = 2\ y_3 = λ(x_1-x_3)-y_1 = 8*(9-2) - 7 = 49 = 10 $$ 最终得到P+Q=(2,10)

椭圆曲线的一些名词

椭圆曲线的阶

我们之前说到每个在有限域上的椭圆曲线都由有限个点组成。那么我们不禁要问：到底是多少个点？

首先，我们要定义一下在一个群有多少个点就叫做这个群的“阶”（order）【在此放上wiki关于order的解释】。

椭圆曲线$E(\mathbb{F}_{13}) : Y ^2 = X^3 + 3X + 8 \;over \;the \;field \;\mathbb{F}_{13}$的阶为9，因为总共有9个点

椭圆曲线的基点

先介绍一下循环子群的概率，方便我们之后基点的讨论

数乘和循环子群

在实数域乘法的定义是： \[ nP=\underbrace{P+P+⋯+P}_{n个p} \] 使用倍加运算可以实现，其时间复杂度我们之后再讨论

举个例子：

已知椭圆曲线 \[ E(\mathbb{F}_p) = \{(x, y) | x, y \in \mathbb{F}_p \; satisfying \; y^2 = x^3 + 3x + 8\}∪ \{O\} \] 取其上一点P(9,7)，其倍点运算为：

P = (9,7)

2P = (9,6)

3P = (1,12)

4P = (9,7)

可以看出，点P运算4次还是点P，这就是循环子群的概念，可以看出点P的运算结果可以构成一个有限的集合，这种群算方式可以看成群的运算，从而构成循环子群。

基点的概念

可以看出基点就是椭圆曲线上的一个点，基点有一个阶n，点加运算n次即可得到再次得到基点。

基点的生成

基点可以构成一个循环子群，sm2椭圆加密就是在这个子群上进行运算。

假设椭圆曲线的阶为N，对于椭圆曲线上的每一个点，我们都有$NP = O$，同时基点G（其阶为n）也能构成一个元素个数的为n的子群，由群论的拉格朗日定理，一定有$n|N$，我们设置一个辅因子$h = N/n$，随机取椭圆曲线上的一点，可以看出点hP循环n次就是无穷远点，如果n为素数，那么点hp生成的子群阶就是n \[ n*h*p = O \] 通过下面的运算步骤，我们可以寻找到阶为n的椭圆曲线下的子群。

计算椭圆曲线的阶 N 。

选择一个阶为 n 的子群。n必须是素数且必须是 N 的因子。

计算辅因子 h=N/n 。

在曲线上选择一个随机的点 P 。

计算 G=hP 。

如果 G 是0，那么回到步骤4。否则我们就已经找到了阶为 n 和辅因子是 h 的子群的基点。

sm2国密算法流程介绍

最详尽的算法流程步骤请参见这个pdf：sm2技术文档，本人在这里仅仅简单介绍一下流程和步骤，并讲解其中原理。其中流程中的一些符号不给予解释，文档里都有。

加密算法

A1：用随机数发生器产生随机数k∈[1,n-1]；
A2：计算椭圆曲线点C1=[k]G=(x1,y1)，按本文本第1部分4.2.8和4.2.4给出的细节，将C1的数据类型转换为比特串；
A3：计算椭圆曲线点S=[h]PB，若S是无穷远点，则报错并退出；
A4：计算椭圆曲线点[k]PB=(x2,y2)，将坐标x2、y2 的数据类型转换为比特串；
A5：计算t=KDF(x2 ∥ y2, klen)，若t为全0比特串，则返回A1；
A6：计算C2 = M ⊕ t；
A7：计算C3 = Hash(x2 ∥ M ∥ y2)；
A8：输出密文C = C1 ∥ C2 ∥ C3。

其中KDF为密钥派生函数，函数原型为KDF(z,klen)可以根据比特串z，和长度klen，从而输出一个长度为klen的比特串。

下面给一个流程图

加密算法原理的一些理解

为什么随机数k范围是[1,n-1]，因为如果k=n，那么[k]G = O，这样的点G是无法用坐标表示出来的，也无法参加后面的运算
真正的密文在C2里，想要获取C2，必须用密钥派生函数KDF得到比特串t的值，而想要得到t的值，又必须计算$[k]P_B$，但是k是随机生成的，想要准确的获取k有两种方法。一种是从$[1,n-1]$遍历，要知道这样的遍历需要计算$n*(n-1)/2$次加法运算，n一般是个很大的数（数量级为$10^{78}$），这种计算量算到宇宙毁灭也算不出来。另一种是解方程$[x]G=(x_1,y_1)$，这就回到椭圆曲线上的困难问题ECDLP，所以想要解密是十分的困难。
C2的长度为klen，即和密文的长度想当
C3是哈希函数，更多的是为验证接受到的密文是否出现改动

解密算法

设klen为密文中C2的比特长度。为了对密文C=C1 ∥ C2 ∥ C3 进行解密，作为解密者的用户B应实现以下运算步骤：

B1：从C中取出比特串C1，将C1的数据类型转换为椭圆曲线上的点，验证C1是否满足椭圆曲线方程，若不满足则报错并退出；
B2：计算椭圆曲线点S=[h]C1，若S是无穷远点，则报错并退出；
B3：计算[dB]C1=(x2,y2)，将坐标x2、y2的数据类型转换为比特串；
B4：计算t=KDF(x2 ∥ y2, klen)，若t为全0比特串，则报错并退出；
B5：从C中取出比特串C2，计算M′ = C2 ⊕ t；
B6：计算u = Hash(x2 ∥ M′ ∥ y2)，从C中取出比特串C3，若$u \neq C3$则报错并退出；
B7：输出明文M′。

还是给出一个流程图

解密算法的一些理解

为什么计算$[d_B]C_1=[d_B*k]G$即可得到加密流程中的(x2,y2)，要知道只有$[k]P_B=(x2,y2)$，理由是：

密钥对的设定是有规则的，即$[d_B]G=P_B$，所以这里可以等价。试想一下，想要根据公钥PB得到私钥$d_B$，又需要解一个ECDLP问题，sm2的安全性就基于此。

Python程序实现

sm2算法的实现有专门的库gmssl，调用其中的库即可实现加解密。

这里仿照gmssl库做出了一些改进，实现了sm2的加解密。

其中椭圆曲线参数的设定和基点的选择参照官方技术文档

基点G的阶:     FFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFF7203DF6B21C6052B53BBF40939D54123
素数p:        FFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF00000000FFFFFFFFFFFFFFFF
基点G横坐标:   32c4ae2c1f1981195f9904466a39c9948fe30bbff2660be1715a4589334c74c7
基点G纵坐标:   bc3736a2f4f6779c59bdcee36b692153d0a9877cc62a474002df32e52139f0a0
椭圆曲线系数a: FFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF00000000FFFFFFFFFFFFFFFC
椭圆曲线系数b: 28E9FA9E9D9F5E344D5A9E4BCF6509A7F39789F515AB8F92DDBCBD414D940E93

代码实现如下

import binascii
from random import choice
from gmssl import sm3,func
import base64
# 选择素域，设置椭圆曲线参数

default_ecc_table = \
{
    'n': 'FFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFF7203DF6B21C6052B53BBF40939D54123',
    'p': 'FFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF00000000FFFFFFFFFFFFFFFF',
    'g_x': '32c4ae2c1f1981195f9904466a39c9948fe30bbff2660be1715a4589334c74c7',
    'g_y': 'bc3736a2f4f6779c59bdcee36b692153d0a9877cc62a474002df32e52139f0a0',
    'a': 'FFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF00000000FFFFFFFFFFFFFFFC',
    'b': '28E9FA9E9D9F5E344D5A9E4BCF6509A7F39789F515AB8F92DDBCBD414D940E93',
}

# default_ecc_table = \
# {
#     'n': '1',
#     'p': 'D',
#     'g_x': '1',
#     'g_y': '1',
#     'a': '3',
#     'b': '8',
# }

class CryptSM2(object):

    def __init__(self, private_key, public_key, ecc_table=default_ecc_table):
        #初始化函数，需要输入公钥和私钥
        self.private_key = private_key
        self.para_len = len(ecc_table['n'])
        self.public_key = self.Str_coordinate_to_jacobian(public_key)
        self.ecc_a3 = (
            int(ecc_table['a'], base=16) + 3) % int(ecc_table['p'], base=16)
        self.ecc_table = ecc_table

    def double_point(self, Point):  # 倍点
        l = len(''.join(Point)) #计算点的16进制长度 
        len_2 = 2 * self.para_len
        x1 = int(Point[0], 16)
        y1 = int(Point[1], 16)
        z1 = int(Point[2], 16)

        p = int(self.ecc_table['p'], base=16)
        #使用计算倍点的公式
        T6 = (z1 * z1) % p
        T2 = (y1 * y1) % p
        T3 = (x1 + T6) % p
        T4 = (x1 - T6) % p
        T1 = (T3 * T4) % p
        T3 = (y1 * z1) % p
        T4 = (T2 * 8) % p
        T5 = (x1 * T4) % p
        T1 = (T1 * 3) % p
        T6 = (T6 * T6) % p
        T6 = (self.ecc_a3 * T6) % p
        T1 = (T1 + T6) % p
        z3 = (T3 + T3) % p
        T3 = (T1 * T1) % p
        T2 = (T2 * T4) % p
        x3 = (T3 - T5) % p

        if (T5 % 2) == 1:
            T4 = (T5 + ((T5 + p) >> 1) - T3) % p
        else:
            T4 = (T5 + (T5 >> 1) - T3) % p

        T1 = (T1 * T4) % p
        y3 = (T1 - T2) % p

        form = '%%0%dx' % self.para_len
        return (form % x3, form % y3, form % z3)

    def add_point(self, P1, P2):  # 点加函数，P2点为仿射坐标即z=1，P1为Jacobian加重射影坐标
        p = int(self.ecc_table['p'], base=16)
        X1 = int(P1[0], 16)
        Y1 = int(P1[1], 16)
        Z1 = int(P1[2],16)

        x2 = int(P2[0], 16)
        y2 = int(P2[1], 16)

        T1 = (Z1 * Z1) % p
        T2 = (y2 * Z1) % p
        T3 = (x2 * T1) % p
        T1 = (T1 * T2) % p
        T2 = (T3 - X1) % p
        T3 = (T3 + X1) % p
        T4 = (T2 * T2) % p
        T1 = (T1 - Y1) % p
        Z3 = (Z1 * T2) % p
        T2 = (T2 * T4) % p
        T3 = (T3 * T4) % p
        T5 = (T1 * T1) % p
        T4 = (X1 * T4) % p
        X3 = (T5 - T3) % p
        T2 = (Y1 * T2) % p
        T3 = (T4 - X3) % p
        T1 = (T1 * T3) % p
        Y3 = (T1 - T2) % p


        form = '%%0%dx' % self.para_len
        
        return (form % X3, form % Y3, form % Z3)

    def kp(self, k_int, Point_xy):  
        # kP运算，即k倍点的运算函数
        Point = Point_xy
        k = k_int
        mask_str = '8'
        for i in range(self.para_len - 1):
            mask_str += '0'
        mask = int(mask_str, 16)
        Temp = Point
        flag = False
        for n in range(self.para_len * 4): #每个16进制4个bit
            if (flag):
                Temp = self.double_point(Temp)
            if (k & mask) != 0: #用与操作判断k的最高位是否为0
                if (flag):
                    Temp = self.add_point(Temp, Point)
                else:
                    flag = True
                    Temp = Point
            k = k << 1
        return self._convert_jacb_to_nor(Temp)
    
    def _convert_jacb_to_nor(self, Point): # Jacobian加重射影坐标转换成仿射坐标
        len_2 = 2 * self.para_len
        x = int(Point[0], 16)
        y = int(Point[1], 16)
        z = int(Point[2], 16)
        p = int(self.ecc_table['p'], base=16)
        z_inv = pow(z, p - 2, p)
        z_invSquar = (z_inv * z_inv) % p
        z_invQube = (z_invSquar * z_inv) % p
        x_new = (x * z_invSquar) % p
        y_new = (y * z_invQube) % p
        z_new = (z * z_inv) % p
        if z_new == 1:
            form = '%%0%dx' % self.para_len
            
            return (form % x_new, form % y_new, form % z_new)
        else:
            return None

    def encrypt(self, data_str):
        # 加密函数，data消息字符串
        data = data_str.encode() #将消息转化为bytes流
        msg = data.hex() # 消息转化为16进制字符串
        G = (self.ecc_table['g_x'],self.ecc_table['g_y'],"1")
        k = 2#int(func.random_hex(self.para_len),16) #1.产生随机数k \in [1,n-1]
        C1 = self.kp(k,G) #2.计算[k]G = (x1,y1)
        C1 = self._convert_jacb_to_nor(C1)
        C1 = C1[0]+C1[1]
        xy = self.kp(k,self.public_key) #3.计算点s = [k]pk
        x2 = xy[0]
        y2 = xy[1]
        m_len = len(msg)
        t = sm3.sm3_kdf((x2+y2).encode('utf8'), m_len/2)#5.计算t = KDF(x2||y2,klen)
        if int(t,16)==0:
            return None
        else:
            form = '%%0%dx' % m_len #两个百分号代表%
            C2 = form % (int(msg, 16) ^ int(t, 16)) #6.计算C2 = M 异或 t，C2的长度理应为消息M的长度
            C3 = sm3.sm3_hash([
                i for i in bytes.fromhex('%s%s%s'% (x2,msg,y2))#7.计算哈希函数C3 = Hash(x2 || M || y2)
            ])
            return bytes.fromhex('%s%s%s' % (C1,C3,C2))

    def decrypt(self, data):
            # 解密函数，data密文（bytes）
            data = data.hex()
            len_2 = 2 * self.para_len
            len_3 = len_2 + 64 
            C1 = self.Str_coordinate_to_jacobian(data[0:len_2]) #1.提取出C1,并转化为坐标点
            C3 = data[len_2:len_3]
            C2 = data[len_3:]
            xyz = self.kp(int(self.private_key,16),C1)#3.计算[sk]C1 = (x2,y2)
            xy = self._convert_jacb_to_nor(xyz)
            # print('xy = %s' % xy)
            x2 = xy[0]
            y2 = xy[1]
            cl = len(C2)
            t = sm3.sm3_kdf((x2+y2).encode('utf8'), cl/2)#4.计算t = KDF(x2||y2,klen)
            if int(t, 16) == 0:
                return None
            else:
                form = '%%0%dx' % cl
                M = form % (int(C2,16) ^ int(t,16))#5.恢复明文M = C2 异或 t
                u = sm3.sm3_hash([
                    i for i in bytes.fromhex('%s%s%s'% (x2,M,y2))
                ])
                return bytes.fromhex(M)
    
    def Str_coordinate_to_jacobian(self,Point_str):
        Point = Point_str
        x = Point[0:self.para_len]
        y = Point[self.para_len:2*self.para_len]
        z = "1"
        return (x,y,z)

if __name__ == "__main__":
    # sm2的公私钥
    SM2_PRIVATE_KEY = '00B9AB0B828FF68872F21A837FC303668428DEA11DCD1B24429D0C99E24EED83D5'
    SM2_PUBLIC_KEY = 'B9C9A6E04E9C91F7BA880429273747D7EF5DDEB0BB2FF6317EB00BEF331A83081A6994B8993F3F5D6EADDDB81872266C87C018FB4162F5AF347B483E24620207'
    operator = CryptSM2(SM2_PRIVATE_KEY,SM2_PUBLIC_KEY)
    c = operator.encrypt("网络空间安全")
    print(base64.b64encode(c))
    result = operator.decrypt(c).decode(encoding="utf-8")
    print(result)

    # n = int(default_ecc_table['n'],16)
    # G = operator.Str_coordinate_to_jacobian(SM2_PRIVATE_KEY)
    # print(operator.kp(n,G))

密码学

#python #椭圆曲线加密算法

sm2椭圆曲线公钥密码算法

http://example.com/2022/12/05/sm2椭圆曲线公钥密码算法/

作者

harper

发布于

2022年12月5日

许可协议

NTT的实现及其加速上一篇

python常用库(密码学) 下一篇