reference python code according https://ppfocus.com/0/ed0cb869b.html, http://zongyue.top:8090/archives/aes%E5%92%8Csm4s%E7%9B%92%E5%A4%8D%E5%90%88%E5%9F%9F%E5%AE%9E%E7%8E%B0%E6%96%B9%E6%B3%95

aes_sm4_sbox.py

from pyfinite import ffield
from pyfinite import genericmatrix

XOR = lambda x,y:x^y
AND = lambda x,y:x&y
DIV = lambda x,y:x

def aes_f():
    gen = 0b100011011
    return ffield.FField(8, gen, useLUT=0)

def sm4_f():
    gen = 0b111110101
    return ffield.FField(8, gen, useLUT=0)    

aesf = aes_f()
sm4f = sm4_f()

def field_pow2(x, F):
    return F.Multiply(x, x)

def field_pow3(x, F):
    return F.Multiply(x, field_pow2(x, F))

def field_pow4(x, F):
    return field_pow2(field_pow2(x, F), F)

def field_pow16(x, F):
    return field_pow4(field_pow4(x, F), F)    

def get_all_WZY(F):
    result_list = []
    for i in range(256):
        if field_pow2(i, F)^i^1 == 0:
            W=i
            W_2 = field_pow2(W, F)
            N = W_2
            for j in range(256):
                if field_pow2(j, F)^j^W_2 == 0:
                    Z = j
                    Z_4 = field_pow4(Z, F)
                    u = F.Multiply(field_pow2(N, F), Z)
                    for k in range(256):
                        if field_pow2(k, F)^k^u == 0:
                            Y = k
                            Y_16 = field_pow16(k, F)
                            result_list.append([W, W_2, Z, Z_4, Y, Y_16])
    return result_list

def gen_X(F, W, W_2, Z, Z_4, Y, Y_16):
    W_2_Z_4_Y_16 = F.Multiply(F.Multiply(W_2, Z_4), Y_16)
    W_Z_4_Y_16 = F.Multiply(F.Multiply(W, Z_4), Y_16)
    W_2_Z_Y_16 = F.Multiply(F.Multiply(W_2, Z), Y_16)
    W_Z_Y_16 = F.Multiply(F.Multiply(W, Z), Y_16)
    W_2_Z_4_Y = F.Multiply(F.Multiply(W_2, Z_4), Y)
    W_Z_4_Y = F.Multiply(F.Multiply(W, Z_4), Y)
    W_2_Z_Y = F.Multiply(F.Multiply(W_2, Z), Y)
    W_Z_Y = F.Multiply(F.Multiply(W, Z), Y)
    return [W_2_Z_4_Y_16, W_Z_4_Y_16, W_2_Z_Y_16, W_Z_Y_16, W_2_Z_4_Y, W_Z_4_Y, W_2_Z_Y, W_Z_Y]

def to_matrix(x):
    m = genericmatrix.GenericMatrix(size=(8,8), zeroElement=0, identityElement=1, add=XOR, mul=AND, sub=XOR, div=DIV)
    m.SetRow(0, [(x[0] & 0x80) >> 7, (x[1] & 0x80) >> 7, (x[2] & 0x80) >> 7, (x[3] & 0x80) >> 7, (x[4] & 0x80) >> 7, (x[5] & 0x80) >> 7, (x[6] & 0x80) >> 7, (x[7] & 0x80) >> 7]) 
    m.SetRow(1, [(x[0] & 0x40) >> 6, (x[1] & 0x40) >> 6, (x[2] & 0x40) >> 6, (x[3] & 0x40) >> 6, (x[4] & 0x40) >> 6, (x[5] & 0x40) >> 6, (x[6] & 0x40) >> 6, (x[7] & 0x40) >> 6]) 
    m.SetRow(2, [(x[0] & 0x20) >> 5, (x[1] & 0x20) >> 5, (x[2] & 0x20) >> 5, (x[3] & 0x20) >> 5, (x[4] & 0x20) >> 5, (x[5] & 0x20) >> 5, (x[6] & 0x20) >> 5, (x[7] & 0x20) >> 5]) 
    m.SetRow(3, [(x[0] & 0x10) >> 4, (x[1] & 0x10) >> 4, (x[2] & 0x10) >> 4, (x[3] & 0x10) >> 4, (x[4] & 0x10) >> 4, (x[5] & 0x10) >> 4, (x[6] & 0x10) >> 4, (x[7] & 0x10) >> 4]) 
    m.SetRow(4, [(x[0] & 0x08) >> 3, (x[1] & 0x08) >> 3, (x[2] & 0x08) >> 3, (x[3] & 0x08) >> 3, (x[4] & 0x08) >> 3, (x[5] & 0x08) >> 3, (x[6] & 0x08) >> 3, (x[7] & 0x08) >> 3]) 
    m.SetRow(5, [(x[0] & 0x04) >> 2, (x[1] & 0x04) >> 2, (x[2] & 0x04) >> 2, (x[3] & 0x04) >> 2, (x[4] & 0x04) >> 2, (x[5] & 0x04) >> 2, (x[6] & 0x04) >> 2, (x[7] & 0x04) >> 2]) 
    m.SetRow(6, [(x[0] & 0x02) >> 1, (x[1] & 0x02) >> 1, (x[2] & 0x02) >> 1, (x[3] & 0x02) >> 1, (x[4] & 0x02) >> 1, (x[5] & 0x02) >> 1, (x[6] & 0x02) >> 1, (x[7] & 0x02) >> 1]) 
    m.SetRow(7, [(x[0] & 0x01) >> 0, (x[1] & 0x01) >> 0, (x[2] & 0x01) >> 0, (x[3] & 0x01) >> 0, (x[4] & 0x01) >> 0, (x[5] & 0x01) >> 0, (x[6] & 0x01) >> 0, (x[7] & 0x01) >> 0]) 
    return m

def matrix_col_byte(c):
    return (c[0] << 7) ^ (c[1] << 6) ^ (c[2] << 5) ^ (c[3] << 4) ^ (c[4] << 3) ^ (c[5] << 2) ^ (c[6] << 1) ^ (c[7] << 0)

def matrix_row_byte(c):
    return (c[0] << 7) ^ (c[1] << 6) ^ (c[2] << 5) ^ (c[3] << 4) ^ (c[4] << 3) ^ (c[5] << 2) ^ (c[6] << 1) ^ (c[7] << 0)

def matrix_cols(m):
    x = []
    for i in range(8):
        c = m.GetColumn(i)
        x.append(matrix_col_byte(c))
    return x

def matrix_rows(m):
    x = []
    for i in range(8):
        r = m.GetRow(i)
        x.append(matrix_row_byte(r))
    return x

def gen_X_inv(x):
    m = to_matrix(x)
    m_inv = m.Inverse()
    return matrix_cols(m_inv)

def G4_mul(x, y):
    '''
    GF(2^2) multiply operator, normal basis is {W^2, W}
    '''
    a = (x & 0x02) >> 1
    b = x & 0x01
    c = (y & 0x02) >> 1
    d = y & 0x01
    e = (a ^ b) & (c ^ d)
    return (((a & c) ^ e) << 1) | ((b & d) ^ e)

def G4_mul_N(x):
    '''
    GF(2^2) multiply N, normal basis is {W^2, W}, N = W^2
    '''
    a = (x & 0x02) >> 1
    b = x & 0x01
    p = b
    q = a ^ b
    return (p << 1) | q

def G4_mul_N2(x):
    '''
    GF(2^2) multiply N^2, normal basis is {W^2, W}, N = W^2
    '''
    a = (x & 0x02) >> 1
    b = x & 0x01
    return ((a ^ b) << 1) | a

def G4_inv(x):
    '''
    GF(2^2) inverse opertor
    '''        
    a = (x & 0x02) >> 1
    b = x & 0x01
    return (b << 1) | a

def G16_mul(x, y):
    '''
    GF(2^4) multiply operator, normal basis is {Z^4, Z}
    '''
    a = (x & 0xc) >> 2
    b = x & 0x03
    c = (y & 0xc) >> 2
    d = y & 0x03
    e = G4_mul(a ^ b, c ^ d)
    e = G4_mul_N(e)
    p = G4_mul(a, c) ^ e
    q = G4_mul(b, d) ^ e
    return (p << 2) | q

def G16_sq_mul_u(x):
    '''
    GF(2^4) x^2 * u operator, u = N^2 Z, N = W^2
    '''    
    a = (x & 0xc) >> 2
    b = x & 0x03
    p = G4_inv(a ^ b)
    q = G4_mul_N2(G4_inv(b))
    return (p << 2) | q

def G16_inv(x):
    '''
    GF(2^4) inverse opertor
    '''
    a = (x & 0xc) >> 2
    b = x & 0x03
    c = G4_mul_N(G4_inv(a ^ b))
    d = G4_mul(a, b)
    e = G4_inv(c ^ d)
    p = G4_mul(e, b)
    q = G4_mul(e, a)
    return (p << 2) | q

def G256_inv(x):
    '''
    GF(2^8) inverse opertor
    '''
    a = (x & 0xf0) >> 4
    b = x & 0x0f
    c = G16_sq_mul_u(a ^ b)
    d = G16_mul(a, b)
    e = G16_inv(c ^ d)
    p = G16_mul(e, b)
    q = G16_mul(e, a)
    return (p << 4) | q

def G256_new_basis(x, b):
    '''
    x presentation under new basis b
    '''
    y = 0
    for i in range(8):
        if x & (1<<((7-i))):
            y ^= b[i]
    return y

AES_A = [0b10001111, 0b11000111, 0b11100011, 0b11110001, 0b11111000, 0b01111100, 0b00111110, 0b00011111]
AES_C = [0, 1, 1, 0, 0, 0, 1, 1]

def AES_SBOX(X, X_inv):
    sbox = []
    for i in range(256):
        t = G256_new_basis(i, X_inv)
        t = G256_inv(t)
        t = G256_new_basis(t, X)
        t = G256_new_basis(t, AES_A)
        sbox.append(t ^ 0x63)
    return sbox

def print_sbox(sbox):
    for i, s in enumerate(sbox):
        print(f'%02x'%s,',', end='')    
        if (i+1) % 16 == 0:
            print()

def print_all_aes_sbox():
    result_list = get_all_WZY(aesf)
    for i, v in enumerate(result_list):
        X = gen_X(aesf, v[0], v[1], v[2], v[3], v[4], v[5])
        X_inv = gen_X_inv(X)
        print_sbox(AES_SBOX(X, X_inv))
        print()

SM4_A = [0b11100101, 0b11110010, 0b01111001, 0b10111100, 0b01011110, 0b00101111, 0b10010111, 0b11001011]
SM4_C = [1, 1, 0, 1, 0, 0, 1, 1]

def SM4_SBOX(X, X_inv):
    sbox = []
    for i in range(256):
        t = G256_new_basis(i, SM4_A)
        t ^= 0xd3
        t = G256_new_basis(t, X_inv)
        t = G256_inv(t)
        t = G256_new_basis(t, X)
        t = G256_new_basis(t, SM4_A)
        sbox.append(t ^ 0xd3)
    return sbox

def print_all_sm4_sbox():
    result_list = get_all_WZY(sm4f)
    for i, v in enumerate(result_list):
        X = gen_X(sm4f, v[0], v[1], v[2], v[3], v[4], v[5])
        X_inv = gen_X_inv(X)
        print_sbox(SM4_SBOX(X, X_inv))
        print()    

def print_m(m):
    for i, s in enumerate(m):
        print(f'0x%02x'%s,',', end='')  

def gen_all_m1_c1_m2_c2():
    aes_result_list = get_all_WZY(aesf)
    sm4_result_list = get_all_WZY(sm4f)
    Aaes = to_matrix(AES_A)
    Aaes_inv = Aaes.Inverse()
    Asm4 = to_matrix(SM4_A)
    Caes = genericmatrix.GenericMatrix(size=(8, 1), zeroElement=0, identityElement=1, add=XOR, mul=AND, sub=XOR, div=DIV)
    for i in range(8):
        Caes.SetRow(i, [AES_C[i]])
    Csm4 = genericmatrix.GenericMatrix(size=(8, 1), zeroElement=0, identityElement=1, add=XOR, mul=AND, sub=XOR, div=DIV)
    for i in range(8):
        Csm4.SetRow(i, [SM4_C[i]])
    for i, v1 in enumerate(aes_result_list):
        Xaes = to_matrix(gen_X(aesf, v1[0], v1[1], v1[2], v1[3], v1[4], v1[5]))
        Xaes_inv = Xaes.Inverse()
        for j, v2 in enumerate(sm4_result_list):
            Xsm4 = to_matrix(gen_X(sm4f, v2[0], v2[1], v2[2], v2[3], v2[4], v2[5]))
            Xsm4_inv = Xsm4.Inverse()
            M1 = Xaes * Xsm4_inv * Asm4
            C1 = Xaes * (Xsm4_inv * Csm4)
            M2 = Asm4 * Xsm4 * Xaes_inv * Aaes_inv
            C2 = M2 * Caes
            print(f'M1=','', end='')
            print_m(matrix_rows(M1))
            print(f' C1=','', end='')
            print(hex(matrix_col_byte(C1.GetColumn(0))))
            print(f'M2=','', end='')
            print_m(matrix_rows(M2))
            print(f' C2=','', end='')
            print(hex(0xd3 ^ matrix_col_byte(C2.GetColumn(0))))
            print()

gen_all_m1_c1_m2_c2()

M1= 0x96 ,0x47 ,0xe9 ,0x3d ,0xde ,0x65 ,0xac ,0xa7
C1= 0x69
M2= 0xfa ,0x64 ,0xb4 ,0x0a ,0x41 ,0xdd ,0x01 ,0xc1
C2= 0x61

M1= 0x52 ,0xbc ,0x2d ,0x02 ,0x9e ,0x25 ,0xac ,0x34
C1= 0x65
M2= 0xcb ,0x9a ,0x0a ,0xb4 ,0xc7 ,0xac ,0x87 ,0x4e
C2= 0x2f

M1= 0x5d ,0x50 ,0x22 ,0x1a ,0xb9 ,0x7d ,0x28 ,0x4c
C1= 0x3e
M2= 0xd3 ,0xba ,0x1d ,0x65 ,0x47 ,0x4c ,0x0e ,0x48
C2= 0x6c

M1= 0xe6 ,0xab ,0x99 ,0x5a ,0x86 ,0x42 ,0x28 ,0x24
C1= 0x8e
M2= 0x2d ,0x8b ,0x65 ,0x1d ,0xc8 ,0xfb ,0x81 ,0xce
C2= 0xe9

M1= 0xd1 ,0x37 ,0xae ,0xce ,0x05 ,0x45 ,0xec ,0xdd
C1= 0x86
M2= 0x50 ,0x16 ,0x5b ,0x2a ,0x53 ,0x92 ,0x62 ,0x33
C2= 0x3c

M1= 0xee ,0xb3 ,0x91 ,0x75 ,0xc1 ,0x81 ,0xec ,0x8a
C1= 0xd6
M2= 0x19 ,0x56 ,0x2a ,0x5b ,0xa4 ,0xea ,0x95 ,0x0b
C2= 0x4d

M1= 0x4d ,0x1f ,0x32 ,0xfe ,0x8e ,0xb1 ,0x17 ,0xd5
C1= 0xce
M2= 0xe8 ,0x28 ,0x74 ,0xc3 ,0xfc ,0x32 ,0x02 ,0x6b
C2= 0x81

M1= 0x0d ,0x9b ,0x72 ,0x3a ,0x35 ,0x0a ,0x17 ,0x06
C1= 0x23
M2= 0xa8 ,0x61 ,0xc3 ,0x74 ,0xc4 ,0x8c ,0x3a ,0x9c
C2= 0x3b