[x,ix,y,iy,z,iz,ierror]=xadd(x,ix,y,iy,z,iz,ierror);

function [x,ix,y,iy,z,iz,ierror]=xadd(x,ix,y,iy,z,iz,ierror);
persistent i i1 i2 igo is j s t ; 

if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(i2), i2=0; end;
if isempty(is), is=0; end;
if isempty(j), j=0; end;
if isempty(igo), igo=0; end;
if isempty(s), s=0; end;
if isempty(t), t=0; end;
%***BEGIN PROLOGUE  XADD
%***PURPOSE  To provide single-precision floating-point arithmetic
%            with an extended exponent range.
%***LIBRARY   SLATEC
%***CATEGORY  A3D
%***TYPE      SINGLE PRECISION (XADD-S, DXADD-D)
%***KEYWORDS  EXTENDED-RANGE SINGLE-PRECISION ARITHMETIC
%***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
%           Smith, John M., (NBS and George Mason University)
%***DESCRIPTION
%     REAL X, Y, Z
%     INTEGER IX, IY, IZ
%
%                  FORMS THE EXTENDED-RANGE SUM  (Z,IZ) =
%                  (X,IX) + (Y,IY).  (Z,IZ) IS ADJUSTED
%                  BEFORE RETURNING. THE INPUT OPERANDS
%                  NEED NOT BE IN ADJUSTED FORM, BUT THEIR
%                  PRINCIPAL PARTS MUST SATISFY
%                  RADIX**(-2L).LE.ABS(X).LE.RADIX**(2L),
%                  RADIX**(-2L).LE.ABS(Y).LE.RADIX**(2L).
%
%***SEE ALSO  XSET
%***REFERENCES  (NONE)
%***ROUTINES CALLED  XADJ
%***COMMON BLOCKS    XBLK2
%***REVISION HISTORY  (YYMMDD)
%   820712  DATE WRITTEN
%   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
%   901019  Revisions to prologue.  (DWL and WRB)
%   901106  Changed all specific intrinsics to generic.  (WRB)
%           Corrected order of sections in prologue and added TYPE
%           section.  (WRB)
%   920127  Revised PURPOSE section of prologue.  (DWL)
%***end PROLOGUE  XADD
global xblk2_1; if isempty(xblk2_1), xblk2_1=0; end;
global xblk2_2; if isempty(xblk2_2), xblk2_2=0; end;
global xblk2_3; if isempty(xblk2_3), xblk2_3=0; end;
global xblk2_4; if isempty(xblk2_4), xblk2_4=0; end;
global xblk2_5; if isempty(xblk2_5), xblk2_5=0; end;
global xblk2_6; if isempty(xblk2_6), xblk2_6=0; end;
global xblk2_7; if isempty(xblk2_7), xblk2_7=0; end;
% common :: ;
%% common /xblk2 / radix , radixl , rad2l , dlg10r , l , l2 , kmax;
%% common /xblk2 / xblk2_1 , xblk2_2 , xblk2_3 , xblk2_4 , xblk2_5 , xblk2_6 , xblk2_7;
% save :: ;
% save ::       ;
%
%
%   THE CONDITIONS IMPOSED ON L AND KMAX BY THIS SUBROUTINE
% ARE
%     (1) 1 .LT. L .LE. 0.5*LOGR(0.5*DZERO)
%
%     (2) NRADPL .LT. L .LE. KMAX/6
%
%     (3) KMAX .LE. (2**NBITS - 4*L - 1)/2
%
% THESE CONDITIONS MUST BE MET BY APPROPRIATE CODING
% IN SUBROUTINE XSET.
%
%***FIRST EXECUTABLE STATEMENT  XADD
ierror = 0;
if( x==0.0 )
z = y;
iz = fix(iy);
[z,iz,ierror]=xadj(z,iz,ierror);
elseif( y==0.0 ) ;
z = x;
iz = fix(ix);
[z,iz,ierror]=xadj(z,iz,ierror);
else;
if( ix<0 || iy<0 )
if( ix>=0 || iy>=0 )
if( abs(ix)>6.*xblk2_5 || abs(iy)>6.*xblk2_5 )
if( ix>=0 )
z = x;
iz = fix(ix);
else;
z = y;
iz = fix(iy);
end;
[z,iz,ierror]=xadj(z,iz,ierror);
return;
end;
end;
end;
i = fix(ix - iy);
if( i<0 )
s = y;
is = fix(iy);
t = x;
elseif( i==0 ) ;
if( abs(x)>1.0 && abs(y)>1.0 )
s = x./xblk2_2;
t = y./xblk2_2;
z = s + t;
iz = fix(ix + xblk2_5);
elseif( abs(x)<1.0 && abs(y)<1.0 ) ;
s = x.*xblk2_2;
t = y.*xblk2_2;
z = s + t;
iz = fix(ix - xblk2_5);
else;
z = x + y;
iz = fix(ix);
end;
[z,iz,ierror]=xadj(z,iz,ierror);
return;
else;
s = x;
is = fix(ix);
t = y;
end;
%
%  AT THIS POINT, THE ONE OF (X,IX) OR (Y,IY) THAT HAS THE
% LARGER AUXILIARY INDEX IS STORED IN (S,IS). THE PRINCIPAL
% PART OF THE OTHER INPUT IS STORED IN T.
%
i1 = fix(abs(i)./xblk2_5);
i2 = fix(rem(abs(i),xblk2_5));
while (1);
if( abs(t)>=xblk2_2 )
j = fix(i1 - 2);
if( j>=0 )
t = t.*xblk2_1.^(-i2)./xblk2_3;
break;
end;
elseif( abs(t)<1.0 ) ;
if( xblk2_2.*abs(t)>=1.0 )
j = fix(i1);
t = t.*xblk2_1.^(-i2);
break;
end;
j = fix(i1 + 1);
t = t.*xblk2_1.^(xblk2_5-i2);
break;
end;
j = fix(i1 - 1);
if( j>=0 )
t = t.*xblk2_1.^(-i2)./xblk2_2;
break;
end;
j = fix(i1);
t = t.*xblk2_1.^(-i2);
break;
end;
%
%  AT THIS POINT, SOME OR ALL OF THE DIFFERENCE IN THE
% AUXILIARY INDICES HAS BEEN USED TO EFFECT A LEFT SHIFT
% OF T.  THE SHIFTED VALUE OF T SATISFIES
%
%       RADIX**(-2*L) .LE. ABS(T) .LE. 1.0
%
% AND, IF J=0, NO FURTHER SHIFTING REMAINS TO BE DONE.
%
if( j~=0 )
igo=0;
while (1);
if( abs(s)<xblk2_2 && j<=3 )
if( abs(s)>=1.0 )
if( j==1 )
s = s.*xblk2_2;
igo=1;
break;
elseif( j==2 || j==3 ) ;
break;
end;
end;
if( xblk2_2.*abs(s)>=1.0 )
if( j==1 )
s = s.*xblk2_2;
igo=1;
break;
elseif( j==2 ) ;
s = s.*xblk2_2;
s = s.*xblk2_2;
igo=1;
break;
elseif( j==3 ) ;
break;
end;
end;
if( j==1 )
s = s.*xblk2_2;
else;
if( j~=2 )
if( j~=3 )
break;
end;
s = s.*xblk2_2;
end;
s = s.*xblk2_2;
s = s.*xblk2_2;
end;
igo=1;
break;
end;
break;
end;
if(igo==0)
z = s;
iz = fix(is);
[z,iz,ierror]=xadj(z,iz,ierror);
return;
end;
end;
%
%   AT THIS POINT, THE REMAINING DIFFERENCE IN THE
% AUXILIARY INDICES HAS BEEN USED TO EFFECT A RIGHT SHIFT
% OF S.  IF THE SHIFTED VALUE OF S WOULD HAVE EXCEEDED
% RADIX**L, THEN (S,IS) IS RETURNED AS THE VALUE OF THE
% SUM.
%
if( abs(s)>1.0 && abs(t)>1.0 )
s = s./xblk2_2;
t = t./xblk2_2;
z = s + t;
iz = fix(is - j.*xblk2_5 + xblk2_5);
elseif( abs(s)<1.0 && abs(t)<1.0 ) ;
s = s.*xblk2_2;
t = t.*xblk2_2;
z = s + t;
iz = fix(is - j.*xblk2_5 - xblk2_5);
else;
z = s + t;
iz = fix(is - j.*xblk2_5);
end;
[z,iz,ierror]=xadj(z,iz,ierror);
end;
end %subroutine xadd
%DECK XADJ