Advanced Encryption Standard (AES) polynomial multiplication in octave / matlab

Question

I am trying to perform polynomial multiplication as explained in the Advanced Encryption Standard (AES) draft.

Here's what I have so far:

function  y  =  multiply_x2( a )

  % a is a 1 x 8 binary vector 
  % y is a 1 x 8 binary vector 
  % t is a 1 x 8 binary vector corresponding to the "AES irreducible polynomial"

  y  =  [ a(2:8) 0 ];     % Multiply byte 'a' by 2 using "left shift" 
  t  =  [ 0 0 0 a(1) a(1) 0 a(1) a(1) ]; %  't' only becomes the "AES irreducible
                                         %  polynomial" if a(1) == 1, otherwise
                                         %  it becomes the "zero" byte array
  y  =  mod( y + t , 2 ) ;   % XOR operation on y and t, as described in AES.                
end

The code above is for

"y = {02} . a"

(where "{}" denotes hexadecimal notation, whose binary representation can be interpreted as the presense of the respective power of x in a polynomial. For example, {02} corresponds to 00000010, which corresponds to the polynomial "x", {06} would correspond to "x²+x", etc, as per the AES docs)

I want to multiply a with 0e , 09 , 0d, and 0b. How will the code be for each of them? i.e. for:

"y= ( {0e} . a )"
"y= ( {09} . a )"
"y= ( {0d} . a )"
"y= ( {02} . a )"

Answer 1

Thank you for your interest. I am trying to implement AES in matlab. And I found the solution in the pages of http://csrc.nist.gov/publications/fips/fips197/fips-197.pdf . Here is the function for the y= ( {09} . a ) multiply operation ;

function  y  =  multiply_x9( a )

% Multiply byte A by 9 over finite field of AES
y2  =  multiply_x2( a ) ;
y4  =  multiply_x2( y2 ) ;
y8  =  multiply_x2( y4 ) ;

y  =  mod( y8 + a , 2 ) ;

end

And also for the any matrix multiplication multiply_xx (a, b, p ) function can be used. Here is the function;

function  y  =  multiply_xx (a, b, p ) 


% Determine input lengths and check if they are equal
n  =  length( a ) ;
if  ( n ~= length(b) ) ,
    error( 'Operand input lengths not equal to each other!' ) ;
elseif  ( n ~= length(p) ) ,
    error( 'Operand input lengths not equal to modulus length!' ) ;
end


% Initialize result to zeros and start iteration row by row
y  =  zeros( 1 , n ) ;
for  i  =  1  :  n ,
    m  =  a(i) * b ;
    y  =  mod( y(1) * p  +  [ y(2:n) 0 ]  +  m , 2 ) ;
end


end

Answer 2

In MATLAB/Octave, conv and deconv are respectively correspond to multiplication and (division/modulo) operations for polynomials.

function out = multiply(A, x)
    mult = mod(conv(A,x), 2); % poynomial multiplication
    [~, modulo] = deconv(mult, [1 0 0 0 1 1 0 1 1]); %modulo operation
    out = mod(modulo(end-7:end) , 2); %extract last 8 bits
end

For example to multiply 0x57 and 0x13

a = [1 0 1 0 1 1 1]; %0x57
b = [1 0 0 1 1]; %0x13
result = multiply(a,b)

result = 
    1   1   1   1   1   1   1   0

that is binary representation of 0xFE

Answer 3

This was an interesting problem. Here is a general solution for multiplication as defined in the AES doc you linked. I use the name xtime for the {02} multiplication, and I implement "addition" (xadd) as an XOR operation (i.e. !=) directly, since this is easier.

Helper functions:

%% in file isByte.m
function Out = isByte (a)
  a = a(:).'; % ensure horizontal vector representation
  if all (a == 0 | a == 1) && length (a) == 8; Out = true; return; end
  Out = false;
end

%% in file byte2hex.m
function Out = byte2hex (a)
  a = a(:).'; % ensure horizontal vector
  assert (isByte (a), 'Input needs to be a valid "byte" array');
  Out = sum (a .* ([2,2,2,2,2,2,2,2] .^ [7:-1:0])); Out = dec2hex (Out);  
end

%% in file hex2byte.m
function Out = hex2byte (h)
  assert (isxdigit (h) && hex2dec (h) < 256, 'Input needs to be a valid hex number below FF');
  Out = dec2bin (hex2dec (h));  Out = sprintf ('%8s', Out); 
  Out = Out - 48             ;  Out(Out == -16) = 0; % convert spaces to zeros
end

Polynomial functions:

%% in file xadd.m
function Out = xadd (a, b)
  a = a(:).'; b = b(:).'; % ensure horizontal vector representations
  assert (isByte (a) && isByte (b), 'Inputs need to be valid "byte" arrays');
  Out = a != b;
end

%% in file xtime.m
function Out = xtime (a)
  a = a(:).'; % ensure horizontal vector
  assert (isByte (a), 'Input needs to be a valid "byte" array');
  Out = [a(2 : 8), 0]; % left shift
  if a(1) == 1
    Out = xadd (Out, [0, 0, 0, 1, 1, 0, 1, 1]);  % subtract (= add) the AES m(x) polynomial
  end  
end

% in file xmultiply.m
function Out = xmultiply (a, b)
  a = a(:)'; b = b(:)'; % ensure horizontal vector representations
  assert (isByte(a) && isByte(b), 'Inputs need to be valid "byte" arrays');

  Out = [0,0,0,0,0,0,0,0];
  if a == [0, 0, 0, 0, 0, 0, 0, 0] || b == [ 0, 0, 0, 0, 0, 0, 0, 0]; return; end
  if b(8) == 1; Out = xadd (Out, a); end % yes this could be done recursively but, why bother.
  if b(7) == 1; Out = xadd (Out, xtime (a)); end
  if b(6) == 1; Out = xadd (Out, xtime (xtime (a))); end
  if b(5) == 1; Out = xadd (Out, xtime (xtime (xtime (a)))); end
  if b(4) == 1; Out = xadd (Out, xtime (xtime (xtime (xtime (a))))); end
  if b(3) == 1; Out = xadd (Out, xtime (xtime (xtime (xtime (xtime (a)))))); end
  if b(2) == 1; Out = xadd (Out, xtime (xtime (xtime (xtime (xtime (xtime (a))))))); end
  if b(1) == 1; Out = xadd (Out, xtime (xtime (xtime (xtime (xtime (xtime (xtime (a)))))))); end  
end

Example use: (same example as in the AES doc)

octave:1> a = hex2byte("57")
  a =
     0   1   0   1   0   1   1   1
octave:2> b = hex2byte("13")
  b =
     0   0   0   1   0   0   1   1
octave:3> c = xmultiply(a, b)
  c =
     1  1  1  1  1  1  1  0
octave:4> byte2hex(c)
  ans = FE