```function [oe,Ao,bo]=veroettprag(A,b,x)
%    VEROETTPRAG    Verified solution of the Oettli-Prager inequality, and the data to it.
%
%    For a rectangular interval matrix A and matching interval vector b and
%    real vector x,
%    [oe,Ao,bo]=veroettprag(A,b,x)
%    either verifies that x solves the Oettli-Prager inequality
%    or verifies that it is not so, or yields no verified result.
%    In the positive case it also finds data Ao, bo to x.
%
%    Possible outputs:
%
%    oe= 1       x is verified to satisfy (1);
%                Ao is a very tight interval matrix which is a part of A and
%                bo is a very tight interval vector which is a part of b
%                such that x is verified to satisfy A1*x=b1 (in exact arithmetic)
%                for some A1 in Ao, b1 in bo,
%    oe= 0       x is verified not to satisfy (1); Ao, bo consist of NaN's,
%    oe=-1       no verified result.
%
%    COMMENT. The Oettli-Prager theorem (Num. Math. 1964) says that a
%    vector x satisfies the inequality (1) if and only if A1*x=b1 holds for
%    some A1 in A, b1 in b. Hence, the main task here is not only to verify
%    (1), but also to find A1, b1. They are enclosed in the computed Ao, bo.
%

%
%    Based on Theorem 2.9 and Proposition 2.10 in
%    M. Fiedler, J. Nedoma, J. Ramik, J. Rohn and K. Zimmermann, Linear
%    Optimization Problems with Inexact Data, Springer-Verlag, New York
%    2006.
%
%    This work was supported by the Czech Republic National Research
%    Program "Information Society", project 1ET400300415.
%
%    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.
%
gr=getround;
setround(0);
[m,n]=size(A);
b=b(:);
oe=-1;
Ao=repmat(infsup(NaN,NaN),m,n);
bo=repmat(infsup(NaN,NaN),m,1);
if nargin<3||m~=length(b)||n~=length(x)||~isreal(A)||~isreal(b)
setround(gr); return
end
if ~isintval(A) % allows for real input
A=infsup(A,A);
end
if ~isintval(b)
b=infsup(b,b);
end
z=ones(n,1); % preallocation
for i=1:n
if x(i)>=0, z(i)=1; else z(i)=-1; end % sign vector of x
end
xi=infsup(x,x);
oeprl=abs(Ac*xi-bc);               % Oettli-Prager inequality, left  side
oeprr=Delta*abs(xi)+delta;         % Oettli-Prager inequality, right side
if all(oeprl.sup<=oeprr.inf)       % Oettli-Prager inequality satisfied
y=(Ac*xi-bc)./oeprr;           % componentvise division
y(find(isnan(y)))=infsup(1,1); % case od both numerator and denominator being zero
Ao=Ac-(diag(y)*Delta)*diag(z); % construction of Ao
bo=bc+diag(y)*delta;           % construction of bo
Ao=intersect(Ao,A);            % Ao made part of A
bo=intersect(bo,b);            % bo made part of b
if ~any(any(isnan(Ao)))&&~any(any(isnan(bo))) % intersections nowhere empty
oe=1;                      % x verified to satisfy the inequality, Ao, bo found
setround(gr); return
else
oe=-1;                     % x verified to satisfy the inequality, Ao, bo not found
Ao=repmat(infsup(NaN,NaN),m,n);
bo=repmat(infsup(NaN,NaN),m,1);
setround(gr); return
end
end
if any(oeprr.sup<oeprl.inf)
oe=0;                          % x verified not to satisfy the inequality
setround(gr); return
end
setround(gr);                      % no verified result
%
% Subfunction
%