function[X,degX,Y,degY]=dioin(A,degA,B,degB,C,degC,ALPHA,SJ,r)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [X,degX,Y,degY]=dioin(A,degA,B,degB,C,degC,ALPHA,SJ,r) %
%                                                        %
% MATLAB function solves the linear Diophantine          %
% equations for polynomial matrices of the form          % 
%       X(s)A(s) + Y(s)B(s) = C(s)                       %
% where  the polynomial matrices A(s), B(s) and C(s)     %
% are given and X(s) and Y(s) are to be found.           %
% A method based on the polynomial matrix interpolation  %
% theory.                                                %
% ALPHA = [alpha1,...,alphaL] is a matrix, where the     %
% columns alpha(j) are nonzero.                          %
% SJ is a vector containing the interpolation points  sj %
% (distinct)  SJ = [s1,s2,...,sL]                        %
% The solution degree is r.                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   used function: polsize, alpha, beta, gama            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   S. Pejchova, May 11th 1994             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %   dimensions check             %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[rA,cA]=polsize(A,degA);
[rB,cB]=polsize(B,degB);
[rC,cC]=polsize(C,degC);

if (rA~=cA) | (cA~=cB) | (cA~=cC)
   error('dioel: Inconsistent dimensions of input matrices');
end

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Creating matrix S = [ A(s) |  %
          %                      | B(s) ]  %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[rd,cd]=size(A);
[rn,cn]=size(B);
if degA > degB
          NUL=zeros(rB,(cd-cn));
          B=[B,NUL];
          degB=degA;
   else
       if degA < degB
          NUL=zeros(rA,(cn-cd));
          A=[A,NUL];
          degA=degB;
       end
end

S=[A;B];
degS=degA;

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Find matrix AL=[S(s1)alpha1,...,S(sL)alphaL] %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

AL=alpha(S,degS,ALPHA,SJ);

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Find matrix BL=[Q(s1)alpha1,...,Q(sL)alphaL] %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

BL=alpha(C,degC,ALPHA,SJ);

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Dimensions of X and Y                        %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rX=rC;
cX=rA;
rY=rC;
cY=rB;
dY=[rY,cY];
Xr=eye(rX);

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Creating matrix BT = [beta(1),...,beta(L)]   %
          %  where beta(j) = bl(j) - [sj^r.Xr,0].al(j)    %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

BT=beta(AL,BL,Xr,dY,SJ,r);

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Matrix SDL calculation                       %
          %  SDL=[SD(s1)al1,...,SD(sL)alL]                %
          %  where SD(s)=[I,...,I.s^r]'                   %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SDL=gama(AL,SJ,r);
[rSDL,cSDL]=size(SDL);
if cSDL>rSDL
        error('dioin: Error in dimensions of SDL');
end

for i=1:cSDL
    if SDL(:,i)==zeros(rSDL,1)
       error('dioin: Matrix SDL has not full rank');
    end
end

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %  Calculate the solution  M = BT/(SDL)         %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=BT/(SDL);

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %   Creating the matrices  X  and  Y            %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degX=r;
degY=r;
      X=M(:,1:cX);
      Y=M(:,(cX+1):(cX+cY));

  if r~=0
        for i=1:r
          Y=[Y,M(:,(i*(cX+cY)+cX+1):((i+1)*(cX+cY)))];
          POM=M(:,(i*(cX+cY)+1):(i*(cX+cY)+cX));
          if i==r

             POM=POM+eye(rX,cX);
          end         
          X=[X,POM];
        end
  end