| [irad,nradpl,dzero,nbits,ierror]=dxset(irad,nradpl,dzero,nbits,ierror); |
function [irad,nradpl,dzero,nbits,ierror]=dxset(irad,nradpl,dzero,nbits,ierror);
persistent dzerox firstCall i ic iflag ii imaxex iminex iradx it j k kk lg102x lgtemp log102 log2r lx nb nbitsx np1 nrdplc ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(ii), ii=0; end;
if isempty(imaxex), imaxex=0; end;
if isempty(iminex), iminex=0; end;
if isempty(iradx), iradx=0; end;
if isempty(it), it=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(lg102x), lg102x=0; end;
if isempty(lgtemp), lgtemp=zeros(1,20); end;
if isempty(log102), log102=zeros(1,20); end;
if isempty(log2r), log2r=0; end;
if isempty(lx), lx=0; end;
if isempty(nb), nb=0; end;
global dxblk1_1; if isempty(dxblk1_1), dxblk1_1=0; end;
if isempty(nbitsx), nbitsx=0; end;
if isempty(np1), np1=0; end;
if isempty(nrdplc), nrdplc=0; end;
%***BEGIN PROLOGUE DXSET
%***PURPOSE To provide double-precision floating-point arithmetic
% with an extended exponent range.
%***LIBRARY SLATEC
%***CATEGORY A3D
%***TYPE doubleprecision (XSET-S, DXSET-D)
%***KEYWORDS EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC
%***AUTHOR Lozier, Daniel W., (National Bureau of Standards)
% Smith, John M., (NBS and George Mason University)
%***DESCRIPTION
%
% subroutine DXSET MUST BE CALLED PRIOR TO CALLING ANY OTHER
% EXTENDED-RANGE SUBROUTINE. IT CALCULATES AND STORES SEVERAL
% MACHINE-DEPENDENT CONSTANTS IN COMMON BLOCKS. THE USER MUST
% SUPPLY FOUR CONSTANTS THAT PERTAIN TO HIS PARTICULAR COMPUTER.
% THE CONSTANTS ARE
%
% IRAD = THE INTERNAL BASE OF DOUBLE-PRECISION
% ARITHMETIC IN THE COMPUTER.
% NRADPL = THE NUMBER OF RADIX PLACES CARRIED IN
% THE DOUBLE-PRECISION REPRESENTATION.
% DZERO = THE SMALLEST OF 1/DMIN, DMAX, DMAXLN WHERE
% DMIN = THE SMALLEST POSITIVE DOUBLE-PRECISION
% NUMBER OR AN UPPER BOUND TO THIS NUMBER,
% DMAX = THE LARGEST DOUBLE-PRECISION NUMBER
% OR A LOWER BOUND TO THIS NUMBER,
% DMAXLN = THE LARGEST DOUBLE-PRECISION NUMBER
% SUCH THAT LOG10(DMAXLN) CAN BE COMPUTED BY THE
% FORTRAN SYSTEM (ON MOST SYSTEMS DMAXLN = DMAX).
% NBITS = THE NUMBER OF BITS (EXCLUSIVE OF SIGN) IN
% AN INTEGER COMPUTER WORD.
%
% ALTERNATIVELY, ANY OR ALL OF THE CONSTANTS CAN BE GIVEN
% THE VALUE 0 (0.0D0 FOR DZERO). IF A CONSTANT IS ZERO, DXSET TRIES
% TO ASSIGN AN APPROPRIATE VALUE BY CALLING I1MACH
% (SEE P.A.FOX, A.D.HALL, N.L.SCHRYER, ALGORITHM 528 FRAMEWORK
% FOR A PORTABLE LIBRARY, ACM TRANSACTIONS ON MATH SOFTWARE,
% V.4, NO.2, JUNE 1978, 177-188).
%
% THIS IS THE SETTING-UP SUBROUTINE FOR A PACKAGE OF SUBROUTINES
% THAT FACILITATE THE USE OF EXTENDED-RANGE ARITHMETIC. EXTENDED-RANGE
% ARITHMETIC ON A PARTICULAR COMPUTER IS DEFINED ON THE SET OF NUMBERS
% OF THE FORM
%
% (X,IX) = X*RADIX**IX
%
% WHERE X IS A DOUBLE-PRECISION NUMBER CALLED THE PRINCIPAL PART,
% IX IS AN INTEGER CALLED THE AUXILIARY INDEX, AND RADIX IS THE
% INTERNAL BASE OF THE DOUBLE-PRECISION ARITHMETIC. OBVIOUSLY,
% EACH REAL NUMBER IS REPRESENTABLE WITHOUT ERROR BY MORE THAN ONE
% EXTENDED-RANGE FORM. CONVERSIONS BETWEEN DIFFERENT FORMS ARE
% ESSENTIAL IN CARRYING OUT ARITHMETIC OPERATIONS. WITH THE CHOICE
% OF RADIX WE HAVE MADE, AND THE SUBROUTINES WE HAVE WRITTEN, THESE
% CONVERSIONS ARE PERFORMED WITHOUT ERROR (AT LEAST ON MOST COMPUTERS).
% (SEE SMITH, J.M., OLVER, F.W.J., AND LOZIER, D.W., EXTENDED-RANGE
% ARITHMETIC AND NORMALIZED LEGENDRE POLYNOMIALS, ACM TRANSACTIONS ON
% MATHEMATICAL SOFTWARE, MARCH 1981).
%
% AN EXTENDED-RANGE NUMBER (X,IX) IS SAID TO BE IN ADJUSTED FORM IF
% X AND IX ARE ZERO OR
%
% RADIX**(-L) .LE. ABS(X) .LT. RADIX**L
%
% IS SATISFIED, WHERE L IS A COMPUTER-DEPENDENT INTEGER DEFINED IN THIS
% subroutine. TWO EXTENDED-RANGE NUMBERS IN ADJUSTED FORM CAN BE ADDED,
% SUBTRACTED, MULTIPLIED OR DIVIDED (IF THE DIVISOR IS NONZERO) WITHOUT
% CAUSING OVERFLOW OR UNDERFLOW IN THE PRINCIPAL PART OF THE RESULT.
% WITH PROPER USE OF THE EXTENDED-RANGE SUBROUTINES, THE ONLY OVERFLOW
% THAT CAN OCCUR IS INTEGER OVERFLOW IN THE AUXILIARY INDEX. IF THIS
% IS DETECTED, THE SOFTWARE CALLS XERROR (A GENERAL ERROR-HANDLING
% FORTRAN SUBROUTINE PACKAGE).
%
% MULTIPLICATION AND DIVISION IS PERFORMED BY SETTING
%
% (X,IX)*(Y,IY) = (X*Y,IX+IY)
% OR
% (X,IX)/(Y,IY) = (X/Y,IX-IY).
%
% PRE-ADJUSTMENT OF THE OPERANDS IS ESSENTIAL TO AVOID
% OVERFLOW OR UNDERFLOW OF THE PRINCIPAL PART. SUBROUTINE
% DXADJ (SEE BELOW) MAY BE CALLED TO TRANSFORM ANY EXTENDED-
% RANGE NUMBER INTO ADJUSTED FORM.
%
% ADDITION AND SUBTRACTION REQUIRE THE USE OF SUBROUTINE DXADD
% (SEE BELOW). THE INPUT OPERANDS NEED NOT BE IN ADJUSTED FORM.
% HOWEVER, THE RESULT OF ADDITION OR SUBTRACTION IS RETURNED
% IN ADJUSTED FORM. THUS, FOR EXAMPLE, IF (X,IX),(Y,IY),
% (U,IU), AND (V,IV) ARE IN ADJUSTED FORM, THEN
%
% (X,IX)*(Y,IY) + (U,IU)*(V,IV)
%
% CAN BE COMPUTED AND STORED IN ADJUSTED FORM WITH NO EXPLICIT
% CALLS TO DXADJ.
%
% WHEN AN EXTENDED-RANGE NUMBER IS TO BE PRINTED, IT MUST BE
% CONVERTED TO AN EXTENDED-RANGE FORM WITH DECIMAL RADIX. SUBROUTINE
% DXCON IS PROVIDED FOR THIS PURPOSE.
%
% THE SUBROUTINES CONTAINED IN THIS PACKAGE ARE
%
% subroutine DXADD
% USAGE
% CALL DXADD(X,IX,Y,IY,Z,IZ,IERROR)
% IF (IERROR.NE.0) RETURN
% DESCRIPTION
% 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).
%
% subroutine DXADJ
% USAGE
% CALL DXADJ(X,IX,IERROR)
% IF (IERROR.NE.0) RETURN
% DESCRIPTION
% TRANSFORMS (X,IX) SO THAT
% RADIX**(-L) .LE. ABS(X) .LT. RADIX**L.
% ON MOST COMPUTERS THIS TRANSFORMATION DOES
% NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS
% THE NUMBER BASE OF DOUBLE-PRECISION ARITHMETIC.
%
% subroutine DXC210
% USAGE
% CALL DXC210(K,Z,J,IERROR)
% IF (IERROR.NE.0) RETURN
% DESCRIPTION
% GIVEN K THIS SUBROUTINE COMPUTES J AND Z
% SUCH THAT RADIX**K = Z*10**J, WHERE Z IS IN
% THE RANGE 1/10 .LE. Z .LT. 1.
% THE VALUE OF Z WILL BE ACCURATE TO FULL
% DOUBLE-PRECISION PROVIDED THE NUMBER
% OF DECIMAL PLACES IN THE LARGEST
% INTEGER PLUS THE NUMBER OF DECIMAL
% PLACES CARRIED IN DOUBLE-PRECISION DOES NOT
% EXCEED 60. DXC210 IS CALLED BY SUBROUTINE
% DXCON WHEN NECESSARY. THE USER SHOULD
% NEVER NEED TO CALL DXC210 DIRECTLY.
%
% subroutine DXCON
% USAGE
% CALL DXCON(X,IX,IERROR)
% IF (IERROR.NE.0) RETURN
% DESCRIPTION
% 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.
%
% subroutine DXRED
% USAGE
% CALL DXRED(X,IX,IERROR)
% IF (IERROR.NE.0) RETURN
% DESCRIPTION
% IF
% RADIX**(-2L) .LE. (ABS(X),IX) .LE. RADIX**(2L)
% THEN DXRED TRANSFORMS (X,IX) SO THAT IX=0.
% IF (X,IX) IS OUTSIDE THE ABOVE RANGE,
% THEN DXRED TAKES NO ACTION.
% THIS SUBROUTINE IS USEFUL IF THE
% RESULTS OF EXTENDED-RANGE CALCULATIONS
% ARE TO BE USED IN SUBSEQUENT ORDINARY
% DOUBLE-PRECISION CALCULATIONS.
%
%***REFERENCES Smith, Olver and Lozier, Extended-Range Arithmetic and
% Normalized Legendre Polynomials, ACM Trans on Math
% Softw, v 7, n 1, March 1981, pp 93--105.
%***ROUTINES CALLED I1MACH, XERMSG
%***COMMON BLOCKS DXBLK1, DXBLK2, DXBLK3
%***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)
% CALLs to XERROR changed to CALLs to XERMSG. (WRB)
% 920127 Revised PURPOSE section of prologue. (DWL)
%***end PROLOGUE DXSET
if isempty(dzerox), dzerox=0; end;
% common :: ;
%% common /dxblk1/ nbitsf;
%% common /dxblk1/ dxblk1_1;
% save :: ;
% save :: ;
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 :: ;
global dxblk3_1; if isempty(dxblk3_1), dxblk3_1=0; end;
global dxblk3_2; if isempty(dxblk3_2), dxblk3_2=0; end;
global dxblk3_3; if isempty(dxblk3_3), dxblk3_3=zeros(1,21); end;
% common :: ;
%% common /dxblk3/ nlg102 , mlg102 , lg102(21);
%% common /dxblk3/ dxblk3_1 , dxblk3_2 , dxblk3_3(21);
% save :: ;
% save :: ;
if isempty(iflag), iflag=0; end;
%
%
% LOG102 CONTAINS THE FIRST 60 DIGITS OF LOG10(2) FOR USE IN
% CONVERSION OF EXTENDED-RANGE NUMBERS TO BASE 10 .
if firstCall, log102=[301,029,995,663,981,195,213,738,894,724,493,026,768,189,881,462,108,541,310,428]; end;
%
% FOLLOWING CODING PREVENTS DXSET FROM BEING EXECUTED MORE THAN ONCE.
% THIS IS IMPORTANT BECAUSE SOME SUBROUTINES (SUCH AS DXNRMP AND
% DXLEGF) CALL DXSET TO MAKE SURE EXTENDED-RANGE ARITHMETIC HAS
% BEEN INITIALIZED. THE USER MAY WANT TO PRE-EMPT THIS CALL, FOR
% EXAMPLE WHEN I1MACH IS NOT AVAILABLE. SEE CODING BELOW.
if firstCall, iflag=[0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DXSET
ierror = 0;
if( iflag~=0 )
return;
end;
iradx = fix(irad);
nrdplc = fix(nradpl);
dzerox = dzero;
iminex = 0;
imaxex = 0;
nbitsx = fix(nbits);
% FOLLOWING 5 STATEMENTS SHOULD BE DELETED IF I1MACH IS
% NOT AVAILABLE OR NOT CONFIGURED TO RETURN THE CORRECT
% MACHINE-DEPENDENT VALUES.
if( iradx==0 )
[ iradx ]=i1mach(10);
end;
if( nrdplc==0 )
[ nrdplc ]=i1mach(14);
end;
if( dzerox==0.0d0 )
[ iminex ]=i1mach(15);
end;
if( dzerox==0.0d0 )
[ imaxex ]=i1mach(16);
end;
if( nbitsx==0 )
[ nbitsx ]=i1mach(8);
end;
if( iradx~=2 )
if( iradx~=4 )
if( iradx~=8 )
if( iradx~=16 )
xermsg('SLATEC','DXSET','IMPROPER VALUE OF IRAD',201,1);
ierror = 201;
return;
end;
end;
end;
end;
log2r = 0;
if( iradx==2 )
log2r = 1;
end;
if( iradx==4 )
log2r = 2;
end;
if( iradx==8 )
log2r = 3;
end;
if( iradx==16 )
log2r = 4;
end;
dxblk1_1 = fix(log2r.*nrdplc);
dxblk2_1 = iradx;
dxblk2_4 = log10(dxblk2_1);
if( dzerox~=0.0d0 )
lx = fix(0.5d0.*log10(dzerox)./dxblk2_4);
% RADIX**(2*L) SHOULD NOT OVERFLOW, BUT REDUCE L BY 1 FOR FURTHER
% PROTECTION.
lx = fix(lx - 1);
else;
lx = fix(min(fix((1-iminex)./2),fix((imaxex-1)./2)));
end;
dxblk2_6 = fix(2.*lx);
if( lx>=4 )
dxblk2_5 = fix(lx);
dxblk2_2 = dxblk2_1.^dxblk2_5;
dxblk2_3 = dxblk2_2.^2;
% IT IS NECESSARY TO RESTRICT NBITS (OR NBITSX) TO BE LESS THAN SOME
% UPPER LIMIT BECAUSE OF BINARY-TO-DECIMAL CONVERSION. SUCH CONVERSION
% IS DONE BY DXC210 AND REQUIRES A CONSTANT THAT IS STORED TO SOME FIXED
% PRECISION. THE STORED CONSTANT (LOG102 IN THIS ROUTINE) PROVIDES
% FOR CONVERSIONS ACCURATE TO THE LAST DECIMAL DIGIT WHEN THE INTEGER
% WORD LENGTH DOES NOT EXCEED 63. A LOWER LIMIT OF 15 BITS IS IMPOSED
% BECAUSE THE SOFTWARE IS DESIGNED TO RUN ON COMPUTERS WITH INTEGER WORD
% LENGTH OF AT LEAST 16 BITS.
if( 15<=nbitsx && nbitsx<=63 )
dxblk2_7 = fix(2.^(nbitsx-1) - dxblk2_6);
nb =fix(fix((nbitsx-1)./2));
dxblk3_2 = fix(2.^nb);
if( 1<=nrdplc.*log2r && nrdplc.*log2r<=120 )
dxblk3_1 = fix(fix((nrdplc.*log2r)./nb) + 3);
np1 = fix(dxblk3_1 + 1);
%
% AFTER COMPLETION OF THE FOLLOWING LOOP, IC CONTAINS
% THE INTEGER PART AND LGTEMP CONTAINS THE FRACTIONAL PART
% OF LOG10(IRADX) IN RADIX 1000.
ic = 0;
for ii = 1 : 20;
i = fix(21 - ii);
it = fix(log2r.*log102(i) + ic);
ic = fix(fix(it./1000));
lgtemp(i) = fix(rem(it,1000));
end; ii = fix(20+1);
%
% AFTER COMPLETION OF THE FOLLOWING LOOP, LG102 CONTAINS
% LOG10(IRADX) IN RADIX MLG102. THE RADIX POINT IS
% BETWEEN LG102(1) AND LG102(2).
dxblk3_3(1) = fix(ic);
for i = 2 : np1;
lg102x = 0;
for j = 1 : nb;
ic = 0;
for kk = 1 : 20;
k = fix(21 - kk);
it = fix(2.*lgtemp(k) + ic);
ic = fix(fix(it./1000));
lgtemp(k) = fix(rem(it,1000));
end; kk = fix(20+1);
lg102x = fix(2.*lg102x + ic);
end; j = fix(nb+1);
dxblk3_3(i) = fix(lg102x);
end; i = fix(np1+1);
%
% CHECK SPECIAL CONDITIONS REQUIRED BY SUBROUTINES...
if( nrdplc>=dxblk2_5 )
xermsg('SLATEC','DXSET','NRADPL .GE. L',205,1);
ierror = 205;
return;
elseif( 6.*dxblk2_5<=dxblk2_7 ) ;
iflag = 1;
else;
xermsg('SLATEC','DXSET','6*L .GT. KMAX',206,1);
ierror = 206;
return;
end;
else;
xermsg('SLATEC','DXSET','IMPROPER VALUE OF NRADPL',204,1);
ierror = 204;
return;
end;
else;
xermsg('SLATEC','DXSET','IMPROPER VALUE OF NBITS',203,1);
ierror = 203;
return;
end;
else;
xermsg('SLATEC','DXSET','IMPROPER VALUE OF DZERO',202,1);
ierror = 202;
return;
end;
end
%DECK DY4
|
|
