% programed by Shin'ichi Oishi, Waseda University (c)
% opened 2000/10/23
function [error2,x]=linear(A,b)
   R=inv(A);
   x=R*b;
   [N,M]=size(A);
   setround(-Inf);
   Gd=abs(R*A-eye(N,N));
   rd=A*x-b;
   setround(Inf);
   Gu=abs(R*A-eye(N,N));
   ru=A*x-b;
   center=(ru+rd)/2;
   radius=center-rd;
   ru=max(rd,ru);
   Gu=max(Gd,Gu);
   R_norm=norm(R,inf);
   G_norm=norm(Gu,inf);
   Rru=abs(R*center+abs(R)*radius);
   setround(-Inf);
   Rrd=abs(R*center-abs(R)*radius);
   D=1-G_norm;
   if D>0
      setround(Inf);
      Rru=max(Rrd,Rru);
      Rr_norm=norm(Rru,inf);
      err=Rr_norm/D;
      kappa=Gu*ones(N,1);
      error2=Rru+(err*kappa);
      nerror=max(error2);
   else
	  fprintf('%d\n',D)
   end