Nimber multiplication algorithm in Matlab

function y = nimprod(x1,x2)
% nimber multiplication of bin objects  ( for the bin class see http://pastebin.com/K02wGsae )



if x1==bin('0')  ||  x2==bin('0')
% if x1 or x2 is 0

    y = bin('0') ;

elseif x1==bin('1')  ||  x2==bin('1')
% if x1 or x2 is 1

    y = x1*x2 ;

elseif length(x1.Expo)==1 && length(x2.Expo)==1 && ~mod( log2(double(x1.Expo)) ,1) && ~mod( log2(double(x2.Expo)) ,1)
% if x1 and x2 are Fermat powers of 2

    if   x1 ~= x2
         y = x1*x2 ;  % Nim product of distinct Fermat powers of 2 is the normal product.
    else
         product = bin('11') * x1 ;     % Nim product of equal Fermat powers of 2 is 3/2 * entry.
         y = bin( product.Expo - 1 ) ;  % Division by two is rightshift.
    end

elseif length(x1.Expo)==1 && length(x2.Expo)==1
% if x1 and x2 are powers of 2 (although at least one is not a Fermat power of 2)

    % The Expo of powers of 2 is a single number, which now has to be converted to a bin object.
    e1 = bin( x1.Expo ,'dec') ;
    e2 = bin( x2.Expo ,'dec') ;
    xor = bitxor(e1,e2) ;
    and = bitand(e1,e2) ;

    XOR = bin('1') ;
    for m = 1 : xor.Weight
        XOR = XOR * bin( 2^(xor.Expo(m)) ) ;
    end

    if and.Weight==1
        % AND = 3/2 * bin( 2^(and.Expo) )
        product = bin('11') * bin( 2^(and.Expo) ) ;
        AND = bin( product.Expo - 1 ) ;  % Division by two is rightshift.
    else
        AND = bin('1') ;
        for m = 1 : and.Weight
            % AND = nimprod(   AND   ,   3/2 * bin( 2^(and.Expo(m)) )   )
            product = bin('11') * bin( 2^(and.Expo(m)) ) ;
            factor = bin( product.Expo - 1 ) ;  % Division by two is rightshift.
            AND = nimprod( AND , factor ) ;
        end
    end

    y = nimprod( XOR , AND ) ;

else  % if x1 or x2 is not a power of 2.

    Y = bin('0') ;

    for m = 1 : x1.Weight
        for n = 1 : x2.Weight
            M = x1.Expo(m) ;  N = x2.Expo(n) ;
            if M==0 || N==0                     % if bin(M)=1 or bin(N)=1 ,  i.e. if M=0 or N=0
                Y = bitxor( Y ,   bin(M)*bin(N)   ) ;
            elseif ~mod( log2(double(M)) ,1) && ~mod( log2(double(N)) ,1)  % if bin(M) and bin(N) are Fermat powers of 2 ,  i.e. if M and N are powers of 2
                if M ~= N
                    Y = bitxor( Y ,   bin(M)*bin(N)   ) ;             % Nim product of distinct Fermat powers of 2 is the normal product.
                else
                    product = bin('11') * bin(M) ;                    % Nim product of   equal  Fermat powers of 2 is 3/2 * entry.
                    Y = bitxor( Y ,   bin( product.Expo - 1 )   ) ;   % Division by two is rightshift.
                end
            else
                binM = bin(M) ;
                binN = bin(N) ;
                Y = bitxor( Y ,     nimprod(  bin(binM.Expo)  ,  bin(binN.Expo)  )     ) ;
            end
        end
    end

    y = Y ;

end



end