Evaluating Algebraic Complexity of a S-box While studying the design and the desirable properties of an AES S-box , I came to know that Algebraic Complexity is also an important property of an S-box which is usually considered while evaluating the properties of an S-box. 

After reading several papers (such as APA), I get to know that Algebraic Complexity is basically the number of terms appearing in the "Algebraic Equation" representing that S-box.

After doing a bit of a research (source) I get to know there are three main methods for finding the polynomial equation to represent S-box which are: 


Lagrangian interpolation
Polynomial linearization
Q-ary polynomial deduction


But honestly I am unable to understand any of the above methods to find the algebraic equation of the given S-box. So that ultimately I can evaluate it's Algebraic Complexity.

So, if someone can explain it in detail (or algorithm for the same), it would be really helpful for me to write the code for it. Or if there is some other way (or codes) to find Algebraic Complexity or Algebraic Equation of an AES S-box please share it.

TL;DR: Given an S-box ($n \times n$) , how can we evaluate its Algebraic Equation or Algebraic Complexity?
The interpolation of Sbox is based on Lagrange algorithm.
the following three steps are used to calculate the interpolation of S-box:


represent the S-box values as polynomial system, where P and C are polynomial representation of S-box input and output.
M (N x N ) is a matrix composed of P power polynomial. N is sbox size (256 for 8-bit sbox)
solve B form the following equation: M.B=C 



# the Equation is M.B=C.
#     1  p0  p0^2  p0^3  p0^4  p0^5 .......  p0^d   
#     1  p1  p1^2  p1^3  p1^4  p1^5 .......  p1^d
#  M= 1  p2  p2^2  p2^3  p2^4  p2^5 .......  p2^d
#     1  p3  p3^2  p3^3  p3^4  p3^5 .......  p3^d   
#       .
#       .
#       .
#     1  pn  pn^2  pn^3  pn^4  pn^5 .......  pn^d 
#  size of M is dXn , d and n are equal for sbox, they are GF (2^8) (8-bit sbox)
# Calculation steps. 
#   1. convert Sbox table to polynomial of input and output, vectors (P is input  and C is output)
#   2. calculate M table  
#   3. B=M.solve_right(C)
#



The following code is implementation of interpolation representation of Rijndael S-box

f=x^8+x^4+x^3+x^1+1  # irreducible polynomial
L.&lt;z&gt;= GF(2^8,modulus=f)
R=PolynomialRing(L,'x')
N=256
M=matrix(L,N,N)
S=[0x63, 0x7C, 0x77, 0x7B, 0xF2, 0x6B, 0x6F, 0xC5, 0x30, 0x01, 0x67, 0x2B, 0xFE, 0xD7, 0xAB, 0x76, 
                             0xCA, 0x82, 0xC9, 0x7D, 0xFA, 0x59, 0x47, 0xF0, 0xAD, 0xD4, 0xA2, 0xAF, 0x9C, 0xA4, 0x72, 0xC0,
                             0xB7, 0xFD, 0x93, 0x26, 0x36, 0x3F, 0xF7, 0xCC, 0x34, 0xA5, 0xE5, 0xF1, 0x71, 0xD8, 0x31, 0x15,
                             0x04, 0xC7, 0x23, 0xC3, 0x18, 0x96, 0x05, 0x9A, 0x07, 0x12, 0x80, 0xE2, 0xEB, 0x27, 0xB2, 0x75,
                             0x09, 0x83, 0x2C, 0x1A, 0x1B, 0x6E, 0x5A, 0xA0, 0x52, 0x3B, 0xD6, 0xB3, 0x29, 0xE3, 0x2F, 0x84,
                             0x53, 0xD1, 0x00, 0xED, 0x20, 0xFC, 0xB1, 0x5B, 0x6A, 0xCB, 0xBE, 0x39, 0x4A, 0x4C, 0x58, 0xCF,
                             0xD0, 0xEF, 0xAA, 0xFB, 0x43, 0x4D, 0x33, 0x85, 0x45, 0xF9, 0x02, 0x7F, 0x50, 0x3C, 0x9F, 0xA8,
                             0x51, 0xA3, 0x40, 0x8F, 0x92, 0x9D, 0x38, 0xF5, 0xBC, 0xB6, 0xDA, 0x21, 0x10, 0xFF, 0xF3, 0xD2,
                             0xCD, 0x0C, 0x13, 0xEC, 0x5F, 0x97, 0x44, 0x17, 0xC4, 0xA7, 0x7E, 0x3D, 0x64, 0x5D, 0x19, 0x73,
                             0x60, 0x81, 0x4F, 0xDC, 0x22, 0x2A, 0x90, 0x88, 0x46, 0xEE, 0xB8, 0x14, 0xDE, 0x5E, 0x0B, 0xDB,
                             0xE0, 0x32, 0x3A, 0x0A, 0x49, 0x06, 0x24, 0x5C, 0xC2, 0xD3, 0xAC, 0x62, 0x91, 0x95, 0xE4, 0x79,
                             0xE7, 0xC8, 0x37, 0x6D, 0x8D, 0xD5, 0x4E, 0xA9, 0x6C, 0x56, 0xF4, 0xEA, 0x65, 0x7A, 0xAE, 0x08,
                             0xBA, 0x78, 0x25, 0x2E, 0x1C, 0xA6, 0xB4, 0xC6, 0xE8, 0xDD, 0x74, 0x1F, 0x4B, 0xBD, 0x8B, 0x8A,
                             0x70, 0x3E, 0xB5, 0x66, 0x48, 0x03, 0xF6, 0x0E, 0x61, 0x35, 0x57, 0xB9, 0x86, 0xC1, 0x1D, 0x9E,
                             0xE1, 0xF8, 0x98, 0x11, 0x69, 0xD9, 0x8E, 0x94, 0x9B, 0x1E, 0x87, 0xE9, 0xCE, 0x55, 0x28, 0xDF,
                             0x8C, 0xA1, 0x89, 0x0D, 0xBF, 0xE6, 0x42, 0x68, 0x41, 0x99, 0x2D, 0x0F, 0xB0, 0x54, 0xBB, 0x16] # rajS-box 

C=vector(L,N)
P=vector(L,N)

# step 1
for i in range(0,N):

    P[i]= (i) + z^1*(i&gt;&gt;1) + z^2*(i&gt;&gt;2) + z^3*(i&gt;&gt;3) + z^4*(i&gt;&gt;4) + z^5*(i&gt;&gt;5) + z^6*(i&gt;&gt;6) + z^7*(i&gt;&gt;7)
    C[i]= S[(i)] + z^1*(S[i]&gt;&gt;1) + z^2*(S[i]&gt;&gt;2) + z^3*(S[i]&gt;&gt;3) + z^4*(S[i]&gt;&gt;4) + z^5*(S[i]&gt;&gt;5) + z^6*(S[i]&gt;&gt;6) + z^7*(S[i]&gt;&gt;7)


#print P
#print C
# step 2
for j in range(0,N):
    for i in range(0,N):
        M[j,i]= P[j]^i





#step 3
B=M.solve_right(C)



print B