function x = spanish(A,b)
%
%   Computing  formal solutions  to  interval  linear systems
%   with "Spanish operator" in Kaucher complete interval arithmetic  
%   by subdifferential Newton method 
  
% Sergey P. Shary, 2009
  
IterLim = 100;  %   restriction on the number of iterations 
Eps = 1.e-7;    %   relative defect of the approximate solution 
tau = 1.;       %   damping factor (a real within (0,1], 1 is recommended)
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
%   inputting data, checking their correctness  
%
  A = intval(A);
  b = intval(b);
  
  m = size(A,1);
  n = size(A,2);
  if ( m~=n ) 
      error('The input matrix must be square!')
  end
  
  mb = size(b,1);
  nb = size(b,2);
  if ( mb~=m || nb>1 ) 
      error('Right-hand side size does no match to that of the matrix!')
  end 
    
  if( Eps<0 )
      error('Incorrect relative defect!') 
  end
  
  if( tau<=0 || tau>1 )
      error('Incorrect damping factor!') 
  end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   forming the "middle" point matrix F of the interval system under solution
% 
   
Dim = n;
F = zeros(2*Dim);
xx = zeros(2*Dim,1);
  
for i=1:Dim
    for j=1:Dim  
        p = mid(A(i,j));
        if( p>=0. )
            F(i,j)=p; F(i+Dim,j+Dim)=p;
        else
            F(i+Dim,j)=p; F(i,j+Dim)=p; 
        end
    end
end
  
for i=1:Dim
    %   dualizing the right-hand side vector
    d(2*i) = sup(b(i));  d(2*i+1) = inf(b(i));
    xx(i) = d(2*i); xx(i+Dim) = d(2*i+1);
end
  
%      solving the embedded "middle" point system
%                we find the starting approximation of the method
  
xx = F\xx;
  
%   launching the iteration of the subdifferential Newton method  
ni=0; r=Inf; q=1;
  
while ( r/q>Eps && ni<IterLim )
    r = 0;      
    ni = ni+1; 
    x = xx;
    F = 0;
    for i=1:Dim 
        s0 = 0;   s1 = 0;
        for j=1:Dim 
            g0 = inf(A(i,j)); g1 = sup(A(i,j));
  
            %   handling diagonal and off-diagonal entries separately  
            if( i~=j ) 
                h0=x(j);  h1=x(j+Dim); 
                %   determining the types of intervals under operation 
                if ( g0*g1>0 ) 
                    if( g0>0 ), l=0; else l=2; end 
                else
                    if( g0<=g1 ), l=1; else l=3; end
                end
                if ( h0*h1>0 ) 
                    if( h0>0 ), m=1; else m=3; end
                else
                    if( h0<=h1 ), m=2; else m=4; end
                end
  
                %   computing multiplication in Kaucher arithmetic,
                %                    forming the subgradient matrix F   
                switch ( 4*l+m ) 
                    case 1
                        t0=g0*h0; t1=g1*h1; F(i,j)=g0; F(i+Dim,j+Dim)=g1; 
                    case 2
                        t0=g1*h0; t1=g1*h1; F(i,j)=g1; F(i+Dim,j+Dim)=g1; 
                    case 3
                        t0=g1*h0; t1=g0*h1; F(i,j)=g1; F(i+Dim,j+Dim)=g0;
                    case 4
                        t0=g0*h0; t1=g0*h1; F(i,j)=g0; F(i+Dim,j+Dim)=g0;
                    case 5
                        t0=g0*h1; t1=g1*h1; F(i,j+Dim)=g0; F(i+Dim,j+Dim)=g1;
                    case 6
                        u0=g0*h1; v0=g1*h0;  u1=g0*h0; v1=g1*h1;
                        if u0<v0  
                            t0=u0; F(i,j+Dim)=g0; 
                        else
                            t0=v0; F(i,j)=g1;
                        end
                        if u1>v1 
                            t1=u1; F(i+Dim,j)=g0; 
                        else
                            t1=v1; F(i+Dim,j+Dim)=g1;
                        end
                    case 7
                        t0=g1*h0; t1=g0*h0; F(i,j)=g1; F(i+Dim,j)=g0;
                    case 8 
                        t0=0; t1=0;
                    case 9
                        t0=g0*h1; t1=g1*h0; F(i,j+Dim)=g0; F(i+Dim,j)=g1;
                    case 10
                        t0=g0*h1; t1=g0*h0; F(i,j+Dim)=g0; F(i+Dim,j)=g0;
                    case 11
                        t0=g1*h1; t1=g0*h0; F(i,j+Dim)=g1; F(i+Dim,j)=g0;
                    case 12
                        t0=g1*h1; t1=g1*h0; F(i,j+Dim)=g1; F(i+Dim,j)=g1;
                    case 13
                        t0=g0*h0; t1=g1*h0; F(i,j)=g0; F(i+Dim,j)=g1;
                    case 14
                        t0=0; t1=0;
                    case 15
                        t0=g1*h1; t1=g0*h1; F(i,j+Dim)=g1; F(i+Dim,j+Dim)=g0;
                    case 16
                        u0=g0*h0; v0=g1*h1;  u1=g0*h1; v1=g1*h0;
                        if u0>v0 
                            t0=u0; F(i,j)=g0; 
                        else
                            t0=v0; F(i,j+Dim)=g1;
                        end
                        if u1<v1 
                            t1=u1; F(i+Dim,j+Dim)=g0;
                        else
                            t1=v1; F(i+Dim,j)=g1;
                        end
                end 
            else
                h0 = x(j+Dim);  h1 = x(j);  
                %   determining the types of intervals under operation 
                if ( g0*g1>0 ) 
                    if(g0>0), l=0; else l=2; end
                else
                    if(g0<=g1), l=1; else l=3; end 
                end
                if ( h0*h1>0 ) 
                    if(h0>0), m=1; else m=3; end
                else
                    if(h0<=h1), m=2; else m=4; end 
                end
  
                %   computing multiplication in Kaucher arithmetic,
                %                   forming the subgradient matrix F   
                switch ( 4*l+m ) 
                    case 1
                        t0=g0*h0; t1=g1*h1; F(i,j+Dim)=g0; F(i+Dim,j)=g1;
                    case 2
                        t0=g1*h0; t1=g1*h1; F(i,j+Dim)=g1; F(i+Dim,j)=g1;
                    case 3
                        t0=g1*h0; t1=g0*h1; F(i,j+Dim)=g1; F(i+Dim,j)=g0;
                    case 4
                        t0=g0*h0; t1=g0*h1; F(i,j+Dim)=g0; F(i+Dim,j)=g0;
                    case 5
                        t0=g0*h1; t1=g1*h1; F(i,j)=g0; F(i+Dim,j)=g1;
                    case 6
                        u0=g0*h1; v0=g1*h0;  u1=g0*h0; v1=g1*h1;
                        if u0<v0 
                            t0=u0; F(i,j)=g0; 
                        else
                            t0=v0; F(i,j+Dim)=g1;
                        end
                        if u1>v1 
                            t1=u1; F(i+Dim,j+Dim)=g0;
                        else
                            t1=v1; F(i+Dim,j)=g1;
                        end
                    case 7
                        t0=g1*h0; t1=g0*h0; F(i,j+Dim)=g1; F(i+Dim,j+Dim)=g0;
                    case 8
                        t0=0; t1=0;
                    case 9
                        t0=g0*h1; t1=g1*h0; F(i,j)=g0; F(i+Dim,j+Dim)=g1;
                    case 10
                        t0=g0*h1; t1=g0*h0; F(i,j)=g0; F(i+Dim,j+Dim)=g0;
                    case 11
                        t0=g1*h1; t1=g0*h0; F(i,j)=g1; F(i+Dim,j+Dim)=g0;
                    case 12
                        t0=g1*h1; t1=g1*h0; F(i,j)=g1; F(i+Dim,j+Dim)=g1;
                    case 13
                        t0=g0*h0; t1=g1*h0; F(i,j+Dim)=g0; F(i+Dim,j+Dim)=g1;
                    case 14
                        t0=0; t1=0;
                    case 15
                        t0=g1*h1; t1=g0*h1; F(i,j)=g1; F(i+Dim,j)=g0;
                    case 16 
                        u0=g0*h0; v0=g1*h1;  u1=g0*h1; v1=g1*h0;
                        if u0>v0 
                            t0=u0; F(i,j+Dim)=g0; 
                        else
                            t0=v0; F(i,j)=g1; 
                        end
                        if u1<v1
                            t1=u1; F(i+Dim,j)=g0; 
                        else
                            t1=v1; F(i+Dim,j+Dim)=g1;
                        end
                end
            end  
                    s0 = s0+t0;       s1 = s1+t1;
        end
  
        t0 = s0 - d(2*i);  t1 = s1 - d(2*i+1);
        xx(i) = t0;  xx(i+Dim) = t1;
        r = r + max(abs(t0),abs(t1));
    end
  
    %    solving the (point) linear algebraic 
    %                        system with the subgradient matrix F 
  
    xx = x - tau*(F\xx);
    q = sum(abs(xx));
    if q==0  
        q=1; 
    end

end
  
  
if( ni >= IterLim )
        display('       The algorithm may diverge:')
        display('       the number of iterations alloted is exhausted!')
end
 
display(sprintf('\n Formal solution to the "Spanish" interval system: \n'));
 
for i=1:Dim
    u0=xx(i);  u1=xx(i+Dim);
    if u0<=u1 
        prop = '    ->';
    else
        prop = '    <-';
    end
    display(sprintf('%d %7s %f %s %f %s %s',i,'.   [',u0,',',u1,']',prop));
end
  
display(sprintf('\n\r Number of iterations = %d',ni));
display(sprintf(' 1-norm of the defect = %f \n',r));
  
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  