```function [sol,x,B,E]=verlinsys(A,b)
%    VERLINSYS       Verified description of all solutions of a system of linear equations.
%
%    This is an INTLAB file. It requires to have INTLAB installed under
%    MATLAB to function properly.
%
%    For an m-by-n real matrix A and a matching real vector b,
%        [sol,x,B,E]=verlinsys(A,b)
%    yields a verified description of ALL solutions of the system of linear
%    equations A*x=b. The number of solutions is given in the variable sol:
%
%    sol = Inf     A*x=b is verified to have infinitely many solutions;
%                  x is an n-by-1 interval vector and B is an n-by-n
%                  interval matrix that are verified to contain a real
%                  vector xo and a real matrix Bo such that the set X of
%                  all solutions of A*x=b is described by
%                      X = { xo+Bo*y | y in R^n };
%                  moreover, x has at least n-r zero entries, where r is
%                  the rank of A ("basic solution"),
%    sol =  1      A*x=b is verified to have exactly one solution xo which
%                  is verified to be contained in the interval vector x;
%                  B is a matrix of interval zeros,
%    sol =  0      A*x=b is verified to possess no solution; x and B
%                  consist of NaN's,
%    sol = -1      no verified result; x and B consist of NaN's.
%
%    The structure E explains reasons for NaN output.
%
%    If sol==0, you may try VERLSQ for verified least squares solution(s).
%

%
%    Based on the description of the set X of all solutions (if nonempty) as
%        X = { xo+(eye(n,n)-pinv(A)*A)*y | y in R^n },
%    where xo in X, with a verified pseudoinverse being computed by VERPINV.
%
%    WARRANTY
%
%    Because the program is licensed free of charge, there is
%    no warranty for the program, to the extent permitted by applicable
%    law. Except when otherwise stated in writing the copyright holder
%    and/or other parties provide the program "as is" without warranty
%    of any kind, either expressed or implied, including, but not
%    limited to, the implied warranties of merchantability and fitness
%    for a particular purpose. The entire risk as to the quality and
%    performance of the program is with you. Should the program prove
%    defective, you assume the cost of all necessary servicing, repair
%    or correction.
%
%    History
%
%    2008-04-14   first version
%    2008-04-20   version for posting
%
gr=getround;
setround(0);
b=b(:);
[m,n]=size(A);
% defaults
sol=-1;
x=repmat(infsup(NaN,NaN),n,1);
B=repmat(infsup(NaN,NaN),n,n);
E.error='verlinsys: none';
E.where='NaN';
E.value='NaN';
% data check
if nargin~=2||nargout>4||m~=length(b)||~isreal(A)||~isreal(b)||isintval(A)||isintval(b)
E.error='verlinsys: wrong data';
setround(gr); return
end
% square case first
if m==n
xx=verifylss(A,b);
if ~isnan(xx.inf(1)) % computed
sol=1;
x=xx;
B=repmat(infsup(0,0),n,n); % B set to zero interval matrix
setround(gr); return
end
% not computed
end
% quantities I (for sol==0 and sol==1)
X=verpinv(A); % pseudoinverse
if isnan(X.inf(1,1))
E.error='verlinsys: verified pseudoinverse of A could not be computed';
setround(gr); return % X needed in all three conditions
end
res=(eye(m,m)-A*X)*b; % should be ~=0 for nonexistence of solution
BB=eye(n,n)-X*A; % should be nonzero for infinitely many solutions
fcrA=verfullcolrank(A); % should be fcrA==1, fcrAb==0 for existence of solution
fcrAb=verfullcolrank([A b]);
% no solution
if ~isnan(X.inf(1,1))&&(any(res.sup<0)||any(res.inf>0)) % verified no solution
sol=0;
setround(gr); return
end
% unique solution
if ~isnan(X.inf(1,1))&&fcrA==1&&fcrAb==0 % verified unique solution
sol=1;
x=X*b;
B=repmat(infsup(0,0),n,n); % B set to zero interval matrix
setround(gr); return
end
% quantities II (for sol==Inf)
[BS,K]=verbasis(A); % basis
if isnan(BS(1,1))
E.error='verlinsys: verified pseudoinverse of the basis could not be computed';
setround(gr); return
end
fcrBSb=verfullcolrank([BS b]); % should be fcrBSb==0 for existence of solution
XBS=verpinv(BS); % used for computation of x
% infinitely many solutions
if ~isnan(X.inf(1,1))&&~isnan(BS(1,1))&&~isnan(XBS.inf(1,1))&&fcrBSb==0&&(any(any(BB.inf>0))||any(any(BB.sup<0)))
sol=Inf;
xK=XBS*b; % xK verified to enclose a solution of BS*X=b
x=infsup(zeros(n,1),zeros(n,1)); x(K)=xK; % x verified to enclose a solution of A*X=b
B=BB;
setround(gr); return
end
E.error='verlinsys: quantities needed could not be computed; no verified result';
setround(gr);
```