devbox
我家小花儿
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

国密算法 SM9 公钥加密 数字签名 密钥交换 基于身份的密码算法（IBC）完整高效的开源python代码

作者：我家小花儿 | 2024-03-04 06:22:20

赞

踩

sm9

上篇文章（发布于2023-09-18）给自己挖了个坑，说是要搞定SM9。从国庆前一周开始，到现在一个月时间，这个坑终于填上了。此前信息安全数学基础太差，理解不了SM9双线性对、扩域计算等等，为此还特意选修了现代密码学和近世代数2门专业课，边写代码边上课，带着问题学确实收获不少。说实话，弯路比预想的多，但实现效果却出乎意料的好。

SM9原理就不赘述了。目前，互联网上开源的基于Python原生实现的、确保正确（输出数据与《GBT 38635.2-2020 信息安全技术 SM9标识密码算法 第2部分：算法》附录A列举的数据完全一致）的SM9貌似是没有（我没找到）。我参考了以下代码： 

  1. GitHub - gongxian-ding/gmssl-Python: a Python crypto for sm2/sm3/sm4/sm9
  2. GitHub - guanzhi/GmSSL: 支持国密SM2/SM3/SM4/SM9/SSL的密码工具箱
  3. GitHub - funfungho/GmSSL: http://gmssl.org/ v2.5.4
  4. GitHub - GmSSL/GmSSL-Python: Python binding to the GmSSL library
  5. GitHub - yaoyuanyylyy/abestudy: Study project for sm9, IBE, ABE implements with JPBC and BC under IDEA

其中，代码①应该是搬运了国外的IBC算法或有限域等运算的Python代码来尝试实现 SM9，虽然其实现并不完全正确，但给了我最初的重要参考。代码②③是C语言实现的GmSSL库，②是新版（还包含js代码实现的版本）包含更多数学上的优化，③是旧版但跟国标描述的更接近。代码④是代码②在Python下的调用，需要先编译②，我还没尝试过，应该是Python环境下最高效的，但没有基于Python原生实现SM9。代码⑤是SM9的Java实现，还有完整的SM9国标PDF。在GitHub上还能搜到不少，没下载的就不列举了。

本次用Python实现SM9国标描述的公钥加密、数字签名、密钥交换算法，除了保证正确性之外，还对实现细节做了若干数学和代码实现上的优化（包括但不限于：用2-NAF优化Miller循环和常数t模幂、最终模幂分解优化、用共轭实现frobenius映射、将Miller循环中的分子分母分别计算以减少模逆、针对Fq4零元较多而设计稀疏乘法、分圆循环子群Gϕ6(Fp2)下平方运算的优化、用Karatsuba 思想方法减少乘法）。不止是把他人的研究成果写进代码，也有自己独创的优化。另有些文章描述的优化方法在理论上是没问题的，但用Python实际实现却效果不佳，这些都要通过大量的代码测试来甄别。最后在不依赖第三方库的情况下，凝练成总计不到700行代码，我相信高效的代码一定不会太冗长。除了SM9国标，还参考了以下论文，在此对作者表示衷心感谢！

[1] 胡芯忆,何德彪,彭聪,等. 一种SM9算法R-ate对的快速实现方法[J]. 密码学报,2022,9(5):936-948. DOI:10.13868/j.cnki.jcr.000559.

[2] 甘植旺,廖方圆. 国密SM9中R-ate双线性对快速计算[J]. 计算机工程,2019,45(6):171-174. DOI:10.19678/j.issn.1000-3428.0054123.

[3] 王明东,何卫国,李军,等. 国密SM9算法R-ate对计算的优化设计[J]. 通信技术,2020,53(9):2241-2244. DOI:10.3969/j.issn.1002-0802.2020.09.025.

[4] 王江涛,樊荣,黄哲. SM9中高次幂运算的快速实现方法[J]. 计算机工程,2023,49(9):118-124,136. DOI:10.19678/j.issn.1000-3428.0065618.

[5] 付柱. R-ate双线性对密码算法的高效实现[D]. 天津:天津大学,2017.

[6] 孙铭玮. SM9标识密码算法关键技术研究[D]. 黑龙江:哈尔滨理工大学,2022.

[7] 李江峰. SM9算法的研究与FPGA实现[D]. 西安:西安电子科技大学,2021.

总之，非常感谢前辈们的工作，没有这些参考，对于一个非密码专业的网安学生，我手上连一本能讲清SM9具体实现的书都没有，想做出自己的实现更是难上加难，事实上这一个月我也的确走了很多弯路。不像SM2、SM3、SM4、ZUC对着国标或教材的算法描述就能实现，SM9涉及的数学基础要厚实得多，把图书馆里椭圆曲线的书都借遍了，它们和国标描述差不多，数学上都是那么几条式子，可这其中暗藏着巨大的知识鸿沟，让我对代码实现一时无从下手。另一个困难，是手上没有已经正确实现的Python代码（写惯了Python再来看Java会感觉代码太散，看C语言就像是当初学了C看汇编语言的感觉。当然，再难也要啃下来）。做其他国密算法实现的时候，更像是在优化，因为前人已经正确实现了，我只是研究如何把性能往上提，只要哪个步骤输出数据跟参考代码不一致，就说明实现错了，方便及时纠正。但做SM9实现，没有Python下的正确参考，双线性对有那么多复杂步骤，任何一步理解错了或一个粗心，都得不到正确结果，还不知道错在哪…………………………………………直到上周五突然输出了正确结果。再经过近一周优化，达到了满意的效果，及时将代码分享出来。 

废话多了，上源码：

  1. from random import randrange
  2. from math import ceil
  3. from .SM3 import digest as sm3
  4.  
  5. # SM9总则(GB_T 38635.1-2020) A.1 系统参数
  6. q = 0XB640000002A3A6F1D603AB4FF58EC74521F2934B1A7AEEDBE56F9B27E351457D  # 基域特征
  7. N = 0XB640000002A3A6F1D603AB4FF58EC74449F2934B18EA8BEEE56EE19CD69ECF25  # 群的阶
  8. # 群G1的生成元 P1=(x_p1, y_p1)
  9. x_p1 = 0X93DE051D62BF718FF5ED0704487D01D6E1E4086909DC3280E8C4E4817C66DDDD
  10. y_p1 = 0X21FE8DDA4F21E607631065125C395BBC1C1C00CBFA6024350C464CD70A3EA616
  11. # 群G2的生成元 P2=(x_p2, y_p2)
  12. x_p2 = (0X85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141,
  13.         0X3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B)
  14. y_p2 = (0X17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96,
  15.         0XA7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7)
  16. HASH_SIZE = 32  # sm3输出256位（32字节）
  17. N_SIZE = 32  # 阶的字节数
  18. KEY_LEN = 128  # 默认密钥位数
  19. K2_len = 256  # MAC函数中密钥K2的位数
  20.  
  21.  
  22. def to_byte(x, size=None):
  23.     if type(x) is int:
  24.         return x.to_bytes(size if size else ceil(x.bit_length() / 8), byteorder='big')
  25.     elif type(x) in (str, bytes):
  26.         x = x.encode() if type(x) is str else x
  27.         return x[:size] if size and len(x) > size else x  # 超过指定长度，则截取左侧字符
  28.     elif type(x) in (tuple, list):
  29.         return b''.join(to_byte(c, size) for c in x)
  30.     return bytes(x)[:size] if size else bytes(x)
  31.  
  32.  
  33. # 将字节转换为int
  34. def to_int(byte):
  35.     return int.from_bytes(byte, byteorder='big')
  36.  
  37.  
  38. # 广义的欧几里得除法求模逆（耗时约为slow/SM2代码内get_inverse函数的43%）
  39. def mod_inv(a, mod=q):
  40.     if a == 0:
  41.         return 0
  42.     lm, low, hm, high = 1, a % mod, 0, mod
  43.     while low > 1:
  44.         r = high // low
  45.         lm, low, hm, high = hm - lm * r, high - low * r, lm, low
  46.     return lm % mod
  47.  
  48.  
  49. class FQ:
  50.     def __init__(self, n):
  51.         self.n = n
  52.  
  53.     def __add__(self, other):
  54.         return FQ(self.n + other.n)
  55.  
  56.     def __sub__(self, other):
  57.         return FQ(self.n - other.n)
  58.  
  59.     def __mul__(self, other):  # 右操作数可为int
  60.         return FQ(self.n * (other.n if type(other) is FQ else other) % q)
  61.  
  62.     def __truediv__(self, other):  # 右操作数可为int
  63.         return FQ(self.n * mod_inv(other.n if type(other) is FQ else other) % q)
  64.  
  65.     def __pow__(self, other):  # 操作数应为int
  66.         return FQ(pow(self.n, other, q) if other else 1)
  67.  
  68.     def __eq__(self, other):  # 右操作数可为int
  69.         return self.n % q == (other.n if type(other) is FQ else other) % q
  70.  
  71.     def __neg__(self):
  72.         return FQ(-self.n)
  73.  
  74.     def __repr__(self):
  75.         return 'FQ(%064X)' % (self.n % q)
  76.  
  77.     def __bytes__(self):
  78.         return to_byte(self.n % q, N_SIZE)
  79.  
  80.     def is_zero(self):
  81.         return self.n % q == 0
  82.  
  83.     def inv(self):
  84.         return FQ(mod_inv(self.n))
  85.  
  86.     def sqr(self):
  87.         return FQ(self.n * self.n % q)
  88.  
  89.     @classmethod
  90.     def one(cls):
  91.         return cls(1)
  92.  
  93.     @classmethod
  94.     def zero(cls):
  95.         return cls(0)
  96.  
  97.  
  98. class FQ2:
  99.     def __init__(self, *coeffs):  # 国标中的表示是高位在前，而此处coeffs是低位在前
  100.         self.coeffs = coeffs
  101.  
  102.     def __add__(self, other):
  103.         (a0, a1), (b0, b1) = self.coeffs, other.coeffs
  104.         return FQ2(a0 + b0, a1 + b1)
  105.  
  106.     def __sub__(self, other):
  107.         (a0, a1), (b0, b1) = self.coeffs, other.coeffs
  108.         return FQ2(a0 - b0, a1 - b1)
  109.  
  110.     def sqr(self):
  111.         a0, a1 = self.coeffs
  112.         return FQ2((a0 * a0 - (a1 * a1 << 1)) % q, (a0 * a1 << 1) % q)  # (a0^2 - 2 * a1^2, 2 * a0 * a1)
  113.  
  114.     def sqr_u(self):
  115.         a0, a1 = self.coeffs
  116.         return FQ2(-(a0 * a1 << 2) % q, (a0 * a0 - (a1 * a1 << 1)) % q)  # (-4 * a0 * a1, a0^2 - 2 * a1^2)
  117.  
  118.     def mul_b_u(self, b):  # 带参数乘法
  119.         (a0, a1), (b0, b1) = self.coeffs, b.coeffs
  120.         return FQ2(-(a0 * b1 + a1 * b0 << 1) % q, (a0 * b0 - (a1 * b1 << 1)) % q)  # (-2*(a0*b1+a1*b0), a0*b0-2*a1*b1)
  121.  
  122.     def __mul__(self, other):
  123.         if type(other) is int:
  124.             a0, a1 = self.coeffs
  125.             return FQ2(a0 << 1, a1 << 1) if other == 2 else FQ2(a0 * other % q, a1 * other % q)
  126.         (a0, a1), (b0, b1) = self.coeffs, other.coeffs
  127.         a0b0, a1b1 = a0 * b0, a1 * b1  # Karatsuba 思想方法（节约一次乘法），实测此处约有5%提升，用在其他地方未见性能提升
  128.         return FQ2((a0b0 - (a1b1 << 1)) % q, ((a0 + a1) * (b0 + b1) - (a0b0 + a1b1)) % q)  # (a0*b0-2*a1*b1,a0*b1+a1*b0)
  129.  
  130.     def __rmul__(self, other):
  131.         return self.__mul__(other)
  132.  
  133.     def __truediv__(self, other):
  134.         if type(other) is int:
  135.             other_inv = mod_inv(other)
  136.             return FQ2([c * other_inv % q for c in self.coeffs])
  137.         return self * other.inv()
  138.  
  139.     def __pow__(self, other):  # 实际运行此函数的对象都是分圆循环子群Gϕ6(Fp2)中的元素
  140.         if other == 0:
  141.             return self.one()
  142.         t = self
  143.         for ri in bin(other)[3:]:
  144.             t = t.sqr2() * self if ri == '1' else t.sqr2()
  145.         return t
  146.  
  147.     def inv(self):
  148.         a0, a1 = self.coeffs
  149.         if a0 == 0:
  150.             return FQ2(0, -mod_inv(a1 << 1))  # (0, -(2 * a1)^-1)
  151.         if a1 == 0:
  152.             return FQ2(mod_inv(a0), 0)  # (a0^-1, 0)
  153.         k = mod_inv(a0 * a0 + (a1 * a1 << 1))  # k = (a0^2 + 2 * a1^2)^-1
  154.         return FQ2(a0 * k % q, -a1 * k % q)  # (a0 * k, -a1 * k)
  155.  
  156.     def conjugate(self):  # 共轭
  157.         a0, a1 = self.coeffs
  158.         return self.__class__(a0, -a1)
  159.  
  160.     def get_fp_list(self):  # 返回所有基域元素（高位在前）
  161.         if type(self) is FQ2:
  162.             return [i % q for i in self[::-1]]
  163.         return [y for x in self[::-1] for y in x.get_fp_list()] if self.coeffs else [0] * 4  # 注意FQ4对象零值的处理
  164.  
  165.     def __repr__(self):
  166.         return '%s(%s)' % (self.__class__.__name__, ', '.join('%064X' % i for i in self.get_fp_list()))
  167.  
  168.     def __bytes__(self):  # 字节串高位在前
  169.         return to_byte(self.get_fp_list(), N_SIZE)
  170.  
  171.     def __eq__(self, other):
  172.         return self.get_fp_list() == other.get_fp_list()
  173.  
  174.     def __neg__(self):
  175.         return self.__class__(*[-c for c in self.coeffs])
  176.  
  177.     def __getitem__(self, item):
  178.         return self.coeffs[item]
  179.  
  180.     def is_zero(self):
  181.         return all(c % q == 0 for c in self.coeffs) if type(self) is FQ2 else all(c.is_zero() for c in self.coeffs)
  182.  
  183.     @classmethod
  184.     def one(cls):
  185.         return FQ2_one if cls is FQ2 else (FQ12_one if cls is FQ12 else FQ4_one)
  186.  
  187.     @classmethod
  188.     def zero(cls):
  189.         return FQ2_zero if cls is FQ2 else ()
  190.  
  191.  
  192. class FQ4(FQ2):  # 零元的coeffs为空，可优化FQ12稀疏乘法运算
  193.     def __add__(self, other):
  194.         if not self.coeffs:
  195.             return other
  196.         if not other.coeffs:
  197.             return self
  198.         (a0, a1), (b0, b1) = self.coeffs, other.coeffs
  199.         return FQ4(a0 + b0, a1 + b1)
  200.  
  201.     def __sub__(self, other):
  202.         if not self.coeffs:
  203.             return -other
  204.         if not other.coeffs:
  205.             return self
  206.         (a0, a1), (b0, b1) = self.coeffs, other.coeffs
  207.         return FQ4(a0 - b0, a1 - b1)
  208.  
  209.     def sqr(self):
  210.         if not self.coeffs:
  211.             return FQ4_zero
  212.         a0, a1 = self.coeffs
  213.         return FQ4(a0.sqr() + a1.sqr_u(), a0 * a1 * 2)  # (a0^2 + a1^2 * u, 2 * a0 * a1)
  214.  
  215.     def sqr_v(self):
  216.         if not self.coeffs:
  217.             return FQ4_zero
  218.         a0, a1 = self.coeffs
  219.         return FQ4(a0.mul_b_u(a1) * 2, a0.sqr() + a1.sqr_u())  # (2 * a0 * a1 * u, a0^2 + a1^2 * u)
  220.  
  221.     def mul_b_v(self, b):  # 带参数乘法
  222.         if not self.coeffs or not b.coeffs:
  223.             return FQ4_zero
  224.         (a0, a1), (b0, b1) = self.coeffs, b.coeffs
  225.         return FQ4(a0.mul_b_u(b1) + a1.mul_b_u(b0), a0 * b0 + a1.mul_b_u(b1))  # (a0*b1*u+a1*b0*u, a0*b0+a1*b1*u)
  226.  
  227.     def __mul__(self, other):
  228.         if not self.coeffs:
  229.             return FQ4_zero
  230.         if type(other) is int:
  231.             a0, a1 = self.coeffs
  232.             return FQ4(a0 * other, a1 * other)
  233.         if not other.coeffs:
  234.             return FQ4_zero
  235.         (a0, a1), (b0, b1) = self.coeffs, other.coeffs
  236.         return FQ4(a0 * b0 + a1.mul_b_u(b1), a0 * b1 + a1 * b0)  # (a0*b0+a1*b1*u, a0*b1+a1*b0)
  237.  
  238.     def inv(self):
  239.         if not self.coeffs:
  240.             return FQ4_zero
  241.         a0, a1 = self.coeffs
  242.         k = (a1.sqr_u() - a0.sqr()).inv()
  243.         return FQ4((-a0 * k), a1 * k)
  244.  
  245.  
  246. class FQ12(FQ2):
  247.     def __add__(self, other):
  248.         (a0, a1, a2), (b0, b1, b2) = self.coeffs, other.coeffs
  249.         return FQ12(a0 + b0, a1 + b1, a2 + b2)
  250.  
  251.     def __sub__(self, other):
  252.         (a0, a1, a2), (b0, b1, b2) = self.coeffs, other.coeffs
  253.         return FQ12(a0 - b0, a1 - b1, a2 - b2)
  254.  
  255.     def sqr(self):
  256.         a0, a1, a2 = self.coeffs
  257.         return FQ12(a0.sqr() + a1.mul_b_v(a2) * 2, a0 * a1 * 2 + a2.sqr_v(), a0 * a2 * 2 + a1.sqr())
  258.  
  259.     def __mul__(self, other):
  260.         (a0, a1, a2), (b0, b1, b2) = self.coeffs, other.coeffs
  261.         return FQ12(a0 * b0 + a1.mul_b_v(b2) + a2.mul_b_v(b1), a0 * b1 + a1 * b0 + a2.mul_b_v(b2),
  262.                     a0 * b2 + a1 * b1 + a2 * b0)
  263.  
  264.     def sqr2(self):  # 分圆循环子群Gϕ6(Fp2)中的元素平方
  265.         a, b, c = self.coeffs
  266.         a2, b2, c2v = a.sqr(), b.sqr(), c.sqr_v()
  267.         return FQ12(a2 + (a2 - a.conjugate()) * 2, c2v + (c2v + b.conjugate()) * 2, b2 + (b2 - c.conjugate()) * 2)
  268.  
  269.     def pow_t(self):  # 只可用于分圆循环子群Gϕ6(Fp2)中的元素，求逆为共轭
  270.         c, _inv = self, self.frobenius6()
  271.         for ti in t_NAF:
  272.             c = c.sqr2()
  273.             if ti == '1':
  274.                 c = c * self
  275.             elif ti == '2':  # 用2代替-1
  276.                 c = c * _inv
  277.         return c
  278.  
  279.     def inv(self):
  280.         a0, a1, a2 = self.coeffs
  281.         a0_2, a1_2 = a0.sqr(), a1.sqr()
  282.         if a2.is_zero():
  283.             k = (a0 * a0_2 + a1.mul_b_v(a1_2)).inv()
  284.             return FQ12(a0_2 * k, (-a0 * a1 * k), a1_2 * k)
  285.         t0, t1, t2 = a1_2 - a0 * a2, a0 * a1 - a2.sqr_v(), a0_2 - a1.mul_b_v(a2)
  286.         t3 = a2 * (t1.sqr() - t0 * t2).inv()
  287.         return FQ12(t2 * t3, (-t1 * t3), t0 * t3)
  288.  
  289.     def frobenius(self):
  290.         (a0, a1), (b0, b1), (c0, c1) = self.coeffs
  291.         a = FQ4(a0.conjugate(), a1.conjugate() * alpha3)
  292.         b = FQ4(b0.conjugate() * alpha1, b1.conjugate() * alpha4)
  293.         c = FQ4(c0.conjugate() * alpha2, c1.conjugate() * alpha5)
  294.         return FQ12(a, b, c)
  295.  
  296.     def frobenius2(self):
  297.         a, b, c = self.coeffs
  298.         return FQ12(a.conjugate(), b.conjugate() * alpha2, c.conjugate() * alpha4)
  299.  
  300.     def frobenius3(self):
  301.         (a0, a1), (b0, b1), (c0, c1) = self.coeffs
  302.         a = FQ4(a0.conjugate(), -a1.conjugate() * alpha3)
  303.         b = FQ4(b0.conjugate() * alpha3, b1.conjugate())
  304.         c = FQ4(-c0.conjugate(), c1.conjugate() * alpha3)
  305.         return FQ12(a, b, c)
  306.  
  307.     def frobenius6(self):
  308.         a, b, c = self.coeffs
  309.         return FQ12(a.conjugate(), -b.conjugate(), c.conjugate())
  310.  
  311.  
  312. class ECC_Point:
  313.     def __init__(self, *pt):  # 采用Jacobian射影坐标计算，输入仿射坐标后会转换为Jacobian射影坐标
  314.         self.pt = pt if len(pt) == 3 else (*pt, pt[0].one())
  315.  
  316.     @classmethod
  317.     def from_byte(cls, byte):  # 输入bytes类型仿射坐标，构建点对象
  318.         fp_num = len(byte) // (N_SIZE << 1)  # 单个坐标包含的域元素个数
  319.         if fp_num in (1, 2) and len(byte) % N_SIZE == 0:
  320.             fp_list = [to_int(byte[i:i + N_SIZE]) for i in range(0, len(byte), N_SIZE)]  # 将bytes转换为域元素列表
  321.             if fp_num == 1:
  322.                 return cls(FQ(fp_list[0]), FQ(fp_list[1]))
  323.             x_list, y_list = fp_list[fp_num - 1::-1], fp_list[:fp_num - 1:-1]  # 从bytes到FQ2对象保存的域元素，需翻转高低位顺序
  324.             return cls(FQ2(*x_list), FQ2(*y_list))
  325.         return False
  326.  
  327.     def is_inf(self):
  328.         return self[2].is_zero()
  329.  
  330.     def is_on_curve(self):  # 检查点是否满足曲线方程 y^2 == x^3 + b
  331.         x, y, z = self.pt
  332.         return y ** 2 == x ** 3 + (_b1 if type(x) is FQ else _b2) * z ** 6
  333.  
  334.     def double(self):
  335.         x, y, z = self.pt
  336.         T1, _y2 = x.sqr() * 3, y.sqr()
  337.         T2, T3 = x * _y2 * 4, _y2.sqr() * 8
  338.         x3 = T1.sqr() - T2 * 2
  339.         y3 = T1 * (T2 - x3) - T3
  340.         z3 = y * z * 2
  341.         return ECC_Point(x3, y3, z3)
  342.  
  343.     def zero(self):
  344.         cls = self[0].__class__
  345.         return ECC_Point(cls.one(), cls.one(), cls.zero())
  346.  
  347.     def __add__(self, p2):
  348.         if self.is_inf():
  349.             return p2
  350.         if p2.is_inf():
  351.             return self
  352.         (x1, y1, z1), (x2, y2, z2) = self.pt, p2.pt
  353.         z1_2, z2_2 = z1.sqr(), z2.sqr()
  354.         T1, T2 = x1 * z2_2, x2 * z1_2
  355.         T3, T4, T5 = T1 - T2, y1 * z2_2 * z2, y2 * z1_2 * z1
  356.         T6, T7, T3_2 = T4 - T5, T1 + T2, T3.sqr()
  357.         T8, T9 = T4 + T5, T7 * T3_2
  358.         x3 = T6.sqr() - T9
  359.         T10 = T9 - x3 * 2
  360.         y3 = (T10 * T6 - T8 * T3_2 * T3) * _2_inv
  361.         z3 = z1 * z2 * T3
  362.         return ECC_Point(x3, y3, z3)
  363.  
  364.     def multiply(self, n):  # 算法一：二进制展开法
  365.         if n in (0, 1):
  366.             return self if n else self.zero()
  367.         Q = self
  368.         for i in bin(n)[3:]:
  369.             Q = Q.double() + self if i == '1' else Q.double()
  370.         return Q
  371.  
  372.     def __mul__(self, n):  # 算法三：滑动窗法
  373.         k = bin(n)[2:]
  374.         l, r = len(k), 5  # 滑动窗口为5效果较好
  375.         if r >= l:  # 如果窗口大于k的二进制位数，则本算法无意义
  376.             return self.multiply(n)
  377.         P_ = {1: self, 2: self.double()}  # 保存P[j]值的字典
  378.         for i in range(1, 1 << (r - 1)):
  379.             P_[(i << 1) + 1] = P_[(i << 1) - 1] + P_[2]
  380.         t = r
  381.         while k[t - 1] != '1':
  382.             t -= 1
  383.         hj = int(k[:t], 2)
  384.         Q, j = P_[hj], t
  385.         while j < l:
  386.             if k[j] == '0':
  387.                 Q = Q.double()
  388.                 j += 1
  389.             else:
  390.                 t = min(r, l - j)
  391.                 while k[j + t - 1] != '1':
  392.                     t -= 1
  393.                 hj = int(k[j:j + t], 2)
  394.                 Q = Q.multiply(1 << t) + P_[hj]
  395.                 j += t
  396.         return Q
  397.  
  398.     def __rmul__(self, n):
  399.         return self.__mul__(n)
  400.  
  401.     def __eq__(self, p2):
  402.         (x1, y1, z1), (x2, y2, z2) = self.pt, p2.pt
  403.         z1_2, z2_2 = z1.sqr(), z2.sqr()
  404.         return x1 * z2_2 == x2 * z1_2 and y1 * z2_2 * z2 == y2 * z1_2 * z1
  405.  
  406.     def __neg__(self):
  407.         x, y, z = self.pt
  408.         return ECC_Point(x, -y, z)
  409.  
  410.     def __getitem__(self, item):
  411.         return self.pt[item]
  412.  
  413.     def __repr__(self):
  414.         return '%s%s' % (self.__class__.__name__, self.normalize())
  415.  
  416.     def __bytes__(self):
  417.         return to_byte(self.normalize())
  418.  
  419.     def normalize(self):
  420.         x, y, z = self.pt
  421.         if not hasattr(self, 'normalize_tuple'):
  422.             if z != z.one():
  423.                 z_inv = z.inv()
  424.                 z_inv_2 = z_inv.sqr()
  425.                 x, y = x * z_inv_2, y * z_inv_2 * z_inv
  426.             self.normalize_tuple = (x.n, y.n) if type(x) is FQ else (x, y)
  427.         return self.normalize_tuple
  428.  
  429.     def frobenius(self):
  430.         x, y, z = self.pt
  431.         return ECC_Point(x.conjugate(), y.conjugate(), z.conjugate() * alpha1)
  432.  
  433.     def frobenius2_neg(self):
  434.         x, y, z = self.pt
  435.         return ECC_Point(x, -y, z * alpha2)
  436.  
  437.  
  438. FQ2_one, FQ2_zero = FQ2(1, 0), FQ2(0, 0)  # FQ2单位元、零元
  439. FQ4_one, FQ4_zero = FQ4(FQ2_one, FQ2_zero), FQ4()  # FQ4单位元、零元
  440. FQ12_one = FQ12(FQ4_one, FQ4_zero, FQ4_zero)  # FQ12单位元
  441. P1 = ECC_Point(FQ(x_p1), FQ(y_p1))  # 群G1的生成元
  442. P2 = ECC_Point(FQ2(*x_p2[::-1]), FQ2(*y_p2[::-1]))  # 群G2的生成元
  443. _b1, _b2 = FQ(5), FQ2(0, 5)  # b2=βb=(1,0)*5
  444. alpha1 = 0X3F23EA58E5720BDB843C6CFA9C08674947C5C86E0DDD04EDA91D8354377B698B  # -2^((q - 1)/12)
  445. alpha2 = 0XF300000002A3A6F2780272354F8B78F4D5FC11967BE65334  # -2^((q - 1)/6)
  446. alpha3 = 0X6C648DE5DC0A3F2CF55ACC93EE0BAF159F9D411806DC5177F5B21FD3DA24D011  # -2^((q - 1)/4)
  447. alpha4 = 0XF300000002A3A6F2780272354F8B78F4D5FC11967BE65333  # -2^((q - 1)/3)
  448. alpha5 = 0X2D40A38CF6983351711E5F99520347CC57D778A9F8FF4C8A4C949C7FA2A96686
  449. _2_inv = 0X5B2000000151D378EB01D5A7FAC763A290F949A58D3D776DF2B7CD93F1A8A2BF  # 1/2
  450. _3div2 = 0X5B2000000151D378EB01D5A7FAC763A290F949A58D3D776DF2B7CD93F1A8A2C0  # 3/2
  451. R_ate_a_NAF = '00100000000000000000000000000000000000010001020200020200101000020'  # a=6t+2的二进制非相邻表示(2-NAF)（去首1）
  452. t_NAF = '10000000000000000000000000000000000000102020010000201020001010'  # t的二进制非相邻表示(2-NAF)（去首1）
  453. hlen = 320  # 8 * ceil(5 * log(N, 2) / 32)
  454.  
  455.  
  456. # 线函数g T,Q(P)，求过点T和Q的直线在P上的值（分母在最终模幂时值为1，可消去）
  457. def g(T, Q, P=None):
  458.     if P:
  459.         (xT, yT, zT), (xQ, yQ, zQ), (xP, yP) = T, Q, P
  460.         zT_2, zQ_2 = zT.sqr(), zQ.sqr()
  461.         zQ_3, t1 = zQ * zQ_2, (xT * zQ_2 - xQ * zT_2) * zT * zQ
  462.         b1, t2 = t1 * zQ_3, (yT * zQ_3 - yQ * zT * zT_2) * zQ
  463.         a0, a4 = t1 * yQ - t2 * xQ, t2 * zQ_2 * xP
  464.     else:  # 当P为空时，g T,T(P)，求过点T的切线在P上的值
  465.         (xT, yT, zT), (xP, yP) = T, Q
  466.         zT_2, t1 = zT.sqr(), xT.sqr() * _3div2
  467.         b1, a0, a4 = zT * zT_2 * yT, yT.sqr() - t1 * xT, t1 * zT_2 * xP
  468.     return FQ12(FQ4(a0, -b1 * yP), FQ4_zero, FQ4(a4, FQ2_zero))
  469.  
  470.  
  471. # BN曲线上R_ate对的计算
  472. def e(P, Q):
  473.     T, nQ, f, P_xy = Q, -Q, FQ12_one, P.normalize()
  474.     for ai in R_ate_a_NAF:
  475.         f, T = f.sqr() * g(T, P_xy), T.double()
  476.         if ai == '1':
  477.             f, T = f * g(T, Q, P_xy), T + Q
  478.         elif ai == '2':  # 用2代替-1
  479.             f, T = f * g(T, nQ, P_xy), T + nQ
  480.     Q1, nQ2 = Q.frobenius(), Q.frobenius2_neg()
  481.     return final_exp(f * g(T, Q1, P_xy) * g(T + Q1, nQ2, P_xy))
  482.  
  483.  
  484. # 最终模幂
  485. def final_exp(f):
  486.     m = f.frobenius6() * f.inv()  # f^(p^6 - 1)
  487.     s = m.frobenius2() * m  # m^(p^2 + 1)
  488.     # 困难部分 s^(p^3 + (6t^2+1)p^2 + (-36t^3-18t^2-12t+1)p + (-36t^3-30t^2-18t-2))
  489.     s_6t = s.pow_t() ** 6
  490.     s_6t2, s_12t = s_6t.pow_t(), s_6t.sqr2()
  491.     s_6t3, s_12t2, s_18t = s_6t2.pow_t(), s_6t2.sqr2(), s_6t * s_12t
  492.     s_36t3, s_18t2, a2 = s_6t3 ** 6, s_6t2 * s_12t2, s_6t2 * s
  493.     s_30t2, a1 = s_12t2 * s_18t2, (s_36t3 * s_18t2 * s_12t).frobenius6() * s
  494.     a0 = (s_36t3 * s_30t2 * s_18t * s.sqr2()).frobenius6()
  495.     return s.frobenius3() * a2.frobenius2() * a1.frobenius() * a0
  496.  
  497.  
  498. # SM9算法(GB_T 38635.2-2020) 5.3.6定义的密钥派生函数
  499. # Z为bytes类型，klen表示输出密钥比特长度（8的倍数）；输出为bytes类型
  500. def KDF(Z, klen=KEY_LEN):
  501.     ksize, K = klen >> 3, bytearray()
  502.     for ct in range(1, ceil(ksize / HASH_SIZE) + 1):
  503.         K.extend(sm3(Z + to_byte(ct, 4)))
  504.     return K[:ksize]
  505.  
  506.  
  507. # SM9算法(GB_T 38635.2-2020) 5.3.2.2和5.3.2.3定义的密码函数
  508. def H(i, Z):
  509.     Ha = to_int(KDF(to_byte(i, 1) + Z, hlen))
  510.     return Ha % (N - 1) + 1
  511.  
  512.  
  513. # SM9算法(GB_T 38635.2-2020) 5.3.5定义的消息认证码函数
  514. def MAC(K2, Z):
  515.     return sm3(Z + K2)
  516.  
  517.  
  518. class SM9:  # SM9算法(GB_T 38635.2-2020)
  519.     def __init__(self, ID='', ks=None, Ppub_s=None, ke=None, Ppub_e=None, hid_s=1, hid_e=3, is_KGC=False):
  520.         self.ID = ID
  521.         if is_KGC:  # 作为密钥生成中心
  522.             self.ks = ks if ks and 0 < ks % N < N else randrange(1, N)  # 未提供签名主密钥则随机生成
  523.             self.Ppub_s = Ppub_s if Ppub_s and Ppub_s == P2 * self.ks else P2 * self.ks  # 确保签名主公钥与签名主私钥匹配
  524.             self.ke = ke if ke and 0 < ke % N < N else randrange(1, N)  # 未提供加密主密钥则随机生成
  525.             self.Ppub_e = Ppub_e if Ppub_e and Ppub_e == P1 * self.ke else P1 * self.ke  # 确保加密主公钥与加密主私钥匹配
  526.         else:  # 作为用户
  527.             self.ds, self.Ppub_s, self.de, self.Ppub_e = ks, Ppub_s, ke, Ppub_e
  528.             self.gs, self.ge = e(P1, Ppub_s), e(Ppub_e, P2)
  529.         self.hid_s_byte, self.hid_e_byte = [hid if type(hid) is bytes else to_byte(hid, 1) for hid in [hid_s, hid_e]]
  530.  
  531.     def KGC_gen_user(self, ID):
  532.         ID_byte = to_byte(ID)
  533.         t1 = (H(1, ID_byte + self.hid_s_byte) + self.ks) % N
  534.         if t1 == 0:  # 需重新产生签名主密钥，并更新所有用户的签名密钥
  535.             return False
  536.         t2 = self.ks * mod_inv(t1, N) % N
  537.         ds = P1 * t2  # 用户签名私钥
  538.         t1 = (H(1, ID_byte + self.hid_e_byte) + self.ke) % N
  539.         if t1 == 0:  # 需重新产生加密主密钥，并更新所有用户的加密密钥
  540.             return False
  541.         t2 = self.ke * mod_inv(t1, N) % N
  542.         de = P2 * t2  # 用户加密私钥
  543.         return SM9(ID, ds, self.Ppub_s, de, self.Ppub_e, self.hid_s_byte, self.hid_e_byte)
  544.  
  545.     # 6.2 数字签名生成算法
  546.     def sign(self, M, r=None, outbytes=True):
  547.         l = 0
  548.         while l == 0:
  549.             r = r if r else randrange(1, N)  # A2
  550.             w = bytes(self.gs ** r)  # A3
  551.             h = H(2, to_byte(M) + w)  # A4
  552.             l = (r - h) % N  # A5
  553.         S = self.ds * l  # A6
  554.         return to_byte([h, S]) if outbytes else (h, S)
  555.  
  556.     # 6.4 数字签名验证算法
  557.     def verify(self, ID, M_, sig):
  558.         h_, S_ = (to_int(sig[:N_SIZE]), ECC_Point.from_byte(sig[N_SIZE:])) if type(sig) is bytes else sig
  559.         if not 0 < h_ < N or not S_ or not S_.is_on_curve():  # B1、B2
  560.             return False
  561.         t = self.gs ** h_  # B4
  562.         h1 = H(1, to_byte(ID) + self.hid_s_byte)  # B5
  563.         P = P2 * h1 + self.Ppub_s  # B6
  564.         u = e(S_, P)  # B7
  565.         w_ = bytes(u * t)  # B8
  566.         h2 = H(2, to_byte(M_) + w_)  # B9
  567.         return h_ == h2
  568.  
  569.     # A 发起协商（也可用作B生成rB、RB；outbytes=True时输出bytes）
  570.     # 7.2 密钥交换协议 A1-A3
  571.     def agreement_initiate(self, IDB, r=None, outbytes=True):
  572.         QB = P1 * H(1, to_byte(IDB) + self.hid_e_byte) + self.Ppub_e  # A1
  573.         rA = r if r else randrange(1, N)  # A2
  574.         RA = QB * rA  # A3
  575.         return rA, bytes(RA) if outbytes else RA
  576.  
  577.     # B 响应协商（option=True时计算选项部分）
  578.     # 7.2 密钥交换协议 B1-B6
  579.     def agreement_response(self, RA, IDA, option=False, rB=None, klen=KEY_LEN, outbytes=True):
  580.         RA = ECC_Point.from_byte(RA) if type(RA) is bytes else RA
  581.         if not RA or not RA.is_on_curve():  # B4
  582.             return False, 'RA不属于椭圆曲线群G1'
  583.         rB, RB = self.agreement_initiate(IDA, rB, outbytes)  # B1-B3
  584.         g1, g2 = e(RA, self.de), bytes(self.ge ** rB)  # B4
  585.         g1, g3 = bytes(g1), bytes(g1 ** rB)  # B4
  586.         tmp_byte = to_byte([IDA, self.ID, RA, RB])
  587.         SKB = KDF(tmp_byte + g1 + g2 + g3, klen)  # B5
  588.         if not option:
  589.             return True, (RB, SKB)
  590.         self.tmp_byte2 = g1 + sm3(g2 + g3 + tmp_byte)
  591.         SB = sm3(to_byte(0x82, 1) + self.tmp_byte2)  # B6（可选部分）
  592.         return True, (RB, SKB, SB)
  593.  
  594.     # A 协商确认
  595.     # 7.2 密钥交换协议 A5-A8
  596.     def agreement_confirm(self, rA, RA, RB, IDB, SB=None, option=False, klen=KEY_LEN):
  597.         RB = ECC_Point.from_byte(RB) if type(RB) is bytes else RB
  598.         if not RB or not RB.is_on_curve():  # A5
  599.             return False, 'RB不属于椭圆曲线群G1'
  600.         g1_, g2_ = bytes(self.ge ** rA), e(RB, self.de)  # A5
  601.         g2_, g3_ = bytes(g2_), bytes(g2_ ** rA)  # A5
  602.         tmp_byte = to_byte([self.ID, IDB, RA, RB])
  603.         if option and SB:  # A6（可选部分）
  604.             tmp_byte2 = g1_ + sm3(g2_ + g3_ + tmp_byte)
  605.             S1 = sm3(to_byte(0x82, 1) + tmp_byte2)
  606.             if S1 != SB:
  607.                 return False, 'S1 != SB'
  608.         SKA = KDF(tmp_byte + g1_ + g2_ + g3_, klen)  # A7
  609.         if not option or not SB:
  610.             return True, SKA
  611.         SA = sm3(to_byte(0x83, 1) + tmp_byte2)  # A8
  612.         return True, (SKA, SA)
  613.  
  614.     # B 协商确认（可选部分）
  615.     # 7.2 密钥交换协议 B8
  616.     def agreement_confirm2(self, SA):
  617.         if not hasattr(self, 'tmp_byte2'):
  618.             return False, 'step error'
  619.         S2 = sm3(to_byte(0x83, 1) + self.tmp_byte2)
  620.         if S2 == SA:
  621.             del self.tmp_byte2
  622.             return True, ''
  623.         return False, 'S2 != SA'
  624.  
  625.     # 8.2 密钥封装算法
  626.     def encaps(self, IDB, klen, r=None, outbytes=True):
  627.         K = bytes()
  628.         while K == bytes(len(K)):
  629.             r, C = self.agreement_initiate(IDB, r, outbytes)  # A1-A3
  630.             w = bytes(self.ge ** r)  # A5
  631.             K = KDF(to_byte([C, w, IDB]), klen)
  632.         return K, C
  633.  
  634.     # 8.4 密钥封装算法
  635.     def decaps(self, C, klen):
  636.         C = ECC_Point.from_byte(C) if type(C) is bytes else C
  637.         if not C or not C.is_on_curve():  # B1
  638.             return False, 'C不属于椭圆曲线群G1'
  639.         w_ = bytes(e(C, self.de))  # B2
  640.         K_ = KDF(to_byte([C, w_, self.ID]), klen)  # B3
  641.         return (True, K_) if K_ != bytes(len(K_)) else (False, 'K为全0比特串')
  642.  
  643.     # 9.2 加密算法
  644.     def encrypt(self, IDB, M, r=None, outbytes=True):
  645.         M = to_byte(M)
  646.         K, C1 = self.encaps(IDB, (len(M) << 3) + K2_len, r, outbytes)  # A1-A6.a.1
  647.         K1, K2 = K[:len(M)], K[len(M):]  # A6.a.1
  648.         C2 = bytes(M[i] ^ K1[i] for i in range(len(M)))  # A6.a.2
  649.         C3 = MAC(K2, C2)  # A7
  650.         return to_byte([C1, C3, C2]) if outbytes else (C1, C3, C2)
  651.  
  652.     # 9.4 解密算法
  653.     def decrypt(self, C):
  654.         C3_start, C3_end = N_SIZE << 1, (N_SIZE << 1) + HASH_SIZE
  655.         C1, C3, C2 = (C[:C3_start], C[C3_start:C3_end], C[C3_end:]) if type(C) is bytes else C
  656.         res, K_ = self.decaps(C1, (len(C2) << 3) + K2_len)  # B1-B3.a.1
  657.         if not res:
  658.             return False, K_.replace('C', 'C1')
  659.         K1_, K2_ = K_[:len(C2)], K_[len(C2):]  # B3.a.1
  660.         if K1_ == bytes(len(K_)):
  661.             return False, 'K1\'为全0比特串'
  662.         u = MAC(K2_, C2)  # B4
  663.         if u != C3:
  664.             return False, 'u != C3'
  665.         return True, bytes(C2[i] ^ K1_[i] for i in range(len(C2)))  # B3.a.2

使用时先创建KGC对象，由它生成用户对象，再由用户对象完成公钥加密、数字签名或密钥交换，例子如下：

  1.     IDA, IDB, message = 'Alice', 'Bob', 'Chinese IBS standard'
  2.     kgc = SM9(ks=0x130E78459D78545CB54C587E02CF480CE0B66340F319F348A1D5B1F2DC5F4,
  3.               ke=0x2E65B0762D042F51F0D23542B13ED8CFA2E9A0E7206361E013A283905E31F, is_KGC=True)
  4.     sm9_A, sm9_B = kgc.KGC_gen_user(IDA), kgc.KGC_gen_user(IDB)
  5.     assert bytes(sm9_A.gs).hex().swapcase().endswith('F0F071D7D284FCFB')
  6.     print("-----------------test sign and verify---------------")
  7.     r = 0x033C8616B06704813203DFD00965022ED15975C662337AED648835DC4B1CBE
  8.     signature = sm9_A.sign(message, r)
  9.     assert signature.hex().swapcase().endswith('827CC2ACED9BAA05')
  10.     assert sm9_B.verify(IDA, message, signature)
  11.     print("success")
  12.  
  13.     print("-----------------test key agreement---------------")
  14.     rA = 0x5879DD1D51E175946F23B1B41E93BA31C584AE59A426EC1046A4D03B06C8
  15.     rA, RA = sm9_A.agreement_initiate(IDB, rA)  # A发起协商
  16.     # A将RA发送给B
  17.     rB = 0x018B98C44BEF9F8537FB7D071B2C928B3BC65BD3D69E1EEE213564905634FE
  18.     res, content = sm9_B.agreement_response(RA, IDA, True, rB)  # B响应协商
  19.     if not res:
  20.         print('B报告协商错误：', content)
  21.         return
  22.     RB, SKB, SB = content
  23.     # B将RB、SB发送给A
  24.     res, content = sm9_A.agreement_confirm(rA, RA, RB, IDB, SB, True)  # A协商确认
  25.     if not res:
  26.         print('A报告协商错误：', content)
  27.         return
  28.     SKA, SA = content
  29.     assert SKA.hex().swapcase() == '68B20D3077EA6E2B825315836FDBC633'
  30.     # A将SA发送给B
  31.     res, content = sm9_B.agreement_confirm2(SA)  # B协商确认
  32.     if not res:
  33.         print('B报告协商错误：', content)
  34.         return
  35.     assert SKA == SKB
  36.     print("success")
  37.  
  38.     print("-----------------test encrypt and decrypt---------------")
  39.     message = 'Chinese IBE standard'
  40.     kgc = SM9(ks=kgc.ks, Ppub_s=kgc.Ppub_s,
  41.               ke=0x01EDEE3778F441F8DEA3D9FA0ACC4E07EE36C93F9A08618AF4AD85CEDE1C22, is_KGC=True)
  42.     sm9_A, sm9_B = kgc.KGC_gen_user(IDA), kgc.KGC_gen_user(IDB)
  43.     C = sm9_A.encrypt(IDB, message, 0xAAC0541779C8FC45E3E2CB25C12B5D2576B2129AE8BB5EE2CBE5EC9E785C)
  44.     assert C.hex().swapcase().endswith('378CDD5DA9513B1C')
  45.     res, content = sm9_B.decrypt(C)
  46.     if not res:
  47.         print('解密错误：', content)
  48.         return
  49.     assert message == content.decode()
  50.     print("success")

虽然和代码①实现的结果不一样，但手头也只有这一个Python的原生实现，但还是比较一下。

此前介绍国密算法的系列文章如下：

三篇SM2：

国密算法 SM2 公钥加密 数字签名 密钥交换 更高效、依赖更少的开源python代码_国密算法 开源-CSDN博客

国密算法 SM2 公钥加密 数字签名 密钥交换 全网最高效的开源python代码_-CSDN博客

国密算法 SM2 公钥加密 非对称加密 数字签名 密钥协商 python实现完整代码_qq_43339242的博客-CSDN博客_python sm2

SM3：国密算法 SM3 消息摘要 杂凑算法 哈希函数 散列函数 python实现完整代码_qq_43339242的博客-CSDN博客_国密sm3

SM4：国密算法 SM4 对称加密 分组密码 python实现完整代码_qq_43339242的博客-CSDN博客_python国密算法库

ZUC：国密算法 ZUC流密码 祖冲之密码 python代码完整实现_qq_43339242的博客-CSDN博客_国密算法代码

对上述几个算法和实现不了解的，建议点进去看看。下面这篇文章是对上述的汇总：

国密算法 SM2公钥密码 SM3杂凑算法 SM4分组密码 python代码完整实现_qq_43339242的博客-CSDN博客_python sm2

上述所有国密算法的完整Python实现代码和测试代码，库名叫hggm，托管在码云：hggm - 国密算法 SM2 SM3 SM4 python实现完整代码: 国密算法 SM2公钥密码 SM3杂凑算法 SM4分组密码 python代码完整实现 效率高于所有公开的python国密算法库 (gitee.com)

至此，公开国密算法的一家子（SM2、SM3、SM4、SM9和ZUC）已经团聚了，在国密算法Python原生实现领域，目前来看，性能应该都做到了网上开源最优。

咱们所处的世界还很不太平，建设网络安全强国的脚步一刻也不会停，任重道远。

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/我家小花儿/article/detail/186700?site
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号