% galois.m.
% Matlab 7.0 functions to find bit-level equations to calculate  
% output of the MDS and RS operations in Twofish algorithm.
% Writes to file, equations, and appends if there is already a file
% named equations.
% John Bryan 2004

function z1
   clear all; clear java; clear classes;
   whos
   syms x;
   x_input='a7*x^7+a6*x^6+a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x^1+a0*x^0';
   x_input_sym=sym(x_input);
   mds = char('5b','ef');
   rs = char('01','a4','55','87','5a','58','db','9e', ...
             '56','82','f3','1e','c6','68','e5','02', ...
             'a1','fc','c1','47','ae','3d','19','03');
   one(mds,x_input,x_input_sym,x,0);
   one(rs,x_input,x_input_sym,x,1);
return

function one(matrix,x_input,x_input_sym,x,mds_rs)
   polynomials(1,:) = blanks(24);
   polynomials_sym(1,:) = sym(hex2polynomial(matrix(1,:)));
   matrix_input_product(1,:)=times(polynomials_sym(1,:),x_input_sym);
   collected_matrix_input_product(1,:)=collect(matrix_input_product(1,:));
   calculate(collected_matrix_input_product(1,:),x,matrix(1,:),mds_rs);
   m=1;
   for n=2:length(matrix),
      m=m+1;
      polynomials_sym(m,:) = sym(hex2polynomial(matrix(m,:)));
      matrix_input_product(m,:)=times(polynomials_sym(m,:),x_input_sym);
      collected_matrix_input_product(m,:)=collect(matrix_input_product(m,:));
      calculate(collected_matrix_input_product(m,:),x,matrix(n,:),mds_rs);
   end
return

function calculate(collected_product,x,rs,mds_rs)
   if (mds_rs==0)
      x8='x^6+x^5+x^3+1';           x8sym=sym(x8);
      x9='x^7+x^6+x^4+x';           x9sym=sym(x9);
      x10='x^7+x^6+x^3+x^2+1';      x10sym=sym(x10);
      x11='x^7+x^6+x^5+x^4+x+1';    x11sym=sym(x11);
      x12='x^7+x^3+x^2+x+1';        x12sym=sym(x12);
      x13='x^6+x^5+x^4+x^2+x+1';    x13sym=sym(x13);
      x14='x^7+x^6+x^5+x^3+x^2+x';  x14sym=sym(x14);
   else
      x8='x^6+x^3+x^2+1';           x8sym=sym(x8);
      x9='x^7+x^4+x^3+x';           x9sym=sym(x9);
      x10='x^6+x^5+x^4+x^3+1';      x10sym=sym(x10);
      x11='x^7+x^6+x^5+x^4+x';      x11sym=sym(x11);
      x12='x^7+x^5+x^3+1';          x12sym=sym(x12);
      x13='x^4+x^3+x^2+x+1';        x13sym=sym(x13);
      x14='x^5+x^4+x^3+x^2+x';      x14sym=sym(x14);
   end
   substitute_for_x8=subs(collected_product,x^8,x8sym);
   substitute_for_x9=subs(substitute_for_x8,x^9,x9sym);
   substitute_for_x10=subs(substitute_for_x9,x^10,x10sym);
   substitute_for_x11=subs(substitute_for_x10,x^11,x11sym);
   substitute_for_x12=subs(substitute_for_x11,x^12,x12sym);
   substitute_for_x13=subs(substitute_for_x12,x^13,x13sym);
   substitute_for_x14=subs(substitute_for_x13,x^14,x14sym);
   substituted_ef_input_product=collect(substitute_for_x14);
   [coefficients,x_power]=coeffs(substituted_ef_input_product,x);
   for i=1:8
      final_odd_multiple_coefficients(i,:)=blanks(56);
   end
   for i=1:8,
      temp=char(coefficients(i));
      for j=1:length(temp),
         odd_multiple_coefficients(i,j)=temp(j);
      end
   end
   [d1,d2]=size(odd_multiple_coefficients);
   for i=1:d1
      for j=1:(d2-3)
         if ((odd_multiple_coefficients(i,j)=='2') || ...
             (odd_multiple_coefficients(i,j)=='4'))
            if (odd_multiple_coefficients(i,j+1)=='*')
               odd_multiple_coefficients(i,j:j+4)=blanks(5);
            end
         end
         if ((odd_multiple_coefficients(i,j)=='3') || ...
            (odd_multiple_coefficients(i,j)=='5'))
            if (odd_multiple_coefficients(i,j+1)=='*')
               odd_multiple_coefficients(i,j:j+1)=blanks(2);
            end
         end
      end
      for j=1:(d2)
         if (odd_multiple_coefficients(i,j)=='+')
               odd_multiple_coefficients(i,j)=blanks(1);
         end
      end
   end
   for i=1:d1
      a_flag=0;
      for j=1:(d2-1)
         if (odd_multiple_coefficients(i,j:j+1)=='a0')
            if (a_flag==1)
               final_odd_multiple_coefficients(i,1:5)=' x.b0';
            else
               final_odd_multiple_coefficients(i,1:5)=' x.b0';
               a_flag=1;
            end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a1')
            if (a_flag==1)
               final_odd_multiple_coefficients(i,6:12)='   x.b1';
            else
               final_odd_multiple_coefficients(i,6:12)='   x.b1';
               a_flag=1;
            end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a2')
            if (a_flag==1)
                final_odd_multiple_coefficients(i,13:19)='   x.b2';
             end
             if (a_flag==0)
                final_odd_multiple_coefficients(i,13:19)='   x.b2';
                a_flag=1;
             end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a3')
             if (a_flag==1)
                final_odd_multiple_coefficients(i,20:26)='   x.b3';
             else
                final_odd_multiple_coefficients(i,20:26)='   x.b3';
                a_flag=1;
             end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a4')
             if a_flag~=0
                final_odd_multiple_coefficients(i,27:33)='   x.b4';
             else
                final_odd_multiple_coefficients(i,27:33)='   x.b4';
                a_flag=1;
             end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a5')
             if a_flag~=0
                final_odd_multiple_coefficients(i,34:40)='   x.b5';
             else
                final_odd_multiple_coefficients(i,34:40)='   x.b5';
                a_flag=1;
             end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a6')
             if a_flag~=0
                final_odd_multiple_coefficients(i,41:47)='   x.b6';
             else
                final_odd_multiple_coefficients(i,41:47)='   x.b6';
                a_flag=1;
             end
         end
         if (odd_multiple_coefficients(i,j:j+1)=='a7')
            if a_flag~=0
               final_odd_multiple_coefficients(i,48:54)='   x.b7';
            else
               final_odd_multiple_coefficients(i,48:54)='   x.b7';
               a_flag=1;
            end
         end
      end
      final_odd_multiple_coefficients(i,55:56)=' ;';
   end
   x_power_char=char(x_power);
   if (x_power_char(10)=='1')
      if (mds_rs==0)
         zz(1,:) = ['mds',rs,'.b0 = '];
         x_superscript(1)=0;
      else
        zz(1,:) = ['rs_',rs,'.b0 = '];
        x_superscript(1)=0;
      end
   end
   if (x_power_char(13)==',')
      if (mds_rs==0)
         zz(2,:) = ['mds',rs,'.b1 = '];
         x_superscript(2)=1;
      else
         zz(2,:) = ['rs_',rs,'.b1 = '];
         x_superscript(2)=1;
      end
   end
   j=3;
   for i=4:9
      switch x_power_char(4*i)
         case '2',
            if (mds_rs==0)
               zz(3,:) = ['mds',rs,'.b2 = '];
               x_superscript(j)=2;
            else
               zz(3,:) = ['rs_',rs,'.b2 = '];
               x_superscript(j)=2;
            end
         case '3',
            if (mds_rs==0)
               zz(4,:) = ['mds',rs,'.b3 = '];
               x_superscript(j)=3;
            else
               zz(4,:) = ['rs_',rs,'.b3 = '];
               x_superscript(j)=3;
            end
         case '4',
            if (mds_rs==0)
               zz(5,:) = ['mds',rs,'.b4 = '];
               x_superscript(j)=4;
            else
               zz(5,:) = ['rs_',rs,'.b4 = '];
               x_superscript(j)=4;
            end
         case '5',
            if (mds_rs==0)
               zz(6,:) = ['mds',rs,'.b5 = '];
               x_superscript(j)=5;
            else
               zz(6,:) = ['rs_',rs,'.b5 = '];
               x_superscript(j)=5;
            end
         case '6',
            if (mds_rs==0)
               zz(7,:) = ['mds',rs,'.b6 = '];
               x_superscript(j)=6;
            else
               zz(7,:) = ['rs_',rs,'.b6 = '];
               x_superscript(j)=6;
            end
         case '7',
            if (mds_rs==0)
               zz(8,:) = ['mds',rs,'.b7 = '];
               x_superscript(j)=7;
            else
               zz(8,:) = ['rs_',rs,'.b7 = '];
               x_superscript(j)=7;
            end
      end
      j=j+1;
   end
   for i=1:8,
      [d3,d4]=size(final_odd_multiple_coefficients);
      flag=0;
      for j=1:d4-3,
         if ((final_odd_multiple_coefficients(i,j) =='x') && (flag == 0)),
            flag=1;
         elseif ((final_odd_multiple_coefficients(i,j) =='x') && (flag == 1)),
                final_odd_multiple_coefficients(i,j-2) ='^';
         end
      end
   end
   for i=1:8
      switch x_superscript(i)
         case 0,
            xs(1,:)=final_odd_multiple_coefficients(i,:);
         case 1,
            xs(2,:)=final_odd_multiple_coefficients(i,:);
         case 2,
            xs(3,:)=final_odd_multiple_coefficients(i,:);
         case 3,
            xs(4,:)=final_odd_multiple_coefficients(i,:);
         case 4,
            xs(5,:)=final_odd_multiple_coefficients(i,:);
         case 5,
            xs(6,:)=final_odd_multiple_coefficients(i,:);
         case 6,
            xs(7,:)=final_odd_multiple_coefficients(i,:);
         case 7,
            xs(8,:)=final_odd_multiple_coefficients(i,:);
      end
   end
   for i=1:8
      zz_char=char(zz);
      xs_char=char(xs);
   end
   ef_eq=strcat(zz_char,xs_char);
   fid=fopen('equations','a+');
   for i=1:8
        for j=1:66
           if ~isspace(ef_eq(i,j))
             fprintf(fid,'%1c',ef_eq(i,j));
           else
             fprintf(fid,' ');
           end
        end
        fprintf(fid,'\n');
   end;
   fprintf(fid,'\n \n');
   fclose(fid);
return

function polynomial = hex2polynomial(polynomial_input)
   bit_value=dec2bin(hex2dec(polynomial_input),8);
   flag0=0;
   j=1;
   if (bit_value(1)=='1')
      polynomial_term(j:j+2) ='x^7';
      j=j+3;
      flag0=1;
   end
   if (bit_value(2)=='1')
      if (flag0==1)
         polynomial_term(j:j+3) ='+x^6';
         j=j+4;
      else
         polynomial_term(j:j+2) ='x^6';
         j=j+3;
         flag0=1;
      end
   end
   if (bit_value(3)=='1')
      if (flag0==1)
         polynomial_term(j:j+3) ='+x^5';
         j=j+4;
      else
         polynomial_term(j:j+2) ='x^5';
         j=j+3;
         flag0=1;
      end
   end
   if (bit_value(4)=='1')
      if (flag0==1)
         polynomial_term(j:j+3) ='+x^4';
         j=j+4;
      else
         polynomial_term(j:j+2) ='x^4';
         j=j+3;
         flag0=1;
      end
   end
   if (bit_value(5)=='1')
      if (flag0==1)
         polynomial_term(j:j+3) ='+x^3';
         j=j+4;
      else
         polynomial_term(j:j+2) ='x^3';
         j=j+3;
         flag0=1;
      end
   end
   if (bit_value(6)=='1')
      if (flag0==1)
         polynomial_term(j:j+3) ='+x^2';
         j=j+4;
      else
         polynomial_term(j:j+2) ='x^2';
         j=j+3;
         flag0=1;
      end
   end
   if (bit_value(7)=='1')
      if (flag0==1)
         polynomial_term(j:j+1) ='+x';
         j=j+2;
      else
         polynomial_term(j) ='x';
         j=j+1;
         flag0=1;
      end
   end
   if (bit_value(8)=='1')
      if (flag0==1)
         polynomial_term(j:j+1) ='+1';
         j=j+2;
      else
         polynomial_term(j) ='1';
         j=j+1;
      end
   end
   polynomial = polynomial_term;
return