function MSpolyR = multiply_MSpoly(MSpolyA,MSpolyB,ech)
%function MSpolyR = multiply_MSpoly(MSpolyA,MSpolyB,ech)
% multiply two polynomials with parameters in multi-simplexes.
%
% inputs: MSpolyA -> polynomial to the left.
%         MSpolyB -> polynomial to the right.
%         ech     -> when exists, the trace of execution is shown on the screen.
%
% outputs: MSpolyR -> resulting polynomial
%
% Example: P(a) = a_1 + a_2;  Q(b)= b_1 + b_2
% Aim: P(a)*Q(b)
% MSpolyA = create_MSpolynomial(1,1,1,1,2,'P')
% MSpolyB = create_MSpolynomial(1,1,2,1,2,'Q')
% MSpolyR = multiply_MSpoly(MSpolyA,MSpolyB,1)
%
% Date: 19/02/2008
% Author: ricfow@dt.fee.unicamp.br

if nargin == 2
    ech = 0;
end
if ech == 1
    fprintf('multiplying... (');
end
t=clock;
MSpolyR.name = strcat(MSpolyA.name,MSpolyB.name);
MSpolyR.nmon = 0;
MSpolyR.nxm = MSpolyA.nxm;
[MSpolyR.var,MSpolyR.N,MSpolyR.deg] = determine_final_N(MSpolyA.var,MSpolyA.N,MSpolyA.deg,MSpolyB.var,MSpolyB.N,MSpolyB.deg);

MSpolyR.mon = [];
MSpolyR.nmon = 1;
for i=1:MSpolyA.nmon
    for j=1:MSpolyB.nmon
        monR = multiply_two_monomials(MSpolyA.mon(:,i),MSpolyA.var,MSpolyA.N,MSpolyB.mon(:,j),MSpolyB.var,MSpolyB.N,MSpolyR.var,MSpolyR.N);        
        MSpolyR.coef{MSpolyR.nmon}.coef_val = MSpolyA.coef{i}.coef_val*MSpolyB.coef{j}.coef_val;
        MSpolyR.mon = [MSpolyR.mon  monR];
        MSpolyR.nmon = MSpolyR.nmon + 1;
    end
end
MSpolyR.nmon = MSpolyR.nmon - 1;
t = etime(clock,t);
if ech == 1
    fprintf('%.2f s) (%d mon)\n',t,MSpolyR.nmon);
    fprintf('simplifying... (',t);
end
t=clock;
MSpolyR = simplify_polynomials(MSpolyR);
t=etime(clock,t);
if ech == 1
    fprintf('%.2f s) (%d mon)\n',t,MSpolyR.nmon);    
end
MSpolyR = reorder_monomials(MSpolyR);


% Various auxiliary routines
%
%________________________________________________________________________________________________________
function [varf,Nf,degf] = determine_final_N(var1,N1,deg1,var2,N2,deg2)

varf = var1;
Nf = N1;
degf = deg1;
for i=1:size(var2,2)
    r = (varf == var2(1,i));
    if sum(r) == 0
        varf = [varf var2(1,i)];
        Nf = [Nf N2(1,i)];
        degf = [degf deg2(1,i)];
    else
        [r,c]=find(varf == var2(1,i));
        degf(1,c) = degf(1,c) + deg2(1,i);
    end
end

%________________________________________________________________________________________________________
function monR = multiply_two_monomials(mon1,var1,N1,mon2,var2,N2,varf,Nf)

monR = zeros(sum(Nf),1);
monR(1:sum(N1),1) = mon1;
pos = 1;
for i=1:size(N2,2)
    tuple_now = mon2(pos:pos+N2(i)-1,:);
    monR = add_term_monomial(monR,varf,Nf,tuple_now,var2(1,i));
    pos = pos + N2(i);
end

%________________________________________________________________________________________________________
function monf = add_term_monomial(monf,varf,Nf,term,varterm)

pos = 1;
for i=1:size(Nf,2)
    tuple_now = monf(pos:pos+Nf(i)-1,:);
    if varf(i) == varterm
        tuple_now = tuple_now + term;
        monf(pos:pos+Nf(i)-1,:) = tuple_now;
        break;
    end
    pos = pos + Nf(i);
end

%________________________________________________________________________________________________________
function MSpolyR = simplify_polynomials(MSpoly)

MSpolyR.name = MSpoly.name;
MSpolyR.deg = MSpoly.deg;
MSpolyR.nxm = MSpoly.nxm;
MSpolyR.N = MSpoly.N;
MSpolyR.var = MSpoly.var;
MSpolyR.nmon = 1;
MSpolyR.mon = [];
out = ones(MSpoly.nmon,1);
for i=1:MSpoly.nmon

    if out(i) == 0
        continue;
    end
    NMON = MSpoly.mon(:,i);    
    coef_val = MSpoly.coef{i}.coef_val;
    for j=i+1:MSpoly.nmon
        if out(j) == 0
            continue;
        end
        if sum(NMON == MSpoly.mon(:,j)) == size(NMON,1)            
            coef_val = coef_val + MSpoly.coef{j}.coef_val;
            out(j) = 0;
        end
    end
    out(i) = 0;
    MSpolyR.mon = [MSpolyR.mon NMON];    
    MSpolyR.coef{MSpolyR.nmon}.coef_val = coef_val;    
    MSpolyR.nmon = MSpolyR.nmon + 1;
end
MSpolyR.nmon = MSpolyR.nmon - 1;


%________________________________________________________________________________________________________
function MSPolyR = reorder_monomials(MSPolyA)

[f,or] = sort(MSPolyA.var);
MSPolyR.var = f;
MSPolyR.mon = [];
MSPolyR.N = [];
for i=1:size(or,2)
    N = MSPolyA.N(1,or(i));
    mon_a = MSPolyA.mon((or(i)-1)*N+1:or(i)*N,:);
    MSPolyR.mon = [MSPolyR.mon;mon_a];
    MSPolyR.N = [MSPolyR.N N];
end
MSPolyR.nxm = MSPolyA.nxm;
MSPolyR.coef = MSPolyA.coef;
MSPolyR.m = size(MSPolyR.var,2);
MSPolyR.nmon = MSPolyA.nmon;
MSPolyR.deg = MSPolyA.deg;
MSPolyR.name = MSPolyA.name;