clear all
warning off
syms s
syms U V matrix

% data
A=[1 2
    1 0];
B=[2 2
    3 4];
C=[0 1
    1 0];
D=[1 2
    2 0];

A=[-1 -2
    0 1];
B=[1 2
    2 1];
C=[1 1
    1 1];
D=[0 0
    0 0];

A=[-1 0
    0 1];
B=[1
    0];
C=[1 0];
D=[1];

G=factor(C*inv(s*eye(size(A))-A)*B+D);

G=[2*(s+1)*(s+2)/s/(s+3)/(s+4) (s+2)/(s+1)/(s+3)];

G=[1/(s+1) (s+3)/(s+1)/(s-2)
    10/(s-2) 5/(s+3)];

G=[1/(s+1)/(s+2) (2*s+1)/(s+1)/(s+2) s/(s+1)/(s+2)
    1/(s+1)^2 (s^2+5*s+3)/(s+1)^2 s/(s+1)^2
    1/(s+1)^2/(s+2) (2*s+1)/(s+1)^2/(s+2) s/(s+1)^2/(s+2)];

G=[(s+1)/(s+2)
    (s+2)/(s+1)];

G=[s/(s+1)^2/(s+2)^2 s/(s+2)^2
    -s/(s+2)^2 -s/(s+2)^2];

G=[1/(s^2+3*s+2) -1/(s^2+3*s+2)
    (s^2+s-4)/(s^2+3*s+2) (2*s^2-s-8)/(s^2+3*s+2)
    (s-2)/(s+1) (2*s-4)/(s+1)];

% construct a left fraction and a right fraction
[p,q]=size(G);
if (p==1 & q==1)
    r=1;
else
    r=eval(rank(G));
end
[nG,dG]=numden(G);
d=1;
for k=1:p
    for l=1:q
        d=maple('lcm',d,dG(k,l));
    end
end
nGr=factor(G*d);
dGr=eye(q,q)*d;
nGl=factor(d*G);
dGl=d*eye(p,p);

% a left coprime fraction, L=gcld(dGl,nGl)
dnGl=[dGl nGl].';
L=factor(maple('hermite',dnGl,s));
L=L(1:p,1:p).';
nGls=factor(inv(L)*nGl);
dGls=factor(inv(L)*dGl);
det(dGls);

% a right coprime fraction, R=gcrd(dGr,nGr)
dnGr=[dGr
    nGr];
R=factor(maple('hermite',dnGr,s));
R=R(1:q,1:q);
nGrs=factor(nGr*inv(R));
dGrs=factor(dGr*inv(R));
det(dGrs);

% McMillan Form
if (p==1 & q==1)
    M=G;
else
    M=factor(maple('smith',simple(d*G),s,U,V)/d);
end
tz=[];
tp=[];
for k=1:r
    tz=[tz solve(M(k,k)).'];
    tp=[tp solve(1/M(k,k)).'];
end
U=maple('eval',U);
V=maple('eval',V);