Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

[pfq]=genHyper(a,b,z,lnpfq,ix,nsigfig);
function [pfq]=genHyper(a,b,z,lnpfq,ix,nsigfig);
% function [pfq]=pfq(a,b,z,lnpfq,ix,nsigfig)
% Description : A numerical evaluator for the generalized hypergeometric
%               function for complex arguments with large magnitudes
%               using a direct summation of the Gauss series.
%               pFq isdefined by (borrowed from Maple):
%   pFq = sum(z^k / k! * product(pochhammer(n[i], k), i=1..p) /
%         product(pochhammer(d[j], k), j=1..q), k=0..infinity )
%
% INPUTS:       a => array containing numerator parameters
%               b => array containing denominator parameters
%               z => complex argument (scalar)
%           lnpfq => (optional) set to 1 if desired result is the natural
%                    log of pfq (default is 0)
%              ix => (optional) maximum number of terms in a,b (see below)
%         nsigfig => number of desired significant figures (default=10)
%
% OUPUT:      pfq => result
%
% EXAMPLES:     a=[1+i,1]; b=[2-i,3,3]; z=1.5;
%               >> genHyper(a,b,z)
%               ans =
%                          1.02992154295955 +     0.106416425916656i
%               or with more precision,
%               >> genHyper(a,b,z,0,0,15)
%               ans =
%                          1.02992154295896 +     0.106416425915575i
%               using the log option,
%               >> genHyper(a,b,z,1,0,15)
%               ans =
%                        0.0347923403326305 +     0.102959427435454i
%               >> exp(ans)
%               ans =
%                          1.02992154295896 +     0.106416425915575i
%
%
% Translated from the original fortran using f2matlab.m
%  by Ben E. Barrowes - barrowes@alum.mit.edu, 7/04.
%  

%% Original fortran documentation
%     ACPAPFQ.  A NUMERICAL EVALUATOR FOR THE GENERALIZED HYPERGEOMETRIC
%
%     1  SERIES.  W.F. PERGER, A. BHALLA, M. NARDIN.
%
%     REF. IN COMP. PHYS. COMMUN. 77 (1993) 249
%
%     ****************************************************************
%     *                                                              *
%     *    SOLUTION TO THE GENERALIZED HYPERGEOMETRIC FUNCTION       *
%     *                                                              *
%     *                           by                                 *
%     *                                                              *
%     *                      W. F. PERGER,                           *
%     *                                                              *
%     *              MARK NARDIN  and ATUL BHALLA                    *
%     *                                                              *
%     *                                                              *
%     *            Electrical Engineering Department                 *
%     *            Michigan Technological University                 *
%     *                  1400 Townsend Drive                         *
%     *                Houghton, MI  49931-1295   USA                *
%     *                     Copyright 1993                           *
%     *                                                              *
%     *               e-mail address: wfp@mtu.edu                    *
%     *                                                              *
%     *  Description : A numerical evaluator for the generalized     *
%     *    hypergeometric function for complex arguments with large  *
%     *    magnitudes using a direct summation of the Gauss series.  *
%     *    The method used allows an accuracy of up to thirteen      *
%     *    decimal places through the use of large integer arrays    *
%     *    and a single final division.                              *
%     *    (original subroutines for the confluent hypergeometric    *
%     *    written by Mark Nardin, 1989; modifications made to cal-  *
%     *    culate the generalized hypergeometric function were       *
%     *    written by W.F. Perger and A. Bhalla, June, 1990)         *
%     *                                                              *
%     *  The evaluation of the pFq series is accomplished by a func- *
%     *  ion call to PFQ, which is a double precision complex func-  *
%     *  tion.  The required input is:                               *
%     *  1. Double precision complex arrays A and B.  These are the  *
%     *     arrays containing the parameters in the numerator and de-*
%     *     nominator, respectively.                                 *
%     *  2. Integers IP and IQ.  These integers indicate the number  *
%     *     of numerator and denominator terms, respectively (these  *
%     *     are p and q in the pFq function).                        *
%     *  3. Double precision complex argument Z.                     *
%     *  4. Integer LNPFQ.  This integer should be set to '1' if the *
%     *     result from PFQ is to be returned as the natural logaritm*
%     *     of the series, or '0' if not.  The user can generally set*
%     *     LNPFQ = '0' and change it if required.                   *
%     *  5. Integer IX.  This integer should be set to '0' if the    *
%     *     user desires the program PFQ to estimate the number of   *
%     *     array terms (in A and B) to be used, or an integer       *
%     *     greater than zero specifying the number of integer pos-  *
%     *     itions to be used.  This input parameter is escpecially  *
%     *     useful as a means to check the results of a given run.   *
%     *     Specificially, if the user obtains a result for a given  *
%     *     set of parameters, then changes IX and re-runs the eval- *
%     *     uator, and if the number of array positions was insuffi- *
%     *     cient, then the two results will likely differ.  The rec-*
%     *     commended would be to generally set IX = '0' and then set*
%     *     it to 100 or so for a second run.  Note that the LENGTH  *
%     *     parameter currently sets the upper limit on IX to 777,   *
%     *     but that can easily be changed (it is a single PARAMETER *
%     *     statement) and the program recompiled.                   *
%     *  6. Integer NSIGFIG.  This integer specifies the requested   *
%     *     number of significant figures in the final result.  If   *
%     *     the user attempts to request more than the number of bits*
%     *     in the mantissa allows, the program will abort with an   *
%     *     appropriate error message.  The recommended value is 10. *
%     *                                                              *
%     *     Note: The variable NOUT is the file to which error mess- *
%     *           ages are written (default is 6).  This can be      *
%     *           changed in the FUNCTION PFQ to accomodate re-      *
%     *           of output to another file                          *
%     *                                                              *
%     *  Subprograms called: HYPER.                                  *
%     *                                                              *
%     ****************************************************************
%
%
%
%

if nargin<6
 nsigfig=10;
elseif isempty(nsigfig)
 nsigfig=10;
end
if nargin<5
 ix=0;
elseif isempty(ix)
 ix=0;
end
if nargin<4
 lnpfq=0;
elseif isempty(lnpfq)
 lnpfq=0;
end
ip=length(a);
iq=length(b);

global  zero   half   one   two   ten   eps;
[zero , half , one , two , ten , eps]=deal(0.0d0,0.5d0,1.0d0,2.0d0,10.0d0,1.0d-10);
global  nout;
%
%
%
%
a1=zeros(2,1);b1=zeros(1,1);gam1=0;gam2=0;gam3=0;gam4=0;gam5=0;gam6=0;gam7=0;hyper1=0;hyper2=0;z1=0;
argi=0;argr=0;diff=0;dnum=0;precis=0;
%
%
i=0;
%
%
%
nout=6;
if ((lnpfq~=0) & (lnpfq~=1)) ;
 ' error in input arguments: lnpfq ~= 0 or 1',
 error('stop encountered in original fortran code');
end;
if ((ip>iq) & (abs(z)>one)) ;
 ip , iq , abs(z),
 %format [,1x,'ip=',1i2,3x,'iq=',1i2,3x,'and abs(z)=',1e12.5,2x,./,' which is greater than one--series does',' not converge');
 error('stop encountered in original fortran code');
end;
if (ip==2 & iq==1 & abs(z)>0.9) ;
 if (lnpfq~=1) ;
  %
  %      Check to see if the Gamma function arguments are o.k.; if not,
  %
  %      then the series will have to be used.
  %
  %
  %
  %      PRECIS - MACHINE PRECISION
  %
  %
  precis=one;
  precis=precis./two;
  dnum=precis+one;
  while (dnum>one);
   precis=precis./two;
   dnum=precis+one;
  end;
  precis=two.*precis;
  for i=1 : 6;
   if (i==1) ;
    argi=imag(b(1));
    argr=real(b(1));
   elseif (i==2);
    argi=imag(b(1)-a(1)-a(2));
    argr=real(b(1)-a(1)-a(2));
   elseif (i==3);
    argi=imag(b(1)-a(1));
    argr=real(b(1)-a(1));
   elseif (i==4);
    argi=imag(a(1)+a(2)-b(1));
    argr=real(a(1)+a(2)-b(1));
   elseif (i==5);
    argi=imag(a(1));
    argr=real(a(1));
   elseif (i==6);
    argi=imag(a(2));
    argr=real(a(2));
   end;
   %
   %       CASES WHERE THE ARGUMENT IS REAL
   %
   %
   if (argi==0.0) ;
    %
    %        CASES WHERE THE ARGUMENT IS REAL AND NEGATIVE
    %
    %
    if (argr<=0.0) ;
     %
     %         USE THE SERIES EXPANSION IF THE ARGUMENT IS TOO NEAR A POLE
     %
     %
     diff=abs(real(round(argr))-argr);
     if (diff<=two.*precis) ;
      pfq=hyper(a,b,ip,iq,z,lnpfq,ix,nsigfig);
      return;
     end;
    end;
   end;
  end;
  gam1=cgamma(b(1),lnpfq);
  gam2=cgamma(b(1)-a(1)-a(2),lnpfq);
  gam3=cgamma(b(1)-a(1),lnpfq);
  gam4=cgamma(b(1)-a(2),lnpfq);
  gam5=cgamma(a(1)+a(2)-b(1),lnpfq);
  gam6=cgamma(a(1),lnpfq);
  gam7=cgamma(a(2),lnpfq);
  a1(1)=a(1);
  a1(2)=a(2);
  b1(1)=a(1)+a(2)-b(1)+one;
  z1=one-z;
  hyper1=hyper(a1,b1,ip,iq,z1,lnpfq,ix,nsigfig);
  a1(1)=b(1)-a(1);
  a1(2)=b(1)-a(2);
  b1(1)=b(1)-a(1)-a(2)+one;
  hyper2=hyper(a1,b1,ip,iq,z1,lnpfq,ix,nsigfig);
  pfq=gam1.*gam2.*hyper1./(gam3.*gam4)+(one-z).^(b(1)-a(1)-a(2)).*gam1.*gam5.*hyper2./(gam6.*gam7);
  return;
 end;
end;
pfq=hyper(a,b,ip,iq,z,lnpfq,ix,nsigfig);
return;


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION BITS                              *
%     *                                                              *
%     *                                                              *
%     *  Description : Determines the number of significant figures  *
%     *    of machine precision to arrive at the size of the array   *
%     *    the numbers must be stored in to get the accuracy of the  *
%     *    solution.                                                 *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [bits]=bits();
%
bit2=0;
%
%
%
bit=1.0;
nnz=0;
nnz=nnz+1;
bit2=bit.*2.0;
bit=bit2+1.0;
while ((bit-bit2)~=0.0);
 nnz=nnz+1;
 bit2=bit.*2.0;
 bit=bit2+1.0;
end;
bits=nnz-3;


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION HYPER                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Function that sums the Gauss series.          *
%     *                                                              *
%     *  Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL,   *
%     *                      IPREMAX.                                *
%     *                                                              *
%     ****************************************************************
%
function [hyper]=hyper(a,b,ip,iq,z,lnpfq,ix,nsigfig);
%
%
% PARAMETER definitions
%
sumr=[];sumi=[];denomr=[];denomi=[];final=[];l=[];rmax=[];ibit=[];temp=[];cr=[];i1=[];ci=[];qr1=[];qi1=[];wk1=[];wk2=[];wk3=[];wk4=[];wk5=[];wk6=[];cr2=[];ci2=[];qr2=[];qi2=[];foo1=[];cnt=[];foo2=[];sigfig=[];numr=[];numi=[];ar=[];ai=[];ar2=[];ai2=[];xr=[];xi=[];xr2=[];xi2=[];bar1=[];bar2=[];
length=0;
length=777;
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
%
accy=0;ai=zeros(10,1);ai2=zeros(10,1);ar=zeros(10,1);ar2=zeros(10,1);ci=zeros(10,1);ci2=zeros(10,1);cnt=0;cr=zeros(10,1);cr2=zeros(10,1);creal=0;denomi=zeros(length+2,1);denomr=zeros(length+2,1);dum1=0;dum2=0;expon=0;log2=0;mx1=0;mx2=0;numi=zeros(length+2,1);numr=zeros(length+2,1);qi1=zeros(length+2,1);qi2=zeros(length+2,1);qr1=zeros(length+2,1);qr2=zeros(length+2,1);ri10=0;rmax=0;rr10=0;sigfig=0;sumi=zeros(length+2,1);sumr=zeros(length+2,1);wk1=zeros(length+2,1);wk2=zeros(length+2,1);wk3=zeros(length+2,1);wk4=zeros(length+2,1);wk5=zeros(length+2,1);wk6=zeros(length+2,1);x=0;xi=0;xi2=0;xl=0;xr=0;xr2=0;
%
cdum1=0;cdum2=0;final=0;oldtemp=0;temp=0;temp1=0;
%
%
%
i=0;i1=0;ibit=0;icount=0;ii10=0;ir10=0;ixcnt=0;l=0;lmax=0;nmach=0;rexp=0;
%
%
%
goon1=0;
foo1=zeros(length+2,1);foo2=zeros(length+2,1);bar1=zeros(length+2,1);bar2=zeros(length+2,1);
%
%
zero=0.0d0;
log2=log10(two);
ibit=fix(bits);
rmax=two.^(fix(ibit./2));
sigfig=two.^(fix(ibit./4));
%
for i1=1 : ip;
 ar2(i1)=real(a(i1)).*sigfig;
 ar(i1)=fix(ar2(i1));
 ar2(i1)=round((ar2(i1)-ar(i1)).*rmax);
 ai2(i1)=imag(a(i1)).*sigfig;
 ai(i1)=fix(ai2(i1));
 ai2(i1)=round((ai2(i1)-ai(i1)).*rmax);
end;
for i1=1 : iq;
 cr2(i1)=real(b(i1)).*sigfig;
 cr(i1)=fix(cr2(i1));
 cr2(i1)=round((cr2(i1)-cr(i1)).*rmax);
 ci2(i1)=imag(b(i1)).*sigfig;
 ci(i1)=fix(ci2(i1));
 ci2(i1)=round((ci2(i1)-ci(i1)).*rmax);
end;
xr2=real(z).*sigfig;
xr=fix(xr2);
xr2=round((xr2-xr).*rmax);
xi2=imag(z).*sigfig;
xi=fix(xi2);
xi2=round((xi2-xi).*rmax);
%
%     WARN THE USER THAT THE INPUT VALUE WAS SO CLOSE TO ZERO THAT IT
%     WAS SET EQUAL TO ZERO.
%
for i1=1 : ip;
 if ((real(a(i1))~=0.0) & (ar(i1)==0.0) & (ar2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - real part of a(',1i2,') was set to zero');
 if ((imag(a(i1))~=0.0) & (ai(i1)==0.0) & (ai2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - imag part of a(',1i2,') was set to zero');
end;
for i1=1 : iq;
 if ((real(b(i1))~=0.0) & (cr(i1)==0.0) & (cr2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - real part of b(',1i2,') was set to zero');
 if ((imag(b(i1))~=0.0) & (ci(i1)==0.0) & (ci2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - imag part of b(',1i2,') was set to zero');
end;
if ((real(z)~=0.0) & (xr==0.0) & (xr2==0.0)) ;
 ' warning - real part of z was set to zero',
 z=complex(0.0,imag(z));
end;
if ((imag(z)~=0.0) & (xi==0.0) & (xi2==0.0)) ;
 ' warning - imag part of z was set to zero',
 z=complex(real(z),0.0);
end;
%
%
%     SCREENING OF NUMERATOR ARGUMENTS FOR NEGATIVE INTEGERS OR ZERO.
%     ICOUNT WILL FORCE THE SERIES TO TERMINATE CORRECTLY.
%
nmach=fix(log10(two.^fix(bits)));
icount=-1;
for i1=1 : ip;
 if ((ar2(i1)==0.0) & (ar(i1)==0.0) & (ai2(i1)==0.0) &(ai(i1)==0.0)) ;
  hyper=complex(one,0.0);
  return;
 end;
 if ((ai(i1)==0.0) & (ai2(i1)==0.0) & (real(a(i1))<0.0));
  if (abs(real(a(i1))-real(round(real(a(i1)))))<ten.^(-nmach)) ;
   if (icount~=-1) ;
    icount=min([icount,-round(real(a(i1)))]);
   else;
    icount=-round(real(a(i1)));
   end;
  end;
 end;
end;
%
%     SCREENING OF DENOMINATOR ARGUMENTS FOR ZEROES OR NEGATIVE INTEGERS
%     .
%
for i1=1 : iq;
 if ((cr(i1)==0.0) & (cr2(i1)==0.0) & (ci(i1)==0.0) &(ci2(i1)==0.0)) ;
  i1,
  %format (1x,'error - argument b(',1i2,') was equal to zero');
  error('stop encountered in original fortran code');
 end;
 if ((ci(i1)==0.0) & (ci2(i1)==0.0) & (real(b(i1))<0.0));
  if ((abs(real(b(i1))-real(round(real(b(i1)))))<ten.^(-nmach)) &(icount>=-round(real(b(i1))) | icount==-1)) ;
   i1,
   %format (1x,'error - argument b(',1i2,') was a negative',' integer');
   error('stop encountered in original fortran code');
  end;
 end;
end;
%
nmach=fix(log10(two.^ibit));
nsigfig=min([nsigfig,fix(log10(two.^ibit))]);
accy=ten.^(-nsigfig);
l=ipremax(a,b,ip,iq,z);
if (l~=1) ;
 %
 %     First, estimate the exponent of the maximum term in the pFq series
 %     .
 %
 expon=0.0;
 xl=real(l);
 for i=1 : ip;
  expon=expon+real(factor(a(i)+xl-one))-real(factor(a(i)-one));
 end;
 for i=1 : iq;
  expon=expon-real(factor(b(i)+xl-one))+real(factor(b(i)-one));
 end;
 expon=expon+xl.*real(log(z))-real(factor(complex(xl,0.0)));
 lmax=fix(log10(exp(one)).*expon);
 l=lmax;
 %
 %     Now, estimate the exponent of where the pFq series will terminate.
 %
 temp1=complex(one,0.0);
 creal=one;
 for i1=1 : ip;
  temp1=temp1.*complex(ar(i1),ai(i1))./sigfig;
 end;
 for i1=1 : iq;
  temp1=temp1./(complex(cr(i1),ci(i1))./sigfig);
  creal=creal.*cr(i1);
 end;
 temp1=temp1.*complex(xr,xi);
 %
 %     Triple it to make sure.
 %
 l=3.*l;
 %
 %     Divide the number of significant figures necessary by the number
 %     of
 %     digits available per array position.
 %
 %
 l=fix((2.*l+nsigfig)./nmach)+2;
end;
%
%     Make sure there are at least 5 array positions used.
%
l=max([l,5]);
l=max([l,ix]);
%     write (6,*) ' Estimated value of L=',L
if ((l<0) | (l>length)) ;
 length,
 %format (1x,['error in fn hyper: l must be < '],1i4);
 error('stop encountered in original fortran code');
end;
if (nsigfig>nmach) ;
 nmach,
 %format (1x,' warning--the number of significant figures requ','ested',./,'is greater than the machine precision--','final answer',./,'will be accurate to only',i3,' digits');
end;
%
sumr(-1+2)=one;
sumi(-1+2)=one;
numr(-1+2)=one;
numi(-1+2)=one;
denomr(-1+2)=one;
denomi(-1+2)=one;
for i=0 : l+1;
 sumr(i+2)=0.0;
 sumi(i+2)=0.0;
 numr(i+2)=0.0;
 numi(i+2)=0.0;
 denomr(i+2)=0.0;
 denomi(i+2)=0.0;
end;
sumr(1+2)=one;
numr(1+2)=one;
denomr(1+2)=one;
cnt=sigfig;
temp=complex(0.0,0.0);
oldtemp=temp;
ixcnt=0;
rexp=fix(ibit./2);
x=rexp.*(sumr(l+1+2)-2);
rr10=x.*log2;
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(sumi(l+1+2)-2);
ri10=x.*log2;
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(sumr(1+2).*rmax.*rmax+sumr(2+2).*rmax+sumr(3+2)).*sign(sumr(-1+2)));
dum2=(abs(sumi(1+2).*rmax.*rmax+sumi(2+2).*rmax+sumi(3+2)).*sign(sumi(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum1=complex(dum1,dum2);
x=rexp.*(denomr(l+1+2)-2);
rr10=x.*log2;
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(denomi(l+1+2)-2);
ri10=x.*log2;
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(denomr(1+2).*rmax.*rmax+denomr(2+2).*rmax+denomr(3+2)).*sign(denomr(-1+2)));
dum2=(abs(denomi(1+2).*rmax.*rmax+denomi(2+2).*rmax+denomi(3+2)).*sign(denomi(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum2=complex(dum1,dum2);
temp=cdum1./cdum2;
%
%     130 IF (IP .GT. 0) THEN
goon1=1;
while (goon1==1);
 goon1=0;
 if (ip<0) ;
  if (sumr(1+2)<half) ;
   mx1=sumi(l+1+2);
  elseif (sumi(1+2)<half);
   mx1=sumr(l+1+2);
  else;
   mx1=max([sumr(l+1+2),sumi(l+1+2)]);
  end;
  if (numr(1+2)<half) ;
   mx2=numi(l+1+2);
  elseif (numi(1+2)<half);
   mx2=numr(l+1+2);
  else;
   mx2=max([numr(l+1+2),numi(l+1+2)]);
  end;
  if (mx1-mx2>2.0) ;
   if (creal>=0.0) ;
    %        write (6,*) ' cdabs(temp1/cnt)=',cdabs(temp1/cnt)
    %
    if (abs(temp1./cnt)<=one) ;
     [sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit);
     hyper=final;
     return;
    end;
   end;
  end;
 else;
  [sumr,sumi,denomr,denomi,temp,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,temp,l,lnpfq,rmax,ibit);
  %
  %      First, estimate the exponent of the maximum term in the pFq
  %      series.
  %
  expon=0.0;
  xl=real(ixcnt);
  for i=1 : ip;
   expon=expon+real(factor(a(i)+xl-one))-real(factor(a(i)-one));
  end;
  for i=1 : iq;
   expon=expon-real(factor(b(i)+xl-one))+real(factor(b(i)-one));
  end;
  expon=expon+xl.*real(log(z))-real(factor(complex(xl,0.0)));
  lmax=fix(log10(exp(one)).*expon);
  if (abs(oldtemp-temp)<abs(temp.*accy)) ;
   [sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit);
   hyper=final;
   return;
  end;
  oldtemp=temp;
 end;
 if (ixcnt~=icount) ;
  ixcnt=ixcnt+1;
  for i1=1 : iq;
   %
   %      TAKE THE CURRENT SUM AND MULTIPLY BY THE DENOMINATOR OF THE NEXT
   %
   %      TERM, FOR BOTH THE MOST SIGNIFICANT HALF (CR,CI) AND THE LEAST
   %
   %      SIGNIFICANT HALF (CR2,CI2).
   %
   %
   [sumr,sumi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(sumr,sumi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   [sumr,sumi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(sumr,sumi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   qr2(l+1+2)=qr2(l+1+2)-1;
   qi2(l+1+2)=qi2(l+1+2)-1;
   %
   %      STORE THIS TEMPORARILY IN THE SUM ARRAYS.
   %
   %
   [qr1,qi1,qr2,qi2,sumr,sumi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,sumr,sumi,wk1,l,rmax);
  end;
  %
  %
  %     MULTIPLY BY THE FACTORIAL TERM.
  %
  foo1=sumr;
  foo2=sumr;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  sumr=foo2;
  foo1=sumi;
  foo2=sumi;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  sumi=foo2;
  %
  %     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT.
  %
  for i1=1 : ip-iq;
   foo1=sumr;
   foo2=sumr;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   sumr=foo2;
   foo1=sumi;
   foo2=sumi;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   sumi=foo2;
  end;
  for i1=1 : iq;
   %
   %      UPDATE THE DENOMINATOR.
   %
   %
   [denomr,denomi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(denomr,denomi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   [denomr,denomi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(denomr,denomi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   qr2(l+1+2)=qr2(l+1+2)-1;
   qi2(l+1+2)=qi2(l+1+2)-1;
   [qr1,qi1,qr2,qi2,denomr,denomi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,denomr,denomi,wk1,l,rmax);
  end;
  %
  %
  %     MULTIPLY BY THE FACTORIAL TERM.
  %
  foo1=denomr;
  foo2=denomr;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  denomr=foo2;
  foo1=denomi;
  foo2=denomi;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  denomi=foo2;
  %
  %     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT.
  %
  for i1=1 : ip-iq;
   foo1=denomr;
   foo2=denomr;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   denomr=foo2;
   foo1=denomi;
   foo2=denomi;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   denomi=foo2;
  end;
  %
  %     FORM THE NEXT NUMERATOR TERM BY MULTIPLYING THE CURRENT
  %     NUMERATOR TERM (AN ARRAY) WITH THE A ARGUMENT (A SCALAR).
  %
  for i1=1 : ip;
   [numr,numi,ar(i1),ai(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,ar(i1),ai(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   [numr,numi,ar2(i1),ai2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,ar2(i1),ai2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   qr2(l+1+2)=qr2(l+1+2)-1;
   qi2(l+1+2)=qi2(l+1+2)-1;
   [qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax);
  end;
  %
  %     FINISH THE NEW NUMERATOR TERM BY MULTIPLYING BY THE Z ARGUMENT.
  %
  [numr,numi,xr,xi,qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,xr,xi,qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
  [numr,numi,xr2,xi2,qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,xr2,xi2,qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
  qr2(l+1+2)=qr2(l+1+2)-1;
  qi2(l+1+2)=qi2(l+1+2)-1;
  [qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax);
  %
  %     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT.
  %
  for i1=1 : iq-ip;
   foo1=numr;
   foo2=numr;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   numr=foo2;
   foo1=numi;
   foo2=numi;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   numi=foo2;
  end;
  %
  %     FINALLY, ADD THE NEW NUMERATOR TERM WITH THE CURRENT RUNNING
  %     SUM OF THE NUMERATOR AND STORE THE NEW RUNNING SUM IN SUMR, SUMI.
  %
  foo1=sumr;
  foo2=sumr;
  bar1=sumi;
  bar2=sumi;
  [foo1,bar1,numr,numi,foo2,bar2,wk1,l,rmax]=cmpadd(foo1,bar1,numr,numi,foo2,bar2,wk1,l,rmax);
  sumi=bar2;
  sumr=foo2;

  %
  %     BECAUSE SIGFIG REPRESENTS "ONE" ON THE NEW SCALE, ADD SIGFIG
  %     TO THE CURRENT COUNT AND, CONSEQUENTLY, TO THE IP ARGUMENTS
  %     IN THE NUMERATOR AND THE IQ ARGUMENTS IN THE DENOMINATOR.
  %
  cnt=cnt+sigfig;
  for i1=1 : ip;
   ar(i1)=ar(i1)+sigfig;
  end;
  for i1=1 : iq;
   cr(i1)=cr(i1)+sigfig;
  end;
  goon1=1;
 end;
end;
[sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit);
%     write (6,*) 'Number of terms=',ixcnt
hyper=final;
return;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARADD                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Accepts two arrays of numbers and returns     *
%     *    the sum of the array.  Each array is holding the value    *
%     *    of one number in the series.  The parameter L is the      *
%     *    size of the array representing the number and RMAX is     *
%     *    the actual number of digits needed to give the numbers    *
%     *    the desired accuracy.                                     *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c,z,l,rmax]=aradd(a,b,c,z,l,rmax);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
ediff=0;i=0;j=0;
%
%
for i=0 : l+1;
 z(i+2)=0.0;
end;
ediff=round(a(l+1+2)-b(l+1+2));
if (abs(a(1+2))<half | ediff<=-l) ;
 for i=-1 : l+1;
  c(i+2)=b(i+2);
 end;
 if (c(1+2)<half) ;
  c(-1+2)=one;
  c(l+1+2)=0.0;
 end;
 return;
else;
 if (abs(b(1+2))<half | ediff>=l) ;
  for i=-1 : l+1;
   c(i+2)=a(i+2);
  end;
  if (c(1+2)<half) ;
   c(-1+2)=one;
   c(l+1+2)=0.0;
  end;
  return;
 else;
  z(-1+2)=a(-1+2);
  goon300=1;
  goon190=1;
  if (abs(a(-1+2)-b(-1+2))>=half) ;
   goon300=0;
   if (ediff>0) ;
    z(l+1+2)=a(l+1+2);
   elseif (ediff<0);
    z(l+1+2)=b(l+1+2);
    z(-1+2)=b(-1+2);
    goon190=0;
   else;
    for i=1 : l;
     if (a(i+2)>b(i+2)) ;
      z(l+1+2)=a(l+1+2);
      break;
     end;
     if (a(i+2)<b(i+2)) ;
      z(l+1+2)=b(l+1+2);
      z(-1+2)=b(-1+2);
      goon190=0;
     end;
    end;
   end;
   %
  elseif (ediff>0);
   z(l+1+2)=a(l+1+2);
   for i=l : -1: 1+ediff ;
    z(i+2)=a(i+2)+b(i-ediff+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   for i=ediff : -1: 1 ;
    z(i+2)=a(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   if (z(0+2)>half) ;
    for i=l : -1: 1 ;
     z(i+2)=z(i-1+2);
    end;
    z(l+1+2)=z(l+1+2)+1;
    z(0+2)=0.0;
   end;
  elseif (ediff<0);
   z(l+1+2)=b(l+1+2);
   for i=l : -1: 1-ediff ;
    z(i+2)=a(i+ediff+2)+b(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   for i=0-ediff : -1: 1 ;
    z(i+2)=b(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   if (z(0+2)>half) ;
    for i=l : -1: 1 ;
     z(i+2)=z(i-1+2);
    end;
    z(l+1+2)=z(l+1+2)+one;
    z(0+2)=0.0;
   end;
  else;
   z(l+1+2)=a(l+1+2);
   for i=l : -1: 1 ;
    z(i+2)=a(i+2)+b(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   if (z(0+2)>half) ;
    for i=l : -1: 1 ;
     z(i+2)=z(i-1+2);
    end;
    z(l+1+2)=z(l+1+2)+one;
    z(0+2)=0.0;
   end;
  end;
  if (goon300==1) ;
   i=i; %here is the line that had a +1 taken from it.
   while (z(i+2)<half & i<l+1);
    i=i+1;
   end;
   if (i==l+1) ;
    z(-1+2)=one;
    z(l+1+2)=0.0;
    for i=-1 : l+1;
     c(i+2)=z(i+2);
    end;
    if (c(1+2)<half) ;
     c(-1+2)=one;
     c(l+1+2)=0.0;
    end;
    return;
   end;
   for j=1 : l+1-i;
    z(j+2)=z(j+i-1+2);
   end;
   for j=l+2-i : l;
    z(j+2)=0.0;
   end;
   z(l+1+2)=z(l+1+2)-i+1;
   for i=-1 : l+1;
    c(i+2)=z(i+2);
   end;
   if (c(1+2)<half) ;
    c(-1+2)=one;
    c(l+1+2)=0.0;
   end;
   return;
  end;
  %
  if (goon190==1) ;
   if (ediff>0) ;
    for i=l : -1: 1+ediff ;
     z(i+2)=a(i+2)-b(i-ediff+2)+z(i+2);
     if (z(i+2)<0.0) ;
      z(i+2)=z(i+2)+rmax;
      z(i-1+2)=-one;
     end;
    end;
    for i=ediff : -1: 1 ;
     z(i+2)=a(i+2)+z(i+2);
     if (z(i+2)<0.0) ;
      z(i+2)=z(i+2)+rmax;
      z(i-1+2)=-one;
     end;
    end;
   else;
    for i=l : -1: 1 ;
     z(i+2)=a(i+2)-b(i+2)+z(i+2);
     if (z(i+2)<0.0) ;
      z(i+2)=z(i+2)+rmax;
      z(i-1+2)=-one;
     end;
    end;
   end;
   if (z(1+2)>half) ;
    for i=-1 : l+1;
     c(i+2)=z(i+2);
    end;
    if (c(1+2)<half) ;
     c(-1+2)=one;
     c(l+1+2)=0.0;
    end;
    return;
   end;
   i=1;
   i=i+1;
   while (z(i+2)<half & i<l+1);
    i=i+1;
   end;
   if (i==l+1) ;
    z(-1+2)=one;
    z(l+1+2)=0.0;
    for i=-1 : l+1;
     c(i+2)=z(i+2);
    end;
    if (c(1+2)<half) ;
     c(-1+2)=one;
     c(l+1+2)=0.0;
    end;
    return;
   end;
   for j=1 : l+1-i;
    z(j+2)=z(j+i-1+2);
   end;
   for j=l+2-i : l;
    z(j+2)=0.0;
   end;
   z(l+1+2)=z(l+1+2)-i+1;
   for i=-1 : l+1;
    c(i+2)=z(i+2);
   end;
   if (c(1+2)<half) ;
    c(-1+2)=one;
    c(l+1+2)=0.0;
   end;
   return;
  end;
 end;
 %
 if (ediff<0) ;
  for i=l : -1: 1-ediff ;
   z(i+2)=b(i+2)-a(i+ediff+2)+z(i+2);
   if (z(i+2)<0.0) ;
    z(i+2)=z(i+2)+rmax;
    z(i-1+2)=-one;
   end;
  end;
  for i=0-ediff : -1: 1 ;
   z(i+2)=b(i+2)+z(i+2);
   if (z(i+2)<0.0) ;
    z(i+2)=z(i+2)+rmax;
    z(i-1+2)=-one;
   end;
  end;
 else;
  for i=l : -1: 1 ;
   z(i+2)=b(i+2)-a(i+2)+z(i+2);
   if (z(i+2)<0.0) ;
    z(i+2)=z(i+2)+rmax;
    z(i-1+2)=-one;
   end;
  end;
 end;
end;
%
if (z(1+2)>half) ;
 for i=-1 : l+1;
  c(i+2)=z(i+2);
 end;
 if (c(1+2)<half) ;
  c(-1+2)=one;
  c(l+1+2)=0.0;
 end;
 return;
end;
i=1;
i=i+1;
while (z(i+2)<half & i<l+1);
 i=i+1;
end;
if (i==l+1) ;
 z(-1+2)=one;
 z(l+1+2)=0.0;
 for i=-1 : l+1;
  c(i+2)=z(i+2);
 end;
 if (c(1+2)<half) ;
  c(-1+2)=one;
  c(l+1+2)=0.0;
 end;
 return;
end;
for j=1 : l+1-i;
 z(j+2)=z(j+i-1+2);
end;
for j=l+2-i : l;
 z(j+2)=0.0;
end;
z(l+1+2)=z(l+1+2)-i+1;
for i=-1 : l+1;
 c(i+2)=z(i+2);
end;
if (c(1+2)<half) ;
 c(-1+2)=one;
 c(l+1+2)=0.0;
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARSUB                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Accepts two arrays and subtracts each element *
%     *    in the second array from the element in the first array   *
%     *    and returns the solution.  The parameters L and RMAX are  *
%     *    the size of the array and the number of digits needed for *
%     *    the accuracy, respectively.                               *
%     *                                                              *
%     *  Subprograms called: ARADD                                   *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c,wk1,wk2,l,rmax]=arsub(a,b,c,wk1,wk2,l,rmax);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
i=0;
%
%
for i=-1 : l+1;
 wk2(i+2)=b(i+2);
end;
wk2(-1+2)=(-one).*wk2(-1+2);
[a,wk2,c,wk1,l,rmax]=aradd(a,wk2,c,wk1,l,rmax);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARMULT                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Accepts two arrays and returns the product.   *
%     *    L and RMAX are the size of the arrays and the number of   *
%     *    digits needed to represent the numbers with the required  *
%     *    accuracy.                                                 *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c,z,l,rmax]=armult(a,b,c,z,l,rmax);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
b2=0;carry=0;
i=0;
%
%
z(-1+2)=(abs(one).*sign(b)).*a(-1+2);
b2=abs(b);
z(l+1+2)=a(l+1+2);
for i=0 : l;
 z(i+2)=0.0;
end;
if (b2<=eps | a(1+2)<=eps) ;
 z(-1+2)=one;
 z(l+1+2)=0.0;
else;
 for i=l : -1: 1 ;
  z(i+2)=a(i+2).*b2+z(i+2);
  if (z(i+2)>=rmax) ;
   carry=fix(z(i+2)./rmax);
   z(i+2)=z(i+2)-carry.*rmax;
   z(i-1+2)=carry;
  end;
 end;
 if (z(0+2)>=half) ;
  for i=l : -1: 1 ;
   z(i+2)=z(i-1+2);
  end;
  z(l+1+2)=z(l+1+2)+one;
  if (z(1+2)>=rmax) ;
   for i=l : -1: 1 ;
    z(i+2)=z(i-1+2);
   end;
   carry=fix(z(1+2)./rmax);
   z(2+2)=z(2+2)-carry.*rmax;
   z(1+2)=carry;
   z(l+1+2)=z(l+1+2)+one;
  end;
  z(0+2)=0.0;
 end;
end;
for i=-1 : l+1;
 c(i+2)=z(i+2);
end;
if (c(1+2)<half) ;
 c(-1+2)=one;
 c(l+1+2)=0.0;
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CMPADD                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes two arrays representing one real and    *
%     *    one imaginary part, and adds two arrays representing      *
%     *    another complex number and returns two array holding the  *
%     *    complex sum.                                              *
%     *              (CR,CI) = (AR+BR, AI+BI)                        *
%     *                                                              *
%     *  Subprograms called: ARADD                                   *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,cr,ci,wk1,l,rmax]=cmpadd(ar,ai,br,bi,cr,ci,wk1,l,rmax);
%
%
%
%
%
%
[ar,br,cr,wk1,l,rmax]=aradd(ar,br,cr,wk1,l,rmax);
[ai,bi,ci,wk1,l,rmax]=aradd(ai,bi,ci,wk1,l,rmax);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CMPSUB                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes two arrays representing one real and    *
%     *    one imaginary part, and subtracts two arrays representing *
%     *    another complex number and returns two array holding the  *
%     *    complex sum.                                              *
%     *              (CR,CI) = (AR+BR, AI+BI)                        *
%     *                                                              *
%     *  Subprograms called: ARADD                                   *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,cr,ci,wk1,wk2,l,rmax]=cmpsub(ar,ai,br,bi,cr,ci,wk1,wk2,l,rmax);
%
%
%
%
%
%
[ar,br,cr,wk1,wk2,l,rmax]=arsub(ar,br,cr,wk1,wk2,l,rmax);
[ai,bi,ci,wk1,wk2,l,rmax]=arsub(ai,bi,ci,wk1,wk2,l,rmax);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CMPMUL                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes two arrays representing one real and    *
%     *    one imaginary part, and multiplies it with two arrays     *
%     *    representing another complex number and returns the       *
%     *    complex product.                                          *
%     *                                                              *
%     *  Subprograms called: ARMULT, ARSUB, ARADD                    *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,cr,ci,wk1,wk2,cr2,d1,d2,wk6,l,rmax]=cmpmul(ar,ai,br,bi,cr,ci,wk1,wk2,cr2,d1,d2,wk6,l,rmax);
%
%
%
%
i=0;
%
%
[ar,br,d1,wk6,l,rmax]=armult(ar,br,d1,wk6,l,rmax);
[ai,bi,d2,wk6,l,rmax]=armult(ai,bi,d2,wk6,l,rmax);
[d1,d2,cr2,wk1,wk2,l,rmax]=arsub(d1,d2,cr2,wk1,wk2,l,rmax);
[ar,bi,d1,wk6,l,rmax]=armult(ar,bi,d1,wk6,l,rmax);
[ai,br,d2,wk6,l,rmax]=armult(ai,br,d2,wk6,l,rmax);
[d1,d2,ci,wk1,l,rmax]=aradd(d1,d2,ci,wk1,l,rmax);
for i=-1 : l+1;
 cr(i+2)=cr2(i+2);
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARYDIV                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Returns the double precision complex number   *
%     *    resulting from the division of four arrays, representing  *
%     *    two complex numbers.  The number returned will be in one  *
%     *    of two different forms:  either standard scientific or as *
%     *    the log (base 10) of the number.                          *
%     *                                                              *
%     *  Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT.    *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,c,l,lnpfq,rmax,ibit]=arydiv(ar,ai,br,bi,c,l,lnpfq,rmax,ibit);
%
%
cdum=[];ae=[];be=[];ce=[];n1=[];e1=[];n2=[];e2=[];n3=[];e3=[];
global  zero   half   one   two   ten   eps;
%
%
%
%
ae=zeros(2,2);be=zeros(2,2);ce=zeros(2,2);dum1=0;dum2=0;e1=0;e2=0;e3=0;n1=0;n2=0;n3=0;phi=0;ri10=0;rr10=0;tenmax=0;x=0;x1=0;x2=0;
cdum=0;
%
dnum=0;
ii10=0;ir10=0;itnmax=0;rexp=0;
%
%
%
rexp=fix(ibit./2);
x=rexp.*(ar(l+1+2)-2);
rr10=x.*log10(two)./log10(ten);
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(ai(l+1+2)-2);
ri10=x.*log10(two)./log10(ten);
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(ar(1+2).*rmax.*rmax+ar(2+2).*rmax+ar(3+2)).*sign(ar(-1+2)));
dum2=(abs(ai(1+2).*rmax.*rmax+ai(2+2).*rmax+ai(3+2)).*sign(ai(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum=complex(dum1,dum2);
[cdum,ae]=conv12(cdum,ae);
ae(1,2)=ae(1,2)+ir10;
ae(2,2)=ae(2,2)+ii10;
x=rexp.*(br(l+1+2)-2);
rr10=x.*log10(two)./log10(ten);
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(bi(l+1+2)-2);
ri10=x.*log10(two)./log10(ten);
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(br(1+2).*rmax.*rmax+br(2+2).*rmax+br(3+2)).*sign(br(-1+2)));
dum2=(abs(bi(1+2).*rmax.*rmax+bi(2+2).*rmax+bi(3+2)).*sign(bi(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum=complex(dum1,dum2);
[cdum,be]=conv12(cdum,be);
be(1,2)=be(1,2)+ir10;
be(2,2)=be(2,2)+ii10;
[ae,be,ce]=ecpdiv(ae,be,ce);
if (lnpfq==0) ;
 [ce,c]=conv21(ce,c);
else;
 [ce(1,1),ce(1,2),ce(1,1),ce(1,2),n1,e1]=emult(ce(1,1),ce(1,2),ce(1,1),ce(1,2),n1,e1);
 [ce(2,1),ce(2,2),ce(2,1),ce(2,2),n2,e2]=emult(ce(2,1),ce(2,2),ce(2,1),ce(2,2),n2,e2);
 [n1,e1,n2,e2,n3,e3]=eadd(n1,e1,n2,e2,n3,e3);
 n1=ce(1,1);
 e1=ce(1,2)-ce(2,2);
 x2=ce(2,1);
 %
 %      TENMAX - MAXIMUM SIZE OF EXPONENT OF 10
 %
 %      THE FOLLOWING CODE CAN BE USED TO DETERMINE TENMAX, BUT IT
 %
 %      WILL LIKELY GENERATE AN IEEE FLOATING POINT UNDERFLOW ERROR
 %
 %      ON A SUN WORKSTATION.  REPLACE TENMAX WITH THE VALUE APPROPRIATE
 %
 %      FOR YOUR MACHINE.
 %
 %
 tenmax=320;
 itnmax=1;
 dnum=0.1d0;
 itnmax=itnmax+1;
 dnum=dnum.*0.1d0;
 while (dnum>0.0);
  itnmax=itnmax+1;
  dnum=dnum.*0.1d0;
 end;
 itnmax=itnmax-1;
 tenmax=real(itnmax);
 %
 if (e1>tenmax) ;
  x1=tenmax;
 elseif (e1<-tenmax);
  x1=0.0;
 else;
  x1=n1.*(ten.^e1);
 end;
 if (x2~=0.0) ;
  phi=atan2(x2,x1);
 else;
  phi=0.0;
 end;
 c=complex(half.*(log(n3)+e3.*log(ten)),phi);
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE EMULT                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes one base and exponent and multiplies it *
%     *    by another numbers base and exponent to give the product  *
%     *    in the form of base and exponent.                         *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=emult(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
nf=n1.*n2;
ef=e1+e2;
if (abs(nf)>=ten) ;
 nf=nf./ten;
 ef=ef+one;
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE EDIV                              *
%     *                                                              *
%     *                                                              *
%     *  Description : returns the solution in the form of base and  *
%     *    exponent of the division of two exponential numbers.      *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=ediv(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
nf=n1./n2;
ef=e1-e2;
if ((abs(nf)<one) & (nf~=zero)) ;
 nf=nf.*ten;
 ef=ef-one;
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE EADD                              *
%     *                                                              *
%     *                                                              *
%     *  Description : Returns the sum of two numbers in the form    *
%     *    of a base and an exponent.                                *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=eadd(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
ediff=0;
%
%
ediff=e1-e2;
if (ediff>36.0d0) ;
 nf=n1;
 ef=e1;
elseif (ediff<-36.0d0);
 nf=n2;
 ef=e2;
else;
 nf=n1.*(ten.^ediff)+n2;
 ef=e2;
 while (1);
  if (abs(nf)<ten) ;
   while ((abs(nf)<one) & (nf~=0.0));
    nf=nf.*ten;
    ef=ef-one;
   end;
   break;
  else;
   nf=nf./ten;
   ef=ef+one;
  end;
 end;
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ESUB                              *
%     *                                                              *
%     *                                                              *
%     *  Description : Returns the solution to the subtraction of    *
%     *    two numbers in the form of base and exponent.             *
%     *                                                              *
%     *  Subprograms called: EADD                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=esub(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
[n1,e1,dumvar3,e2,nf,ef]=eadd(n1,e1,n2.*(-one),e2,nf,ef);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CONV12                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Converts a number from complex notation to a  *
%     *    form of a 2x2 real array.                                 *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [cn,cae]=conv12(cn,cae);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
%
%
cae(1,1)=real(cn);
cae(1,2)=0.0;
while (1);
 if (abs(cae(1,1))<ten) ;
  while (1);
   if ((abs(cae(1,1))>=one) | (cae(1,1)==0.0)) ;
    cae(2,1)=imag(cn);
    cae(2,2)=0.0;
    while (1);
     if (abs(cae(2,1))<ten) ;
      while ((abs(cae(2,1))<one) & (cae(2,1)~=0.0));
       cae(2,1)=cae(2,1).*ten;
       cae(2,2)=cae(2,2)-one;
      end;
      break;
     else;
      cae(2,1)=cae(2,1)./ten;
      cae(2,2)=cae(2,2)+one;
     end;
    end;
    break;
   else;
    cae(1,1)=cae(1,1).*ten;
    cae(1,2)=cae(1,2)-one;
   end;
  end;
  break;
 else;
  cae(1,1)=cae(1,1)./ten;
  cae(1,2)=cae(1,2)+one;
 end;
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CONV21                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Converts a number represented in a 2x2 real   *
%     *    array to the form of a complex number.                    *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [cae,cn]=conv21(cae,cn);
%
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
dnum=0;tenmax=0;
itnmax=0;
%
%
%     TENMAX - MAXIMUM SIZE OF EXPONENT OF 10
%
itnmax=1;
dnum=0.1d0;
itnmax=itnmax+1;
dnum=dnum.*0.1d0;
while  (dnum>0.0);
 itnmax=itnmax+1;
 dnum=dnum.*0.1d0;
end;
itnmax=itnmax-2;
tenmax=real(itnmax);
%
if (cae(1,2)>tenmax | cae(2,2)>tenmax) ;
 %      CN=CMPLX(TENMAX,TENMAX)
 %
 itnmax,
 %format (' error - value of exponent required for summation',' was larger',./,' than the maximum machine exponent ',1i3,./,[' suggestions:'],./,' 1) re-run using lnpfq=1.',./,' 2) if you are using a vax, try using the',' fortran./g_floating option');
 error('stop encountered in original fortran code');
elseif (cae(2,2)<-tenmax);
 cn=complex(cae(1,1).*(10.^cae(1,2)),0.0);
else;
 cn=complex(cae(1,1).*(10.^cae(1,2)),cae(2,1).*(10.^cae(2,2)));
end;
return;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ECPMUL                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Multiplies two numbers which are each         *
%     *    represented in the form of a two by two array and returns *
%     *    the solution in the same form.                            *
%     *                                                              *
%     *  Subprograms called: EMULT, ESUB, EADD                       *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c]=ecpmul(a,b,c);
%
%
n1=[];e1=[];n2=[];e2=[];c2=[];
c2=zeros(2,2);e1=0;e2=0;n1=0;n2=0;
%
%
[a(1,1),a(1,2),b(1,1),b(1,2),n1,e1]=emult(a(1,1),a(1,2),b(1,1),b(1,2),n1,e1);
[a(2,1),a(2,2),b(2,1),b(2,2),n2,e2]=emult(a(2,1),a(2,2),b(2,1),b(2,2),n2,e2);
[n1,e1,n2,e2,c2(1,1),c2(1,2)]=esub(n1,e1,n2,e2,c2(1,1),c2(1,2));
[a(1,1),a(1,2),b(2,1),b(2,2),n1,e1]=emult(a(1,1),a(1,2),b(2,1),b(2,2),n1,e1);
[a(2,1),a(2,2),b(1,1),b(1,2),n2,e2]=emult(a(2,1),a(2,2),b(1,1),b(1,2),n2,e2);
[n1,e1,n2,e2,c(2,1),c(2,2)]=eadd(n1,e1,n2,e2,c(2,1),c(2,2));
c(1,1)=c2(1,1);
c(1,2)=c2(1,2);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ECPDIV                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Divides two numbers and returns the solution. *
%     *    All numbers are represented by a 2x2 array.               *
%     *                                                              *
%     *  Subprograms called: EADD, ECPMUL, EDIV, EMULT               *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c]=ecpdiv(a,b,c);
%
%
b2=[];c2=[];n1=[];e1=[];n2=[];e2=[];n3=[];e3=[];
global  zero   half   one   two   ten   eps;
%
b2=zeros(2,2);c2=zeros(2,2);e1=0;e2=0;e3=0;n1=0;n2=0;n3=0;
%
%
b2(1,1)=b(1,1);
b2(1,2)=b(1,2);
b2(2,1)=-one.*b(2,1);
b2(2,2)=b(2,2);
[a,b2,c2]=ecpmul(a,b2,c2);
[b(1,1),b(1,2),b(1,1),b(1,2),n1,e1]=emult(b(1,1),b(1,2),b(1,1),b(1,2),n1,e1);
[b(2,1),b(2,2),b(2,1),b(2,2),n2,e2]=emult(b(2,1),b(2,2),b(2,1),b(2,2),n2,e2);
[n1,e1,n2,e2,n3,e3]=eadd(n1,e1,n2,e2,n3,e3);
[c2(1,1),c2(1,2),n3,e3,c(1,1),c(1,2)]=ediv(c2(1,1),c2(1,2),n3,e3,c(1,1),c(1,2));
[c2(2,1),c2(2,2),n3,e3,c(2,1),c(2,2)]=ediv(c2(2,1),c2(2,2),n3,e3,c(2,1),c(2,2));


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION IPREMAX                           *
%     *                                                              *
%     *                                                              *
%     *  Description : Predicts the maximum term in the pFq series   *
%     *    via a simple scanning of arguments.                       *
%     *                                                              *
%     *  Subprograms called: none.                                   *
%     *                                                              *
%     ****************************************************************
%
function [ipremax]=ipremax(a,b,ip,iq,z);
%
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
%
%
expon=0;xl=0;xmax=0;xterm=0;
%
i=0;j=0;
%
xterm=0;
for j=1 : 100000;
 %
 %      Estimate the exponent of the maximum term in the pFq series.
 %
 %
 expon=zero;
 xl=real(j);
 for i=1 : ip;
  expon=expon+real(factor(a(i)+xl-one))-real(factor(a(i)-one));
 end;
 for i=1 : iq;
  expon=expon-real(factor(b(i)+xl-one))+real(factor(b(i)-one));
 end;
 expon=expon+xl.*real(log(z))-real(factor(complex(xl,zero)));
 xmax=log10(exp(one)).*expon;
 if ((xmax<xterm) & (j>2)) ;
  ipremax=j;
  return;
 end;
 xterm=max([xmax,xterm]);
end;
' error in ipremax--did not find maximum exponent',
error('stop encountered in original fortran code');


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION FACTOR                            *
%     *                                                              *
%     *                                                              *
%     *  Description : This function is the log of the factorial.    *
%     *                                                              *
%     *  Subprograms called: none.                                   *
%     *                                                              *
%     ****************************************************************
%
function [factor]=factor(z);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
pi=0;
%
if (((real(z)==one) & (imag(z)==zero)) | (abs(z)==zero)) ;
 factor=complex(zero,zero);
 return;
end;
pi=two.*two.*atan(one);
factor=half.*log(two.*pi)+(z+half).*log(z)-z+(one./(12.0d0.*z)).*(one-(one./(30.d0.*z.*z)).*(one-(two./(7.0d0.*z.*z))));


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION CGAMMA                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Calculates the complex gamma function.  Based *
%     *     on a program written by F.A. Parpia published in Computer*
%     *     Physics Communications as the `GRASP2' program (public   *
%     *     domain).                                                 *
%     *                                                              *
%     *                                                              *
%     *  Subprograms called: none.                                   *
%     *                                                              *
%     ****************************************************************
function [cgamma]=cgamma(arg,lnpfq);
%
%
%
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
argi=0;argr=0;argui=0;argui2=0;argum=0;argur=0;argur2=0;clngi=0;clngr=0;diff=0;dnum=0;expmax=0;fac=0;facneg=0;fd=zeros(7,1);fn=zeros(7,1);hlntpi=0;obasq=0;obasqi=0;obasqr=0;ovlfac=0;ovlfi=0;ovlfr=0;pi=0;precis=0;resi=0;resr=0;tenmax=0;tenth=0;termi=0;termr=0;twoi=0;zfaci=0;zfacr=0;
%
first=0;negarg=0;
i=0;itnmax=0;
%
%
%
%----------------------------------------------------------------------*
%     *
%     THESE ARE THE BERNOULLI NUMBERS B02, B04, ..., B14, EXPRESSED AS *
%     RATIONAL NUMBERS. FROM ABRAMOWITZ AND STEGUN, P. 810.            *
%     *
fn=[1.0d00,-1.0d00,1.0d00,-1.0d00,5.0d00,-691.0d00,7.0d00];
fd=[6.0d00,30.0d00,42.0d00,30.0d00,66.0d00,2730.0d00,6.0d00];
%
%----------------------------------------------------------------------*
%
hlntpi=[1.0d00];
%
first=[true];
%
tenth=[0.1d00];
%
argr=real(arg);
argi=imag(arg);
%
%     ON THE FIRST ENTRY TO THIS ROUTINE, SET UP THE CONSTANTS REQUIRED
%     FOR THE REFLECTION FORMULA (CF. ABRAMOWITZ AND STEGUN 6.1.17) AND
%     STIRLING'S APPROXIMATION (CF. ABRAMOWITZ AND STEGUN 6.1.40).
%
if (first) ;
 pi=4.0d0.*atan(one);
 %
 %      SET THE MACHINE-DEPENDENT PARAMETERS:
 %
 %
 %      TENMAX - MAXIMUM SIZE OF EXPONENT OF 10
 %
 %
 itnmax=1;
 dnum=tenth;
 itnmax=itnmax+1;
 dnum=dnum.*tenth;
 while (dnum>0.0);
  itnmax=itnmax+1;
  dnum=dnum.*tenth;
 end;
 itnmax=itnmax-1;
 tenmax=real(itnmax);
 %
 %      EXPMAX - MAXIMUM SIZE OF EXPONENT OF E
 %
 %
 dnum=tenth.^itnmax;
 expmax=-log(dnum);
 %
 %      PRECIS - MACHINE PRECISION
 %
 %
 precis=one;
 precis=precis./two;
 dnum=precis+one;
 while (dnum>one);
  precis=precis./two;
  dnum=precis+one;
 end;
 precis=two.*precis;
 %
 hlntpi=half.*log(two.*pi);
 %
 for i=1 : 7;
  fn(i)=fn(i)./fd(i);
  twoi=two.*real(i);
  fn(i)=fn(i)./(twoi.*(twoi-one));
 end;
 %
 first=false;
 %
end;
%
%     CASES WHERE THE ARGUMENT IS REAL
%
if (argi==0.0) ;
 %
 %      CASES WHERE THE ARGUMENT IS REAL AND NEGATIVE
 %
 %
 if (argr<=0.0) ;
  %
  %       STOP WITH AN ERROR MESSAGE IF THE ARGUMENT IS TOO NEAR A POLE
  %
  %
  diff=abs(real(round(argr))-argr);
  if (diff<=two.*precis) ;
   ,
   argr , argi,
   %format (' argument (',1p,1d14.7,',',1d14.7,') too close to a',' pole.');
   error('stop encountered in original fortran code');
  else;
   %
   %        OTHERWISE USE THE REFLECTION FORMULA (ABRAMOWITZ AND STEGUN 6.1
   %        .17)
   %        TO ENSURE THAT THE ARGUMENT IS SUITABLE FOR STIRLING'S
   %
   %        FORMULA
   %
   %
   argum=pi./(-argr.*sin(pi.*argr));
   if (argum<0.0) ;
    argum=-argum;
    clngi=pi;
   else;
    clngi=0.0;
   end;
   facneg=log(argum);
   argur=-argr;
   negarg=true;
   %
  end;
  %
  %       CASES WHERE THE ARGUMENT IS REAL AND POSITIVE
  %
  %
 else;
  %
  clngi=0.0;
  argur=argr;
  negarg=false;
  %
 end;
 %
 %      USE ABRAMOWITZ AND STEGUN FORMULA 6.1.15 TO ENSURE THAT
 %
 %      THE ARGUMENT IN STIRLING'S FORMULA IS GREATER THAN 10
 %
 %
 ovlfac=one;
 while (argur<ten);
  ovlfac=ovlfac.*argur;
  argur=argur+one;
 end;
 %
 %      NOW USE STIRLING'S FORMULA TO COMPUTE LOG (GAMMA (ARGUM))
 %
 %
 clngr=(argur-half).*log(argur)-argur+hlntpi;
 fac=argur;
 obasq=one./(argur.*argur);
 for i=1 : 7;
  fac=fac.*obasq;
  clngr=clngr+fn(i).*fac;
 end;
 %
 %      INCLUDE THE CONTRIBUTIONS FROM THE RECURRENCE AND REFLECTION
 %
 %      FORMULAE
 %
 %
 clngr=clngr-log(ovlfac);
 if (negarg) ;
  clngr=facneg-clngr;
 end;
 %
else;
 %
 %      CASES WHERE THE ARGUMENT IS COMPLEX
 %
 %
 argur=argr;
 argui=argi;
 argui2=argui.*argui;
 %
 %      USE THE RECURRENCE FORMULA (ABRAMOWITZ AND STEGUN 6.1.15)
 %
 %      TO ENSURE THAT THE MAGNITUDE OF THE ARGUMENT IN STIRLING'S
 %
 %      FORMULA IS GREATER THAN 10
 %
 %
 ovlfr=one;
 ovlfi=0.0;
 argum=sqrt(argur.*argur+argui2);
 while (argum<ten);
  termr=ovlfr.*argur-ovlfi.*argui;
  termi=ovlfr.*argui+ovlfi.*argur;
  ovlfr=termr;
  ovlfi=termi;
  argur=argur+one;
  argum=sqrt(argur.*argur+argui2);
 end;
 %
 %      NOW USE STIRLING'S FORMULA TO COMPUTE LOG (GAMMA (ARGUM))
 %
 %
 argur2=argur.*argur;
 termr=half.*log(argur2+argui2);
 termi=atan2(argui,argur);
 clngr=(argur-half).*termr-argui.*termi-argur+hlntpi;
 clngi=(argur-half).*termi+argui.*termr-argui;
 fac=(argur2+argui2).^(-2);
 obasqr=(argur2-argui2).*fac;
 obasqi=-two.*argur.*argui.*fac;
 zfacr=argur;
 zfaci=argui;
 for i=1 : 7;
  termr=zfacr.*obasqr-zfaci.*obasqi;
  termi=zfacr.*obasqi+zfaci.*obasqr;
  fac=fn(i);
  clngr=clngr+termr.*fac;
  clngi=clngi+termi.*fac;
  zfacr=termr;
  zfaci=termi;
 end;
 %
 %      ADD IN THE RELEVANT PIECES FROM THE RECURRENCE FORMULA
 %
 %
 clngr=clngr-half.*log(ovlfr.*ovlfr+ovlfi.*ovlfi);
 clngi=clngi-atan2(ovlfi,ovlfr);
 %
end;
if (lnpfq==1) ;
 cgamma=complex(clngr,clngi);
 return;
end;
%
%     NOW EXPONENTIATE THE COMPLEX LOG GAMMA FUNCTION TO GET
%     THE COMPLEX GAMMA FUNCTION
%
if ((clngr<=expmax) & (clngr>=-expmax)) ;
 fac=exp(clngr);
else;
 ,
 clngr,
 %format (' argument to exponential function (',1p,1d14.7,') out of range.');
 error('stop encountered in original fortran code');
end;
resr=fac.*cos(clngi);
resi=fac.*sin(clngi);
cgamma=complex(resr,resi);
%
return;

Contact us