| [rs,re,z,trunc]=mpnzr(rs,re,z,trunc); |
function [rs,re,z,trunc]=mpnzr(rs,re,z,trunc);
persistent b2 gt i i2 i2m i2p ii is it j k ;
if isempty(i), i=0; end;
if isempty(i2), i2=0; end;
if isempty(i2m), i2m=0; end;
if isempty(i2p), i2p=0; end;
if isempty(is), is=0; end;
if isempty(it), it=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
global mpcom_4; if isempty(mpcom_4), mpcom_4=0; end;
global mpcom_3; if isempty(mpcom_3), mpcom_3=0; end;
global mpcom_5; if isempty(mpcom_5), mpcom_5=0; end;
if isempty(ii), ii=0; end;
%***BEGIN PROLOGUE MPNZR
%***SUBSIDIARY
%***PURPOSE Subsidiary to DQDOTA and DQDOTI
%***LIBRARY SLATEC
%***TYPE ALL (MPNZR-A)
%***AUTHOR (UNKNOWN)
%***DESCRIPTION
%
% Modified for use with BLAS. Blank COMMON changed to named COMMON.
% Assumes long (i.e. (mpcom_2+4)-DIGIT) fraction in R, sign = RS, exponent
% = RE. Normalizes, and returns 'mp' result in Z. Integer arguments
% RS and RE are not preserved. R*-rounding is used if TRUNC.EQ.0
%
% The argument Z(*) and the variable R in COMMON are INTEGER arrays
% of size 30. See the comments in the routine MPBLAS for the reason
% for this choice.
%
%***SEE ALSO DQDOTA, DQDOTI, MPBLAS
%***ROUTINES CALLED MPERR, MPOVFL, MPUNFL
%***COMMON BLOCKS MPCOM
%***REVISION HISTORY (YYMMDD)
% 791001 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900402 Added TYPE section. (WRB)
% 930124 Increased Array size in MPCON for SUN -r8. (RWC)
%***end PROLOGUE MPNZR
% common :: ;
global mpcom_1; if isempty(mpcom_1), mpcom_1=0; end;
global mpcom_2; if isempty(mpcom_2), mpcom_2=0; end;
global mpcom_6; if isempty(mpcom_6), mpcom_6=zeros(1,30); end;
%% common /mpcom / b , t , m , lun , mxr , r(30);
%% common /mpcom / mpcom_1 , mpcom_2 , mpcom_3 , mpcom_4 , mpcom_5 , mpcom_6(30);
z_shape=size(z);z=reshape(z,1,[]);
if isempty(b2), b2=0; end;
if isempty(gt), gt=zeros(1,5); end;
%***FIRST EXECUTABLE STATEMENT MPNZR
i2 = fix(mpcom_2 + 4);
gt(:)=0;
if( rs~=0 )
% CHECK THAT SIGN = +-1
if( abs(rs)<=1 )
% LOOK FOR FIRST NONZERO DIGIT
for i = 1 : i2;
is = fix(i - 1);
if( mpcom_6(i)>0 )
if( is~=0 )
% NORMALIZE
re = fix(re - is);
i2m = fix(i2 - is);
for j = 1 : i2m;
k = fix(j + is);
mpcom_6(j) = fix(mpcom_6(k));
end; j = fix(i2m+1);
i2p = fix(i2m + 1);
for j = i2p : i2;
mpcom_6(j) = 0;
end; j = fix(i2+1);
end;
% CHECK TO SEE IF TRUNCATION IS DESIRED
if( trunc==0 )
% SEE IF ROUNDING NECESSARY
% TREAT EVEN AND ODD BASES DIFFERENTLY
b2 = fix(fix(mpcom_1./2));
gt(:)=0;
while (1);
if((2.*b2)~=mpcom_1 )
% ODD BASE, ROUND IF R(T+1)... .GT. 1/2
for ii = 1 : 4;
it = fix(mpcom_2 + ii);
if( mpcom_6(it)<b2 )
break;
end;
if( mpcom_6(it)~=b2 )
gt(2)=1;
break;
end;
end;
if(gt(2)==1)
break;
end;
gt(5)=1;
break;
else;
% B EVEN. ROUND IF R(T+1).GE.B2 UNLESS R(T) ODD AND ALL ZEROS
% AFTER R(T+2).
if( mpcom_6(mpcom_2+1)<b2 )
gt(5)=1;
break;
end;
if( mpcom_6(mpcom_2+1)==b2 )
if( rem(mpcom_6(mpcom_2),2)~=0 )
if((mpcom_6(mpcom_2+2)+mpcom_6(mpcom_2+3)+mpcom_6(mpcom_2+4))==0 )
gt(5)=1;
break;
end;
end;
end;
end;
break;
end;
gt(2)=0;
% ROUND
if(gt(5)==0)
for j = 1 : mpcom_2;
ii = fix(mpcom_2 + 1 - j);
mpcom_6(ii) = fix(mpcom_6(ii) + 1);
if( mpcom_6(ii)<mpcom_1 )
gt(5)=1;
break;
end;
mpcom_6(ii) = 0;
end;
if(gt(5)==0)
% EXCEPTIONAL CASE, ROUNDED UP TO .10000...
re = fix(re + 1);
mpcom_6(1) = 1;
end;
end;
gt(5)=0;
end;
% CHECK FOR OVERFLOW
if( re>mpcom_3 )
writef(mpcom_4,[' *** OVERFLOW OCCURRED IN MPNZR ***' ' \n']);
%format (' *** OVERFLOW OCCURRED IN MPNZR ***');
[z]=mpovfl(z);
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
% CHECK FOR UNDERFLOW
elseif( re<(-mpcom_3) ) ;
% UNDERFLOW HERE
[z]=mpunfl(z);
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
else;
% STORE RESULT IN Z
z(1) = fix(rs);
z(2) = fix(re);
for ii = 1 : mpcom_2;
z(ii+2) = fix(mpcom_6(ii));
end; ii = fix(mpcom_2+1);
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
end;
end;
% FRACTION ZERO
end;
else;
writef(mpcom_4,[' *** SIGN NOT 0, +1 OR -1 IN CALL TO MPNZR,',' POSSIBLE OVERWRITING PROBLEM ***' ' \n']);
%format (' *** SIGN NOT 0, +1 OR -1 IN CALL TO MPNZR,',' POSSIBLE OVERWRITING PROBLEM ***');
mperr;
end;
end;
% STORE ZERO IN Z
z(1) = 0;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
end %subroutine mpnzr
%DECK MPOVFL
|
|