[x,ix,ierror]=dxcon(x,ix,ierror);

function [x,ix,ierror]=dxcon(x,ix,ierror);
persistent a b firstCall gt i i1 icase ispace itemp j j1 j2 z ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(icase), icase=0; end;
if isempty(ispace), ispace=0; end;
if isempty(itemp), itemp=0; end;
if isempty(j), j=0; end;
if isempty(j1), j1=0; end;
if isempty(j2), j2=0; end;
%***BEGIN PROLOGUE  DXCON
%***PURPOSE  To provide double-precision floating-point arithmetic
%            with an extended exponent range.
%***LIBRARY   SLATEC
%***CATEGORY  A3D
%***TYPE      doubleprecision (XCON-S, DXCON-D)
%***KEYWORDS  EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC
%***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
%           Smith, John M., (NBS and George Mason University)
%***DESCRIPTION
%     doubleprecision X
%     INTEGER IX
%
%                  CONVERTS (X,IX) = X*RADIX**IX
%                  TO DECIMAL FORM IN PREPARATION FOR
%                  PRINTING, SO THAT (X,IX) = X*10**IX
%                  WHERE 1/10 .LE. ABS(X) .LT. 1
%                  IS RETURNED, EXCEPT THAT IF
%                  (ABS(X),IX) IS BETWEEN RADIX**(-2L)
%                  AND RADIX**(2L) THEN THE REDUCED
%                  FORM WITH IX = 0 IS RETURNED.
%
%***SEE ALSO  DXSET
%***REFERENCES  (NONE)
%***ROUTINES CALLED  DXADJ, DXC210, DXRED
%***COMMON BLOCKS    DXBLK2
%***REVISION HISTORY  (YYMMDD)
%   820712  DATE WRITTEN
%   890126  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  DXCON
if isempty(gt), gt=0; end;
%
%   THE CONDITIONS IMPOSED ON L AND KMAX BY THIS SUBROUTINE
% ARE
%    (1) 4 .LE. L .LE. 2**NBITS - 1 - KMAX
%
%    (2) KMAX .LE. ((2**NBITS)-2)/LOG10R - L
%
% THESE CONDITIONS MUST BE MET BY APPROPRIATE CODING
% IN SUBROUTINE DXSET.
%
global dxblk2_1; if isempty(dxblk2_1), dxblk2_1=0; end;
global dxblk2_2; if isempty(dxblk2_2), dxblk2_2=0; end;
global dxblk2_3; if isempty(dxblk2_3), dxblk2_3=0; end;
global dxblk2_4; if isempty(dxblk2_4), dxblk2_4=0; end;
global dxblk2_5; if isempty(dxblk2_5), dxblk2_5=0; end;
global dxblk2_6; if isempty(dxblk2_6), dxblk2_6=0; end;
global dxblk2_7; if isempty(dxblk2_7), dxblk2_7=0; end;
% common :: ;
%% common /dxblk2/ radix , radixl , rad2l , dlg10r , l , l2 , kmax;
%% common /dxblk2/ dxblk2_1 , dxblk2_2 , dxblk2_3 , dxblk2_4 , dxblk2_5 , dxblk2_6 , dxblk2_7;
% save :: ;
% save ::         ;
%
if isempty(a), a=0; end;
if isempty(b), b=0; end;
if isempty(z), z=0; end;
%
if firstCall,   ispace=[1];  end;
firstCall=0;
%   THE PARAMETER ISPACE IS THE INCREMENT USED IN FORM-
% ING THE AUXILIARY INDEX OF THE DECIMAL EXTENDED-RANGE
% FORM. THE RETURNED VALUE OF IX WILL BE AN INTEGER MULT-
% IPLE OF ISPACE. ISPACE MUST SATISFY 1 .LE. ISPACE .LE.
% L/2. IF A VALUE GREATER THAN 1 IS TAKEN, THE RETURNED
% VALUE OF X WILL SATISFY 10**(-ISPACE) .LE. ABS(X) .LE. 1
% WHEN (ABS(X),IX) .LT. RADIX**(-2L), AND 1/10 .LE. ABS(X)
% .LT. 10**(ISPACE-1) WHEN (ABS(X),IX) .GT. RADIX**(2L).
%
%***FIRST EXECUTABLE STATEMENT  DXCON
ierror = 0;
[x,ix,ierror]=dxred(x,ix,ierror);
if( ierror==0 )
if( ix~=0 )
[x,ix,ierror]=dxadj(x,ix,ierror);
if( ierror==0 )
%
% CASE 1 IS WHEN (X,IX) IS LESS THAN RADIX**(-2L) IN MAGNITUDE,
% CASE 2 IS WHEN (X,IX) IS GREATER THAN RADIX**(2L) IN MAGNITUDE.
itemp = 1;
icase =fix(fix((3+(abs(itemp).*sign(ix)))./2));
if( icase==2 )
if( abs(x)<1.0d0 )
x = x.*dxblk2_2;
ix = fix(ix - dxblk2_5);
end;
elseif( abs(x)>=1.0d0 ) ;
x = x./dxblk2_2;
ix = fix(ix + dxblk2_5);
end;
%
% AT THIS POINT, RADIX**(-L) .LE. ABS(X) .LT. 1.0D0     IN CASE 1,
%                      1.0D0 .LE. ABS(X) .LT. RADIX**L  IN CASE 2.
i = fix(log10(abs(x))./dxblk2_4);
a = dxblk2_1.^i;
if( icase==2 )
while( a>abs(x) );
i = fix(i - 1);
a = a./dxblk2_1;
end;
while( abs(x)>=dxblk2_1.*a );
i = fix(i + 1);
a = a.*dxblk2_1;
end;
else;
while( a>dxblk2_1.*abs(x) );
i = fix(i - 1);
a = a./dxblk2_1;
end;
while( abs(x)>=a );
i = fix(i + 1);
a = a.*dxblk2_1;
end;
end;
%
% AT THIS POINT I IS SUCH THAT
% RADIX**(I-1) .LE. ABS(X) .LT. RADIX**I      IN CASE 1,
%     RADIX**I .LE. ABS(X) .LT. RADIX**(I+1)  IN CASE 2.
itemp = fix(ispace./dxblk2_4);
a = dxblk2_1.^itemp;
b = 10.0d0.^ispace;
while( a>b );
itemp = fix(itemp - 1);
a = a./dxblk2_1;
end;
while( b>=a.*dxblk2_1 );
itemp = fix(itemp + 1);
a = a.*dxblk2_1;
end;
%
% AT THIS POINT ITEMP IS SUCH THAT
% RADIX**ITEMP .LE. 10**ISPACE .LT. RADIX**(ITEMP+1).
gt=0;
while (1);
if( itemp<=0 )
% ITEMP = 0 IF, AND ONLY IF, ISPACE = 1 AND RADIX = 16.0D0
x = x.*dxblk2_1.^(-i);
ix = fix(ix + i);
[ix,z,j,ierror]=dxc210(ix,z,j,ierror);
if( ierror~=0 )
return;
end;
x = x.*z;
ix = fix(j);
if( icase==1 )
gt=1;
break;
end;
if( icase==2 )
break;
end;
end;
i1 = fix(fix(i./itemp));
x = x.*dxblk2_1.^(-i1.*itemp);
ix = fix(ix + i1.*itemp);
%
% AT THIS POINT,
% RADIX**(-ITEMP) .LE. ABS(X) .LT. 1.0D0        IN CASE 1,
%           1.0D0 .LE. ABS(X) .LT. RADIX**ITEMP IN CASE 2.
[ix,z,j,ierror]=dxc210(ix,z,j,ierror);
if( ierror~=0 )
return;
end;
j1 = fix(fix(j./ispace));
j2 = fix(j - j1.*ispace);
x = x.*z.*10.0d0.^j2;
ix = fix(j1.*ispace);
%
% AT THIS POINT,
%  10.0D0**(-2*ISPACE) .LE. ABS(X) .LT. 1.0D0                IN CASE 1,
%           10.0D0**-1 .LE. ABS(X) .LT. 10.0D0**(2*ISPACE-1) IN CASE 2.
if( icase~=2 )
gt=1;
break;
end;
break;
end;
if(gt==0)
while( 10.0d0.*abs(x)>=b );
x = x./b;
ix = fix(ix + ispace);
end;
return;
end;
while( b.*abs(x)<1.0d0 );
x = x.*b;
ix = fix(ix - ispace);
end;
end;
end;
end;
end %subroutine dxcon
%DECK DX