[x,iy,z]=mpdivi(x,iy,z);

function [x,iy,z]=mpdivi(x,iy,z);
persistent b2 c c2 gt i i2 iq iqj ir j j1 j11 j2 k kh r1 re rs ; 

if isempty(i), i=0; end;
if isempty(i2), i2=0; end;
if isempty(iq), iq=0; end;
if isempty(iqj), iqj=0; end;
if isempty(ir), ir=0; end;
if isempty(j), j=0; end;
if isempty(j1), j1=0; end;
if isempty(j11), j11=0; end;
if isempty(j2), j2=0; end;
if isempty(k), k=0; end;
if isempty(kh), kh=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;
%***BEGIN PROLOGUE  MPDIVI
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DQDOTA and DQDOTI
%***LIBRARY   SLATEC
%***TYPE      ALL (MPDIVI-A)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%  Divides 'mp' X by the single-precision integer IY giving 'mp' Z.
%  This is much faster than division by an 'mp' number.
%
%  The arguments X(*) and Z(*), and the variable R in COMMON are all
%  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  MPCHK, MPERR, MPNZR, MPSTR, MPUNFL
%***COMMON BLOCKS    MPCOM
%***REVISION HISTORY  (YYMMDD)
%   791001  DATE WRITTEN
%   ??????  Modified for use with BLAS.  Blank COMMON changed to named
%           COMMON.  R given dimension 12.
%   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  MPDIVI
% 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);
x_shape=size(x);x=reshape(x,1,[]);
z_shape=size(z);z=reshape(z,1,[]);
if isempty(rs), rs=0; end;
if isempty(re), re=0; end;
if isempty(r1), r1=0; end;
if isempty(c), c=0; end;
if isempty(c2), c2=0; end;
if isempty(b2), b2=0; end;
if isempty(gt), gt=zeros(1,3); end;
%***FIRST EXECUTABLE STATEMENT  MPDIVI
rs = fix(x(1));
j = fix(iy);
gt(:)=0;
if( j<0 )
j = fix(-j);
rs = fix(-rs);
elseif( j==0 ) ;
writef(mpcom_4,[' *** ATTEMPTED DIVISION BY ZERO IN CALL TO MPDIVI ***' ' \n']);
%format (' *** ATTEMPTED DIVISION BY ZERO IN CALL TO MPDIVI ***');
gt(3)=1;
end;
if(gt(3)==0)
re = fix(x(2));
% CHECK FOR ZERO DIVIDEND
if( rs~=0 )
% CHECK FOR DIVISION BY B
if( j==mpcom_1 )
[x,z]=mpstr(x,z);
if( re<=(-mpcom_3) )
% CHECK FOR DIVISION BY 1 OR -1
% UNDERFLOW HERE
[z]=mpunfl(z);
else;
z(1) = fix(rs);
z(2) = fix(re - 1);
end;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
elseif( j==1 ) ;
[x,z]=mpstr(x,z);
z(1) = fix(rs);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
else;
while (1);
if(gt(1)==0)
c = 0;
i2 = fix(mpcom_2 + 4);
i = 0;
% IF J*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
% LONG DIVISION.   ASSUME AT LEAST 16-BIT WORD.
b2 = fix(max((8.*mpcom_1),fix(32767./mpcom_1)));
if( j>=b2 )
% HERE NEED SIMULATED DOUBLE-PRECISION DIVISION
c2 = 0;
j1 = fix(fix(j./mpcom_1));
j2 = fix(j - j1.*mpcom_1);
j11 = fix(j1 + 1);
% LOOK FOR FIRST NONZERO DIGIT
while( true );
i = fix(i + 1);
c = fix(mpcom_1.*c + c2);
c2 = 0;
if( i<=mpcom_2 )
c2 = fix(x(i+2));
end;
if( c>=j1 )
if( c~=j1 )
break;
end;
if( c2>=j2 )
break;
end;
end;
end;
% COMPUTE T+4 QUOTIENT DIGITS
re = fix(re + 1 - i);
k = 1;
% GET APPROXIMATE QUOTIENT FIRST
while( true );
ir = fix(fix(c./j11));
% NOW REDUCE SO OVERFLOW DOES NOT OCCUR
iq = fix(c - ir.*j1);
if( iq>=b2 )
% HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR
ir = fix(ir + 1);
iq = fix(iq - j1);
end;
iq = fix(iq.*mpcom_1 - ir.*j2);
if( iq<0 )
% HERE IQ NEGATIVE SO IR WAS TOO LARGE
ir = fix(ir - 1);
iq = fix(iq + j);
end;
if( i<=mpcom_2 )
iq = fix(iq + x(i+2));
end;
iqj = fix(fix(iq./j));
% R(K) = QUOTIENT, C = REMAINDER
mpcom_6(k) = fix(iqj + ir);
c = fix(iq - j.*iqj);
if( c<0 )
break;
end;
% MAIN LOOP FOR LARGE ABS(IY) CASE
k = fix(k + 1);
if( k>i2 )
gt(2)=1;
break;
end;
i = fix(i + 1);
end;
if(gt(2)==1)
break;
end;
else;
% LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT
while( true );
i = fix(i + 1);
c = fix(mpcom_1.*c);
if( i<=mpcom_2 )
c = fix(c + x(i+2));
end;
r1 = fix(fix(c./j));
if( r1<0 )
gt(1)=1;
break;
end;
if( r1~=0 )
break;
end;
end;
if(gt(1)==1)
continue;
end;
% ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT
re = fix(re + 1 - i);
mpcom_6(1) = fix(r1);
c = fix(mpcom_1.*(c-j.*r1));
kh = 2;
if( i<mpcom_2 )
kh = fix(1 + mpcom_2 - i);
for k = 2 : kh;
i = fix(i + 1);
c = fix(c + x(i+2));
mpcom_6(k) = fix(fix(c./j));
c = fix(mpcom_1.*(c-j.*mpcom_6(k)));
end; k = fix(kh+1);
if( c<0 )
gt(1)=1;
continue;
end;
kh = fix(kh + 1);
end;
for k = kh : i2;
mpcom_6(k) = fix(fix(c./j));
c = fix(mpcom_1.*(c-j.*mpcom_6(k)));
end; k = fix(i2+1);
if( c>=0 )
break;
end;
end;
% CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED
end;
gt(1)=0;
mpchk(1,4);
writef(mpcom_4,[' *** INTEGER OVERFLOW IN MPDIVI, B TOO LARGE ***' ' \n']);
%format (' *** INTEGER OVERFLOW IN MPDIVI, B TOO LARGE ***');
gt(3)=1;
break;
end;
end;
end;
% NORMALIZE AND ROUND RESULT
if(gt(3)==0)
[rs,re,z]=mpnzr(rs,re,z,0);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
end;
end;
mperr();
z(1) = 0;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine mpdivi
%DECK MPERR