function [F_L,Puy,G,F,E]=poly2gpcpm(A,B,L)
% POLY2GPCPM
% Syntax: [F_L,Puy,G,F,E]=poly2gpcpm(A,B,L)
% Purpose: Given the polynomial matrices A and B. This function solves
% the Diophantine equations for the GPC prediction problem and form
% the matrices F_L and Puy in the prediction model
% y_f=F_L*Delta u_f + p_L                   (1)
% p_L=Puy*[past inputs;past outputs].
% On input:
% A  - polynomial of order, na. A=1  +a_1*z^(-1)+...+a_na*z^(-na)
% B  - polynomial of order, nb. B=b_0+b_1*z^(-1)+...+b_nb*z^(-nb)
% This corresponds to a model: A y_k=B u_(k-1) + C e_k
% where C=1/(1-z^(-1)).
% On output:
% F_L, Puy - The matrices in the prediction model (1)
% G,F,E    - The polynomials satisfying the Diophantine equation.
%
% Author: DDIR, 25.03.2001

%
% Solving for the F_j, E_j and G_j polynomials. Only the coefficients in
% the F_j and G_j polynomials are used in the PM.
% NOTE: In the following we are using: z:=z^{-1} for simplicity of notation.
na=length(A)-1; 
nb=length(B)-1;
At=conv([1 -1],A);                 % At=(1-z)*A

%%% The F_j polynomials, j=1,2,...L where F_j=f_j0+f_j1*z+...f_jna*z^na
%F1=-At(2:na+2)                    % F1=z*(1-At)=f10+f11*z
F=zeros(4,na+1);
F(1,:)=-At(2:na+2);
for j=1:L-1
R=F(j,1);
for i=1:na+1
   if i<na+1
      F(j+1,i)=F(j,i+1)-R*At(i+1);
   else
      F(j+1,i)=-R*At(i+1);
   end     
end
end
%%% Solve for the E_j polynomials
E=zeros(L,L);
E(1,1)=1;
for j=2:L
   fj0=F(j-1,1);
   E(j,1:j)=[E(j-1,1:j-1) fj0];
end
%%% Solve for the G_j polynomials
G=zeros(L,L+1);
for j=1:L
   G(j,1:j+nb)=conv(E(j,1:j),B);
end
%%% Form the PM model matrices
F_L=zeros(L,L);                       % The F_L^delta matrix.
k=1;
for j=1:L
   F_L=F_L+diag(G(k:L,j),-j+1);
   k=k+1;
end
Puy=[diag(G,1)*[-1 1] fliplr(F)];     % The matrix taking care of the past.
