function [xl,xu,S]=intlineqhull(Ac,Delta,bc,delta) % INTerval LINear EQuations, HULL % TO BE POSTED
%
%  A=[Ac-Delta,Ac+Delta] is a square interval matrix, b=[bc-delta,bc+delta] is a matching interval vector,
% ~isempty(xl): [xl,xu] is the interval hull of the solution set of A*x=b, and S=[],
% ~isempty(S) : S is a singular matrix in A, and xl=[], xu=[].
%
%  MIT License for INTLINEQHULL:
%
%  Copyright 2005-2017 Jiri Rohn (rohn@cs.cas.cz). 
%
%  Permission is hereby granted, free of charge, to any person obtaining a copy
%  of this software and associated documentation files (the "Software"), to deal
%  in the Software without restriction, including without limitation the rights
%  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%  copies of the Software, and to permit persons to whom the Software is
%  furnished to do so, subject to the following conditions:
% 
%  The above copyright notice and this permission notice shall be included in all
%  copies or substantial portions of the Software.
% 
%  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
%  SOFTWARE.
%
n=length(bc); 
ep=max(n,10)*max([norm(Ac-Delta,1) norm(Ac+Delta,1) norm(bc-delta,1) norm(bc+delta,1)])*eps;
if rank(Ac)<n, xl=[]; xu=[]; S=Ac; return, end
xl=Ac\bc; xu=xl;
z=ones(1,n);
for j=1:n, if xl(j)<0, z(j)=-1; end, end
Z=z; D=[]; 
while ~isempty(Z)
    p=size(Z,1); z=Z(p,:); Z=Z(1:p-1,:); D=[D; z];
    [Q,S]=qzmatrix(Ac,Delta,-z);
    if ~isempty(S), xl=[]; xu=[]; return, end
    xlt=Q*bc-abs(Q)*delta;
    [Q,S]=qzmatrix(Ac,Delta,z);
    if ~isempty(S), xl=[]; xu=[]; return, end
    xut=Q*bc+abs(Q)*delta;
    if all(xlt<=xut)
        xl=min(xl,xlt);
        xu=max(xu,xut);
        for j=1:n
            zt=z; zt(j)=-zt(j);
            if (xl(j)*xu(j)<=ep)&&(~ismember(zt,Z,'rows'))&&(~ismember(zt,D,'rows'))
                Z=[Z; zt];
            end
        end
    end
end
S=[];
%
function [Q,S]=qzmatrix(Ac,Delta,z) % $Q_z$ MATRIX
%
%  A=[Ac-Delta,Ac+Delta]: square interval matrix, z a plus/minus-one vector,
% ~isempty(Q): Q solves Q*Ac-abs(Q)*Delta*diag(z)=eye(size(A)), and S=[],
% ~isempty(S): S is a singular matrix in A, and Q=[].
%
n=size(Ac,1); I=eye(n,n); 
Q=zeros(n,n); S=[]; 
for i=1:n
    [x,S]=absvaleqn(Ac',-diag(z)*Delta',I(:,i));
    if ~isempty(S), Q=[]; S=S'; return, end
    Q(i,:)=x';
end
%
function [x,S,iter]=absvaleqn(A,B,b) % ABSolute VALue EQuatioN 
%
% ~isempty(x): x solves A*x+B*abs(x)=b (A,B square), S is empty,
% ~isempty(S): S is a singular matrix satisfying abs(S-A)<=abs(B), x is empty.
%  iter: number of iterations (it may be zero).
%
%  Copyright 2005-2015 Jiri Rohn (rohn@cs.cas.cz)
%
b=b(:); n=length(b); I=eye(n,n);
ep=n*(max([norm(A,inf) norm(B,inf) norm(b,inf)]))*eps;
x=[]; S=[]; iter=0;
if rank(A)<n, S=A; return, end
x=A\b; z=sgn(x);
if rank(A+B*diag(z))<n, S=A+B*diag(z); x=[]; return, end
x=(A+B*diag(z))\b;
C=-inv(A+B*diag(z))*B;
X=zeros(n,n);
r=zeros(1,n);
while any(z.*x<-ep)
    k=find(z.*x<-ep,1);
    iter=iter+1;
    if 1+2*z(k)*C(k,k)<=0
        S=A+B*(diag(z)+(1/C(k,k))*I(:,k)*I(k,:)); x=[]; return
    end
    if ((k<n)&&(all(r(k)>r(k+1:n))))||((k==n)&&(r(k)>0))
        x=x-X(:,k); z=sgn(x);
        ct=A*x;
        jm=abs(B)*abs(x);
        y=zeros(1,n);
        for i=1:n
            if jm(i)>ep, y(i)=ct(i)/jm(i); else y(i)=1; end
        end
        S=A-diag(y)*abs(B)*diag(z); x=[]; return
    end
    X(:,k)=x;
    r(k)=iter;
    z(k)=-z(k);
    alpha=2*z(k)/(1-2*z(k)*C(k,k));
    x=x+alpha*x(k)*C(:,k);
    C=C+alpha*C(:,k)*C(k,:);
end 
%
function z=sgn(x) % SiGN vector of x (column)
n=length(x); 
z=zeros(n,1);
for j=1:n
    if x(j)>=0, z(j)=1; else z(j)=-1; end 
end