function [U,degU,V,degV] = rmod(A,degA,B,degB,mes)

% polynomial matrix right modulo division 
%                 B = U*A + V                               
%                                                         
%  A,degA ........ nonsingular column reduced polynomial matrix              
%  B,degB ........ polynomial matrix              
%  mes (=0,1,2)... message and interrupt level    
%  U, degU     ... polynomial matrix   
%  V, degV     ... polynomial matrix: V*A^{-1} strictly proper     
%
%  call: [U,degU,V,degV] = rmod(A,degA,B,degB,mes)      

% functions used: polsize, polcl, colred, transp; 

% disp('.............Start of right division..............');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting of EPS                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

EPS=1e-8;

[rA,cA]=polsize(A,degA);
if rA ~= cA
  error('rmod: A not square')
end
DegA=coldeg(A,degA);

[rB,cB]=polsize(B,degB);
if cB ~= rA
  error('rmod: A,B not compatible')
end
DegB=coldeg(B,degB);

U=zeros(rB,cB);
degU=0;

if DegB(:)<DegA(:)
   V=B;
   degV=degB;
   if mes>0, U, degU, V, degV, keyboard, end
   return;
end   

% B*A^{-1} not strictly proper, modulo division follows

DegE=degA-DegA;       % powers to equalize column degrees of A

DegBE=DegB+DegE;
degBE=max(DegBE);

degU=degBE-degA;
for i=1:degU
   U=[U zeros(rB,cB)];
end

BE=B;
for j=1:degBE-degB
   BE=[BE zeros(rB,cB)];
end

Al=A;
for j=1:cA
   Ahc(:,j)=A(:,j+cA*DegA(j));
   Al(:,j+cA*DegA(j))=zeros(rA,1);
    
   if DegE(j)~=0
      for k=cB*degBE+j:-cB:cB*DegE(j)+j
         BE(:,k)=BE(:,k-cB*DegE(j));
      end
      for k=j:cB:cB*(DegE(j)-1)+j
         BE(:,k)=zeros(rB,1);
      end 
   end 
end

if Al(:,cA*degA+1:cA*(degA+1))==zeros(cA,cA)
   Al(:,cA*degA+1:cA*(degA+1))=[];
   degAl=degA-1;
   DegAl=DegA-1;
else
   degAl=degA;
   DegAl=DegA;
   error('modulo: failed to pick up the lower-terms');
   disp(' matrix of A');
   break;
end

Ahc=inv(Ahc);

if mes>1, disp('end of preliminaries'), keyboard, end

% Main loop                                             

for ner=1:degU+1

   % split BE                                    

   BEh=BE(:,cB*degBE+1:cB*(degBE+1));
   BEl=BE(:,1:cB*degBE);
   degBEl=degBE-1;

   % update U                                         

   BEh=BEh*Ahc;              % multiplication without EPS
   U(:,cB*(degBE-degA)+1:cB*(degBE-degA+1))=BEh;

   % Computation of new BE                              

   BEh=BEh*Al;
   degBEh=degAl;                    % multiplication without EPS
   for j=1:degBE-degAl-1            % expanding dimension 
      BEh=[BEh zeros(rB,cA)];
   end   
   degBEh=degBE-1;                            

   DegS=degBE-DegA;
   for j=1:cA
      if DegS(j)~=0
         for k=cA*degBEh+j:-cA:cA*DegS(j)+j
            BEh(:,k)=BEh(:,k-cA*DegS(j));
         end
         for k=j:cA:cA*(DegS(j)-1)+j
            BEh(:,k)=zeros(rB,1);
         end
      end
   end  
 
   [BE,degBE]=poladd(BEl,degBEl,-BEh,degBEh);          

if mes>1, disp('another step over'), keyboard, end

   % Earlier termination, if possible                   

   if degBE<degA
      break;
   end
end

% recover V 

for j=1:cA
   if DegE(j)~=0
      for k=j:cB:cB*(DegE(j)-1)+j
         if BE(:,k)~=0
            V=[];
            degV=[];
            error('rmod: failed to recover V');
         end
      end
      for k=j:cB:cB*(degBE-DegE(j))+j
         BE(:,k)=BE(:,k+cB*DegE(j));
      end
      for k=cB*(degBE-DegE(j)+1)+j:cB:cB*(degBE+1)
         BE(:,k)=zeros(rB,1);
      end
   end
end

[V,degV]=polcl(BE,degBE);

if mes>0
  disp('..............End of right division...............');
  disp('U,degU,V,degV available')
  keyboard 
end