Chapter 23: AES Design and Implementation (SageMath)#
Python Version Available
This notebook contains the SageMath implementation of AES. A pure Python/NumPy version is available in Chapter 23: AES Design and Implementation.
SageMath’s native GF(2^8) replaces all manual polynomial arithmetic, giving direct access to
field inversion, matrix algebra over finite fields, and polynomial ring operations.
23.1 Historical Context: The NIST Competition#
By the mid-1990s the Data Encryption Standard (DES), adopted in 1977, was showing its age. Its 56-bit key was vulnerable to brute-force search (the EFF Deep Crack machine broke DES in 22 hours in 1998), and triple-DES, while secure, was slow.
In January 1997 the U.S. National Institute of Standards and Technology (NIST) issued a public call for proposals for a new block cipher standard to replace DES. The requirements were:
128-bit block size
Support for 128-, 192-, and 256-bit keys
Security, efficiency, and simplicity
Fifteen candidates were submitted. After two rounds of public evaluation the field was narrowed to five finalists: MARS, RC6, Rijndael, Serpent, and Twofish.
In October 2000 NIST selected Rijndael, designed by the Belgian cryptographers Joan Daemen and Vincent Rijmen, as the AES. The standard was published as FIPS 197 in November 2001.
Why Rijndael?
Rijndael combined strong security margins with exceptional performance across a wide range of platforms, from 8-bit smart cards to 64-bit workstations. Its mathematical elegance (based on operations in \(\mathrm{GF}(2^8)\)) and the wide trail strategy for provable resistance to differential and linear cryptanalysis were key factors.
The Wide Trail Strategy#
Daemen and Rijmen’s design philosophy, the wide trail strategy, ensures that after a small number of rounds any non-trivial differential or linear trail must activate many S-boxes. This is achieved through careful choice of the linear layer (MixColumns combined with ShiftRows) so that the branch number of the linear map is maximal.
23.2 Overview of AES#
AES operates on a \(4 \times 4\) array of bytes called the state. For AES-128 (the variant we implement here) there are 10 rounds. Each round (except the last) applies four transformations in sequence:
Step |
Transformation |
Purpose |
|---|---|---|
1 |
SubBytes |
Non-linear substitution via S-box |
2 |
ShiftRows |
Cyclic row shifts for inter-column diffusion |
3 |
MixColumns |
Matrix multiplication over \(\mathrm{GF}(2^8)\) for intra-column diffusion |
4 |
AddRoundKey |
XOR with the round key |
The final (10th) round omits MixColumns. Before the first round an initial AddRoundKey is applied.
State layout
The 16 plaintext bytes \(b_0, b_1, \ldots, b_{15}\) are arranged column-major:
23.3 Arithmetic in \(\mathrm{GF}(2^8)\)#
All AES operations work over the finite field \(\mathrm{GF}(2^8) \cong \mathbb{F}_2[x]/(x^8+x^4+x^3+x+1)\).
In the Python version we implemented GF(2^8) multiplication and inversion manually using
bit-level operations and lookup tables. SageMath provides all of this natively through
GF(2^8) – we simply define the field with the AES irreducible polynomial and get exact
arithmetic for free.
SageMath advantage
With GF(2^8, 'a', modulus=x^8+x^4+x^3+x+1), SageMath handles all polynomial reduction,
multiplication, and inversion internally. There is no need for lookup tables or manual
bit manipulation.
# --- GF(2^8) setup using SageMath's native finite field ---
R.<x> = PolynomialRing(GF(2))
F256.<a> = GF(2**8, name='a', modulus=x^8+x^4+x^3+x+1)
print(f"Field: {F256}")
print(f"Generator: {a}")
print(f"Modulus: {F256.modulus()}")
print(f"Field order: {F256.order()}")
Field: Finite Field in a of size 2^8
Generator: a
Modulus: x^8 + x^4 + x^3 + x + 1
Field order: 256
Hexadecimal representation of elements in \(\mathbb{F}_{2^8}\)#
We define utility functions to convert between hexadecimal strings and elements of \(\mathbb{F}_{2^8}\). These are essential for matching AES test vectors.
def HexToF256(n):
"""Convert a hex string to an element of F256."""
xl=list(map(int,bin(int(n,16))[2:]))
xl.reverse()
return F256(R(xl))
def F256ToHex(el):
"""Convert an element of F256 to a hex string."""
al=R(el).list()
if len(al)==0:
return '00'
al.reverse()
return hex(int("".join(map(str,al)),2))[2:].upper()
# Short aliases for convenience
def h(x):
return HexToF256(x)
def g(x):
return F256ToHex(x)
def HexToBits(h):
return leftpadding(list(map(int,bin(int(h,16))[2:])),8)
def b(h):
return HexToBits(h)
# Quick verification of GF(2^8) arithmetic
# Compare with Python chapter: gf_mul(0x57, 0x83) = 0xC1
a_val = h('57')
b_val = h('83')
print(f"h('57') * h('83') = {g(a_val * b_val)}")
# Inversion: no need for a^254 trick, SageMath does it directly
c_val = h('53')
print(f"h('53')^(-1) = {g(c_val^(-1))}")
print(f"h('53') * h('53')^(-1) = {g(c_val * c_val^(-1))}")
print(f"\nAll 255 non-zero elements are invertible in F256: "
f"{all((a^k)^(-1) * (a^k) == F256(1) for k in range(255))}")
h('57') * h('83') = C1
h('53')^(-1) = CA
h('53') * h('53')^(-1) = 1
All 255 non-zero elements are invertible in F256: True
23.4 AES Specification and State Setup#
We define the AES parameters for AES-128. The state is a \(4 \times 4\) array of \(\mathbb{F}_{2^8}\) elements – in SageMath, these are genuine finite field elements, not just integers.
#AES specification
Nb=4
AEStype='128' #choose '128','192' or '256'
if AEStype=='128':
Nk=4
Nr=10
if AEStype=='192':
Nk=6
Nr=12
if AEStype=='256':
Nk=8
Nr=14
#state is a 4 by 4 array of F256 elements, e.g.
state=[[F256(0) for i in range(0,4)] for j in range(0,4)]
print(f"AES-{AEStype}: Nk={Nk}, Nr={Nr}, Nb={Nb}")
state
AES-128: Nk=4, Nr=10, Nb=4
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
23.5 Building the AES S-Box from First Principles#
The AES S-box is the only non-linear component in the cipher. It is constructed in two steps:
Multiplicative inversion in \(\mathrm{GF}(2^8)\): \(b \mapsto b^{-1}\) (with \(0 \mapsto 0\)).
Affine transformation over \(\mathrm{GF}(2)\):
where \(c = \texttt{0x63}\) (binary 01100011).
This construction guarantees high algebraic degree, high non-linearity, and low differential uniformity.
SageMath approach
In the Python version we had to build the affine matrix and inverse table manually.
Here, SageMath’s GF(2^8) handles inversion natively via el^(-1), and we use
Sage’s matrix() and vector() for the affine step. The circulant matrix is
constructed from the row \([1,0,0,0,1,1,1,1]\) using cyclic rotations.
Helper Functions#
#cyclic shifts of lists
def CyclicRot(li,n):
l=len(li)
return li[n%l:]+li[:-((l-n)%l)]
#(right)padding of the list
def padding(li,n):
return li+[0]*(n-len(li))
#left padding of the list
def leftpadding(li,n):
return [0]*(n-len(li))+li
#x is meant to be an element of F256, we lift it to R ring representation and turn into coefficients list
def l(x,p):
li=R(x).list()
return padding(li,p)
S-Box and SubBytes#
The affine matrix is the \(8 \times 8\) circulant matrix over \(\mathrm{GF}(2)\) built
from the row vector \([1,0,0,0,1,1,1,1]\). After computing \(b^{-1}\) in \(\mathrm{GF}(2^8)\),
we multiply by this matrix and add the constant 0x63.
#important matrix used in the SBox design
m=matrix([CyclicRot([1,0,0,0,1,1,1,1],-k) for k in range(0,8)])
def SBox(el):
if not(el==0):
inv=el**(-1)
else:
inv=F256(el)
return F256(list(m*vector(l(inv,8))))+h('63')
def SubBytes(state):
out=[]
for row in state:
li=list(map(SBox,row))
out.append(li)
return out
# Verify against known values
print(f"SBox(0x00) = {g(SBox(h('00')))} (expected 63)")
print(f"SBox(0x01) = {g(SBox(h('01')))} (expected 7C)")
print(f"SBox(0x53) = {g(SBox(h('53')))} (expected ED)")
print(f"SBox(0xFF) = {g(SBox(h('FF')))} (expected 16)")
SBox(0x00) = 63 (expected 63)
SBox(0x01) = 7C (expected 7C)
SBox(0x53) = ED (expected ED)
SBox(0xFF) = 16 (expected 16)
Visualising the S-Box#
The S-box is traditionally displayed as a \(16 \times 16\) lookup table indexed by the high and low nibbles of the input byte.
# Build the full S-Box table as a 16x16 matrix
sbox_values = []
for i in range(16):
row = []
for j in range(16):
val = i * 16 + j
hex_str = f'{val:02X}'
sbox_out = SBox(h(hex_str))
row.append(int(g(sbox_out), 16))
sbox_values.append(row)
# Display as a formatted table
M_sbox = matrix(ZZ, sbox_values)
print("AES S-Box (16 x 16 lookup table):")
print(" " + " ".join(f" {j:X}" for j in range(16)))
print(" " + "-" * 64)
for i in range(16):
row_str = " ".join(f"{sbox_values[i][j]:02X}" for j in range(16))
print(f" {i:X} | {row_str}")
# SageMath matrix plot
P = matrix_plot(M_sbox, cmap='viridis', colorbar=True)
show(P, figsize=8, title='AES S-Box Values')
AES S-Box (16 x 16 lookup table):
0 1 2 3 4 5 6 7 8 9 A B C D E F
----------------------------------------------------------------
0 | 63 7C 77 7B F2 6B 6F C5 30 01 67 2B FE D7 AB 76
1 | CA 82 C9 7D FA 59 47 F0 AD D4 A2 AF 9C A4 72 C0
2 | B7 FD 93 26 36 3F F7 CC 34 A5 E5 F1 71 D8 31 15
3 | 04 C7 23 C3 18 96 05 9A 07 12 80 E2 EB 27 B2 75
4 | 09 83 2C 1A 1B 6E 5A A0 52 3B D6 B3 29 E3 2F 84
5 | 53 D1 00 ED 20 FC B1 5B 6A CB BE 39 4A 4C 58 CF
6 | D0 EF AA FB 43 4D 33 85 45 F9 02 7F 50 3C 9F A8
7 | 51 A3 40 8F 92 9D 38 F5 BC B6 DA 21 10 FF F3 D2
8 | CD 0C 13 EC 5F 97 44 17 C4 A7 7E 3D 64 5D 19 73
9 | 60 81 4F DC 22 2A 90 88 46 EE B8 14 DE 5E 0B DB
A | E0 32 3A 0A 49 06 24 5C C2 D3 AC 62 91 95 E4 79
B | E7 C8 37 6D 8D D5 4E A9 6C 56 F4 EA 65 7A AE 08
C | BA 78 25 2E 1C A6 B4 C6 E8 DD 74 1F 4B BD 8B 8A
D | 70 3E B5 66 48 03 F6 0E 61 35 57 B9 86 C1 1D 9E
E | E1 F8 98 11 69 D9 8E 94 9B 1E 87 E9 CE 55 28 DF
F | 8C A1 89 0D BF E6 42 68 41 99 2D 0F B0 54 BB 16
23.6 The Four AES Transformations#
We now implement each transformation individually before assembling the full cipher.
23.6.1 SubBytes#
SubBytes applies the S-box independently to each of the 16 bytes of the state. It is the only non-linear transformation in AES. (Already defined above.)
23.6.2 ShiftRows#
ShiftRows cyclically shifts each row of the state to the left:
Row 0: no shift
Row 1: shift left by 1
Row 2: shift left by 2
Row 3: shift left by 3
This provides inter-column diffusion: after ShiftRows, each column of the state contains bytes from four different original columns.
#shifts initialization designed in AES
if Nb==4:
C0=0
C1=1
C2=2
C3=3
def ShiftRows(state):
row0,row1,row2,row3=state
return [CyclicRot(row0,C0),CyclicRot(row1,C1),CyclicRot(row2,C2),CyclicRot(row3,C3)]
# Demonstration
demo = [[h(x) for x in row] for row in
[['D4','E0','B8','1E'],
['27','BF','B4','41'],
['11','98','5D','52'],
['AE','F1','E5','30']]]
shifted = ShiftRows(demo)
print("Before ShiftRows:")
for row in demo:
print(' '.join(g(x) for x in row))
print("\nAfter ShiftRows:")
for row in shifted:
print(' '.join(g(x) for x in row))
Before ShiftRows:
D4 E0 B8 1E
27 BF B4 41
11 98 5D 52
AE F1 E5 30
After ShiftRows:
D4 E0 B8 1E
BF B4 41 27
5D 52 11 98
30 AE F1 E5
23.6.3 MixColumns#
MixColumns multiplies each column of the state by a fixed polynomial over \(\mathrm{GF}(2^8)\):
modulo \(x^4 + 1\). This corresponds to multiplying each column by the circulant matrix:
The inverse uses the polynomial \(c^{-1}(x) = \texttt{0B} \cdot x^3 + \texttt{0D} \cdot x^2 + \texttt{09} \cdot x + \texttt{0E}\).
SageMath approach
Instead of building multiplication tables and doing manual GF(2^8) matrix-vector
products, we define a polynomial ring PolynomialRing(F256) and perform the
multiplication modulo \(T^4 + 1\) directly. SageMath handles all the field arithmetic.
#MixColumns setup
R2.<T>=PolynomialRing(F256)
cpol=h('03')*T^3+h('01')*T^2+h('01')*T+h('02')
cpolinv=(h('0B')*T^3+h('0D')*T^2+h('09')*T+h('0E'))
print(f"c(T) = {cpol}")
print(f"c_inv(T) = {cpolinv}")
c(T) = (a + 1)*T^3 + T^2 + T + a
c_inv(T) = (a^3 + a + 1)*T^3 + (a^3 + a^2 + 1)*T^2 + (a^3 + 1)*T + a^3 + a^2 + a
#inverse check: c(T) * c_inv(T) mod (T^4+1) should be 1
product = (cpol*cpolinv)%(T^4+1)
print(f"c(T) * c_inv(T) mod (T^4+1) = {product}")
print(f"Is identity: {product == 1}")
c(T) * c_inv(T) mod (T^4+1) = 1
Is identity: True
def MixColumns(state):
"""Apply MixColumns to the state using polynomial multiplication over F256."""
m=matrix(state)
mt=transpose(m)
c0,c1,c2,c3=list(mt)
vpol=vector([T^3,T^2,T,1])
newstate=[]
for c in [c0,c1,c2,c3]:
p=((vector(R2,padding(list(reversed((c).list())),4-len(c.list())))*vpol)*cpol)%(T^4+1)
newstate+=[padding(p.list(),4)]
return list(transpose(matrix(newstate)))
# Demonstration
demo = [[h(x) for x in row] for row in
[['D4','E0','B8','1E'],
['BF','B4','41','27'],
['5D','52','11','98'],
['30','AE','F1','E5']]]
mixed = MixColumns(demo)
print("Before MixColumns:")
for row in demo:
print(' '.join(g(x) for x in row))
print("\nAfter MixColumns:")
for row in mixed:
print(' '.join(g(x) for x in row))
Before MixColumns:
D4 E0 B8 1E
BF B4 41 27
5D 52 11 98
30 AE F1 E5
After MixColumns:
4 E0 48 28
66 CB F8 6
81 19 D3 26
E5 9A 7A 4C
23.6.4 AddRoundKey#
AddRoundKey adds (XOR in \(\mathrm{GF}(2^8)\) is simply + in the field) the state
with the 128-bit round key. Since addition in \(\mathrm{GF}(2^8)\) is its own inverse,
the same function serves for both encryption and decryption.
SageMath note
In \(\mathrm{GF}(2^8)\), addition + is equivalent to bitwise XOR. SageMath’s
matrix(state) + transpose(matrix(roundkey)) performs element-wise field addition
directly.
def AddRoundKey(state,roundkey):
"""Add (XOR) the state with the round key using GF(2^8) addition."""
return list(matrix(state)+transpose(matrix(roundkey)))
# Demonstration with NIST test values
s = [[h(x) for x in row] for row in
[['32','88','31','E0'],
['43','5A','31','37'],
['F6','30','98','07'],
['A8','8D','A2','34']]]
k_words = [[h(x) for x in word] for word in
[['2B','7E','15','16'],
['28','AE','D2','A6'],
['AB','F7','15','88'],
['09','CF','4F','3C']]]
result = AddRoundKey(s, k_words)
print("State + Key:")
for row in result:
print(' '.join(g(x) for x in row))
State + Key:
19 A0 9A E9
3D F4 C6 F8
E3 E2 8D 48
BE 2B 2A 8
23.7 AES Key Expansion#
AES-128 expands the 16-byte (128-bit) cipher key into 11 round keys (one for the initial AddRoundKey plus one for each of the 10 rounds), giving 44 four-byte words total.
The key schedule uses:
RotWord: cyclic left rotation of a 4-byte word
SubWord: apply the S-box to each byte
Rcon: round constants \(\text{Rcon}[i] = (\texttt{rc}_i, 0, 0, 0)\) where \(\texttt{rc}_1 = 1\), \(\texttt{rc}_i = 2 \cdot \texttt{rc}_{i-1}\) in \(\mathrm{GF}(2^8)\)
SageMath note
The round constants \(\texttt{rc}_i = a^{i-1}\) where \(a\) is the generator of \(\mathrm{GF}(2^8)\). SageMath computes these powers directly – no need to precompute an RCON table.
def KeyExpansion(Key):
"""
AES key expansion.
Key: list of 16 elements from F256.
Returns: list of Nb*(Nr+1) = 44 four-byte words.
"""
W=[0]*(Nb*(Nr+1))
if Nk<=6:
for i in range(0,Nk):
W[i]=[F256(Key[4*i+j]) for j in range(0,4)]
for i in range(Nk,Nb*(Nr+1)):
temp=W[i-1]
if i%Nk==0:
temp=list(map(SBox,CyclicRot(temp,1)))
temp[0]=temp[0]+a**(i//Nk-1) # Rcon uses generator a
W[i]=list(vector(F256,W[i-Nk])+vector(F256,temp))
return W
else: # Nk > 6 (AES-256)
for i in range(0, Nk):
W[i] = [F256(Key[4*i+j]) for j in range(0, 4)]
for i in range(Nk, Nb*(Nr+1)):
temp = W[i-1]
if i % Nk == 0:
temp = list(map(SBox, CyclicRot(temp, 1)))
temp[0] = temp[0] + a**(i//Nk - 1) # Rcon
elif i % Nk == 4:
temp = list(map(SBox, temp))
W[i] = list(vector(F256, W[i-Nk]) + vector(F256, temp))
return W
def RoundKey(expkey,i):
"""Extract round key i from the expanded key."""
return expkey[Nb*i:Nb*(i+1)]
def PrintState(state):
"""Display state as hex strings."""
return [[g(x) for x in r] for r in state]
# Test with NIST key
test_key = [h(x) for x in
['2B','7E','15','16','28','AE','D2','A6',
'AB','F7','15','88','09','CF','4F','3C']]
expkey = KeyExpansion(test_key)
print("Round key 0 (original key):")
print(PrintState(RoundKey(expkey, 0)))
print("\nRound key 1:")
print(PrintState(RoundKey(expkey, 1)))
print("\nRound key 10 (last):")
print(PrintState(RoundKey(expkey, 10)))
Round key 0 (original key):
[['2B', '7E', '15', '16'], ['28', 'AE', 'D2', 'A6'], ['AB', 'F7', '15', '88'], ['9', 'CF', '4F', '3C']]
Round key 1:
[['A0', 'FA', 'FE', '17'], ['88', '54', '2C', 'B1'], ['23', 'A3', '39', '39'], ['2A', '6C', '76', '5']]
Round key 10 (last):
[['D0', '14', 'F9', 'A8'], ['C9', 'EE', '25', '89'], ['E1', '3F', 'C', 'C8'], ['B6', '63', 'C', 'A6']]
23.8 The Complete AES Cipher#
We now assemble all components into the full AES encryption and decryption functions.
State representation
The state is a list of 4 lists (rows), each containing 4 elements of \(\mathbb{F}_{2^8}\). The 16-byte input fills the state column-major to match FIPS 197.
def AESEncryptionVerbose(inp,key):
"""AES encryption with verbose output of every intermediate state."""
expkey=KeyExpansion(key)
state=inp
print("Input")
print(PrintState(state))
state=AddRoundKey(state,RoundKey(expkey,0))
print("Round Key Value")
print(PrintState(RoundKey(expkey,0)))
for rnd in range(1,Nr):
print("\n Start Of Round {}".format(rnd))
print(PrintState(state))
state=SubBytes(state)
print("After SubBytes")
print(PrintState(state))
state=ShiftRows(state)
print("After ShiftRows")
print(PrintState(state))
state=MixColumns(state)
print("After MixColumns")
print(PrintState(state))
state=AddRoundKey(state,RoundKey(expkey,rnd))
print("Round Key Value")
print(PrintState(RoundKey(expkey,rnd)))
print("\n Start Of the Final Round {}".format(10))
print(PrintState(state))
state=SubBytes(state)
print("After SubBytes")
print(PrintState(state))
state=ShiftRows(state)
print("After ShiftRows")
print(PrintState(state))
state=AddRoundKey(state,RoundKey(expkey,Nr))
print("\n Output")
print(PrintState(state))
return [list(x) for x in state]
def AESEncryption(inp,key):
"""AES encryption (compact version without debug output)."""
expkey=KeyExpansion(key)
state=inp
state=AddRoundKey(state,RoundKey(expkey,0))
for rnd in range(1,Nr):
state=SubBytes(state)
state=ShiftRows(state)
state=MixColumns(state)
state=AddRoundKey(state,RoundKey(expkey,rnd))
state=SubBytes(state)
state=ShiftRows(state)
state=AddRoundKey(state,RoundKey(expkey,Nr))
return [list(x) for x in state]
NIST FIPS 197 Test Vector#
We verify our implementation against the official NIST test vector from Appendix B of FIPS 197.
# NIST FIPS 197 Appendix B test vector
inp=matrix(4,4,[h(x) for x in ['32','43','f6','a8','88','5a','30','8d','31','31','98','a2','e0','37','07','34']])
inp=list(transpose(inp))
key=[h(x) for x in ['2b','7e','15','16','28','ae','d2','a6','ab','f7','15','88','09','cf','4f','3c']]
print("=== NIST FIPS 197 Test Vector ===")
print()
cipher=AESEncryptionVerbose(inp,key)
# Expected ciphertext: 39 25 84 1D 02 DC 09 FB DC 11 85 97 19 6A 0B 32
print("\n=== Expected ciphertext: 39 25 84 1D 02 DC 09 FB DC 11 85 97 19 6A 0B 32 ===")
=== NIST FIPS 197 Test Vector ===
Input
[['32', '88', '31', 'E0'], ['43', '5A', '31', '37'], ['F6', '30', '98', '7'], ['A8', '8D', 'A2', '34']]
Round Key Value
[['2B', '7E', '15', '16'], ['28', 'AE', 'D2', 'A6'], ['AB', 'F7', '15', '88'], ['9', 'CF', '4F', '3C']]
Start Of Round 1
[['19', 'A0', '9A', 'E9'], ['3D', 'F4', 'C6', 'F8'], ['E3', 'E2', '8D', '48'], ['BE', '2B', '2A', '8']]
After SubBytes
[['D4', 'E0', 'B8', '1E'], ['27', 'BF', 'B4', '41'], ['11', '98', '5D', '52'], ['AE', 'F1', 'E5', '30']]
After ShiftRows
[['D4', 'E0', 'B8', '1E'], ['BF', 'B4', '41', '27'], ['5D', '52', '11', '98'], ['30', 'AE', 'F1', 'E5']]
After MixColumns
[['4', 'E0', '48', '28'], ['66', 'CB', 'F8', '6'], ['81', '19', 'D3', '26'], ['E5', '9A', '7A', '4C']]
Round Key Value
[['A0', 'FA', 'FE', '17'], ['88', '54', '2C', 'B1'], ['23', 'A3', '39', '39'], ['2A', '6C', '76', '5']]
Start Of Round 2
[['A4', '68', '6B', '2'], ['9C', '9F', '5B', '6A'], ['7F', '35', 'EA', '50'], ['F2', '2B', '43', '49']]
After SubBytes
[['49', '45', '7F', '77'], ['DE', 'DB', '39', '2'], ['D2', '96', '87', '53'], ['89', 'F1', '1A', '3B']]
After ShiftRows
[['49', '45', '7F', '77'], ['DB', '39', '2', 'DE'], ['87', '53', 'D2', '96'], ['3B', '89', 'F1', '1A']]
After MixColumns
[['58', '1B', 'DB', '1B'], ['4D', '4B', 'E7', '6B'], ['CA', '5A', 'CA', 'B0'], ['F1', 'AC', 'A8', 'E5']]
Round Key Value
[['F2', 'C2', '95', 'F2'], ['7A', '96', 'B9', '43'], ['59', '35', '80', '7A'], ['73', '59', 'F6', '7F']]
Start Of Round 3
[['AA', '61', '82', '68'], ['8F', 'DD', 'D2', '32'], ['5F', 'E3', '4A', '46'], ['3', 'EF', 'D2', '9A']]
After SubBytes
[['AC', 'EF', '13', '45'], ['73', 'C1', 'B5', '23'], ['CF', '11', 'D6', '5A'], ['7B', 'DF', 'B5', 'B8']]
After ShiftRows
[['AC', 'EF', '13', '45'], ['C1', 'B5', '23', '73'], ['D6', '5A', 'CF', '11'], ['B8', '7B', 'DF', 'B5']]
After MixColumns
[['75', '20', '53', 'BB'], ['EC', 'B', 'C0', '25'], ['9', '63', 'CF', 'D0'], ['93', '33', '7C', 'DC']]
Round Key Value
[['3D', '80', '47', '7D'], ['47', '16', 'FE', '3E'], ['1E', '23', '7E', '44'], ['6D', '7A', '88', '3B']]
Start Of Round 4
[['48', '67', '4D', 'D6'], ['6C', '1D', 'E3', '5F'], ['4E', '9D', 'B1', '58'], ['EE', 'D', '38', 'E7']]
After SubBytes
[['52', '85', 'E3', 'F6'], ['50', 'A4', '11', 'CF'], ['2F', '5E', 'C8', '6A'], ['28', 'D7', '7', '94']]
After ShiftRows
[['52', '85', 'E3', 'F6'], ['A4', '11', 'CF', '50'], ['C8', '6A', '2F', '5E'], ['94', '28', 'D7', '7']]
After MixColumns
[['F', '60', '6F', '5E'], ['D6', '31', 'C0', 'B3'], ['DA', '38', '10', '13'], ['A9', 'BF', '6B', '1']]
Round Key Value
[['EF', '44', 'A5', '41'], ['A8', '52', '5B', '7F'], ['B6', '71', '25', '3B'], ['DB', 'B', 'AD', '00']]
Start Of Round 5
[['E0', 'C8', 'D9', '85'], ['92', '63', 'B1', 'B8'], ['7F', '63', '35', 'BE'], ['E8', 'C0', '50', '1']]
After SubBytes
[['E1', 'E8', '35', '97'], ['4F', 'FB', 'C8', '6C'], ['D2', 'FB', '96', 'AE'], ['9B', 'BA', '53', '7C']]
After ShiftRows
[['E1', 'E8', '35', '97'], ['FB', 'C8', '6C', '4F'], ['96', 'AE', 'D2', 'FB'], ['7C', '9B', 'BA', '53']]
After MixColumns
[['25', 'BD', 'B6', '4C'], ['D1', '11', '3A', '4C'], ['A9', 'D1', '33', 'C0'], ['AD', '68', '8E', 'B0']]
Round Key Value
[['D4', 'D1', 'C6', 'F8'], ['7C', '83', '9D', '87'], ['CA', 'F2', 'B8', 'BC'], ['11', 'F9', '15', 'BC']]
Start Of Round 6
[['F1', 'C1', '7C', '5D'], ['00', '92', 'C8', 'B5'], ['6F', '4C', '8B', 'D5'], ['55', 'EF', '32', 'C']]
After SubBytes
[['A1', '78', '10', '4C'], ['63', '4F', 'E8', 'D5'], ['A8', '29', '3D', '3'], ['FC', 'DF', '23', 'FE']]
After ShiftRows
[['A1', '78', '10', '4C'], ['4F', 'E8', 'D5', '63'], ['3D', '3', 'A8', '29'], ['FE', 'FC', 'DF', '23']]
After MixColumns
[['4B', '2C', '33', '37'], ['86', '4A', '9D', 'D2'], ['8D', '89', 'F4', '18'], ['6D', '80', 'E8', 'D8']]
Round Key Value
[['6D', '88', 'A3', '7A'], ['11', 'B', '3E', 'FD'], ['DB', 'F9', '86', '41'], ['CA', '00', '93', 'FD']]
Start Of Round 7
[['26', '3D', 'E8', 'FD'], ['E', '41', '64', 'D2'], ['2E', 'B7', '72', '8B'], ['17', '7D', 'A9', '25']]
After SubBytes
[['F7', '27', '9B', '54'], ['AB', '83', '43', 'B5'], ['31', 'A9', '40', '3D'], ['F0', 'FF', 'D3', '3F']]
After ShiftRows
[['F7', '27', '9B', '54'], ['83', '43', 'B5', 'AB'], ['40', '3D', '31', 'A9'], ['3F', 'F0', 'FF', 'D3']]
After MixColumns
[['14', '46', '27', '34'], ['15', '16', '46', '2A'], ['B5', '15', '56', 'D8'], ['BF', 'EC', 'D7', '43']]
Round Key Value
[['4E', '54', 'F7', 'E'], ['5F', '5F', 'C9', 'F3'], ['84', 'A6', '4F', 'B2'], ['4E', 'A6', 'DC', '4F']]
Start Of Round 8
[['5A', '19', 'A3', '7A'], ['41', '49', 'E0', '8C'], ['42', 'DC', '19', '4'], ['B1', '1F', '65', 'C']]
After SubBytes
[['BE', 'D4', 'A', 'DA'], ['83', '3B', 'E1', '64'], ['2C', '86', 'D4', 'F2'], ['C8', 'C0', '4D', 'FE']]
After ShiftRows
[['BE', 'D4', 'A', 'DA'], ['3B', 'E1', '64', '83'], ['D4', 'F2', '2C', '86'], ['FE', 'C8', 'C0', '4D']]
After MixColumns
[['00', 'B1', '54', 'FA'], ['51', 'C8', '76', '1B'], ['2F', '89', '6D', '99'], ['D1', 'FF', 'CD', 'EA']]
Round Key Value
[['EA', 'D2', '73', '21'], ['B5', '8D', 'BA', 'D2'], ['31', '2B', 'F5', '60'], ['7F', '8D', '29', '2F']]
Start Of Round 9
[['EA', '4', '65', '85'], ['83', '45', '5D', '96'], ['5C', '33', '98', 'B0'], ['F0', '2D', 'AD', 'C5']]
After SubBytes
[['87', 'F2', '4D', '97'], ['EC', '6E', '4C', '90'], ['4A', 'C3', '46', 'E7'], ['8C', 'D8', '95', 'A6']]
After ShiftRows
[['87', 'F2', '4D', '97'], ['6E', '4C', '90', 'EC'], ['46', 'E7', '4A', 'C3'], ['A6', '8C', 'D8', '95']]
After MixColumns
[['47', '40', 'A3', '4C'], ['37', 'D4', '70', '9F'], ['94', 'E4', '3A', '42'], ['ED', 'A5', 'A6', 'BC']]
Round Key Value
[['AC', '77', '66', 'F3'], ['19', 'FA', 'DC', '21'], ['28', 'D1', '29', '41'], ['57', '5C', '00', '6E']]
Start Of the Final Round 10
[['EB', '59', '8B', '1B'], ['40', '2E', 'A1', 'C3'], ['F2', '38', '13', '42'], ['1E', '84', 'E7', 'D2']]
After SubBytes
[['E9', 'CB', '3D', 'AF'], ['9', '31', '32', '2E'], ['89', '7', '7D', '2C'], ['72', '5F', '94', 'B5']]
After ShiftRows
[['E9', 'CB', '3D', 'AF'], ['31', '32', '2E', '9'], ['7D', '2C', '89', '7'], ['B5', '72', '5F', '94']]
Output
[['39', '2', 'DC', '19'], ['25', 'DC', '11', '6A'], ['84', '9', '85', 'B'], ['1D', 'FB', '97', '32']]
=== Expected ciphertext: 39 25 84 1D 02 DC 09 FB DC 11 85 97 19 6A 0B 32 ===
23.9 AES Decryption#
Each AES transformation has an inverse. Decryption applies them in reverse order.
Inv#
SubBytes
#important matrix used in the SBox design
m=matrix(GF(2),[CyclicRot([1,0,0,0,1,1,1,1],-k) for k in range(0,8)])
minv=m**(-1)
def SBoxinv(el):
"""Compute the inverse S-Box for a single element of F256."""
val=F256((list(minv*vector(l(el-h('63'),8)))))
if not(val==F256(0)):
return val**(-1)
else:
return val
def InvSubBytes(state):
"""Apply inverse S-Box to every element of the state."""
out=[]
for row in state:
li=list(map(SBoxinv,row))
out.append(li)
return out
# Verify round-trip
test_val = h('53')
print(f"SBox(0x53) = {g(SBox(test_val))}")
print(f"SBoxinv(SBox(0x53)) = {g(SBoxinv(SBox(test_val)))}")
print(f"Round-trip check: {SBoxinv(SBox(test_val)) == test_val}")
SBox(0x53) = ED
SBoxinv(SBox(0x53)) = 53
Round-trip check: True
InvMixColumns#
def InvMixColumns(state):
"""Inverse MixColumns using the inverse polynomial cpolinv."""
m=matrix(state)
mt=transpose(m)
c0,c1,c2,c3=list(mt)
vpol=vector([T^3,T^2,T,1])
newstate=[]
for c in [c0,c1,c2,c3]:
p=((vector(R2,padding(list(reversed((c).list())),4-len(c.list())))*vpol)*cpolinv)%(T^4+1)
newstate+=[padding(p.list(),4)]
return list(transpose(matrix(newstate)))
InvShiftRows#
def InvShiftRows(state):
"""Inverse ShiftRows: shift rows to the right."""
row0,row1,row2,row3=[list(x) for x in state]
return [CyclicRot(row0,4-C0),CyclicRot(row1,4-C1),CyclicRot(row2,4-C2),CyclicRot(row3,4-C3)]
Complete AES Decryption#
def AESDecryption(inp,key):
"""AES-128 decryption."""
expkey=KeyExpansion(key)
state=inp
state=AddRoundKey(state,RoundKey(expkey,Nr))
for rnd in range(Nr-1,0,-1):
state=InvShiftRows(state)
state=InvSubBytes(state)
state=AddRoundKey(state,RoundKey(expkey,rnd))
state=InvMixColumns(state)
state=InvShiftRows(state)
state=InvSubBytes(state)
state=AddRoundKey(state,RoundKey(expkey,0))
return [list(x) for x in state]
Encryption/Decryption Round-Trip Verification#
# Verify encrypt-then-decrypt round-trip with NIST test vector
inp=matrix(4,4,[h(x) for x in ['32','43','f6','a8','88','5a','30','8d','31','31','98','a2','e0','37','07','34']])
inp=transpose(inp)
key=[h(x) for x in ['2b','7e','15','16','28','ae','d2','a6','ab','f7','15','88','09','cf','4f','3c']]
cipher=AESEncryption(inp,key)
plain=AESDecryption(matrix(cipher),key)
print("Ciphertext:")
print([[g(x) for x in row] for row in cipher])
print("\nDecrypted plaintext:")
print([[g(x) for x in row] for row in plain])
print("\nMatches original:", [[g(x) for x in row] for row in plain] ==
[[g(x) for x in row] for row in inp])
Ciphertext:
[['39', '2', 'DC', '19'], ['25', 'DC', '11', '6A'], ['84', '9', '85', 'B'], ['1D', 'FB', '97', '32']]
Decrypted plaintext:
[['32', '88', '31', 'E0'], ['43', '5A', '31', '37'], ['F6', '30', '98', '7'], ['A8', '8D', 'A2', '34']]
Matches original: True
23.10 Experiments#
23.10.1 Algebraic Structure of the S-Box in SageMath#
A unique advantage of SageMath is the ability to study the AES S-box algebraically. We can express the S-box as a polynomial over \(\mathrm{GF}(2^8)\) using Lagrange interpolation.
# Build the S-box as a lookup dictionary over F256
sbox_dict = {}
for i in range(256):
hex_in = f'{i:02X}'
sbox_dict[h(hex_in)] = SBox(h(hex_in))
# Lagrange interpolation: express S-box as a polynomial over GF(2^8)
R_interp.<Y> = PolynomialRing(F256)
points = [(el, sbox_dict[el]) for el in sorted(sbox_dict.keys(),
key=lambda e: int(g(e), 16))]
sbox_poly = R_interp.lagrange_polynomial(points)
print(f"Degree of S-box polynomial: {sbox_poly.degree()}")
print(f"Number of non-zero terms: {len(sbox_poly.coefficients())}")
print(f"\nNote: The S-box has algebraic degree 7 (degree 254 = 2^8 - 2),")
print("which provides high resistance to algebraic attacks.")
# Verify that the polynomial matches the S-box at a few points
for test in ['00', '01', '53', 'AB', 'FF']:
poly_val = sbox_poly(h(test))
sbox_val = SBox(h(test))
print(f" poly({test}) = {g(poly_val)}, SBox({test}) = {g(sbox_val)}, match: {poly_val == sbox_val}")
Degree of S-box polynomial: 254
Number of non-zero terms: 9
Note: The S-box has algebraic degree 7 (degree 254 = 2^8 - 2),
which provides high resistance to algebraic attacks.
poly(00) = 63, SBox(00) = 63, match: True
poly(01) = 7C, SBox(01) = 7C, match: True
poly(53) = ED, SBox(53) = ED, match: True
poly(AB) = 62, SBox(AB) = 62, match: True
poly(FF) = 16, SBox(FF) = 16, match: True
23.10.2 S-Box Nonlinearity Analysis#
The nonlinearity of a Boolean function \(f: \mathbb{F}_2^n \to \mathbb{F}_2\) is the minimum Hamming distance between \(f\) and all affine functions. It can be computed via the Walsh-Hadamard transform (WHT).
For the AES S-box (which maps 8 bits to 8 bits), we analyse each of the 8 component Boolean functions and compute their nonlinearity.
The maximum possible nonlinearity for an 8-bit function is \(2^7 - 2^3 = 112\). The AES S-box achieves this maximum, making it optimally resistant to linear cryptanalysis.
def walsh_hadamard_transform(f):
"""
Compute the Walsh-Hadamard transform of a Boolean function.
f: list of length 2^n with values in {0, 1}
Returns: WHT coefficients as a list.
"""
n = len(f)
# Convert {0,1} -> {1,-1}
h_arr = [1 - 2 * int(fi) for fi in f]
# Fast Walsh-Hadamard
step = 1
while step < n:
for i in range(0, n, 2 * step):
for j in range(step):
u = h_arr[i + j]
v = h_arr[i + j + step]
h_arr[i + j] = u + v
h_arr[i + j + step] = u - v
step *= 2
return h_arr
def nonlinearity_bool(f):
"""Compute the nonlinearity of a Boolean function."""
wht = walsh_hadamard_transform(f)
return (len(f) - max(abs(w) for w in wht)) // 2
# Analyse each component function of the AES S-box
n_bits = 8
nonlinearities = []
for bit in range(n_bits):
f = []
for x_val in range(256):
hex_str = f'{x_val:02X}'
sbox_out = int(g(SBox(h(hex_str))), 16)
f.append((sbox_out >> bit) & 1)
nl = nonlinearity_bool(f)
nonlinearities.append(nl)
print("Component function nonlinearities:")
for bit, nl in enumerate(nonlinearities):
print(f" Bit {bit}: NL = {nl} (max possible = 112)")
# SageMath bar chart
P_nl = bar_chart(nonlinearities, color='teal')
P_nl += line([(0, 112), (7, 112)], color='red', linestyle='--',
legend_label='Maximum (112)')
show(P_nl, figsize=6, title='AES S-Box Component Nonlinearities',
axes_labels=['Output bit', 'Nonlinearity'])
Component function nonlinearities:
Bit 0: NL = 112 (max possible = 112)
Bit 1: NL = 112 (max possible = 112)
Bit 2: NL = 112 (max possible = 112)
Bit 3: NL = 112 (max possible = 112)
Bit 4: NL = 112 (max possible = 112)
Bit 5: NL = 112 (max possible = 112)
Bit 6: NL = 112 (max possible = 112)
Bit 7: NL = 112 (max possible = 112)
23.10.3 Differential Uniformity of the S-Box#
The differential uniformity \(\delta\) of a function \(S\) is:
For the AES S-box, \(\delta = 4\), which is the theoretical minimum for any power mapping over \(\mathrm{GF}(2^8)\). This means no differential characteristic through a single S-box can have probability greater than \(4/256 = 1/64\).
# Compute the Difference Distribution Table (DDT)
ddt = [[0] * 256 for _ in range(256)]
for x_val in range(256):
for dx in range(256):
x_hex = f'{x_val:02X}'
xdx_hex = f'{(x_val ^^ dx):02X}'
sy = int(g(SBox(h(x_hex))), 16)
sxdx = int(g(SBox(h(xdx_hex))), 16)
dy = sy ^^ sxdx
ddt[dx][dy] += 1
# Find differential uniformity (max entry excluding row dx=0)
delta = max(ddt[dx][dy] for dx in range(1, 256) for dy in range(256))
print(f"Differential uniformity of AES S-box: delta = {delta}")
print(f"Maximum differential probability: {delta}/256 = 1/{256 // delta}")
# Distribution of DDT entries (excluding row dx=0)
from collections import Counter
entries = [ddt[dx][dy] for dx in range(1, 256) for dy in range(256)]
dist = Counter(entries)
print("\nDistribution of DDT entries (excluding row dx=0):")
for val in sorted(dist.keys()):
print(f" Value {val}: appears {dist[val]} times")
Differential uniformity of AES S-box: delta = 4
Maximum differential probability: 4/256 = 1/64
Distribution of DDT entries (excluding row dx=0):
Value 0: appears 32895 times
Value 2: appears 32130 times
Value 4: appears 255 times
Security implications
The AES S-box has nonlinearity 112 (the maximum) and differential uniformity 4 (the minimum for power mappings). These two properties together make AES provably resistant to both linear cryptanalysis and differential cryptanalysis after a sufficient number of rounds.
# Visualise a portion of the DDT
sub_ddt = matrix(ZZ, [[ddt[dx][dy] for dy in range(32)] for dx in range(32)])
P_ddt = matrix_plot(sub_ddt, cmap='YlOrRd', colorbar=True)
show(P_ddt, figsize=8, title='AES S-Box DDT (upper-left 32x32 block)',
axes_labels=['Output difference (dy)', 'Input difference (dx)'])
23.10.4 Avalanche Effect in AES#
The avalanche effect states that flipping a single input bit should change approximately half of the output bits. We measure this for AES-128 by:
Encrypting a plaintext \(P\) to get ciphertext \(C\).
Flipping each of the 128 bits of \(P\) one at a time to get \(P'\).
Encrypting \(P'\) to get \(C'\).
Counting the Hamming distance \(d(C, C')\).
The ideal avalanche means \(d(C, C') \approx 64\) for each single-bit flip.
def state_to_int_list(state):
"""Convert a state (list of F256 rows) to a flat list of integers."""
result = []
for row in state:
for el in row:
result.append(int(g(el), 16))
return result
def int_list_to_state(ints):
"""Convert a flat list of 16 integers to a state (as column-major matrix)."""
m_inp = matrix(4, 4, [h(f'{v:02X}') for v in ints])
return list(transpose(m_inp))
def hamming_distance_bytes(a_list, b_list):
"""Hamming distance (in bits) between two byte sequences."""
dist = 0
for x_val, y_val in zip(a_list, b_list):
dist += bin(x_val ^^ y_val).count('1')
return dist
# Use a fixed test plaintext and key
set_random_seed(42)
pt_ints = [ZZ.random_element(0, 256) for _ in range(16)]
key_ints = [ZZ.random_element(0, 256) for _ in range(16)]
key_f256 = [h(f'{v:02X}') for v in key_ints]
pt_state = int_list_to_state(pt_ints)
ct_original = AESEncryption(pt_state, key_f256)
ct_orig_ints = state_to_int_list(ct_original)
# Flip each bit of the plaintext and measure avalanche
distances = []
for byte_idx in range(16):
for bit_idx in range(8):
pt_flipped = list(pt_ints) # copy
pt_flipped[byte_idx] ^^= (1 << bit_idx)
pt_flip_state = int_list_to_state(pt_flipped)
ct_flipped = AESEncryption(pt_flip_state, key_f256)
ct_flip_ints = state_to_int_list(ct_flipped)
distances.append(hamming_distance_bytes(ct_orig_ints, ct_flip_ints))
mean_d = sum(distances) / len(distances)
print(f"Mean Hamming distance: {float(mean_d):.2f} / 128 bits")
print(f"Min / Max: {min(distances)} / {max(distances)}")
# Plot
P_aval = bar_chart(distances, color='steelblue', width=1)
P_aval += line([(0, 64), (127, 64)], color='red', linestyle='--',
legend_label='Ideal (64 bits)')
show(P_aval, figsize=(12, 4),
title='AES-128 Avalanche: Single-Bit Flips',
axes_labels=['Flipped plaintext bit index', 'Hamming distance (bits)'])
Mean Hamming distance: 62.89 / 128 bits
Min / Max: 50 / 76
23.10.5 Round-by-Round Avalanche Build-up#
How quickly does AES achieve full diffusion? We measure the avalanche effect after each round by performing partial encryption.
Observation
After just 2 rounds, the avalanche effect is already close to 50% of bits. By round 4–5, full diffusion is essentially achieved. This rapid diffusion is a direct consequence of the wide trail strategy built into the ShiftRows + MixColumns combination.
def AESEncryptNRounds(inp, key, n_rounds):
"""Encrypt using only the first n_rounds of AES (for avalanche analysis)."""
expkey = KeyExpansion(key)
state = inp
state = AddRoundKey(state, RoundKey(expkey, 0))
for rnd in range(1, min(n_rounds + 1, Nr)):
state = SubBytes(state)
state = ShiftRows(state)
state = MixColumns(state)
state = AddRoundKey(state, RoundKey(expkey, rnd))
if n_rounds >= Nr:
# Final round (no MixColumns)
state = SubBytes(state)
state = ShiftRows(state)
state = AddRoundKey(state, RoundKey(expkey, Nr))
return [list(x) for x in state]
# Measure avalanche build-up round by round
set_random_seed(123)
pt_ints = [ZZ.random_element(0, 256) for _ in range(16)]
key_ints = [ZZ.random_element(0, 256) for _ in range(16)]
key_f256 = [h(f'{v:02X}') for v in key_ints]
pt_state = int_list_to_state(pt_ints)
# Flip bit 0 of byte 0
pt_flipped_ints = list(pt_ints)
pt_flipped_ints[0] ^^= 1
pt_flip_state = int_list_to_state(pt_flipped_ints)
round_distances = []
for n in range(1, Nr + 1):
ct1 = AESEncryptNRounds(pt_state, key_f256, n)
ct2 = AESEncryptNRounds(pt_flip_state, key_f256, n)
ct1_ints = state_to_int_list(ct1)
ct2_ints = state_to_int_list(ct2)
d = hamming_distance_bytes(ct1_ints, ct2_ints)
round_distances.append(d)
# Plot
P_rounds = list_plot(list(enumerate(round_distances, 1)), plotjoined=True,
color='darkgreen', marker='o', markersize=40,
legend_label='Hamming distance')
P_rounds += line([(1, 64), (Nr, 64)], color='red', linestyle='--',
legend_label='Ideal (64 bits)')
show(P_rounds, figsize=(8, 5),
title='Avalanche Build-up: Single Bit Flip Through AES Rounds',
axes_labels=['Number of rounds', 'Hamming distance (bits)'],
ymin=0, ymax=130)
for i, d in enumerate(round_distances, 1):
print(f"After round {int(i):2d}: {int(d):3d} / 128 bits changed ({float(100*d/128):.1f}%)")
After round 1: 21 / 128 bits changed (16.4%)
After round 2: 56 / 128 bits changed (43.8%)
After round 3: 67 / 128 bits changed (52.3%)
After round 4: 58 / 128 bits changed (45.3%)
After round 5: 61 / 128 bits changed (47.7%)
After round 6: 65 / 128 bits changed (50.8%)
After round 7: 52 / 128 bits changed (40.6%)
After round 8: 70 / 128 bits changed (54.7%)
After round 9: 52 / 128 bits changed (40.6%)
After round 10: 64 / 128 bits changed (50.0%)
23.11 Exercises#
Exercise 23.1: Verify the S-Box Algebraically
Compute the AES S-box entry for input 0xAB step by step using SageMath:
Compute \(\texttt{0xAB}^{-1}\) in \(\mathrm{GF}(2^8)\) using
h('AB')^(-1).Apply the affine transformation using the matrix
mand constanth('63').Verify that your result matches
SBox(h('AB')).
Hint
In SageMath you can compute the inverse directly: h('AB')^(-1).
Then extract the coefficient vector with l(inv, 8), multiply by the
affine matrix m, convert back to F256, and add h('63').
Exercise 23.2: MixColumns by Hand
For the column vector \([\texttt{0xDB}, \texttt{0x13}, \texttt{0x53}, \texttt{0x45}]^T\), compute the result of MixColumns using the polynomial representation.
The expected output is \([\texttt{0x8E}, \texttt{0x4D}, \texttt{0xA1}, \texttt{0xBC}]^T\).
Hint
Construct the column as col = vector(F256, [h('DB'), h('13'), h('53'), h('45')]),
form the polynomial p = col[0]*T^3 + col[1]*T^2 + col[2]*T + col[3],
and compute (p * cpol) % (T^4 + 1).
Exercise 23.3: Key Schedule Properties
Using KeyExpansion, verify that the AES key schedule
is not invertible from a single intermediate round key:
Show that knowing round key \(K_5\) (say) does not trivially give you \(K_0\).
However, show that knowing any round key \(K_i\) allows you to recover all round keys (and hence the master key) by running the key schedule both forward and backward.
Hint
To invert the key schedule: from \(W[i]\) you can recover \(W[i-N_k]\) using
\(W[i-N_k] = W[i] + T_i\) (addition in \(\mathrm{GF}(2^8)\) is its own inverse).
The S-box is invertible via SBoxinv, so RotWord and SubWord can be undone.
Exercise 23.4: Avalanche After Reduced Rounds
Modify the encryption function to perform only 1, 2, 3, or 4 rounds. For each variant, measure the avalanche effect over 500 random plaintext pairs (single-bit difference). At how many rounds does the mean Hamming distance first exceed 60 bits?
Hint
Use AESEncryptNRounds defined above. Loop over \(n \in \{1, 2, 3, 4\}\),
generate 500 random plaintexts, flip one random bit, encrypt both, and
compute the mean Hamming distance.
Exercise 23.5: Fixed Points of the S-Box
A fixed point of the S-box is a value \(x\) such that \(S(x) = x\).
How many fixed points does the AES S-box have?
How many values \(x\) satisfy \(S(x) = \overline{x}\) (complement)?
Compare with the expected number for a random permutation of 256 elements.
Hint
Simply iterate over all 256 values and check. For a random permutation of \(n\) elements, the expected number of fixed points is 1 (by the derangement formula: the probability of at least one fixed point approaches \(1 - 1/e \approx 0.632\)).
Exercise 23.6: S-Box as a Polynomial (SageMath Bonus)
Using the Lagrange interpolation polynomial of the S-box (computed in Experiment 23.10.1 above):
What is the algebraic degree of the S-box polynomial?
How many non-zero coefficients does it have?
Evaluate the polynomial at \(a^{10}\) and verify it matches \(S(a^{10})\).
This exercise highlights SageMath’s unique ability to study AES algebraically.
Hint
The Lagrange polynomial was stored in sbox_poly. Use sbox_poly.degree()
and len(sbox_poly.coefficients()). Evaluate with sbox_poly(a^10).
# Starter code for Exercise 23.5: Fixed Points of the S-Box
fixed_points = [i for i in range(256)
if SBox(h(f'{i:02X}')) == h(f'{i:02X}')]
complement_fixed = [i for i in range(256)
if SBox(h(f'{i:02X}')) == h(f'{(i ^^ 0xFF):02X}')]
print(f"Number of fixed points S(x) = x: {len(fixed_points)}")
if fixed_points:
print(f" Values: {[hex(x_val) for x_val in fixed_points]}")
print(f"Number of complement fixed S(x) = ~x: {len(complement_fixed)}")
if complement_fixed:
print(f" Values: {[hex(x_val) for x_val in complement_fixed]}")
print(f"\nExpected fixed points for random permutation of 256: ~1")
Number of fixed points S(x) = x: 0
Number of complement fixed S(x) = ~x: 0
Expected fixed points for random permutation of 256: ~1
23.12 Summary#
In this chapter we built a complete AES-128 implementation using SageMath’s native finite field and polynomial ring infrastructure:
GF(2^8) arithmetic – using
GF(2^8, 'a', modulus=x^8+x^4+x^3+x+1), which replaces all manual polynomial arithmetic with exact algebraic operations.S-box construction from first principles: multiplicative inversion via
el^(-1)in the finite field, followed by an affine transformation using Sage’smatrix()andvector(). We verified maximum nonlinearity (112) and minimum differential uniformity (4).Four transformations – SubBytes, ShiftRows, MixColumns, AddRoundKey – and their inverses. MixColumns uses polynomial multiplication in
PolynomialRing(F256)modulo \(T^4 + 1\), directly mirroring the mathematical specification.Key expansion for AES-128, using the field generator
afor round constants.Algebraic analysis unique to SageMath: Lagrange interpolation of the S-box as a polynomial over \(\mathrm{GF}(2^8)\), enabling direct study of algebraic degree.
Experimental analysis showing excellent avalanche properties and rapid diffusion.
Key Takeaways#
AES is a substitution-permutation network (SPN), fundamentally different from the Feistel structure of DES.
SageMath’s
GF(2^8)makes the algebraic structure transparent – operations that require manual bit-level coding in Python are single-line expressions in Sage.The wide trail strategy guarantees provable resistance to differential and linear cryptanalysis.
AES remains unbroken after more than two decades; the best known attacks are only marginally better than brute force.
References#
Daemen, J., & Rijmen, V. (2002). The Design of Rijndael: AES – The Advanced Encryption Standard. Springer-Verlag.
NIST (2001). FIPS 197: Advanced Encryption Standard (AES). https://csrc.nist.gov/publications/detail/fips/197/final
Daemen, J., & Rijmen, V. (1999). AES Proposal: Rijndael. NIST AES submission.
Heys, H. M. (2002). A Tutorial on Linear and Differential Cryptanalysis. Cryptologia, 26(3), 189–221.
Stinson, D. R., & Paterson, M. (2018). Cryptography: Theory and Practice, 4th edition. CRC Press. Chapter 4.
SageMath Documentation. Finite Fields. https://doc.sagemath.org/html/en/reference/finite_rings/