// Updated
// Implements the method of Ahlen and Sternad

function [R,degR,S,degS,X,degY] = ...
lqg_as(A,degA,B,degB,C,degC,rho,k,V,degV,W,degW,F,degF)

getf specfac.sci;
[r,b,degb] = specfac(A,degA,B,degB,rho,V,degV,W,degW,F,degF);

WFA = flip(convol(A,convol(F,W)));
dWFA = degW + degF + degA;

[rhs1,drhs1] = polmul(W,degW,WFA,dWFA);
[rhs1,drhs1] = polmul(rhs1,drhs1,C,degC);
rhs1 = rho * rhs1;

rhs2 = convol(C,convol(V,flip(convol(B,V))));
drhs2 = degC + 2*degV + degB;
for i = 1:degb-degB-degV,
  rhs2 = convol(rhs2,[0,1]);
end
drhs2 = drhs2 + degb-degB-degV;
C1 = zeros(1,2);

[C1,degC1] = putin(C1,0,rhs1,drhs1,1,1);
[C1,degC1] = putin(C1,degC1,rhs2,drhs2,1,2);

rbf = r * flip(b);
D1 = zeros(2,2);
[D1,degD1] = putin(D1,0,rbf,degb,1,1);
for i = 1:k,
	  rbf = convol(rbf,[0 1]);
end
[D1,degD1] = putin(D1,degD1,rbf,degb+k,2,2);

N = zeros(1,2);
[N,degN] = putin(N,0,-B,degB,1,1);
[AF,dAF] = polmul(A,degA,F,degF);
[N,degN] = putin(N,degN,AF,dAF,1,2);

[Y,degY,X,degX] = ...
xdync(N,degN,D1,degD1,C1,degC1);
[R,degR] = ext(X,degX,1,1);
[S,degS] = ext(X,degX,1,2);
X = flip(Y);
X = clean(X);
endfunction;