SM4 AESNI指令集优化

一、SM4简介

SM4 算法于 2012 年被国家密码管理局确定为国家密码行业标准，最初主要用于 WAPI (WLAN Authentication and Privacy Infrastructure) 无线网络中。SM4 算法的出现为将我国商用产品上的密码算法由国际标准替换为国家标准提供了强有力的支撑。随后，SM4 算法被广泛应用于政府办公、公安、银行、税务、电力等信息系统中，其在我国密码行业中占据着及其重要的位置。类似于 DES、AES 算法，SM4 算法也是一种分组密码算法。

SM4的官方文档：点击此处跳转

SM4 算法的原理可以参考该博文：点击此处跳转

二、AES简介

高级加密标准（Advanced Encryption Standard，AES）又称 Rijindael 算法，由比利时著名密码学家 Joan Daemen 和 Vincent Rijimen 设计，是美国联邦政府采用的一种分组加密标准，用来替代之前的 DES 算法，已被多方分析且使用。高级加密标准由美国国家标准与技术研究所在 2001 年 11 月 26 日发布于 FIPS PUB 197，并在 2002 年 5 月 26 日成为有效标准，目前为对称密钥加密算法中最流行的算法之一。

AES 算法的整体结构为 替代——置换网络（Substitution-Permutation Network，SPN）结构，算法的明文长度为 128 位，密钥长度为 128，192 或 256 位，分别被称为 AES-128，AES-192 和 AES-256。

由于 AES 算法网上已经有很多详细的介绍了（比如说知乎上），此处就不过多叙述。
这里提供一个 FIPS PUB 197 的链接，里面是 AES 的官方说明： 点击此处跳转

三、SIMD指令简介

SIMD(Single Instruction Multiple Data)即单指令流多数据流，是一种采用一个控制器来控制多个处理器，同时对一组数据（又称“数据向量”）中的每一个分别执行相同的操作从而实现空间上的并行性的技术。简单来说就是一个指令能够同时处理多个数据。

intel 的 SIMD 指令集参考文档：点击此处跳转

arm 的 SIMD 指令集参考文档：点击此处跳转

本优化的核心指令是AESNI指令集中的_mm_aesenclast_si128指令

四、主要思想

基本思想是利用SM4与AES中S盒结构的相似性，借助intel的AESNI指令完成S盒操作

AES的S盒结构形如 $S_a(x)=A_a.I_a(x)+C_a$

SM4的S盒结构形如 $S_s(x)=A_s.I_s(A_{s}x+C_s)+C_s$

AES使用的不可约多项式是 $x^8+x^4+x^3+x+1$

SM4使用的不可约多项式是 $x^8+x^7+x^6+x^5+x^4+x^2+1$

其中，最耗时的部分为有限域求逆部分 $I(x)$，借助有限域同构的知识，可以构造出同构映射 $T$

（关于同构映射构造原理可参考博文），将SM4对应的有限域元素映射到AES对应的有限域元素中，再借助指令求逆，最后再逆映射不过实现时要特别注意8bit数与矩阵的对应关系，对于$s=b_0{b}_1\cdots b_7$，在SM4中对应的矩阵为 $[b_0,b_1,\cdots ,b_7{]}^T$，但在AES中为 $[b_7,b_6,\cdots ,b_0{]}^T$
关于这点，修改AES矩阵的行列位置即可抵消

4.1 AES与SM4 S盒映射

记 $y=A_{s}x+C_s$
，由SM4有限域映射到AES有限域的函数为 $T(x_s)\to x_a$那么有 $I_s(y)=T^{-1}\cdot I_a(T(y))$

即有

I s ( y ) = T − 1 A a − 1 . { S a ( T ( y ) ) + C a } = T − 1 A a − 1 ⋅ S a ( T ( y ) ) + T − 1 A a − 1 C a S s ( x ) = A s T − 1 A a − 1 ⋅ S a ( T ( y ) ) + A s T − 1 A a − 1 C a + C s = { A s T − 1 A a − 1 } ⋅ S a ( { T A s } x + { T C s } ) + { A s T − 1 A a − 1 C a + C s } \begin{aligned} I_s(y)&=T^{−1}A_a^{−1}.\{S_a(T(y))+C_a\}\\ &=T^{−1}A^{−1}_a⋅S_a(T(y))+T^{−1}A^{−1}_aC_a\\ S_s(x)&=A_sT^{−1}A^{−1}_a⋅S_a(T(y))+A_sT^{−1}A^{−1}_aC_a+C_s\\ &=\{A_sT^{−1}A^{−1}_a\}⋅S_a(\{TA_s\}x+\{TC_s\})+\{A_sT^{−1}A^{−1}_aC_a+C_s\} \end{aligned} Is​(y)Ss​(x)​=T−1Aa−1​.{Sa​(T(y))+Ca​}=T−1Aa−1​⋅Sa​(T(y))+T−1Aa−1​Ca​=As​T−1Aa−1​⋅Sa​(T(y))+As​T−1Aa−1​Ca​+Cs​={As​T−1Aa−1​}⋅Sa​({TAs​}x+{TCs​})+{As​T−1Aa−1​Ca​+Cs​}​

对应矩阵的值可以预计算，在此给出一个示例
{ A s T − 1 A a − 1 } = [ 1 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 0 0 1 1 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 0 1 1 1 0 1 0 1 0 0 1 1 1 0 0 ] , { A s T − 1 A a − 1 C a + C s } = [ 0 0 1 1 1 0 1 1 ] { T A s } = [ 0 0 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 1 1 1 0 0 1 0 0 0 1 1 1 0 1 0 0 0 1 1 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 1 1 0 0 0 0 0 1 1 0 ] , { T C s } = [ 0 0 1 0 0 0 1 1 ] \{A_sT^{−1}A^{−1}_a\}= \begin{bmatrix} 1&0&1&0&1&0&0&0\\ 0&1&1&0&0&0&0&1\\ 1&1&0&0&0&0&1&1\\ 0&1&1&1&0&1&0&0\\ 1&1&0&0&0&1&0&0\\ 1&0&0&0&1&1&0&0\\ 0&0&1&1&1&0&1&0\\ 1&0&0&1&1&1&0&0\\ \end{bmatrix},\{A_sT^{−1}A^{−1}_aC_a+C_s\}= \begin{bmatrix} 0\\ 0\\ 1\\ 1\\ 1\\ 0\\ 1\\ 1\\ \end{bmatrix}\\ \{TA_s\}= \begin{bmatrix} 0&0&0&0&1&1&0&1\\ 1&0&0&1&1&0&1&1\\ 0&1&1&1&0&0&1&0\\ 0&0&1&1&1&0&1&0\\ 0&0&1&1&0&1&0&1\\ 0&0&0&0&1&0&1&0\\ 0&0&0&1&0&1&1&1\\ 0&0&0&0&0&1&1&0\\ \end{bmatrix},\{TC_s\}= \begin{bmatrix} 0\\ 0\\ 1\\ 0\\ 0\\ 0\\ 1\\ 1\\ \end{bmatrix} {As​T−1Aa−1​}=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​10101101​01111000​11010010​00010011​10000111​00011101​00100010​01100000​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​,{As​T−1Aa−1​Ca​+Cs​}=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​00111011​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​{TAs​}=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​01000000​00100000​00111000​01111010​11010100​10001011​01110111​11001010​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​,{TCs​}=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​00100011​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​

其中 $x=b_7{b}_6\cdots b_0$对应向量为 $[b_7,b_6,\cdots ,b_0{]}^T$

4.2 矩阵乘法优化

利用vpshufb指令可以实现4bit的快速查表运算，可以将矩阵乘法拆分成两个4bit的查表
$$A\times x=\left[\begin{array}{cc}{A}_{0,0}& A_{0,1}\\ A_{1,0}& A_{1,1}\end{array}\right]\times \left[\begin{array}{c}{b}_7\\ b_6\\ b_5\\ b_4\\ b_3\\ b_2\\ b_1\\ b_0\end{array}\right]=\left[\begin{array}{cc}{A}_{0,0}& A_{0,1}\\ A_{1,0}& A_{1,1}\end{array}\right]\times \left(\left[\begin{array}{c}{b}_7\\ b_6\\ b_5\\ b_4\\ 0\\ 0\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ b_3\\ b_2\\ b_1\\ b_0\end{array}\right]\right)$$
即，高4bit对应一个 4→8查找表，低4bit对应一个 4→8
查找表。两次查表结果异或就是最终乘法的结果

4.3 算法流程

因为AESNI指令操作的是128bit的数据，SM4一组消息每轮需要查表的数据仅有32bit，故将4组消息打包至一块，使得处理数据内容达到128bit

1. 输入消息$M_0,M_1,M_2,M_3$
2. $X_0,X_1,X_2,X_3\leftarrow$分组打包$(M_0,M_1,M_2,M_3)$
3. $fori=0\to 31$
4. $S\leftarrow X_1\oplus X_2\oplus X_3\oplus K_i$
5. $S\leftarrow TA\times S+TC$
6. $S\leftarrow AESS盒(S)$
7. $S\leftarrow ATA\times S+ATAC$
8. $S\leftarrow X_0\oplus L(S)$
9. $X_0,X_1,X_2,X_3\leftarrow X_1,X_2,X_3,S$
10. $M_0,M_1,M_2,M_3\leftarrow 分组解包(X_3,X_2,X_1,X_0)$
11. $输出M_0,M_1,M_2,M_3$

第 5-7 步利用AESNI指令完成了SM4的S盒查表操作
在实现时，使用的是_mm_aesenclast_si128指令，其包括行移位、字节代替（S盒）以及轮密钥加三步。为了消除行移位和轮密钥加的影响，使用逆行移位和全零轮密钥来抵消

五、代码实现

代码使用的是AESNI和SSE指令集，向量寄存器宽度为128位。主要为了演示思想，若对效率要求较高，可使用AVX指令集，其支持的向量寄存器宽度为256位

如果采用命令行编译，需要添加命令行参数-msse4和-maes
编译命令为gcc -o main.exe main.c sm4_aesni_x4.c sm4.c -Og -msse4 -maes

在sm4.h中定义了SM4密钥生成算法，sm4_aesni_x4.h中定义了利用AESNI指令集实现的4分组消息加解密，main.c为测试样例，sm4.c和sm4_aesni_x4.c中是对应函数的实现

sm4.h

#ifndef SM4_H
#define SM4_H

#include <stdint.h>

/**
 * @brief SM4 密钥
 */
typedef struct _SM4_Key {
    uint32_t rk[32];//32轮密钥
} SM4_Key;

/**
 * @brief 初始化 SM4 轮密钥
 * @param key 128bit长度密钥
 * @param sm4_key SM4 密钥
 */
void SM4_KeyInit(uint8_t* key, SM4_Key* sm4_key);

#endif // !SM4_H

sm4_aesni_x4.h

#ifndef SM4_AESNI_X4_H
#define SM4_AESNI_X4_H

#include"sm4.h"

void SM4_AESNI_Encrypt_x4(uint8_t* plaintext, uint8_t* ciphertext, SM4_Key* sm4_key);

void SM4_AESNI_Decrypt_x4(uint8_t* ciphertext, uint8_t* plaintext, SM4_Key* sm4_key);

#endif // !SM4_AESNI_X4_H

main.c

#include <stdio.h>
#include <stdlib.h>

#include "sm4_aesni_x4.h"

int main() {
    // 01 23 45 67 89 ab cd ef fe dc ba 98 76 54 32 10
    unsigned char key[16 * 8] = {0x01, 0x23, 0x45, 0x67, 0x89, 0xab,
                                 0xcd, 0xef, 0xfe, 0xdc, 0xba, 0x98,
                                 0x76, 0x54, 0x32, 0x10};
    // 01 23 45 67 89 ab cd ef fe dc ba 98 76 54 32 10
    // 00 00 ... 00
    unsigned char in[16 * 4] = {0x01, 0x23, 0x45, 0x67, 0x89, 0xab, 0xcd, 0xef,
                                0xfe, 0xdc, 0xba, 0x98, 0x76, 0x54, 0x32, 0x10};
    SM4_Key sm4_key;
    SM4_KeyInit(key, &sm4_key);
    SM4_AESNI_Encrypt_x4(in, in, &sm4_key);
    // 68 1e df 34 d2 06 96 5e 86 b3 e9 4f 53 6e 42 46
    // 26 77 f4 6b 09 c1 22 cc 97 55 33 10 5b d4 a2 2a
    // 26 ...
    printf("C:\n");
    for (int j = 0; j < 4; j++) {
        printf("\t");
        for (int i = 0; i < 16; i++) {
            printf("%02x ", in[i + 16 * j]);
        }
        printf("\n");
    }

    printf("P:\n");
    SM4_AESNI_Decrypt_x4(in, in, &sm4_key);
    // 01 23 45 67 89 ab cd ef fe dc ba 98 76 54 32 10
    // 00 00 ... 00
    for (int j = 0; j < 4; j++) {
        printf("\t");
        for (int i = 0; i < 16; i++) {
            printf("%02x ", in[i + 16 * j]);
        }
        printf("\n");
    }
    system("pause");
    return 0;
}

sm4.c

#pragma once
#include"sm4.h"

static uint32_t FK[4] = { 0xa3b1bac6, 0x56aa3350, 0x677d9197, 0xb27022dc };
static uint32_t CK[32] = {
    0x00070e15, 0x1c232a31, 0x383f464d, 0x545b6269, 0x70777e85, 0x8c939aa1,
    0xa8afb6bd, 0xc4cbd2d9, 0xe0e7eef5, 0xfc030a11, 0x181f262d, 0x343b4249,
    0x50575e65, 0x6c737a81, 0x888f969d, 0xa4abb2b9, 0xc0c7ced5, 0xdce3eaf1,
    0xf8ff060d, 0x141b2229, 0x30373e45, 0x4c535a61, 0x686f767d, 0x848b9299,
    0xa0a7aeb5, 0xbcc3cad1, 0xd8dfe6ed, 0xf4fb0209, 0x10171e25, 0x2c333a41,
    0x484f565d, 0x646b7279 };
static uint8_t SBox[256] = {
    0xD6, 0x90, 0xE9, 0xFE, 0xCC, 0xE1, 0x3D, 0xB7, 0x16, 0xB6, 0x14, 0xC2,
    0x28, 0xFB, 0x2C, 0x05, 0x2B, 0x67, 0x9A, 0x76, 0x2A, 0xBE, 0x04, 0xC3,
    0xAA, 0x44, 0x13, 0x26, 0x49, 0x86, 0x06, 0x99, 0x9C, 0x42, 0x50, 0xF4,
    0x91, 0xEF, 0x98, 0x7A, 0x33, 0x54, 0x0B, 0x43, 0xED, 0xCF, 0xAC, 0x62,
    0xE4, 0xB3, 0x1C, 0xA9, 0xC9, 0x08, 0xE8, 0x95, 0x80, 0xDF, 0x94, 0xFA,
    0x75, 0x8F, 0x3F, 0xA6, 0x47, 0x07, 0xA7, 0xFC, 0xF3, 0x73, 0x17, 0xBA,
    0x83, 0x59, 0x3C, 0x19, 0xE6, 0x85, 0x4F, 0xA8, 0x68, 0x6B, 0x81, 0xB2,
    0x71, 0x64, 0xDA, 0x8B, 0xF8, 0xEB, 0x0F, 0x4B, 0x70, 0x56, 0x9D, 0x35,
    0x1E, 0x24, 0x0E, 0x5E, 0x63, 0x58, 0xD1, 0xA2, 0x25, 0x22, 0x7C, 0x3B,
    0x01, 0x21, 0x78, 0x87, 0xD4, 0x00, 0x46, 0x57, 0x9F, 0xD3, 0x27, 0x52,
    0x4C, 0x36, 0x02, 0xE7, 0xA0, 0xC4, 0xC8, 0x9E, 0xEA, 0xBF, 0x8A, 0xD2,
    0x40, 0xC7, 0x38, 0xB5, 0xA3, 0xF7, 0xF2, 0xCE, 0xF9, 0x61, 0x15, 0xA1,
    0xE0, 0xAE, 0x5D, 0xA4, 0x9B, 0x34, 0x1A, 0x55, 0xAD, 0x93, 0x32, 0x30,
    0xF5, 0x8C, 0xB1, 0xE3, 0x1D, 0xF6, 0xE2, 0x2E, 0x82, 0x66, 0xCA, 0x60,
    0xC0, 0x29, 0x23, 0xAB, 0x0D, 0x53, 0x4E, 0x6F, 0xD5, 0xDB, 0x37, 0x45,
    0xDE, 0xFD, 0x8E, 0x2F, 0x03, 0xFF, 0x6A, 0x72, 0x6D, 0x6C, 0x5B, 0x51,
    0x8D, 0x1B, 0xAF, 0x92, 0xBB, 0xDD, 0xBC, 0x7F, 0x11, 0xD9, 0x5C, 0x41,
    0x1F, 0x10, 0x5A, 0xD8, 0x0A, 0xC1, 0x31, 0x88, 0xA5, 0xCD, 0x7B, 0xBD,
    0x2D, 0x74, 0xD0, 0x12, 0xB8, 0xE5, 0xB4, 0xB0, 0x89, 0x69, 0x97, 0x4A,
    0x0C, 0x96, 0x77, 0x7E, 0x65, 0xB9, 0xF1, 0x09, 0xC5, 0x6E, 0xC6, 0x84,
    0x18, 0xF0, 0x7D, 0xEC, 0x3A, 0xDC, 0x4D, 0x20, 0x79, 0xEE, 0x5F, 0x3E,
    0xD7, 0xCB, 0x39, 0x48 };

#define rotl32(value, shift) ((value << shift) | value >> (32 - shift))
void SM4_KeyInit(uint8_t* key, SM4_Key* sm4_key) {
    uint32_t k[4];
    uint32_t tmp;
    uint8_t* tmp_ptr8 = (uint8_t*)&tmp;
    // 初始化密钥
    for (int i = 0; i < 4; i++) {
        int j = 4 * i;
        k[i] = (key[j + 0] << 24) | (key[j + 1] << 16) | (key[j + 2] << 8) |
            (key[j + 3]);
        k[i] = k[i] ^ FK[i];
    }
    // 32轮变换
    for (int i = 0; i < 32; i++) {
        tmp = k[1] ^ k[2] ^ k[3] ^ CK[i];
        // SBox 盒变换
        for (int j = 0; j < 4; j++) {
            tmp_ptr8[j] = SBox[tmp_ptr8[j]];
        }
        // 线性变换
        sm4_key->rk[i] = k[0] ^ tmp ^ rotl32(tmp, 13) ^ rotl32(tmp, 23);
        // 移位
        k[0] = k[1];
        k[1] = k[2];
        k[2] = k[3];
        k[3] = sm4_key->rk[i];
    }
}

sm4_aesni_x4.c

#pragma once
#include "sm4_aesni_x4.h"

#include <immintrin.h>

static void SM4_AESNI_do(uint8_t* in, uint8_t* out, SM4_Key* sm4_key, int enc);

void SM4_AESNI_Encrypt_x4(uint8_t* plaintext, uint8_t* ciphertext,
                       SM4_Key* sm4_key) {
    SM4_AESNI_do(plaintext, ciphertext, sm4_key, 0);
}

void SM4_AESNI_Decrypt_x4(uint8_t* ciphertext, uint8_t* plaintext,
                       SM4_Key* sm4_key) {
    SM4_AESNI_do(ciphertext, plaintext, sm4_key, 1);
}

#define MM_PACK0_EPI32(a, b, c, d) \
    _mm_unpacklo_epi64(_mm_unpacklo_epi32(a, b), _mm_unpacklo_epi32(c, d))
#define MM_PACK1_EPI32(a, b, c, d) \
    _mm_unpackhi_epi64(_mm_unpacklo_epi32(a, b), _mm_unpacklo_epi32(c, d))
#define MM_PACK2_EPI32(a, b, c, d) \
    _mm_unpacklo_epi64(_mm_unpackhi_epi32(a, b), _mm_unpackhi_epi32(c, d))
#define MM_PACK3_EPI32(a, b, c, d) \
    _mm_unpackhi_epi64(_mm_unpackhi_epi32(a, b), _mm_unpackhi_epi32(c, d))

#define MM_XOR2(a, b) _mm_xor_si128(a, b)
#define MM_XOR3(a, b, c) MM_XOR2(a, MM_XOR2(b, c))
#define MM_XOR4(a, b, c, d) MM_XOR2(a, MM_XOR3(b, c, d))
#define MM_XOR5(a, b, c, d, e) MM_XOR2(a, MM_XOR4(b, c, d, e))
#define MM_XOR6(a, b, c, d, e, f) MM_XOR2(a, MM_XOR5(b, c, d, e, f))
#define MM_ROTL_EPI32(a, n) \
    MM_XOR2(_mm_slli_epi32(a, n), _mm_srli_epi32(a, 32 - n))

static __m128i SM4_SBox(__m128i x);

static void SM4_AESNI_do(uint8_t* in, uint8_t* out, SM4_Key* sm4_key, int enc) {
    __m128i X[4], Tmp[4];
    __m128i vindex;
    // Load Data
    Tmp[0] = _mm_loadu_si128((const __m128i*)in + 0);
    Tmp[1] = _mm_loadu_si128((const __m128i*)in + 1);
    Tmp[2] = _mm_loadu_si128((const __m128i*)in + 2);
    Tmp[3] = _mm_loadu_si128((const __m128i*)in + 3);
    vindex =
        _mm_setr_epi8(3, 2, 1, 0, 7, 6, 5, 4, 11, 10, 9, 8, 15, 14, 13, 12);
    // Pack Data
    X[0] = MM_PACK0_EPI32(Tmp[0], Tmp[1], Tmp[2], Tmp[3]);
    X[1] = MM_PACK1_EPI32(Tmp[0], Tmp[1], Tmp[2], Tmp[3]);
    X[2] = MM_PACK2_EPI32(Tmp[0], Tmp[1], Tmp[2], Tmp[3]);
    X[3] = MM_PACK3_EPI32(Tmp[0], Tmp[1], Tmp[2], Tmp[3]);
    // Shuffle Endian
    X[0] = _mm_shuffle_epi8(X[0], vindex);
    X[1] = _mm_shuffle_epi8(X[1], vindex);
    X[2] = _mm_shuffle_epi8(X[2], vindex);
    X[3] = _mm_shuffle_epi8(X[3], vindex);
    // Loop
    for (int i = 0; i < 32; i++) {
        __m128i k =
            _mm_set1_epi32((enc == 0) ? sm4_key->rk[i] : sm4_key->rk[31 - i]);
        Tmp[0] = MM_XOR4(X[1], X[2], X[3], k);
        // SBox
        Tmp[0] = SM4_SBox(Tmp[0]);
        // L
        Tmp[0] = MM_XOR6(X[0], Tmp[0], MM_ROTL_EPI32(Tmp[0], 2),
                         MM_ROTL_EPI32(Tmp[0], 10), MM_ROTL_EPI32(Tmp[0], 18),
                         MM_ROTL_EPI32(Tmp[0], 24));
        //
        X[0] = X[1];
        X[1] = X[2];
        X[2] = X[3];
        X[3] = Tmp[0];
    }
    // Shuffle Endian
    X[0] = _mm_shuffle_epi8(X[0], vindex);
    X[1] = _mm_shuffle_epi8(X[1], vindex);
    X[2] = _mm_shuffle_epi8(X[2], vindex);
    X[3] = _mm_shuffle_epi8(X[3], vindex);
    // Pack and Store
    _mm_storeu_si128((__m128i*)out + 0, MM_PACK0_EPI32(X[3], X[2], X[1], X[0]));
    _mm_storeu_si128((__m128i*)out + 1, MM_PACK1_EPI32(X[3], X[2], X[1], X[0]));
    _mm_storeu_si128((__m128i*)out + 2, MM_PACK2_EPI32(X[3], X[2], X[1], X[0]));
    _mm_storeu_si128((__m128i*)out + 3, MM_PACK3_EPI32(X[3], X[2], X[1], X[0]));
}

static __m128i MulMatrix(__m128i x, __m128i higherMask, __m128i lowerMask) {
    __m128i tmp1, tmp2;
    __m128i andMask = _mm_set1_epi32(0x0f0f0f0f);
    tmp2 = _mm_srli_epi16(x, 4);
    tmp1 = _mm_and_si128(x, andMask);
    tmp2 = _mm_and_si128(tmp2, andMask);
    tmp1 = _mm_shuffle_epi8(lowerMask, tmp1);
    tmp2 = _mm_shuffle_epi8(higherMask, tmp2);
    tmp1 = _mm_xor_si128(tmp1, tmp2);
    return tmp1;
}

static __m128i MulMatrixATA(__m128i x) {
    __m128i higherMask =
        _mm_set_epi8(0x14, 0x07, 0xc6, 0xd5, 0x6c, 0x7f, 0xbe, 0xad, 0xb9, 0xaa,
                     0x6b, 0x78, 0xc1, 0xd2, 0x13, 0x00);
    __m128i lowerMask =
        _mm_set_epi8(0xd8, 0xb8, 0xfa, 0x9a, 0xc5, 0xa5, 0xe7, 0x87, 0x5f, 0x3f,
                     0x7d, 0x1d, 0x42, 0x22, 0x60, 0x00);
    return MulMatrix(x, higherMask, lowerMask);
}

static __m128i MulMatrixTA(__m128i x) {
    __m128i higherMask =
        _mm_set_epi8(0x22, 0x58, 0x1a, 0x60, 0x02, 0x78, 0x3a, 0x40, 0x62, 0x18,
                     0x5a, 0x20, 0x42, 0x38, 0x7a, 0x00);
    __m128i lowerMask =
        _mm_set_epi8(0xe2, 0x28, 0x95, 0x5f, 0x69, 0xa3, 0x1e, 0xd4, 0x36, 0xfc,
                     0x41, 0x8b, 0xbd, 0x77, 0xca, 0x00);
    return MulMatrix(x, higherMask, lowerMask);
}

static __m128i AddTC(__m128i x) {
    __m128i TC = _mm_set1_epi8(0b00100011);
    return _mm_xor_si128(x, TC);
}

static __m128i AddATAC(__m128i x) {
    __m128i ATAC = _mm_set1_epi8(0b00111011);
    return _mm_xor_si128(x, ATAC);
}

static __m128i SM4_SBox(__m128i x) {
    __m128i MASK = _mm_set_epi8(0x03, 0x06, 0x09, 0x0c, 0x0f, 0x02, 0x05, 0x08,
                                0x0b, 0x0e, 0x01, 0x04, 0x07, 0x0a, 0x0d, 0x00);
    x = _mm_shuffle_epi8(x, MASK);  //逆行移位
    x = AddTC(MulMatrixTA(x));
    x = _mm_aesenclast_si128(x, _mm_setzero_si128());
    return AddATAC(MulMatrixATA(x));
}