emmansun/aes_sm4_sbox_poly.py

## aes_sm4_sbox_poly.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, type):
    result_list = []
    for i in range(256):
        if field_pow2(i, F)^i^1 == 0:
            W=i
            N = W
            N_2 = field_pow2(N, F)
            for j in range(256):
                if field_pow2(j, F)^j^N == 0:
                    Z = j
                    if type == 0:
                        u = F.Multiply(N, Z)
                    elif type == 1:
                        u = F.Multiply(N_2, Z)
                    elif type == 2:
                        u = F.Multiply(N, Z) ^ N
                    elif type == 3:
                        u = F.Multiply(N_2, Z) ^ N_2
                    elif type == 4:
                        u = F.Multiply(N, Z) ^ 1
                    elif type == 5:
                        u = F.Multiply(N_2, Z) ^ N
                    elif type == 6:
                        u = F.Multiply(N, Z) ^ N_2
                    elif type == 7:
                        u = F.Multiply(N_2, Z) ^ 1
                    else:
                        return result_list

                    for k in range(256):
                        tmp = field_pow2(k, F)^k^u
                        if tmp == 0:
                            Y = k
                            result_list.append([1, W, 1, Z, 1, Y])
                            print_m([1, W, 1, Z, 1, Y])
                            print()
    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, 1}
    '''
    a = (x & 0x02) >> 1
    b = x & 0x01
    c = (y & 0x02) >> 1
    d = y & 0x01
    e = (a ^ b) & (c ^ d)
    return (((b & d) ^ e) << 1) | ((a & c) ^ (b & d))

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

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

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

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

def G16_sq_mul_u(type, x):
    '''
    GF(2^4) x^2 * u operator, N = W
    '''
    a = (x & 0xc) >> 2
    b = x & 0x03
    if type == 0:
        p = G4_inv(a) ^ G4_mul_N(G4_inv(b))
        q = G4_mul_N2(G4_inv(a))
    elif type == 1:
        p = G4_mul_N(G4_inv(a)) ^ G4_mul_N2(G4_inv(b))
        q = G4_inv(a)
    elif type == 2:
        p = G4_mul_N2(G4_inv(a)) ^ G4_mul_N(G4_inv(b))
        q = G4_mul_N(G4_inv(b))
    elif type == 3:
        p = G4_inv(a) ^ G4_mul_N2(G4_inv(b))
        q = G4_mul_N2(G4_inv(b))
    elif type == 4:
        p = G4_mul_N(G4_inv(b))
        q = G4_inv(a ^ b)
    elif type == 5:
        p = G4_mul_N2(G4_inv(b))
        q = G4_mul_N(G4_inv(a ^ b))
    elif type == 6:
        p = G4_mul_N(G4_inv(a ^ b))
        q = G4_mul_N(G4_inv(a ^ b)) ^ G4_inv(b)
    elif type == 7:
        p = G4_mul_N2(G4_inv(a ^ b))
        q = G4_mul_N2(G4_inv(a ^ b)) ^ G4_mul_N(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 ))
    d = G4_mul(a, b)
    e = G4_inv(b)
    e = G4_inv(c ^ d ^ e)
    p = G4_mul(e, a)
    q = G4_mul(e, a ^ b)
    return (p << 2) | q

def G256_inv(type, x):
    '''
    GF(2^8) inverse opertor
    '''
    a = (x & 0xf0) >> 4
    b = x & 0x0f

    c = G16_sq_mul_u(type, a)
    d = G16_mul(a, b)
    e = G16_mul(b, b)
    e = G16_inv(c ^ d ^ e)
    p = G16_mul(e, a)
    q = G16_mul(e, a ^ b)
    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(type, X, X_inv):
    sbox = []
    for i in range(256):
        t = G256_new_basis(i, X_inv)
        t = G256_inv(type, t)
        t = G256_new_basis(t, X)
        t = G256_new_basis(t, AES_A)
        sbox.append(t ^ 0x63)
    return sbox

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

def SM4_SBOX(type, 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(type, t)
        t = G256_new_basis(t, X)
        t = G256_new_basis(t, SM4_A)
        sbox.append(t ^ 0xd3)
    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():
    for type in range(8):
        result_list = get_all_WZY(aesf, type)
        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(type, X, X_inv))
            print()

def print_all_sm4_sbox():
    for type in range(8):
        result_list = get_all_WZY(sm4f, type)
        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(type, 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():
    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 type in range(8):
        aes_result_list = get_all_WZY(aesf, type)
        sm4_result_list = get_all_WZY(sm4f, type)
        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()
	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, type):
	result_list = []
	for i in range(256):
	if field_pow2(i, F)^i^1 == 0:
	W=i
	N = W
	N_2 = field_pow2(N, F)
	for j in range(256):
	if field_pow2(j, F)^j^N == 0:
	Z = j
	if type == 0:
	u = F.Multiply(N, Z)
	elif type == 1:
	u = F.Multiply(N_2, Z)
	elif type == 2:
	u = F.Multiply(N, Z) ^ N
	elif type == 3:
	u = F.Multiply(N_2, Z) ^ N_2
	elif type == 4:
	u = F.Multiply(N, Z) ^ 1
	elif type == 5:
	u = F.Multiply(N_2, Z) ^ N
	elif type == 6:
	u = F.Multiply(N, Z) ^ N_2
	elif type == 7:
	u = F.Multiply(N_2, Z) ^ 1
	else:
	return result_list

	for k in range(256):
	tmp = field_pow2(k, F)^k^u
	if tmp == 0:
	Y = k
	result_list.append([1, W, 1, Z, 1, Y])
	print_m([1, W, 1, Z, 1, Y])
	print()
	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, 1}
	'''
	a = (x & 0x02) >> 1
	b = x & 0x01
	c = (y & 0x02) >> 1
	d = y & 0x01
	e = (a ^ b) & (c ^ d)
	return (((b & d) ^ e) << 1) \| ((a & c) ^ (b & d))

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

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

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

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

	def G16_sq_mul_u(type, x):
	'''
	GF(2^4) x^2 * u operator, N = W
	'''
	a = (x & 0xc) >> 2
	b = x & 0x03
	if type == 0:
	p = G4_inv(a) ^ G4_mul_N(G4_inv(b))
	q = G4_mul_N2(G4_inv(a))
	elif type == 1:
	p = G4_mul_N(G4_inv(a)) ^ G4_mul_N2(G4_inv(b))
	q = G4_inv(a)
	elif type == 2:
	p = G4_mul_N2(G4_inv(a)) ^ G4_mul_N(G4_inv(b))
	q = G4_mul_N(G4_inv(b))
	elif type == 3:
	p = G4_inv(a) ^ G4_mul_N2(G4_inv(b))
	q = G4_mul_N2(G4_inv(b))
	elif type == 4:
	p = G4_mul_N(G4_inv(b))
	q = G4_inv(a ^ b)
	elif type == 5:
	p = G4_mul_N2(G4_inv(b))
	q = G4_mul_N(G4_inv(a ^ b))
	elif type == 6:
	p = G4_mul_N(G4_inv(a ^ b))
	q = G4_mul_N(G4_inv(a ^ b)) ^ G4_inv(b)
	elif type == 7:
	p = G4_mul_N2(G4_inv(a ^ b))
	q = G4_mul_N2(G4_inv(a ^ b)) ^ G4_mul_N(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 ))
	d = G4_mul(a, b)
	e = G4_inv(b)
	e = G4_inv(c ^ d ^ e)
	p = G4_mul(e, a)
	q = G4_mul(e, a ^ b)
	return (p << 2) \| q

	def G256_inv(type, x):
	'''
	GF(2^8) inverse opertor
	'''
	a = (x & 0xf0) >> 4
	b = x & 0x0f

	c = G16_sq_mul_u(type, a)
	d = G16_mul(a, b)
	e = G16_mul(b, b)
	e = G16_inv(c ^ d ^ e)
	p = G16_mul(e, a)
	q = G16_mul(e, a ^ b)
	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(type, X, X_inv):
	sbox = []
	for i in range(256):
	t = G256_new_basis(i, X_inv)
	t = G256_inv(type, t)
	t = G256_new_basis(t, X)
	t = G256_new_basis(t, AES_A)
	sbox.append(t ^ 0x63)
	return sbox

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

	def SM4_SBOX(type, 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(type, t)
	t = G256_new_basis(t, X)
	t = G256_new_basis(t, SM4_A)
	sbox.append(t ^ 0xd3)
	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():
	for type in range(8):
	result_list = get_all_WZY(aesf, type)
	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(type, X, X_inv))
	print()

	def print_all_sm4_sbox():
	for type in range(8):
	result_list = get_all_WZY(sm4f, type)
	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(type, 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():
	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 type in range(8):
	aes_result_list = get_all_WZY(aesf, type)
	sm4_result_list = get_all_WZY(sm4f, type)
	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()