function [Kinv,Kinvd]= Kinvers(Rn,epsn,Nper1,Nhom1,Nhet12,Nper2,Nhom2,Nhet21,...
	iiper1,iihom1,iihet12,iiper2,iihom2,iihet21)
%
% Upper bounds for ||K^(-1)|| and ||K^(-1)d|| are computed. See Lemma 3 
% and Lemma 4 of the paper 
% 
% input:
%	Rn(N,4)		sequence of triangular 2x2-matrices Rn(i)
%	epsn(N)		magnitude of d_n: |d_n| <= epsn(n)
%	Nper1,Nhom1,		
%	Nhet12,Nper2,		
%	Nhom2,Nhet21	length of the orbits u^1, u^2, u^3, u^4, u^5, u^6 
%	iiper1,iihom1,		
%	iihet12,iiper2,		
%	iihom2,iihet21	indices of the beginning of the  
%			orbits u^1, u^2, u^3, u^4, u^5, u^6
%
% output: 
% 	Kinv: 		interval variable. Upper bound for ||K^(-1)||
% 	Kinvd: 		interval variable. Upper bound for ||K^(-1)d||
%
Rn=abs(Rn);
setround(1)

eta1max=0;
eta2max=0;
qq1max=0;
qq2max=0;

istposeta=1;
istposqq=1;
while ~isempty(istposeta) | ~isempty(istposqq),
%  anzposeta=length(istposeta)
%  anzposq=length(istposqq)
  eta(1)=eta1max;
  qq(1)=qq1max;
  for i=iiper1:iiper1+Nper1-1,
    eta(i+1)=1+Rn(i,4)*eta(i);
    qq(i+1)=epsn(i)+Rn(i,4)*qq(i);
  end
  etaper1=eta(iiper1+Nper1);
  qqper1=qq(iiper1+Nper1);

  eta(iihom1)=eta1max;
  qq(iihom1)=qq1max;
  for i=iihom1:iihom1+Nhom1-1,
    eta(i+1)=1+Rn(i,4)*eta(i);
    qq(i+1)=epsn(i)+Rn(i,4)*qq(i);
  end
  etahom1=eta(iihom1+Nhom1);
  qqhom1=qq(iihom1+Nhom1);

  eta(iihet12)=eta1max;
  qq(iihet12)=qq1max;
  for i=iihet12:iihet12+Nhet12-1,
    eta(i+1)=1+Rn(i,4)*eta(i);
    qq(i+1)=epsn(i)+Rn(i,4)*qq(i);
  end
  etahet12=eta(iihet12+Nhet12);
  qqhet12=qq(iihet12+Nhet12);

  eta(iiper2)=eta2max;
  qq(iiper2)=qq2max;
  for i=iiper2:iiper2+Nper2-1,
    eta(i+1)=1+Rn(i,4)*eta(i);
    qq(i+1)=epsn(i)+Rn(i,4)*qq(i);
  end
  etaper2=eta(iiper2+Nper2);
  qqper2=qq(iiper2+Nper2);

  eta(iihom2)=eta2max;
  qq(iihom2)=qq2max;
  for i=iihom2:iihom2+Nhom2-1,
    eta(i+1)=1+Rn(i,4)*eta(i);
    qq(i+1)=epsn(i)+Rn(i,4)*qq(i);
  end
  etahom2=eta(iihom2+Nhom2);
  qqhom2=qq(iihom2+Nhom2);

  eta(iihet21)=eta2max;
  qq(iihet21)=qq2max;
  for i=iihet21:iihet21+Nhet21-1,
    eta(i+1)=1+Rn(i,4)*eta(i);
    qq(i+1)=epsn(i)+Rn(i,4)*qq(i);
  end
  etahet21=eta(iihet21+Nhet21);
  qqhet21=qq(iihet21+Nhet21);


  istposeta=find([eta1max eta1max eta1max eta2max eta2max eta2max] - [etaper1 etahom1 etahet21 etaper2 etahom2 etahet12] < 0);

  istposqq=find([qq1max qq1max qq1max qq2max qq2max qq2max] - [qqper1 qqhom1 qqhet21 qqper2 qqhom2 qqhet12] < 0);

  eta1max=max([etaper1 etahom1 etahet21]);
  eta2max=max([etaper2 etahom2 etahet12]);
  qq1max=max([qqper1 qqhom1 qqhet21]);
  qq2max=max([qqper2 qqhom2 qqhet12]);
end



xi1max=0;
xi2max=0;
pp1max=0;
pp2max=0;

istposxi=1;
istpospp=1;

while ~isempty(istposxi) | ~isempty(istpospp),
%  anzposxi=length(istposxi)
%  anzposp=length(istpospp)
  xi(iihet21+Nhet21)=xi1max;
  pp(iihet21+Nhet21)=pp1max;
  for i=iihet21+Nhet21-1:-1:iihet21,
    xi(i)=(1+Rn(i,3)*eta(i)+xi(i+1))/Rn(i,1);
    pp(i)=(epsn(i)+Rn(i,3)*qq(i)+pp(i+1))/Rn(i,1);
  end
  xihet21=xi(iihet21);
  pphet21=pp(iihet21);

  xi(iihom2+Nhom2)=xi2max;
  pp(iihom2+Nhom2)=pp2max;
  for i=iihom2+Nhom2-1:-1:iihom2,
    xi(i)=(1+Rn(i,3)*eta(i)+xi(i+1))/Rn(i,1);
    pp(i)=(epsn(i)+Rn(i,3)*qq(i)+pp(i+1))/Rn(i,1);
  end
  xihom2=xi(iihom2);
  pphom2=pp(iihom2);

  xi(iiper2+Nper2)=xi2max;
  pp(iiper2+Nper2)=pp2max;
  for i=iiper2+Nper2-1:-1:iiper2,
    xi(i)=(1+Rn(i,3)*eta(i)+xi(i+1))/Rn(i,1);
    pp(i)=(epsn(i)+Rn(i,3)*qq(i)+pp(i+1))/Rn(i,1);
  end
  xiper2=xi(iiper2);
  ppper2=pp(iiper2);

  xi(iihet12+Nhet12)=xi2max;
  pp(iihet12+Nhet12)=pp2max;
  for i=iihet12+Nhet12-1:-1:iihet12,
    xi(i)=(1+Rn(i,3)*eta(i)+xi(i+1))/Rn(i,1);
    pp(i)=(epsn(i)+Rn(i,3)*qq(i)+pp(i+1))/Rn(i,1);
  end
  xihet12=xi(iihet12);
  pphet12=pp(iihet12);

  xi(iihom1+Nhom1)=xi1max;
  pp(iihom1+Nhom1)=pp1max;
  for i=iihom1+Nhom1-1:-1:iihom1,
    xi(i)=(1+Rn(i,3)*eta(i)+xi(i+1))/Rn(i,1);
    pp(i)=(epsn(i)+Rn(i,3)*qq(i)+pp(i+1))/Rn(i,1);
  end
  xihom1=xi(iihom1);
  pphom1=pp(iihom1);

  xi(iiper1+Nper1)=xi1max;
  pp(iiper1+Nper1)=pp1max;
  for i=iiper1+Nper1-1:-1:iiper1,
    xi(i)=(1+Rn(i,3)*eta(i)+xi(i+1))/Rn(i,1);
    pp(i)=(epsn(i)+Rn(i,3)*qq(i)+pp(i+1))/Rn(i,1);
  end
  xiper1=xi(iiper1);
  ppper1=pp(iiper1);




  istposxi=find([xi1max xi1max xi1max xi2max xi2max xi2max] - [xiper1 xihom1 xihet21 xiper2 xihom2 xihet12] < 0);

  istpospp=find([pp1max pp1max pp1max pp2max pp2max pp2max] - [ppper1 pphom1 pphet21 ppper2 pphom2 pphet12] < 0);

  xi1max=max([xiper1 xihom1 xihet21]);
  xi2max=max([xiper2 xihom2 xihet12]);
  pp1max=max([ppper1 pphom1 pphet21]);
  pp2max=max([ppper2 pphom2 pphet12]);
end

xi(1)=xi1max;
pp(1)=pp1max;

Kinv=sqrt(max(eta.^2+xi.^2));
Kinvd=sqrt(max(qq.^2+pp.^2));

setround(0)

return

%
% written   Sept/10/02   D.Stoffer
%


