[x,y,z]=mpmul(x,y,z);

function [x,y,z]=mpmul(x,y,z);
persistent c gt i i2 i2p j j1 re ri rs xi ; 

if isempty(i), i=0; end;
if isempty(i2), i2=0; end;
if isempty(i2p), i2p=0; end;
if isempty(j), j=0; end;
if isempty(j1), j1=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  MPMUL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DQDOTA and DQDOTI
%***LIBRARY   SLATEC
%***TYPE      ALL (MPMUL-A)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%  Multiplies X and Y, returning result in Z, for 'mp' X, Y and Z.
%  The simple o(mpcom_2**2) algorithm is used, with four guard digits and
%  R*-rounding. Advantage is taken of zero digits in X, but not in Y.
%  Asymptotically faster algorithms are known (see Knuth, VOL. 2),
%  but are difficult to implement in FORTRAN in an efficient and
%  machine-independent manner. In comments to other 'mp' routines,
%  M(mpcom_2) is the time to perform mpcom_2-digit 'mp' multiplication. Thus
%  M(mpcom_2) = o(mpcom_2**2) with the present version of MPMUL, but
%  M(mpcom_2) = o(mpcom_2.log(mpcom_2).log(log(mpcom_2))) is theoretically possible.
%
%  The arguments X(*), Y(*), 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, MPMLP, MPNZR
%***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  MPMUL
% 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,[]);
y_shape=size(y);y=reshape(y,1,[]);
z_shape=size(z);z=reshape(z,1,[]);
if isempty(rs), rs=0; end;
if isempty(re), re=0; end;
if isempty(xi), xi=0; end;
if isempty(c), c=0; end;
if isempty(ri), ri=0; end;
if isempty(gt), gt=zeros(1,2); end;
%***FIRST EXECUTABLE STATEMENT  MPMUL
mpchk(1,4);
i2 = fix(mpcom_2 + 4);
i2p = fix(i2 + 1);
% FORM SIGN OF PRODUCT
rs = fix(x(1).*y(1));
gt(:)=0;
while (1);
if( rs~=0 )
% FORM EXPONENT OF PRODUCT
re = fix(x(2) + y(2));
% CLEAR ACCUMULATOR
for i = 1 : i2;
mpcom_6(i) = 0;
end; i = fix(i2+1);
% PERFORM MULTIPLICATION
c = 8;
for i = 1 : mpcom_2;
xi = fix(x(i+2));
% FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST
if( xi~=0 )
[mpcom_6(sub2ind(size(mpcom_6),max(i+1,1)):end),y(sub2ind(size(y),max(3,1)):end),xi]=mpmlp(mpcom_6(sub2ind(size(mpcom_6),max(i+1,1)):end),y(sub2ind(size(y),max(3,1)):end),xi,min(mpcom_2,i2-i));
c = fix(c - 1);
if( c<=0 )
% CHECK FOR LEGAL BASE B DIGIT
if((xi<0) ||(xi>=mpcom_1) )
gt(2)=1;
break;
end;
% PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
% FASTER THAN DOING IT EVERY TIME.
for j = 1 : i2;
j1 = fix(i2p - j);
ri = fix(mpcom_6(j1) + c);
if( ri<0 )
gt(1)=1;
break;
end;
c = fix(fix(ri./mpcom_1));
mpcom_6(j1) = fix(ri - mpcom_1.*c);
end;
if(gt(1)==1)
break;
end;
if( c~=0 )
gt(2)=1;
break;
end;
c = 8;
end;
end;
end;
if(any(gt==1))
break;
end;
if( c~=8 )
if((xi<0) ||(xi>=mpcom_1) )
gt(2)=1;
break;
end;
c = 0;
for j = 1 : i2;
j1 = fix(i2p - j);
ri = fix(mpcom_6(j1) + c);
if( ri<0 )
gt(1)=1;
break;
end;
c = fix(fix(ri./mpcom_1));
mpcom_6(j1) = fix(ri - mpcom_1.*c);
end;
if(gt(1)==1)
break;
end;
if( c~=0 )
gt(2)=1;
break;
end;
end;
% NORMALIZE AND ROUND RESULT
[rs,re,z]=mpnzr(rs,re,z,0);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
else;
% SET RESULT TO ZERO
z(1) = 0;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
end;
break;
end;
if(gt(2)==0)
writef(mpcom_4,[' *** INTEGER OVERFLOW IN MPMUL, B TOO LARGE ***' ' \n']);
%format (' *** INTEGER OVERFLOW IN MPMUL, B TOO LARGE ***');
end;
if(gt(1)==0)
writef(mpcom_4,[' *** ILLEGAL BASE B DIGIT IN CALL TO MPMUL,',' POSSIBLE OVERWRITING PROBLEM ***' ' \n']);
%format (' *** ILLEGAL BASE B DIGIT IN CALL TO MPMUL,',' POSSIBLE OVERWRITING PROBLEM ***');
end;
mperr;
z(1) = 0;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine mpmul
%DECK MPMULI