function LmzC= AcommB(LmzA,LmzB);
%AcommB  Bivariate polynomial for the commutator of two random matrices
%   LmzC = AcommB(LmzA,LmzB) returns the bivariate polynomial for 
%   C = i*[A, Q*B*Q'] where Q is Haar unitary and [X,Y] denotes X*Y-Y*X
%
%   Input : LmzA, LmzB
%   Output: LmzC
%   
%   Notation: Lmz is the bivariate polynomial satisfied by the Stieltjes
%   transform (see definiton in reference below)
%   
%   Example 1:
%       syms m z
%       LmzA = wignerpol;
%       LmzB = invA(LmzA);
%       LmzC = AcommB(LmzA,LmzB);
%       (should return the cauchy distribution)
%
%   Example 2:
%       syms m z 
%       LmzA = wignerpol;
%       LmzB = wishartpol(c);
%       LmzC = AcommB(LmzA,LmzB);
%       pdfinfo = Lmz2pdf(LmzC);
%       plot(pdfinfo.range,pdfinfo.density,'red','LineWidth',2);
%
%
%
%   Reference: N. Raj Rao and Alan Edelman, "The polynomial method for
%   the commutator of random matrices".
%   
%   N. Raj Rao, Wednesday, March 15, 2006.

syms m z r g s y 

warning off
LrgA = Lmz2Lrg(LmzA);
LrgAe = Luv2Lsq(LrgA,r,g);
LmzAe = Lrg2Lmz(LrgAe);
LsyAe = Lmz2Lsy(LmzAe);

if(LmzA==LmzB),
    LsyBe = LsyAe;
else
    LrgB = Lmz2Lrg(LmzB);
    LrgBe = Luv2Lsq(LrgB,r,g);
    LmzBe = Lrg2Lmz(LrgBe);
    LsyBe = Lmz2Lsy(LmzBe);
end

LsyT = L1timesL2(LsyAe,LsyBe,s);
LsyT = subs(LsyT,s,s*(1+y));
LmzT = Lsy2Lmz(LsyT);
LmzC = AplusB(LmzT,LmzT);
LrgC = Lmz2Lrg(LmzC);
LrgC = subs(LrgC,g,g^2);
LrgC = subs(LrgC,r,r/g);
LmzC = Lrg2Lmz(LrgC);

warning on

LmzC = irreducLuv(LmzC,m,z);